Congested shallow water model: roof modelling in free surface flow

. We are interested in the modeling and the numerical approximation of ﬂows in the presence of a roof, for example ﬂows in sewers or under an ice ﬂoe. A shallow water model with a supplementary congestion constraint describing the roof is derived from the Navier-Stokes equations. The congestion constraint is a challenging problem for the numerical resolution of hyperbolic equations. To overcome this diﬃculty, we follow a pseudo-compressibility relaxation approach. Eventually, a numerical scheme based on a Finite Volume method is proposed. The well-balanced property and the dissipation of the mechanical energy, acting as a mathematical entropy, are ensured under a non-restrictive condition on the time step in spite of the large celerity of the potential waves in the congested areas. Simulations in one dimension for tran-scritical steady ﬂow are carried out and numerical solutions are compared to several analytical (stationary and non-stationary) solutions for validation.


Introduction
We deal with the derivation and the numerical resolution of a model describing shallow water flows in the presence of a roof. The roof denotes an impermeable surface above the flow which constrains the water surface in several parts of the domain. We do not assume that the roof is fixed, but the roof has a prescribed motion and the surface pressure is seen as an unilateral constraint. This assumption is not really essential to the analysis but greatly simplifies the reading. In particular the numerical procedure is not affected by the roof dynamics. The description of the fluid/structure interaction is studied in the literature for the Navier-Stokes equations [43]. However, for vertical-integrated models such as the shallow water model, the interaction with a roof has been considered only recently. The main applications of these models concern the prediction of urban flooding taking into account the flows in sewers or the hydraulic network in small hydroelectric power plants [13,24]. Moreover, the production of sustainable energy from wave power using floating buoys can also be realized considering a dynamical roof.

Formal derivation of the congested shallow water model
A shallow water type model with an additional congestion constraint modeling a roof, i.e. an impermeable boundary above the flow, is derived from the Navier-Stokes equations using classical thin layer arguments.
Let us consider a flow contained between two given surfaces respectively called bottom and roof. We assume that the two surfaces can be parametrized by two mono-valued smooth enough functions B (x, t) (for the bottom) and R (x, t) (for the roof) which satisfy B (x, t) ≤ R (x, t). We consider the domain where Ω x ⊂ R d−1 with d = 2 or d = 3, see Figure 1. The space variable is split into its horizontal component x ∈ Ω x and its vertical component z ∈ [B (x, t) , R (x, t)].

Navier-Stokes model with roof
A Navier-Stokes model describing the dynamics of the fluid is proposed. The bottom (resp. the roof) moves according to its own prescribed velocity with U B (x, t) ∈ R d−1 (resp. U R ) the horizontal component and W B (x, t) ∈ R (resp. W R ) the vertical component of the given velocity. More precisely for Γ = B or Γ = R, it satisfies The fluid is assumed homogeneous, with a density set to 1 for readability. The part of the domain occupied by the fluid, denoted Ω f (t), is marked by the color function ψ (x, z, t) ∈ {0, 1} such that ψ (x, z, t) = 1 iff (x, z) ∈ Ω f (t). More precisely we assume that the fluid surface is regular enough and that there exists a mono-valued function η (x, t) parametrizing the free surface such that ψ writes The velocity of the fluid satisfies the incompressible Navier-Stokes equations with gravity, i.e. for any (x, z, t) ∈ Ω f (t) × R +    ∇ x · u + ∂ z w = 0, ∂ t u + (u · ∇ x ) u + w∂ z u + ∇ x p − ∇ x · σ xx − ∂ z σ xz = 0, ∂ t w + u · ∇ x w + w∂ z w + ∂ z (p + gz) − ∇ x · σ zx − ∂ z σ zz = 0, (2.2) with u (x, z, t) ∈ R d−1 the horizontal component and w (x, z, t) ∈ R the vertical component of the velocity, p (x, z, t) ∈ R the pressure, g the gravity acceleration and σ (u, w) the viscosity tensor. We denote by n η the outward normal to the fluid surface defined by The fluid is advected by the flow and the color function describing the part of the domain occupied by the fluid satisfies the following relation From now on, a Newtonian fluid is considered, i.e. the viscosity tensor is given by where D x (u) = 1 2 ∇ x u + (∇ x u) t is the symmetric gradient and µ > 0 the viscosity coefficient. For each surface Γ = B or Γ = R the outward unit normal, see Figure 1, denoted by n Γ , is defined by At the bottom and at the roof, a no-penetration condition is considered, i.e.
The inequality at the roof is required to allow the fluid to detach itself from the surface. In addition, the friction between the fluid and the surface Γ = B or Γ = R is modeled by the following Navier condition with κ Γ > 0 a friction coefficient at the surface Γ and Π Γ the projection on the tangential plane to the surface. At the free surface, the continuity of the stresses reads with P (x, z, t) a given atmospheric pressure. At this stage of modeling, when the pressure in the fluid reaches the atmospheric pressure, i.e. p = P , the fluid can either create vacuum, i.e. ψ = 0, or the pressure can decrease below the atmospheric pressure p < P , see [10] for more details. In [10] the authors propose a numerical strategy able to forbid the creation of vacuum but it is not clearly described at the continuous level. In the current work, only non-negative relative pressure is considered, i.e. p |z=η ≥ P |z=η . (2.7) The drawback of this assumption is that areas without water trapped under the roof can appear, which is not physically relevant. Some perspectives to overcome this drawback are discussed in the conclusion section.
The system (2.1)-(2.7) has to be completed with compatible initial data. More precisely the fluid domain at the initial time reads Ω f (0) = Ω 0 f with Ω 0 f given and η (x, z, 0) = η 0 (x). Similarly, for any (x, z) ∈ Ω 0 f , we assume u (x, z, 0) = u 0 (x, z) and w (x, z, 0) = w 0 (x, z) where u 0 and w 0 are given. Note that to be physically relevant, the initial fluid has to be contained in the domain, i.e. Ω 0 f ⊂ Ω (0) and the velocity should satisfy the incompressibility condition, i.e. first equation of (2.2).
The following domain invariance holds.
Proof. By integration of (2.3) in R d \ Ω (t) and using a Leibniz integral rule, we get Using the non-penetration inequality (2.4) and integrating in time, it yields The assumption on the initial condition allows to conclude since ψ (x, This property ensures that the fluid stays confined between the bottom and the roof.

Congested shallow water model
Let us introduce the congested shallow water model we are interested in. The unknowns of the model are the water depth h, the vertical-averaged horizontal velocity u and the surface pressure p η which satisfy with the following given parameters: κ the efficient friction parameter, U the efficient velocity of forcing and µ the viscosity coefficient. In addition, the water depth has to satisfy the congestion constraint with the opening between the roof and the bottom H (see Fig. 1) defined by Eventually, the surface pressure has to satisfy 1 − 1 h≥H (p η − P η ) = 0 and p η ≥ P η , (2.10) with P η (x, t) the given atmospheric pressure. The shallow water model (2.8)-(2.10) is completed with the initial conditions The congested shallow water model (2.8)-(2.10) can be derived from the Navier-Stokes model (2.1)-(2.7) following the strategy proposed in [25]. More precisely, considering the two characteristic space dimensions along the vertical and the horizontal dimensions, respectively H and L, and the characteristic time of observation T , the dimensionless variables readû = T L u, For readability, the hatˆis dropped even if from now on the variables are dimensionless. We claim the following proposition.
Proposition 2.2. For ε = H L small enough, assume that the parameters satisfy the following scaling Then (2.8)-(2.10) with the initial conditions is formally an approximation in O ε 2 of (2.1)-(2.7). More precisely Proof. As it was done in [25], the vertical integration of (2.2) yields the mass conservation (first equation of (2.8)) and the following momentum balance (2.12) Thanks to Proposition 2.1, ψ |z=R (x, t) = 1 h≥H (x, t). Let us start by formally identifying a first order approximation. As in [25], the main term of the third equation of (2.2) reads p = p |z=η + g (η − z) + O (ε). The main term of the second equation of (2.2) and the Navier boundary conditions (2.5) lead to u = u + O (ε). Injected in (2.12), the momentum balance of the congested shallow water model reads at first order In order to get a better approximation, the first perturbation terms are taken into account. The vertical momentum equation of (2.2) leads to where the surface pressure reads p η = p |z=η − 2εµ 0 ∂ z w |z=η . Considering the horizontal momentum equation of (2.2) and the Navier conditions at the bottom and the roof (2.5), the horizontal velocity of (2.2) satisfies By taking into account the definition of u, one can check that the horizontal velocity of (2.2) is under the form Inserted in (2.12), the shallow water model (2.8) follows, provided we define U and κ by (2.11). It remains to establish (2.9) and (2.10). A direct consequence of Proposition 2.1 leads to B ≤ η ≤ R and thus to the congestion constraint (2.9). Similarly, the non-negativity of the relative pressure (2.7) implies p η ≥ P |z=η . Introducing the first equation of (2.6) into the second equation of (2.6) gives We get the expected constraint (2.10).
The two cases in the definition of κ B and κ R are chosen to get positive parameters whatever the water depth. More precisely, both correspond to the same approximation at order O ε 2 . It follows that the effective friction parameter κ is always positive.

Preliminary analysis of the congested shallow water model
Let us mention some mathematical properties of the congested shallow water model. The objective is to highlight the physical relevance of the model (2.8)-(2.10). For a deeper analysis, we refer to [33]. In particular, the dissipation of mechanical energy is a crucial point for many applications such as renewable energy production using buoys. From the mathematical point of view, the energy balance, acting as a mathematical entropy, is an argument for the existence of long time solutions. In addition the conservation of the steady state, in particular the lake at rest, is relevant for many applications because the flow is close to this state. To characterize these two properties, the potential of the conservative forces defined by is a useful quantity.
with the mechanical energy E = K + P, the kinetic energy and the potential energy respectively given by K = 1 2 h u 2 and P = gh B + h 2 + hP η and the energy flux Proof. The kinetic energy balance law is obtained by multiplying the momentum balance of (2.8) by the mean velocity u which gives with φ defined by (2.13). Furthermore, the potential energy balance law is obtained by multiplying the mass conservation of (2.8) by the potential φ. It yields Summing the two energy balance laws, it remains to estimate the source term Eventually, considering the constraint (2.10), we get The congested shallow water model (2.8)-(2.10) admits the following lake at rest solution, i.e. a steady solution with a vanishing velocity. To provide an expression of it, we use the positive part function defined by (ψ) + = |ψ| + ψ 2 and we consider a uniform potential of conservative forces φ 0 .
Then for any φ 0 ∈ R, the lake at rest defined by is a steady solution of (2.8)-(2.10) and the surface pressure reads Proof. The lake at rest solution defined by (2.14) satisfies In addition, since u = 0, we conclude that ∂ t (hu) = 0 and ∂ t h = 0. Eventually, it is clear that the lake at rest (2.14) ensures the congestion constraint (2.9) and (2.10).
The lake at rest is not the only steady solution of (2.8)-(2.10). However, it is the simplest one. Some works in the literature consider more complex steady states [14,29,38,39].

Relaxation of the congestion constraint
Due to the congestion constraint (2.9), the congested shallow water model is difficult to handle analytically as well as numerically. For this reason, we propose new formulations of the congestion constraint. In the context of incompressible fluid dynamics, several strategies were proposed in the literature to deal with a divergence free constraint. These strategies require that the constraint is written under the form of a divergence. Thus let us first write the unilateral constraint of the congested shallow water model under the form of a divergence. In such a context, the surface pressure p η acts as a Lagrange multiplier in the congested areas.
Then the following constraints are equivalent for any solution of (2.8) and (2.10).
• i)⇒ii): Assume i). Multiplying the mass conservation equation of (2.8) by the indicator function for h ≥ H, the constraint ii) is recovered since in the congested area we have h = H. • ii)⇒iii): Assume ii). Using the mass conservation, the constraint becomes Remark 2.6. Once the congestion constraint is under the divergence form ii) of Proposition 2.5, one can remark that in the congested area, i.e. h = H, the surface pressure satisfies In order to recover the surface pressure from (2.15), boundary conditions are required. In [33], the continuity of the surface pressure is assumed. However, this strategy assumes a smooth enough flow, which is not always satisfied since hyperbolic models can lead to discontinuous solutions, see Section 3.2.3. In addition, this strategy requires the description of the interface position which is a challenging issue since the interface is strongly linked to the flow dynamics.
The relaxation method is a strategy frequently used to approximate the divergence free constraint, see [17,46]. We introduce λ > 0 the relaxation parameter, characterizing the amplitude of the penetration in the roof. Let us write the approximate model completed with the initial conditions The unknowns of (2.16) are h λ the approximated water depth which does not exactly satisfy (2.9) and u λ the approximated vertical-averaged horizontal velocity. The unknown surface pressure p η in (2.8) is replaced by p λ which is taken as an implicit function of h λ . More precisely, we assume We assume that the relaxation operator L : Λ → L (Λ) is invertible. However since choosing, for example, L = −∆ x needs to specify a boundary value problem, we will not discuss the domain of definition and the domain of invertibility of the operator L in the general case and we start by discussing the choice of L. Classically the regularization operator L is chosen among [45] where these operators are presented in the context of Navier-Stokes equations. However, this list is not exhaustive. Since the congestion constraint can be written under the form of a time derivative, see iii) of Proposition 2.5, it seems interesting to choose L under the form with L an operator chosen among I d or −∆ x . Then L is invertible and L −1 is a strictly positive map. We will however skip the theoretical details and do not precise the spaces on which this property holds when L = I d .
Note that with this choice, an additional initial condition on the surface pressure p λ (x, 0) = p 0 λ (x) is required and is not provided by (2.8). However, the equation (2.17) can be integrated in time. More precisely choosing the initial condition on p λ such that This choice corresponds to a pseudo-compressible approximation, see [45]. In the context of incompressible fluid dynamics, the convergence of the approximate model to the constrained model when λ goes to 0 using the classical regularization operators cited above, can be mathematically proven and the error introduced by the relaxation can be estimated, see for example [44].
The regularization operator L should be chosen in order to preserve the lake at rest solution, see Proposition 2.4 as well as the positivity of the relative pressure (2.10). In addition, the regularization operator should ensure the dissipation of the mechanical energy given in Proposition 2.3. In the following the latter properties are proven in the case of pseudo-compressible approximation (2.18). To characterize these properties, the potential of the conservative forces which now writes (see (2.13)) is again introduced. Let us start by a result on the energy balance. We skip the proof which is similar to the proof of Proposition 2.3.
Proposition 2.7. For smooth enough solutions of (2.16) and (2.18), the following dissipation law holds with the mechanical energy E λ = K λ + P λ , the kinetic energy and the potential energy respectively given by the potential energy relative to the regularization R (ψ) such that and the flux of energy Let us now prove the positivity of the relative pressure (2.10).
Proposition 2.8. The solution of (2.16) and (2.18) satisfies Proof. The positivity of the relative pressure is easy to obtain from the expression (2.18) of the surface pressure. Now multiplying (2.18) by the indicator function of the areas where the water does not reach the roof yields We conclude since L −1 is a strictly positive map.
Eventually, the conservation of the lake at rest is highlighted.
Then for any φ 0 ∈ R, the lake at rest defined by u λ = 0 and is a steady solution of (2.16) and (2.18) and Proof. Let us first consider the case where φ 0 > g H + B + P η . Starting from h λ + λ 2 L (h λ ) using (2.19), we get Using (2.18) it leads to and composing by L −1 we get p λ = φ 0 − g (B + h λ ) then φ λ = φ 0 .
At the limit λ → 0, one can check that the state (2.19) tends to the state (2.14). It is a (very weak) argument to claim that (2.16) and (2.18) is a consistent approximation of (2.8)-(2.10). It is consistent (as λ → 0) at least for some steady solutions.
To conclude, the pseudo-compressibility method seems a suitable approximation of the congested shallow water model (2.8)-(2.10). In particular, the simplest pseudo-compressibility method leading to an approximation satisfying the properties of the congested shallow water model (2.8)-(2.10) is given by L = I d . In addition, in this case, (2.18) gives an explicit function for the surface pressure Substituting it in (2.16) leads to One can note the similarity between (2.22) and the classical Euler model at the low-Froude regime [27]. The convergence to the "incompressible" limit (λ → 0) is studied in [31,32] in the context of compressible fluid flow.
Besides, the approximate model (2.22) is close to the Partially Free Surface model (PFS) proposed in [10]. Both models are able to switch between the congested areas and the free surface domain without following the interface and the surface pressure is an explicit function of the water depth. In [10], the authors link the parameter λ to a physical parameter such as the speed of sound in the water whereas it is a mathematical artifice here. As for the PFS model, one can show that the approximate model (2.22) is hyperbolic.
More precisely, denoting u N = u λ · N with a unit vector N ∈ R d−1 , the eigenvalues in the direction N are Proof. The proof is similar to the classical shallow water case. The approximate model (2.22) can be writen as The eigenvalues in the direction N being the eigenvalues of the matrix J x N x + J y N y , we conclude by simple computations.

Numerical resolution
The current section is devoted to the numerical approximation of the congested shallow water model (2.8)-(2.10). Following the strategy proposed in Section 2.2, the numerical approach is based on a discretization of the hyperbolic model (2.22) at the regime λ 1. We recall that the model (2.22) came from a particular choice of L, which leads, in addition of the mandatory characteristics of the approximate model , i.e. an energy estimation, the positiveness of the relative pressure and the steady state preservation, to an explicit formula for the pressure (2.21). Another choice of L would be less efficient numerically since Newton's fixed point explained hereafter would become less straightforward. From now on, the dissipative terms are neglected, i.e. κ = µ = 0. More precisely, the numerical strategy is based on an operator splitting between the hyperbolic part of the equation and the dissipative terms. We only detail the discretization of the hyperbolic part of the model since the dissipative step might be solved using an usual implicit scheme.
The numerical strategy used to solve the approximate model (2.22) has to be consistent with the "incompressible" regime λ 1. More precisely, in the congested areas, this corresponds to high potential regime since the second term in (2.21) becomes large. By analogy with gas dynamics, this corresponds to the low-Mach regime, see [31,32]. Classical Godunov solvers are known to be too much dissipative in this case and to require a very restrictive CFL condition, see [19,27]. Several schemes in the literature are designed to recover the "incompressible" regime, see [15,28,40].
In addition, a particular attention is given to the preservation of the lake at rest solution. Even more important than in the case of the classical shallow water model, the conservation of the lake at rest at the numerical level, using so-called well-balanced schemes [2,4,26], is a required property for the simulation of the model (2.22).
More precisely, if a scheme is not well-balanced and does not preserve the lake at rest, it will create oscillations, of order of the space step, at the free surface. However, the scheme is still consistent with the classical shallow water model. In the case of the congested shallow water model, small oscillations can artificially separate the surface from the roof and lead to strong variations of the surface pressure.
Most of the schemes designed to be consistent in the "incompressible" regime do not care about steady states at rest. More precisely, this regime is usually analyzed for the low-Mach regime of the Euler equations without source term or with a source term constant in space such as gravity. However, the Centered-Potential Regularization (CPR) scheme introduced in [40] satisfies the two requirements: consistency with the "incompressible" regime and the well-balanced property. The numerical strategy presented in the following is an adaptation of the CPR scheme to take into account the congestion constraint.

The CPR scheme adapted to the congested shallow water model
Let T be a mesh of Ω x composed of star-shaped control volumes, see Figure 2. We denote k ∈ T a control volume, F k the set of its faces, |k| its volume and |∂k| its surface area. Furthermore for a face f , its area is denoted |f | and k f is the neighbor control volume of k such that k ∩ k f = f . The space step is defined by l k = |k| |∂k| and the normalized face size by ν k f = |f | |∂k| . The unit normal to f outward to k is denoted as N k f . The set of the faces at the boundary of the mesh is designed by F ∂T . In addition, the time step is denoted by δ n+1 t , i.e. t n+1 = t n + δ n+1 t .

Description of the CPR scheme
The CPR scheme is an IMplicit-EXplicit Finite Volume method. The numerical unknowns are approximations of the average in a control volume k of the variables. The mean atmospheric pressure, the mean bottom level, the mean roof level and the mean opening in the control volume k and at time t n are respectively denoted by P n k , B n k , R n k and H n k = R n k − B n k . One particularity of the CPR scheme is that the numerical unknowns are not the conserved variables, they are chosen as the potential φ n k and the horizontal velocity u n k . This choice of variables is particularly well adapted to the congested shallow water model since the potential can be used as a parametrization of the water depth h n k and the surface pressure p n k . More precisely, we recall that the potential reads φ λ = g (h λ + B) + p λ and we set For readability, we set h n k = h (k, n, φ n k ) and p n k = p (k, n, φ n k ). Using an implicit finite volume method, the mass conservation, first equation of (2.22), reads with F n+1 f an approximation of the mean mass flux h λ u λ at the face f between the times t n and t n+1 . More precisely if f ∈ F ∂T , we set with a regularization parameter γ ≥ 0 characteristic of the CPR scheme, called CPR-parameter in the sequel, see [40]. We have also introduced the following notation at the face If f ∈ F ∂T , the numerical scheme requires boundary conditions which depend on the regime of the flow. In the current work, we do not detail the boundary conditions treatment (for a hyperbolic system it is known to be a difficult point, see [22] for more details). The scheme (3.3) is an implicit and non-linear scheme on the potential, thus it cannot be solved directly. An iterative Newton process is used to solve it at each time step. More precisely, writing the scheme as S k φ n+1 = 0 with the Newton method consists in constructing a sequence (φ n,q ) q≥0 defined by φ n,0 = φ n and where J denotes the Jacobian matrix of S (φ) = (S k (φ)) k∈T . Furthermore φ n+1 is defined as lim q→∞ φ n,q . A similar Newton method is analyzed in [11]. Numerically, a quadratic convergence of the Newton method is observed.
Once the potential φ n+1 k is recovered, the water depth h n+1 k is known and the momentum balance is computed using the following explicit upwind scheme in order to get u n+1 Eventually, the discrete surface pressure can be estimated a posteriori for output using (3.2).
Remark that thanks to the definition of the variables (3.1)-(3.2), in particular the choice of φ instead of h, the scheme (3.3) and (3.4) can be used in practice with λ set to 0. In other words, there is no longer penetration in the roof and the congestion constraint (2.9) is exactly satisfied. In the next section, we first propose an analysis with λ > 0 then we investigate the case λ = 0 numerically in the simulations.

Numerical analysis
In the following section, the numerical counterpart of the main properties of the congested shallow water model is established. In particular, the dissipation of the mechanical energy in Proposition 2.7, acts as a mathematical entropy and proves the numerical stability of the scheme. Firstly, let us recall the main numerical properties inherited from the CPR scheme, before showing the properties specific to the congested shallow water model.
Let us start by a consistency result valid on a regular cartesian grid. We denote by δ x the uniform space step, i.e. the distance between the center of two neighboring cells. In this particular case the coordinates of the center of each cell are denoted by X k .
Proof. The proof is done in Proposition 2.2 of [40]. i) (Positivity of water depth) Assume h 0 k > 0 for any k ∈ T and assume that the following implicit CFL condition is satisfied for any n ∈ N Then h n k > 0 for any k ∈ T and n ∈ N. ii) (Conservation of water volume) For any n ∈ N, iii) (Well-balanced scheme) For any φ 0 ∈ R, let the initial condition be One can check that the discrete lake at rest iii) defined in Proposition 3.2 is a discrete version of the lake at rest of Proposition 2.9.
Note that the CFL condition (3.5) cannot be satisfied at a wet-dry interface. In practice the simulation of the latter can be realised by taking artificially some small time step but the scheme will not satisfy the dissipation law given below in Proposition 3.3.
Let us now show some properties in the case of the congested shallow water model.
with the discrete mechanical energy E n k = K n k + P n k , the kinetic energy and the potential energy respectively given by the flux of the potential energy Proof. The following proof is an adaptation of Theorem 2.3 from [40]. The estimation of the discrete kinetic energy dissipation law under the CFL condition (3.5) established in Lemma 2.4 of [40] still holds in our case since it is obtained for any potential. However, the discrete potential energy dissipation law ( [40], Lem. 2.5) is obtained for a differentiable potential, which is not our case due to the positive part function of the regularization. Let us establish the discrete potential energy balance in our case. Using the two equalities b and taking into account that Q n+1 k ≥ 0 gives the discrete potential energy dissipation law with the rest coming from the space discretization and the work of the conservative forces respectively given by Finally, the proof of Theorem 2.3 from [40] allows to conclude.
Note that S n+1 k δ n+1 t is consistent with the terms depending on a time derivative in the right-hand side of the energy balance law given in Proposition 2.7. ii) (Upper bound on the surface pressure) Assume that for any f ∈ F ∂T , we have F n f = 0. Then there exists K ∈ R * + such that for any λ > 0, we have p n k ≤ K for any k ∈ T and n ∈ N. Proof. The proof of i) is a consequence of Proposition 2.8.
Let us now prove ii). The mass conservation (3.3) at time t n+1 reads −γδ n+1 Let us consider the latter equality for all k ∈ T and denote φ n+1 the vector with the components φ n+1 k for k ∈ T. We get a system of the form M n+1 φ n+1 = Φ with M n+1 a non-singular M-matrix. Indeed the entries of the main diagonal of M n+1 are non-negative, the off-diagonal entries are non-positive and the matrix is strictly diagonally dominant. The system leads to an inequality of the form In fact, thanks to the water volume conservation ii) of Proposition 3.2, the water depth is bounded by 0 ≤ h n k ≤ h ∞ with h ∞ = k∈T |k|h 0 k min k∈T |k| and the bound on the right hand side follows.
As a consequence, we get a bound on φ n k and the bound on the surface pressure follows. More precisely, ∀n ≥ 0, max k∈T p n k ≤ max k∈T φ n k + g (h ∞ + max k∈T B n k ).
Thanks to ii) from Proposition 3.4, we are able to characterize the penetration with respect to λ. More precisely considering (2.21), we have It shows that (2.9) is verified as λ goes to zero. Furthermore, it shows that the mechanical energy in Proposition 3.3 remains bounded when λ tends to 0.

Simulation and numerical validation
The numerical solution of the congested shallow water model is illustrated in a one dimensional framework, i.e. x ∈ R. Unless indicated otherwise, for all the test cases g = 9.81, P η = 0, γ = 1 and a channel of length l = 1 is considered, i.e. x ∈ [0, 1] where the entry is located at x = 0 and the exit at x = 1.

The lake at rest
The discrete lake at rest iii) of Proposition 3.2 is tested and validated in various configurations with bottom and roof which may be continuous or not. Note that in such configurations, the CFL condition (3.5) is satisfied for any positive time step. Setting λ to zero, the discrete lake at rest iii) of Proposition 3.2 satisfies the constraint (2.9) and thus is a discrete version of Proposition (2.14). This setting is tested and validated as well. The water depth at the exit is fixed to h (l, t) = 8 and the discharge at the entry is set to hu (0, t) = 40. The physical parameters are chosen such that the flow is subcritical (fluvial) at the boundaries but becomes supercritical (torrential) because of the roof. The numerical solutions with the parameter λ = 0 and λ 2 = δ x together with the analytical solution with regular space step δ x = 10 −3 are plotted in Figure 3.

Transcritical steady flow with free hydraulic jump
First of all, let us discuss the results with respect to the value of λ. The simulation was performed with several parameters λ from 1 to 0 and for δ x ∈ {1/300, 1/1000, 1/3000, 1/10 000}. In Figure 4 the L 2 -errors of the water height, the velocity and the surface pressure are plotted versus λ 2 for the different space steps. We notice that for λ < 1, the water depth and the velocity are well approximated. The surface pressure is more difficult to obtain. The best agreement with the analytical solution is obtained for δ x ≤ λ 2 ≤ (δ x ) 1 2 . For λ 2 < δ x , the surface pressure is not in good agreement with the analytical solution. Indeed the numerical error of the scheme on the congestion constraint is fully reflected to the surface pressure when λ tends to 0 whereas it is split between the water level and the surface pressure in the case λ > 0. Figure 5 shows the convergence rates in L 2 -norm to the analytical solution for the water depth, the velocity as well as the surface pressure for λ = 0 and λ 2 = δ x . In both cases, λ = 0 and λ 2 = δ x , the water depth and the velocity converge at first order. In the case λ 2 = δ x , the surface pressure converges at approximatively order 1.5, which is faster than expected. Indeed, as the water height converges at first order, the surface pressure is expected to converge at first order as well. In the case λ = 0, the surface pressure seems to converge but chaotically.

Transcritical steady flow with constrained hydraulic jump
The solution presented in Appendix A is not valid when the top of the hydraulic jump reaches the roof. In this case, the non-conservative product h∇ x p η is not clearly defined. The following simulation illustrates the solution of the numerical scheme in such a configuration.   The water depth at the exit is fixed to h (l, t) = 8 and the discharge at the entry is set to hu (0, t) = 25. The physical parameters are chosen such that the flow is subcritical (fluvial) at the boundaries but becomes supercritical (torrential) because of the nonflat bottom and the hydraulic jump reaches the roof. The simulated solution with δ x = 10 −3 and λ 2 = δ x as well as λ = 0 is plotted in Figure 6. Furthermore a reference solution is plotted for comparison. More precisely the analytical solution of the water height without considering the roof and the corresponding velocity are plotted, see (A.2) as well as the analytical value of the surface pressure considering the roof, see (A.3). As expected, the water depth and the velocity are similar in the areas where the flow does not reach the roof and only the localization of the hydraulic jump is not recovered. At the hydraulic jump, the water depth as well as the surface pressure are discontinuous which makes the product h∇ x p η really non-conservative and the position of the shock not easy to compute. The peak of the pressure at the top of the hydraulic jump is due to oscillation effects of the CPR scheme at discontinuities, see [40]. By increasing the CPR-parameter γ in order to regularize the solution and decrease the oscillation effects, the approximation of the water height gets worse. However, the hydraulic jump is clearly not located at the same position with and without the roof. We conclude that the surface pressure has an impact on the jump condition, at least at the discrete level.

Time-dependent emptying tank
Assuming a fixed and flat bottom, an analytical solution of the congested shallow water model can be obtained for any fixed bottom as long as the interface between the congested and the free surface flow is fixed. More precisely, we claim the following statement.  x r , i.e.
with an arbitrary opening H 0 > 0. Then for any α ≥ 0, h 0 > H 0 and any time t ≤ h 0 −H0 αH0 , the following functions  Proof. The proof relies on simple computations after injecting the functions in the congested shallow water equations (2.8)-(2.10). The non-conservative product is well-defined since the pressure at the bottom gh + p η is continuous. Figure 7 shows the analytical solution of the Proposition 3.5 at time t = 0.5 with the parameters h 0 = 10, H 0 = 5.8, x l = 0.3 and x r = 0.7, together with the numerical solution with δ x = 10 −3 and λ 2 = δ x as well as λ = 0. The numerical solution is in very good agreement with the analytical solution. Figure 8 shows the convergence rates in L 2 -norm to the analytical solution for the water depth, the velocity and the surface pressure for λ = 0 and λ 2 = δ x . The three unknowns converge approximately at order 1.36 for any λ, which is better than the expected first order. In this case the pressure depends linearly on the water depth, which explains the same order of convergence for the pressure for both values of λ compared to the stationary case.

Time-dependent filling tank
The current simulation is a numerical investigation of the jump condition when the top of the hydraulic jump reaches the roof. More precisely, a fixed flat bottom and roof and a symmetrical initial condition given by with H (x, 0) = 5 for all x ∈ [0, l] and C = 8 are considered. The initial condition is chosen such that the flow is supercritical (torrential). Figure 9 shows the numerical result at time t = 0.5 with δ x = 10 −4 , λ 2 = δ x and γ = 10. Due to the shock waves, the CPR-parameter γ is increased to get slightly more regularisation, see [40]. Otherwise oscillations of the pressure in the vicinity of the shock due to oscillation effects of the CPR scheme are observed. However, the oscillations increase when λ becomes smaller and more regularization is needed. As a consequence, we cannot pass λ = 0 in the simulation, as either the oscillations become too large or the regularization effect is too high. As expected, the solution is composed by tree constant states separated by two shock waves. By symmetry, it is clear that the velocity in the middle state should vanish. In addition, the conservation of the mass imposes the velocity of the shock wave. More precisely, it leads to X = 2U with X (t) the position of the shock. However, the computation of the pressure of the middle state is inaccessible since the water depth and the surface pressure are discontinuous at the shock.

Conclusion
A shallow water type model with a constraint of congestion is derived from the Navier-Stokes equations for the modeling of a roof in geophysical applications. The simulation of the congested shallow water model is challenging due to the unknown localization of the interface between the congested areas and the free surface domain. To overcome this difficulty, an approximate model which is hyperbolic in the free surface domain as well as in the congested areas is proposed based on a pseudo-compressible approximation. We show that two properties, mainly well-balanced and low-Froude stability, are needed at the discrete level for the simulation of the approximate model. Eventually, a numerical scheme adapted from the Centered-Potential Regularization scheme is proposed and confronted with analytical solutions to illustrate the efficiency of the scheme.
Concerning the theoretical analysis, several issues have not yet been adressed. Firstly, to confirm the relevance of the pseudo-compressible approximation, the convergence of the approximate model (2.16) and (2.18) as λ goes to zero to the initial congested shallow water model (2.8)-(2.10) seems a natural question but not at hand for the moment. The second open problem in the analysis of the congested shallow water model (2.8)-(2.10) is the definition of the non-conservative product h∇ x p η . More precisely the water height as well as the surface pressure can be discontinuous in the same point as shown in Figure 6. Since the surface pressure plays the role of a Lagrange multiplier, the model (2.8)-(2.10) is not a hyperbolic conservation law. As a consequence a classical analysis of non-conservative products in hyperbolic theory, for example [18], cannot be directly applied here.
To improve the modeling, several ways can be investigated. First of all, the modeling of areas with negative relative pressure, i.e. p η < P η is quite challenging. This issue seems more generally linked to the modeling of trapped air pockets under the roof. Secondly, the congestion constraint can be applied to more complex models of free surface flow such as the layerwise discretized models [1,3], the non-hydrostatic models [12,34] or the combination of the two [23]. Moreover, the bilateral interaction between the roof and the fluid seems interesting for applications, in particular for object buoyancy such as boats or buoys. Last but not least, the modeling of an immersed object, using layerwise discretized models, is still quite far but can be a challenging objective.
surface and the surface pressure is assumed continuous we have p j = P j . To determine this point, we start from the right point where the opening intersects the critical water depth. At the left of this point, the flow is clearly supercritical. Then, we proceed successively node by node going to the left as follows. Assuming the current node is the characteristic node, the hydraulic head is computed using (A.1). Then, the pressure at the left neighbor node is computed using (A.3). If the relative pressure at the left is negative, we conclude that the current node is not the characteristic node and we pass to the next one. Once the characteristic node is determined, the water depth at the right is computed using (A.2) until the hydraulic jump. At the left of the characteristic node, the surface pressure is computed using (A.3) until the pressure becomes negative, which corresponds to the entrance of the water under the roof. Eventually, the velocity can be deduced dividing in each point the mass flux by the water depth. The numerical implementation of the analytical solution is given in Algorithm 1.