On the relevance of the dam break problem in the context of nonlinear shallow water equations

The classical dam break problem has become the de facto standard in validating the Nonlinear Shallow Water Equations (NSWE) solvers. Moreover, the NSWE are widely used for flooding simulations. While applied mathematics community is essentially focused on developing new numerical schemes, we tried to examine the validity of the mathematical model under consideration. The main purpose of this study is to check the pertinence of the NSWE for flooding processes. From the mathematical point of view, the answer is not obvious since all derivation procedures assumes the total water depth positivity. We performed a comparison between the two-fluid Navier-Stokes simulations and the NSWE solved analytically and numerically. Several conclusions are drawn out and perspectives for future research are outlined.


1.
Introduction. During the last century there were more than 200 failures of dams greater than 15 m high [99,118]. They have caused a loss of more than 8000 lives and millions of dollars worth of damage. Consequently, dam break flows have become an important practical problem in civil engineering. Numerical models have become essential as a predictive tool in evaluating the risks associated with the failure of the hydraulic structures. That is why, the number of numerical studies has drastically increased during past decades.
To our knowledge, the dam break problem was studied analytically for the first time in the PhD thesis of Pohle (1950), [89], who used a lagrangian description to solve this problem. The classical analytical solution for the dam break problem in the context of the NSWE can be found in the book of Stoker (1957), [102]. Later, this solution was generalized to the constant slope case by Mangeney et al. (2000), [79]. Note, that Hunt (1982), [58], also considered the sloping channel case and he obtained a closed-form solution using a kinematic wave approximation. Among classical works on this topic, we have to mention the prominent paper by Benjamin (1968), [12]. Recently, Korobkin & Yilmaz (2008), [68], studied the initial stages of the dam break flow in the framework of potential free surface flows.
It is interesting, however, to recall some other known analytical solutions to NSWE even if they are not directly related to the dam break problem. Wave runup on a sloping beach was investigated by Carrier & Greenspan (1958), [22,26], using a hodograph transformation 1 . This solution is extensively used in the tsunami waves community to validate the run-up algorithm of various NSWE solvers [60,105,109,95,94,62,35]. The transform of Carrier & Greenspan was employed later by Synolakis and his collaborators to study analytically tsunami run up on a sloping beach, cf. e.g. [103,107,108,69]. There is also an analytical solution by Liu et al. (2003), [72] of the linearized shallow water equations on a sloping beach, where they used a forcing term to model an underwater landslide. This solution is also currently used to test numerical codes, [35].
On the other hand, the dam break problem and various lock exchange flows were extensively studied experimentally, cf. e.g. [81,111,65,55,96,10,11,97]. In this study, we do not directly appeal to them, since our main concern is to study the validity of NSWE as an approximation to more complex mathematical models in some extreme situations.
The authors decided to perform this study because there is an apparent contradiction between the mathematical origins of the NSWE and some applications of this model. When we look carefully at any derivation procedure of NSWE, we will see that an implicit assumption of water depth positivity is adopted. Moreover, these equations are designed to model infinitely long waves. That is why, strictly speaking, these equations can be valid only in fluid regions. However, using various numerical techniques (sometimes ad-hoc, semiempirical) this model is routinely used for wetting/drying (run-up/run-down) simulations, cf. e.g. [105,109]. This process is considerably more complex and the validity of the NSWE is not obvious a priori. Recall, that the shoreline can be considered as a triple point: water, air and solid (soil) meet their. Of course, this situation is simplified for mathematical modelling.
We choose a Direct Numerical Simulation (DNS) by the two-fluid Navier-Stokes equations [104,90] as the reference solution. This system contains all the necessary physical effects ranging from viscosity to the surface tension. Moreover, the ambient fluid (air) is resolved. In the absence of experimental data, these simulations can be assimilated to an idealized experiment. Up to graphical resolution, our numerical results are very similar to the experiments of J. Martin and W. Moyce [81] and we remain clearly in the laminar régime. We also underline that we consider a realistic density ratio 1:1000 as for the air/water interface (see Table 1). The results of the DNS are compared with several solutions to the NSWE. Namely, the analytical solution of Stoker [102] (see Section 3.3) was used in our comparison. Numerical solutions to the NSWE were obtained using the VOLNA code, cf. [37,88,35].
The present study is organized as follows. In Section 2 we present two mathematical models which are used in this study. In the same section we also discuss several mathematical properties and extensions of the NSWE. In Section 3 we review some known analytical solutions to the NSWE of the dam break problem. After discussing briefly the numerical techniques, (Section 4), we present and discuss our numerical results in Section 5. Conclusions are outlined in Section 6.
2. Mathematical models. In this section we briefly present two mathematical models which are used in the sequel. The first model is the well-known Nonlinear Shallow Water Equations (NSWE) which were derived for the first time by Saint-Venant (1871), cf. [36]. The second model is the two-fluid Navier-Stokes equations written under the assumption of fluids immiscibility. These equations are much more complete from physical and mathematical points of view. That is why, the two-fluid model is supposed to provide us reliable results.
where H( x, t) is the total water depth and u( x, t) : R 2 × R + → R 2 is the depthaveraged horizontal velocity. Traditionally, g denotes the acceleration due to the gravity, h( x, t) is the bathymetry function and I is the identity tensor. We do not provide here the derivation of these equations since it is more than classical and can be found in various sources [102,78]. Remark 1. The bathymetry function h( x, t) can be time-dependent. It is especially important for tsunami generation problems by submarine earthquakes, landslides, etc. The coupling with seismology is usually done through this function. Namely, various earthquake models, cf. e.g. [28,30,29,66] give us the seabed displacements which are then transmitted to the ocean layer. Obviously, in this study we consider the fluid propagation over the flat bottom in view of applying analytical techniques.
Governing equations (1), (2) form the system of balance laws (conservation laws, if the bottom is even h = const). Moreover, this system is strictly hyperbolic provided that H > 0. This property is extensively used in the construction of various numerical schemes and, in particular, in the Characteristic Flux approach, cf. [43,41,42,44,35], which is also implemented in the code VOLNA .
Let us discuss the eigensystem of the advective flux. First, we introduce the so-called conservative variables and rewrite the governing equations as a system of conservation laws: where we introduced the following notations: w( x, t) : R 2 × R + → R 3 , w = (w 1 , w 2 , w 3 ) = (H, Hu, Hv),

DENYS DUTYKH AND DIMITRIOS MITSOTAKIS
After projecting the flux F ( w) on a unit normal direction n = (n x , n y ), | n| = 1, one can compute the Jacobian matrix A n . Its expression in the physical variables has the following form: where u n = u · n is the velocity vector projected on n. The Jacobian matrix A n has three distinct eigenvalues: where c = √ gH is the gravity wave speed in infinite wavelength limit. This quantity plays the same rôle as the speed of sound in compressible fluid mechanics. The hyperbolicity condition for the system (1), (2) follows immediately from (4) and the definition of c. The eigenstructure of the Jacobian matrix A n is fundamental for constructing numerical flux function, cf. [35], and thus, upwinding the discrete solution.
2.1.1. Properties. Nonlinear Shallow Water Equations have many other interesting properties. Some of them will be briefly recalled here. To reveal these properties, we shall take the water wave theory point of view.
Let us recast equations (1), (2) in the following nonconservative form in one space dimension: These equations possess a (non-canonical) Hamiltonian structure [93,117,91]: where the Hamiltonian H is defined as Moreover, the pair of equations (5), (6) possesses an infinity of conservation laws [13,80]. Equations (5), (6) can be also derived from Luke's Lagrangian variational principle [75] if we introduce the velocity potential function φ( x, t) such that u = ∇φ. In this case, the Lagrangian reads Governing equations (5), (6) are obtained by varying L with respect to η and φ.
Recently, a generalized variational principle was proposed by Clamond & Dutykh, cf. [21]. Their approach allows for more flexibility and can be used to derive various generalized shallow and deep water approximations. Figure 1. Two immiscible fuids separated by an interface. (1), (2) arise after a series of approximations applied to complete set of equations. Strictly speaking, they model the propagation and transformation of infinitely long water waves. That is why, their validity for run-up and flooding simulations is not so obvious a priori.
2.2. Two-fluid Navier-Stokes equations. Let us consider two immiscible and incompressible 2 fluids (water and air, for example) occupying domain Ω = Ω + ∪ Ω − , where they are separated by an interface S. This situation is schematically depicted in Figure 1. We note that we do not make any assumption on the interface complexity and topology. In what follows we will denote by superscripts ± all quantities related to the heavy and light fluids respectively.
In each fluid we can write mass and momentum balance equations: The latter may be written in conservative form: where u is the fluid velocity, ρ ± are the fluids densities, µ ± are the fluids dynamic viscosities, D = 1 2 (∂ i u j +∂ j u i ) is the rate of deformation tensor. The surface tension term is a force concentrated at the interface, σ is the surface tension coefficient, κ is the curvature of the interface, n is the unit normal to the interface and δ S is the distribution (Dirac mass function) concentrated on the interface S.
Governing equations (7), (8) have to be completed by the following jump conditions across the interface: • Velocity continuity • Tangential stress condition • Normal stress condition where t is a tangent vector ( t · n = 0) to the interface and notation [·] S represents the jump of a quantity across the surface S. However, for numerical computations it is advantageous to introduce a characteristic function φ, (cf. [61,106,44,33,32,31]) defined as: Thus, φ and n are related by the formula ∇φ = nδ S . In the absence of phase change, φ is simply advected by the fluid motion: In order to write a unique formulation for the entire domain, we express the density and the viscosity as functions of φ: Thus, we have the following momentum balance equation: Along with the mass conservation equation (7) and the volume fraction advection equation (12), it forms the two-fluid Navier-Stokes equations with an interface, which are solved numerically below.
Remark 2. We can recover jump conditions (9) -(11) if we investigate the governing equations (7), (12), (13) in the neighborhood of the surface S and making use of the formula ∇φ = nδ S .
3. Analytical solutions. In this section we review known analytical solutions related to the dam break problem that we use in the comparison with the numerical results.  (5), (6): In some situations, it is advantageous to eliminate the velocity variable u to obtain The Initial Value Problem (IVP) for (14) corresponding to the dam break takes the following form: where H(x) is the Heaviside function. This IVP can be easily solved using the Fourier transform: The sketch of this solution is presented in Figure 2. Namely, it consists of two waves propagating in opposite directions with velocities ±c 0 . Hence, the front speed is equal to −c 0 . Of course, this result is nonphysical as it will be shown below. Similar solutions can be constructed considering the linearized Euler equations for either one or two fluids separated by an interface.

3.2.
Small time asymptotics. Several small time asymptotics were proposed to solve the dam break problem. One of the first solutions was derived by Pohle (1950), [89]. Such methods generally require the use of lagrangian description. The prominent book by Stoker, [102], also contains such a solution: where (X, Y ) are new coordinates of the particle (a, b) at time t. Recently this solution was generalized by Korobkin & Oguz (2008) [68]. We tried to compare the solution (15), (16) with our numerical results and found that its validity time is too short for any practical use. That is why this solution is not plotted bellow. Figure 3. Sketch of the initial condition for the shallow water computations.

DENYS DUTYKH AND DIMITRIOS MITSOTAKIS
Note, that expressions (15) and (16) Schematically it is depicted on Figure 3. Then, by considering the Riemann invariants and using the method of characteristics, [70,51,52,23], one can derive the following solution: where c 0 := √ gh 0 is the gravity wave speed in the undisturbed region. The front position is given by the characteristic outgoing from the fluid region: Recall that recently this solution was generalized to the constant slope case by Mangeney et al. (2000), [79]. 4. Numerical methods. The main purpose of this study is to draw out some conclusions on the validity of NSWE for wetting (flooding) process simulations. That is why, we do not provide here any details about numerical methods used to compute solutions. The interested reader can consult references given below to get technical details.
In order to solve numerically the two-fluid Navier-Stokes equations (7), (13) and (12), we applied the finite volumes method, cf. e.g. [63,92,87]. Namely, a freely available solver interDyMFoam of the OpenFOAM CFD Toolbox [87] was used. The interface between two fluids is reconstructed from the volume fraction φ distribution using the VOF method, cf. [57,104,90]. Let us underline that all twofluid computations presented in this study are 3D with only one cell in z-direction. Everywhere we impose the classical no-slip boundary condition.
Nonlinear Shallow Water Equations are solved with our operational numerical code VOLNA , cf. [37,35]. This code was developed in close collaboration with R. Poncet and F. Dias when the first author was at CMLA, ENS de Cachan. The VOLNA code uses unstructured triangular meshes and is able to run in arbitrary complex coastal regions. The numerical method is a second-order finite volumes MUSCL-TVD scheme along with the SSP-RK3(4) method for the discretization in time, cf. [101]. Details on adopted discretization procedure can be found in [37,88,35]. All the computations we performed are in 2D and only one-dimensional cross sections are presented below. On the lateral boundaries we impose the wall boundary condition u · n = 0. This choice is consistent with two-fluid computation and allows us to have an insight into the impact process.

5.
Comparison results and discussion. In this section we perform a comparison between a two-fluid simulation (DNS), the analytical solution by Stoker (1957) and numerical solutions to the NSWE by the VOLNA code. The initial set-up for the VOLNA code is shown in Figure 3. Sketch of the initial condition for the DNS is depicted on Figure 4. The simulation time and propagation distance is chosen so that the right boundary do not influence obtained results. All parameters used in computations are given in Table 1. These parameters are chosen suitable to simulate the air/water interaction.
The snapshots of our simulations are given on Figures 5 -11. On the left image (a) we represent the volume fraction φ distribution provided by the DNS. On the right image (b) we plot together the analytical solution (17) (red dotted line) and simulation results by the VOLNA code (black solid line) for the free surface elevation η. The analytical solution is almost superposed with our numerical simulation as it is expected. This result can be considered as one more validation test of the wetting/drying algorithm used in the VOLNA code.
In the beginning of the simulation, the water column is slightly deformed due to the gravity force ( Figure 5). Only a small time interval is needed for the heavy fluid to acquire the kinetic energy and to enter into the propagation régime depicted on  Table 1. Parameters used in numerical simulations.    (17) prescribes a parabolic form of the interface. However, the DNS shows somehow different shape. Lower fluid layers undergo stronger acceleration than in NSWE and thus propagate faster. Nonuniform distribution of the velocity field along the heavy fluid is illustrated on Figure 12. This creates a strong distortion of the interface which is elongated near the bottom (it can be easily seen in Figures 10 -11). This effect is not present in NSWE simulations since, the vertical flow structure is not resolved by this approximate model.   Consequently, in NSWE we obtain a piecewise linear distribution of the velocity field as it follows from analytical solution (18). Let us notice another one fact. In Figure 11 one can observe that the NSWE solution has already reached the left vertical wall. From this moment, the analytical solution is not valid anymore. However, the two-fluid simulation has not yet reached the left boundary. This discrepancy comes from the time lag due to initial acceleration stage, on one hand, and slightly different front propagation speeds, on    Table 2. Front speed predicted by four different approaches. stage (t ≤ 0.25 s) which is present in the two-fluid model. Another explanation consists in the front velocity which can be determined by measuring the slope to the curve t → x f (t). Determined in this way front speeds (in permanent régime, t ≥ 0.25 s) are given in Table 2. The LSWE give completely wrong results showing again that the nonlinearity plays the crucial rôle in this process. It is also worth to note that the numerical front speed by the VOLNA code is closer to the DNS. This positive fact can be attributed to the effect of the numerical diffusion on unstructured meshes. 5.1. Impact process. In the previous section we presented results concerning the initial and propagation stages of the dam break problem. However, we continued the computations until the interaction with the left wall and even slightly beyond. The goal is to test again the validity of the NSWE in such extreme conditions. For this kind of situations there is no analytical solution, and thus, we compare only DNS and VOLNA code results. In two-fluid simulation we use the same noslip condition for all solid boundaries. The implementation of the wall boundary condition of VOLNA solver can be found in [35], and it is based on considerations of incoming characteristics. The general methodology is presented in works of J.-M. Ghidaglia and F. Pascal [49,48,50].
The comparison results are presented in Figures 14 and 15. For instance, the wave amplitude on Figure 15     From these results it is obvious that the wave impact process is not correctly modelled by NSWE. It is possible to foresee this conclusion if one recalls two main constutive assumptions behind NSWE: • The pressure is hydrostatic • Vertical velocity and acceleration are neglected Notice, that for infinitely long waves, these two hypotheses are equivalent. Since, the dynamic pressure dominates in the impact process, we get qualitatively wrong results (it is especially clear from Figure 15). 6. Conclusions and perspectives. In the present study we tried to examine the validity of the NSWE for wetting (flooding) process modelling. As a test-case we chose the classical dam break problem which is the de facto standard in this field. This problem was solved in the context of two completely different models in terms of physical precision and, by consequence, of different complexity. The two-fluid DNS was chosen as the reference solution since all necessary physical effects are included in it.
Comparison results presented above show good overall performance of the NSWE. In order to appreciate more these results, one should take into account also the computational cost of the DNS and relatively inexpensive shallow water simulations.
However, we revealed several drawbacks of the depth-integrated model. Namely, the free surface shape differs from the parabolic profile predicted by NSWE. This discrepancy is attributed to the non-piecewise linear distribution of the velocity field inside the water column. To compare, see Figure 12 for the DNS result and formula (18) for the analytical prediction by NSWE. The experimental and theoretical study of P.K. Stansby et al (1998) [96] also revealed some differences during the initial stages. However, their objection concerned essentially some new jet-like phenomena just after release. For later times, they found relatively close agreement with NSWE.
We went beyond the initial purpose of this study and continued our simulations until the impact process with the left wall. It was shown that the NSWE strongly underestimate the wave height. This discrepancy has its origins in the hydrostatic pressure assumption. Actually, the dynamic pressure becomes dominant during the impact process. Its excess is responsible of spectacular splashes that we may have a chance to observe in nature.
Concerning the front propagation speed, we obtained slightly different values between the DNS and the NSWE solutions. We have to note that in a physical experiment this quantity strognly depends on the soil conditions. The standard no-slip boundary condition is clearly insufficient to describe all kinds of soils. We believe that future research activities will focus on developing wall function laws and realistic boundary conditions for Navier-Stokes equations (two-fluid or with free surface). On the other hand, NSWE can also be improved. To produce physically correct results, these equations should be completed by friction laws (Chézy, Manning, Darcy-Weisbach and other laws) with properly adjusted coefficients. Recently, bottom boundary layer effects on long waves were studied [39]. Another direction consists in extending NSWE to account for bed material transport as it was recently proposed by Fraccarollo & Capart (2002) [40].