{'fraction_non_alphanumeric': 0.05098849038841143, 'fraction_numerical': 0.015518236205168696, 'mean_word_length': 3.9200657403257133, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the rigid-lid approximation of shallow water Bingham

This paper discusses the well posedness of an initial value problem describing the motion of a Bingham fluid in a basin with a degenerate bottom topography. A physical interpretation of such motion is discussed. The system governing such motion is obtained from the Shallow Water-Bingham models in the regime where the Froude number degenerates, i.e taking the limit of such equations as the Froude number tends to zero. Since we are considering equations with degenerate coefficients, then we shall work with weighted Sobolev spaces in order to establish the existence of a weak solution. In order to overcome the difficulty of the discontinuity in Bingham's constitutive law, we follow a similar approach to that introduced in [G. DUVAUT and J.-L. LIONS, Springer-Verlag, 1976]. We study also the behavior of this solution when the yield limit vanishes. Finally, a numerical scheme for the system in 1D is furnished.

1. Introduction. Bingham fluids constitute a crucial topic of study on which many applicable researches are conducted. As an interesting example we suggest the snow that appears in some important natural phenomena such as avalanches. Technically speaking, the avalanche dynamics can be described through different perspectives, such as the center of mass consideration [27] and the density one. In addition, many approaches consider the avalanche to be a deformable body whereas others describe it as a granular material. Not only the dynamics of snow, but also its constitutive behavior can be described according to different points of views such as Newtonian fluid and Bingham fluid. In many literatures, snow is considered to be a non-Newtonian fluid. One of the reasons is that Newtonian fluids adapt rapidly to deform themselves where they reach a negligible depth in an unbounded space, on the other hand, snow will pile in such a case and thus will rest with a finite depth. This can be explained by the fact that upon being in rest state, snow achieves a yield value with a non zero shear stress, and thus a threshold stress value should be exceeded to start deformation. Another reason is the experimental results of the analysis of the avalanche's velocity profile along its depth, which reveals the viscosity dependence of the shear [24] [26]. These two reasons triggered many authors in the literature to treat snow as a Bingham fluid which is characterized mainly by a non constant stress that depends on viscosity and that may differ in the same body of mass. Hence, one would notice some portions flowing while others moving in a bulk motion as solids. More precisely, at low shear stress, Bingham fluids possess high viscosity, and thus they behave as a shear thinning fluid at low or zero speed. As the shear stress reaches a certain limit-denoted yield stress-one notices a sudden drop in the viscosity, while above the yield stress, snow behaves like a low viscosity liquid.
The need to use reduced models arises from the difficulty that the three dimensional Cauchy momentum equations induce in the analytical and numerical studies, especially if we consider a fee boundary case. In this sequel, the reduced model adopted is the Shallow water-Bingham model that is derived using depth integration of the Navier-Stokes-Bingham system with free surface, taking into account that the horizontal length scale is much greater than the vertical one. Since there is no well posedness result on the initial system and the approximated one, then any mathematical justification of this procedure is still as far as we know an interesting open problem due to the wide range of applicability especially in numerical applications. For instance, in [9], [16], the authors proposed several models of a Shallow water type system for Bingham fluid. The resulting system is given by ∂ t h + div(hv) = 0 ∂ t (hv) + div(hv ⊗ v) − div(σ) + 1 2 Fr 2 ∇h 2 = hf, where h(t, x) is the water height, v(t, x) is the horizontal velocity. We denote by hv ⊗ v the matrix with component hv i v j , Fr is the Froude number and σ is the stress tensor given by where D(v) represents the symmetric part of the velocity gradient given by As mentioned above, the Shallow water Bingham model is still for today an open problem. We remark that, in general, the usual strategy of the existence theory is based on two steps: the first one consists of finding u n , a solution of a suitably chosen approximating problem (a Galerkin approximation for example). The second one amounts to establish uniform estimates in proper spaces and one can then obtains u, a solution to System (3)-(4), as a "weak" limit of a suitably chosen subsequence of u n . The central difficulty is the passage to the limit in the nonlinear terms and in the discontinuous constitutive law. Actually, the ideas developed in studying the well posedness issue of Shallow water equations (or compressible Navier-Stokes equations) for Newtonian fluids (see [22], [8], [30], [15]) don't work for non Newtonian ones. As a forward step in this approach, we study here the rigid lid approximation of such system, that we shall call Bingham lake equations.
Bingham lake equations characterizing lake equations decode many natural as well as industrial phenomena. The most interesting example is that of mixed flows of Bingham type fluids (such as petroleum) in closed supply pipes. Two types of flows are exhibited inhere characterizing a transition phase from free surface (i.e when only part of the section of pipe is filled and the pressure is being known: atmospheric pressure), to pressurized flow (i.e when the section of the pipe is full, and the pressure is an unknown). The authors in [5] have studied such kind of flows occurring for Newtonian fluids thoroughly. They presented a new model called PFS model, based on coupling the free surface part equations derived from the incompressible Navier Stokes or Euler systems, and the pressurized part equations derived from the Compressible Euler equations. A finite volume discretization have been studied in [5], and a kinetic formulation of such models is presented in [6]. Mathematically, this system can be obtained from (1)-(2) system by passing to the limit (Fr → 0), where the initial height converges to a non constant function b(x) depending on the space variables only. The obtained model is given by: The shear stress σ satisfies the special constitutive law of a Bingham fluid (see [3]): Here, u(t, x) denotes the velocity vector, p the pressure, f (t, x) the known external force, µ and λ the Lamé viscosity coefficients, g the yield limit. We couple the system (3) with the so called Lions' boundary conditions given by and with the initial data (defined in a weak sence, see Theorem 4.2) In (5), n and τ denote respectively the unit normal and tangential vectors to the boundary, whereas κ(x) is the curvature of ∂Ω. Lions' boundary conditions, the particular case of Navier boundary condition that were first used by Navier in 1872, regard that there is a stagnant layer of fluid close to the wall allowing a fluid to slip, and the slip velocity is proportional to the shear stress.
Since the singular term D(u)/| D(u)| is not always defined, numerical and mathematical obstacles appear, which forced many authors to develop new formulations of the problem in order to tackle the difficulty. One of the approaches used was introduced in [10] for an unsteady flow of an incompressible fluid having its Cauchy stress given implicitly and relating both the symmetric part of the velocity gradient D(u) and the deviatoric part of the Cauchy stress. The authors regard the stress tensor σ as a new variable (along with the density and the velocity vector), and this results in a dissipation rate in the form of a Young function depending on D(u) and σ t . The second approach, that we will adopt in the sequel, was initiated by G. Duvaut and J. L. Lions [12] who replaced the original system by a variational inequality to get rid of the singularity.
In this paper, we will first adopt the methodology of variational inequality used in [12] to prove the existence of a weak solution of the incompressible Bingham fluid confined to a shallow basin with a varying bottom topography but with some changes in the nature of spaces used. Again recalling that we are dealing with a degenerate bathymetry b(x), meaning that b(x) may vanish on the boundary, we will prove that the solution exists in some weighted Sobolev space where the weight is assumed to be a Muckenhoupt type. For more details about these spaces, we refer the reader to [23]. Then, we will discuss the behavior of solution when g tends to zero. It is good to mention that the vanishing viscosity limit (µ tends to zero) of Bingham fluids is studied by J.-L. Lions [21] in a bounded open set of R 2 without bottom topography. However, an improvement of weak regularity on the solution is required in order to achieve the convergence on the nonlinear term. Yet, as long as we can't improve the regularity of our solution due to the degeneracy of our equations, then we are not able to study such limit in our case.
The paper consists of 6 parts. In the next part (Section 2), we will introduce the spaces that we prove the solution in and some preliminary results related to the functional analysis of the problem. In addition, we will accomplish the equivalence between the initial system and the variational inequality introduced in Section 3. Then comes Section 4 where we exhibit the existence of solution via a Galerkin method. In Sections 5, we study the behavior of the solution when g tends to zero. Finally, in Section 6, a numerical scheme is implemented in 1D verifying our theoretical results.
2. Functional spaces. Before we start the analysis, and since we are going to work with weighted Sobolev spaces, let us first give a brief definition of some of the spaces and preliminaries that we shall need in the sequel, especially the constraints on the weight b and the domain Ω. Domain: (I) For an integer m, 1 ≤ m ≤ 2, we set Q m = (0, 1) m . We assume that there exists a bi-Lipschitz mapping B : Q 2 → Ω, such that B(Q 1 ) = ∂Ω.
Weight: (II) We define a space function b(x), locally integrable and belonging to the Muckenhoupt class A q . Generally speaking, for a weight of Muckenhoupt type, the definition of trace operator is well defined. More precisely, one can check that if b ∈ A q , we have u ∈ W 1,q b (Ω) → W 1,q loc (Ω) and hence there is a linear trace operator γ 0,b : W 1,q loc (Ω) → L 1 loc (Ω). Though it is well defined, yet we lack characterization of such trace regarding a general Muckenhoupt weight. That's why we restrict ourselves in what follows to a more specific weight that provides a characterization of the boundary terms. Its expression is given in a neighborhood of the boundary In this situation, the definition of the trace is more accurate. For more details about examples of weights satisfying Muckenhoupt condition, we refer the reader to [14], [18], [29].
Weighted Sobolev spaces: As mentioned in the introduction, we will prove our solution in Weighted Sobolev space V b . We introduce the following weighted spaces: For more details about Muckenhoupt classes and Weighted Sobolev spaces we refer the reader to [23] and [2].
Trace: As mentioned above in (7), we will adopt a special choice of Muckenhopt type weights that provide a good characterization of the trace. In fact, this result is proved by A. Nekvinda via the following theorem.
such that T 1,q b ( ∂Ω)(u) = u | ∂Ω . As a consequence of this theorem, we remark that if f ∈ H 1 ρ α (Ω), 0 < α < 1 the trace of f is well defined and belong to L 2 ( ∂Ω) → L 2 b ( ∂Ω). This helps us define the boundary integrals coming from the Navier boundary conditions. Throughout this paper, we assume that the domain Ω satisfies Hypothesis (I), and b satisfies Hypothesis (II), and any other restrictions would be specified.
Notation Remark: Let us fix some notations which will be used throughout the sequel: • We say that u is a b-divergence free or u satisfies the b-incompressibility condition if div(bu) = 0. Remark 1. We should remark that in the sequel of the study, we will need to have q = 3 2 , i.e the weight b should belong to A 3 2 . In fact, this choice is for the sake of bounding the nonlinear term (u · ∇)u in lemma 4.1. Notice that in studying the Stokes problem, we just need to take q = 2 which is necessary to have a weighted version of Poincaré (or Korn's inequality). The reader can refer to lemma 4.1 in [2] for a detailed explanation of such choice.
Without loss of generality, C will denote a positive generic constant in the following literature.
3. Variational inequality. Following Duvaut and Lions in [12], we will derive in this section a variational inequality and prove (at least formally) that the resolution of this variational inequality is equivalent to solve the problem (3)-(4) in a weak sense (solution of variational problem). Indeed, let us establish first the variational formulation of system (3)- (4). For this sake, we consider a set V of functions v : Ω → R 2 of enough regularity and such that div(bv) = 0 in Ω and bv · n = 0 in ∂Ω. Then, we multiply (3) by (v − u) and integrate in space. We get Accounting for the Navier boundary conditions and the relation div(bv) = 0, we can write Thus The physical understanding of this term indicates that it should be interpreted when D(u) = 0 as "any trace-free symmetric matrix with norm less or equal to one". In the mathematical language such quantity is called "multivalued". This fundamental difficulty makes the approximation of the problem (8) already a complex challenge , and has motivated a large literature, see [12], [10] and references therein. Below, we follow the procedure in [12] to derive a "formal" equivalent form of the variational formulation (8). Indeed, regarding the Cauchy-Schwartz inequality Hence we get By now we have proved that if (u, p) is a formal solution of (3), then it would satisfy (9). Reciprocally, assuming that u is a solution of the variational inequality (9) such that D(u) = 0 a.e. in Ω. Letting θw = v − u in the previous inequality, one gets Now, as we divide by θ we get: Concerning the term S, we have in fact: On the other hand, we have By changing w into −w, we find the opposite inequality of (10), and so we get Thus, we establish formally the equivalence between (9) and (3).
Let us now precise our definition of weak solution of System (3)-(4).
Definition 3.1. We say that u is a solution of system (3)-(4) equipped with boundary and initial conditions given in (5)-(6) if u satisfies the following regularity , and the following variational inequality holds for all Remark 2. The well definition of the right hand side will be discussed in the sequel as we take u and v ∈ V b (mainly in lemma 4.1). It is good to mention here that for left hand side term in (12), one should pay attention for the regularity of bf, hence assuming the least regularity possible, i.e in the dual of V b , we can replace the integral form by the dual representation: bf, v − u V b ,V b . Though, this in fact doesn't affect the computations. Without loss of generality, we will assume in the sequel that bf is in ( 4. Main results. We state in this section the existence result of problem (12). The presence of b in the diffusion operator and the "b-incompressibility condition" make the weighted Sobolev spaces the ambient ones to prove existence within. In what follows, we introduce a priori estimates concerning the non linear term in lemma 4.1, and which will later on serve the well definition of the integral forms. Next, we exhibit in theorem 4.2 the existence and uniqueness results. We will rely in the proof on a Galerkin approximation technique. Several operators will be used in the latter method for which different properties will be given in lemma 4.1.
, and there exists a positive real number C independent of u and v such that We define the following linear operators A and B such that where A and B are given by Proof. Here we just want to prove that Au ∈ L 2 (0, T ; V b ) and Bu ∈ L 2 (0, T ; V b ). The other properties are proved in [2]. Indeed, For u ∈ L ∞ (0, and v ∈ V b , we estimate using Hölder's inequality Since u ∈ L 2 (0, T ; V b ), thus by duality we get As for the operator B, we estimate the nonlinear term using Hölder and Gagliardo-Nirenberg inequalities as follows Remark 3. Notice that, in [Proposition 3.1, [2]] , the first author proved a weighted version of Korn's inequality. However, some constraints on the weight b and the domain Ω were needed. More precisely, the author proved that if we exclude that case when the domain Ω is a disc, b is radial and not identically zero on the boundary, then we have Nevertheless, he showed also that these assumptions are not required when studying the evolution problem. That is why in our case we won't suppose such assumptions. For more details, the reader is referred to [2].
The main result of this paper is given below. We suppose that f and u 0 are the applied force and the initial datum given such that f lies in L 2 (0, T ; L 2 b (Ω)) and u 0 belongs to H b . Assume that b satisfies Hypothesis (II), and that κ(x) is in L ∞ (Ω), then there exists a unique vector field u such that and the initial condition is defined in a weak sense: Proof. In the same spirit of Lions and Duvaut's approach, our proof will be composed of fourth main steps. First, let us denote by the operator j the integral form: Step1: Regularizing j. In attempt to approximate the problem, we start first by regularizing the operator j which is derived from Bingham's singular term.
Let ε ∈ (0, 1). We approximate j by a differentiable functional In fact, j (ψ) is well defined since as we deal with a bounded domain and the fact that D(ψ) is in L 2 b (Ω), we infer that j (ψ)<∞. The Gateaux differential of j (·) along v is given by (11).
In fact, we claim that We adopt the notation in [12] ( It is good to mention here that j (w) ∈ L 2 (0, We have the following property Proof. We have first Using Hölder inequality, we get So the above inequality becomes 1+ε is a constant independent of time. Thus, we get Again applying Hölder inequality to the above integral we get Take lim inf ε of the above inequality, we get lim inf As v −→ j(v) is lower semi continuous for the weak topology of L 2 (0, T ; V b ) (from the definition of weak convergence), then Therefore, combining the above two inequalities we get our result We render the problem into a new approximated one corresponding to j such that: where the representation (j (·), ·) stands for (15).
Step 2: Seeking for a sequence (u m ) via Galerkin approximation. Consider the canonical isomorphism ∧ : V b −→ V b with {w 1 , ..., w m , ..} being the set of unit eigenfunctions of the operator ∧, Define now the space V m = Span{w 1 , . . . w m }, where {w 1 , . . . , w m } is a free and total family in H b (it is permissible to choose the set as that since V b is separable). In fact this isomorphism is proven for weighted Sobolev spaces by F röhlich in bounded domains, check for instance [17]. Upon projecting our system on V m , and according to the Cauchy theory in a finite space, we can construct a solution u m := u εm in V m for the approximated variational inequality (17) with u m (0) being the projection of u 0 on V m , and u m being defined in a space sum over j, we get: From Lemma 4.1, we deduce that the second term is equal to zero, and since Notice that using Hölder's and Young's inequalities, we can estimate the term on the right hand of equation (20) as follows , for some arbitrary η > 0. Now, after choosing η sufficiently small, we can deduce using weighted Korn's inequality from Remark 3 that Hence, we deduce Since we have P m ( ∂ t (bu m )) = ∂ t (bu m ), then we can write Therefore, we get Since k m is bounded in L 2 (0, T ; V b ), then what remains to show is that Since ∧ is an isomorphism from V b to V b , then we can deduce that λ 1 2 w j constitute an orthogonal basis of V b for the norm Finally, we have our aimed result bu m ∈ L 2 (0, T ; V b ).
As a result of (22) and (23), we can deduce that up to a sequence (denoted again u m ), we have Now, in order to use compactness result we must establish a fractional estimate in time. Indeed, we have the following Lemma.
where we denote by τ h u m (t) = u m (t + h).
Proof. The proof of lemma 4.4 follows the same line as in proposition 4.1. in [2]. For completeness sake, we will provide the proof again. Let h > 0, indeed, integrating We multiply equation (26) by u mj (t + h) − u mj (t), then we add the resulting equations over j, as j goes from 1 to m, we obtain that Estimating each term on the left hand side of the above equation using lemma 4.1, estimate 16 and trace theorem, we obtain that By Fubini theorem, the second member becomes Since |s − s − h| ≤ |s − h − s| = h, then Hölder inequality and stimate (22) give us Using again stimate (22), we obtain that Combining the above estimates with the one in (28) yields and hence the proof is finished.
Thus we are able to use the compactness result and deduce that u m → u ε strongly in L 2 (0, T ; H b ).
Step 3. Passage to the limit in m. The concern now is in the convergence of terms in (18). We will prove the convergence of the operators as m → ∞.
• Convergence of the linear operator A. As for the operator A given by (Ω)), and hence Before we treat the convergence on the boundary term, let us prove the following lemma.
Recall that κ is the curvature of ∂Ω where we have ∂n ∂τ But v is parallel to τ , so (31) follows by linearity. Now back to prove the convergence of the boundary term, we notice first using (31) and the following identity We easily remark that the first term on the right hand side of the above equality converges to zero because of the weak convergence of u m to u in V b . For the second term, we need to use a Hardy's type inequality to show that Indeed for b = ρ α , we have ∇b/b ∼ 1/ρ, then (see [13]) Thus we get the convergence of the boundary term.
• Convergence of the initial condition. The space V b can be viewed as Now, we choose the orthonormal base of V b such that for all j, we have w j ∈ D b (Ω). In this case, since (bu m · ∇)u m is bounded in L 2 (0, T ; (L 1 (Ω)) 2 ). Moreover we have a(u m , w j ) + b(u m , u m , w j ) ∈ L 2 (0, T ).
Thus Equation (18) shows that and hence Ω u m .w j b dx is bounded in H 1 (0, T ). On the other hand, we know that H 1 (0, T ) is compact in C u (0, T ) (see Theorem III.2.34 in [7]) which yields The above convergence holds also in L 2 (0, T ) (due to strong convergence of u m in L 2 (0, T ; H b )). In particular Notice that, by definition, we have u m (0) = P m (u 0 ), so we can write Passing to the limit in Equation (33) yields to On other hand, by virtue of Equation (17), we have for all Since u ∈ L ∞ (0, T ; H b ) and V b is dense in H b , we infer that (see Lemma 1.4 in [28]).
It remains to show that ξ = j (u ). By now we have gotten We choose a function φ such that φ in L 2 (0, T, V b ), bφ in L 2 (0, T, V b ) then introduce the following quantity One can show obviously that X m is positive. For instance, we can rewrite X m such that Substituting (18) in the above expression Again, take v = u in (35), and then substitute this latter in X ε , we get As X m ≥ 0, then X ε ≥ 0. Take now φ = u − θψ with θ ∈ R in the expression of X ε . This yields Dividing by θ we get This is true for all ψ, hence: Therefore in this step we have proved the existence of a sequence (u ) ε such that: Step 4: Passage to the limit in . Thanks to the energy estimate (21), we've got that the bounds obtained in (36) are uniform in ε. Therefore, we conclude that there exists a subsequence of u , that is still denoted the same, such that Again, by compactness theory we get u → u strongly in L 2 (0, T ; H b ). We introduce the following quantity Taking v = u in (36) (which is allowed since u ∈ L 2 (0, T ; V b )), and then substituting in Z ε , we get In fact, since j is coercive, then we get Z ε ≥ 0. Hence Due to the fact that all terms on the right hand side are lower semi continuous (norms and linear continuous maps are lower semi continuous), and using lemma 4.3, we obtain Finally, we get The last step is to prove that the above inequality not only holds in integral form but also almost everywhere in (0, T ), in particular, we have to prove Proving this in fact allows us to prove the existence of solution for system (12). Let us fix a test function v ∈ V b , and t 0 ∈ (0, T ). Define We introducev Hence taking v =v, inequality (37) becomes Using the Lebesgue theorem in the differentiation of set functions, we have But t 0 is arbitrary in (0,T), thus we get for all t ∈ (0, T ), for all v ∈ L 2 (0, T ; V b ). By this we end the proof of existence.
Uniqueness of solution. Let u 1 , u 2 be two solutions of the variational inequality (39). We take v = u 2 (t) (resp v = u 1 (t)) as a test function in the variational inequality satisfied by u 1 (resp. u 2 ). Let U = u 1 − u 2 . Now, if we add both inequalities we get Therefore, using Remark 3, and the fact that where C 1 is the constant resulting from Korn's inequality. We estimate the term on the right hand side of Inequality (40) using Weighted Gagliardo-Nirenberg inequality and Young's inequality as follows From Inequality (40) and estimate (41), we deduce that . which readily ensures that U = 0, whence u 1 − u 2 = 0 by applying a Gronwall's type inequality. The proof of uniqueness of solution is finished. 5. Newtonian fluids as a limit of Non-Newtonian fluids. Viscous lake systems for Newtonian fluids are asymptotically derived systems according to two successive approximations that are characterized by the smallness of non dimensional parameters. The first is the rigid lid approximation assuming that the typical deviation of the top of the fluid's surface from the mean level is much smaller than the typical depth. This smallness can be regarded also as a Froude number due to the dynamics of the physical case at hand. The second approximation is the shallow water approximation characterizing the ansatz of the smallness of the typical depth compared to the typical horizontal length, see for instance [20]. Nevertheless, in this section, we will prove the existence of a weak solution of the viscous lake system by passing to the zero limit of g (yield stress limit) in the variational inequality satisfied by the weak solution of Bingham system (3). Roughly speaking, one can naively conclude that the fluid will behave as a viscous liquid once g vanishes, yet for the mathematical justification we need to set on a rigorous proof. So, we will prove the following theorem.
Let Ω a bounded Lipschitz domain and b satisfying Hypothesis (II). Consider the Bingham model introduced in section 3. We suppose that all the parameters therein are fixed except for the the yield stress limit g which is assumed to be independent of other parameters and rendered to vary (consequently taken to zero). Denote by u g the constructed solution of system (12) with initial data u| t=0 = u g (0) proved in section 4. Then, if u g (0) converges to u 0 in L 2 b (Ω) as g tends to zero, then up to a sub-sequence, we have where u is a weak solution of the viscous Lake system and thus (Ω) → 0 when g tends to 0, and the solution u of the viscous lake equations enjoys the following property Proof. Let's recall the definition of the weak solution of viscous lake equations: , we find that the weak solution u g satisfies Thus, u g remains in a bounded set with respect to g in L ∞ (0, T ; H b ) ∩ L 2 (0, T ; V b ). In addition, from previous section, we have bu εg bu g ∈ L 2 (0, T ; V b ), so, bu g is also bounded in L 2 (0, T ; V b ) (the bound is uniform in g). Hence, there exists u ∈ L 2 (0, bu g bu ∈ L 2 (0, T ; V b ). Following the same compactness technique used in previous section, we get u g −→ u ∈ L 2 (0, T ; H b ), and the convergence of the corresponding operators follows as in section 4. Hence u g satisfies (Ω) + a(u g , u g ). Integrating in time gives Using the fact that RHS is lower semi continuous , and since | D(u g )| and | D(v)| are bounded in L 1 (0, T ; H b ), then Therefore, for all t in (0, T ), we have lim inf Due to the convergence of the LHS terms and the last term in the RHS of (43) (following the same strategy as in the proof of theorem (4.2) in section 4), and due to the lower semi-continuity of the norm operator for the topology of L 2 b (Ω), we obtain Thus, for all v ∈ L 2 (0, T ; V b ), and for all t in (0, T ), we get Again, using Lebesgue theory in the differentiation of set functions, we obtain for all v ∈ L 2 (0, which means that u satisfies the weak formulation of the viscous Lake system. Thus we end the proof. 6. Numerical scheme. In this section, we propose in the one dimensional case a numerical scheme for the approximation of the studied model i.e. the system (3)-(4).
6.1. Semi-discrete scheme. The system (3)-(4) can be rewritten as follows Notice that in the 1d case, the definition of σ given by (4) reduces to For the time discretisation, we denote t n = k≤n ∆t k where the time steps ∆t k will be precised later though a CFL condition. Following [11], we use an operator splitting technique resulting in a two step scheme where the quantity X n+1 satisfies the divergence free constraint The system (46)-(48) has to be completed with suitable boundary conditions that will be precised later, see paragraph 6.4. More precisely, the prediction step (46) consists in the resolution of advection diffusion equation i.e.

Discrete scheme.
To approximate the solution (X, p) T of the system (3)-(4), we use a combined finite volume/finite element framework. We assume that the computational domain is discretized with I nodes x i , i = 1, . . . , I. We denote C i the cell ( the approximate solution at time t n on the cell C i . The pressure p is discretized on a staggered grid (in fact the dual mesh if we consider the 2d case) Hence, the discrete (in space and time) version of equations (49)-(50),(51) is given by and F n i+1/2 is a numerical flux accounting for the advection term bu 2 and any classical flux e.g. full upwind, Rusanov, etc can be used. In the simulation results given at the end of this section, an upwind formula has been used and having the form with u n i+1/2 = (u n i+1 + u n i )/2. The quantities τ n+1/4 i+1/2− and τ n+1/4 i−1/2+ correspond to discretisation using the apparent topography technique (see [4]) of the source term f and we have τ n i+1/2− = (x i+1/2 − x i )f n i+1/2 , τ n i−1/2+ = (x i − x i−1/2 )f n i−1/2 .
6.4. Boundary conditions. Boundary conditions have to be defined for Eqs. (52) and (55) at both side of the domain, and we have to face two difficulties • equation (52) contains a hyperbolic part and a parabolic part whereas equation (55) is an elliptic equation. • The boundary conditions applied to (55) have to be consistent with those applied to (52).
The proposed solution has been adapted from [1] but notice that other solutions can be investigated since the coupling of the boundary conditions between a hyperbolic step and a parabolic/elliptic step is far from being obvious. Let us consider the boundary at the entry of the domain i.e. at abscissa x 1/2 . We assume the inflow is prescribed typically b 0 u n 0 = q in (t n ) where q in (t n ) is a given quantity then the definition (57) can be used to define F n 1/2 under the form F n 1/2 = u n 1/2 q in (t n )1 u n 1/2 ≥0 + X n 1 1 u n 1/2 ≤0 , with u n 1/2 = (u n 1 + q in (t n )/b 1 )/2 and assuming b 0 = b 1 > 0. For the parabolic part of equation (52) we use Dirichlet boundary conditions defined by q in (t n ).
Since the inflow q in (t n ) is prescribed, it is natural to assume Neumann boundary at the entry for p in equation (55) i.e. ∂ x p| 0 (t n+1 ) = 0. Now we consider the boundary at the exit of the domain i.e. at abscissa x I+1/2 . Assuming Neumann boundary conditions, we define F n I+1/2 = u n I X n I , and ∂ x X| I+1/2 (t n+1/4 ) = 0 with b I+1 = b I for the parabolic part. In this context, it is convenient to assume Dirichlet boundary at the exit for p in equation (55) i.e. p| I+1 (t n+1 ) = 0. For a more complete justification of the choices for the boundary conditions, the reader can refer to [1]. 6.5. Simulation results. We present now some simulations results for the model (3)  The simulated velocity profile at time t = T /2 is depicted over Fig. 2-(a) whereas the variations of the simulated pressure is given over Fig. 2-(b). The variations of σ appear over Fig.3-(b). Near the boundaries, the gradient of the velocity is very small ∂ x u ≈ 0, see Fig. 3-(b).

7.
Acknowledgment. This work was realized during the secondment of the first author in the Ange Inria team as a post-doc position.