Efficient solvers for shallow-water Saint-Venant equations and debris transportation-deposition models
EEfficient solvers for shallow-water Saint-Venant equations anddebris transportation-deposition models
Florian De Vuyst * Université de Technologie de Compiègne, Alliance Sorbonne Université, Laboratoire deMathématiques Appliquées de Compiègne (LMAC), Compiègne, France.
15 February 2021
Abstract
This research is aimed at achieving an efficient digital infrastructure for evaluating risks and dam-ages caused by tsunami flooding. This research has been mainly focused on the suitable modelingof debris dynamics for a simple (but accurate enough) assessment of damages. For different rea-sons including computational performance and Big Data management issues, we focus our researchon Eulerian debris flow modeling. Rather than using complex multiphase debris models, we ratheruse an empirical transportation and deposition model that takes into account the interaction with themain water flow, friction/contact with the ground but also debris interaction. In particular, for de-bris interaction, we have used ideas coming from vehicular traffic flow modeling. We introduce avelocity regularization term similar to the so-called “anticipation term” in traffic flow modeling thattakes into account the local flow between neighboring debris and makes the problem mathematicallywell-posed. It prevents from the generation of “Dirac measures of debris” at shock waves. As a re-sult, the model is able to capture emerging phenomenons like debris aggregation and accumulations,and possibly to react on the main flow by creating hills of debris and make the main stream deviate.We also discuss the way to derive quantities of interest (QoI), especially “damage functions” fromthe debris density and momentum fields. We believe that this original unexplored debris approachcan lead to a valuable analysis of tsunami flooding damage assessment with Physics-based damagefunctions. Numerical experiments show the nice behaviour of the numerical solvers, including thesolution of Saint-Venant’s shallow water equations and debris dynamics equations.
Keywords.
Tsunami flooding, risk assessment, uncertainty quantification, Saint-Venant shal-low water equations, debris dynamics, quantity of interest, damage function, design of computerexperiment, Big Data, data analytics, database, datawarehouse * email: [email protected] a r X i v : . [ c s . C E ] F e b Executive summary
This work is dedicated to the construction of an efficient numerical infrastructure for evaluating risks anddamages caused by tsunami flooding.In this document, we mainly focus on the suitable modeling of debris dynamics and the derivation ofdamage functions. For different reasons including computational performance and Big Data managementissues, we focus our research on Eulerian debris flow modeling. Rather than using complex multiphasedebris models (see for example [13, 25, 26, 27]) or complex non-Newtonian models, we rather use anempirical transportation and deposition model that takes into account the interaction with the main waterflow, friction/contact with the ground but also debris interaction. In particular, for debris interaction, wehave used ideas coming from vehicular traffic flow modeling ([1, 2]). We introduce a velocity regulariza-tion term similar to the so-called “anticipation term” in traffic flow modeling that takes into account thelocal flow between neighboring debris and makes the problem mathematically well-posed. It preventsfrom the generation of “Dirac measures of debris” at shock waves. As a result, the model is able tocapture emerging phenomenons like debris aggregation and accumulations, and possibly to react on themain flow by creating hills of debris and make the main stream deviate.We also discuss the way to derive quantities of interest (QoI), especially “damage functions” fromthe debris density and momentum fields. We believe that this original unexplored debris approach canlead to a valuable analysis of tsunami flooding damage assessment with Physics-based damage functions.Numerical experiments show the nice behavior of the numerical solvers, including the solution of Saint-Venant’s shallow water equations and debris dynamics equations.The next step is to introduce this debris model into a high-performance code like VOLNA-OP2 [3,4, 5], then use a computer design of experiment (DoCE) with damage analysis, sensitivity analysis anduncertainty quantification. There is clearly a Big Data issue in this computational context and we will tryto propose original numerical methodologies but also suitable parallel software environments to processthe data and extract knowledge in this risk assessment framework.
Among the natural disasters, tsunamis generated by earthquakes or landslides may cause the death ofthousands of people and damage important urban infrastructures. Recent events like the 2004 Indonesiaand 2011 Japan tsunamis have caused dramatic damages by severe flooding (see the figure 1 below).More than a pure flooding, tsunami tidal waves drags millions of debris of all kinds: cars, tree trucks,but also building materials, boats, etc. It is important to prevent coastal population and cities from thesemajor risks and to provide forecast tools for that. Among the means of forecast, numerical simulationappears to be powerful tool to assess damages and costs, achieve sensitivity analysis and uncertaintyquantification (UQ) and so on.Acknowledged tsunami and general flooding computational codes are based on the so-called shallowwater equations, or Saint-Venant equations. They are able to return rather good estimates of the tsunamiwave propagation. The GPU-VOLNA code for example is a GPU-accelerated implementation of the2igure 1: Case of debris accumulation during the 2011 Japan tsunami flooding (courtesy Youtube video).VOLNA code initially developed at CMLA ENS CACHAN by Dutykh et al [3]. It has be ported onlarge-scale GPU supercomputers by means of the OP2 platform [4] that allow large-scale computationsand designs of computer experiments (DoCE) for statistical analysis. Runup elevation and runup dis-tances are the commonly damages estimates computed from shallow water codes, and DoCE are able toreturn high-fidelity metamodels (or emulators) of runup. There could be many other candidate damagefunctions like for example instantaneous or accumulated flow rates that act as loadings on mechanicalstructures and as the main factor of fatigue breaking and generation of debris.However, Indonesia and Japan tsunami events have shown us that debris flows are another importantfeature of tsunami flooding and, as such, have to be taken into account into the general flooding flow.This need is originally the reason of the present research.One can observe that damage and risk assessment of tsunami flooding is a multicriteria/multiobjectiveproblem subject to a high-dimensional configuration space, where each evaluation requires heavy com-putations of high-performance parallel supercomputers. There is clearly a “Big Data” dimension of theproblem where we need to mine high-dimensional data (as computational results), extract data sum-maries, descriptors, indicators and risk/damage functions and provide emulators. The issue to performthis data analysis in a parallel efficient way is a current active field of Research (map-reduce type algo-rithms, ...). Another issue is the way to perform a design of computer experiment that can be done inan incremental manner in order to reduce metamodeling errors or uncertainty. This is also a Big Datadatawarehouse dimension that allows for easy analysis and visualization.At our research developmental stage, this work is mainly dedicated to the suitable modeling of debrisdynamics that allows a rather good estimation of debris flow while being able to return damage functions.The models will also be searched as simple as possible to allow “fast” evaluations and a full explorationof the configuration space for statistical purposes. We also discuss the suitable design Saint-Venantsolvers for tidal wave and coast flooding applications where the capture of dry-wet phase transition hasto be done in an efficient and accurate manner. For that we introduce a class of finite volume schemeswe recently developed for multimaterial compressible flow applications: the so-called Lagrange-fluxschemes. 3
Introducing the Saint-Venant shallow-water equations
Saint-Venant shallow water equations are a simplified two-dimensional model of a truly three-dimensionalflow of an incompressible fluid in contact with air and subject to gravity forces. These equations aremainly used for tsunami modeling as well as river flooding and for the analysis of dam break events.In what follows, we will denote g the gravity constant, z = z ( x, y ) will be the bathymetry ( z < ) orlandscape topography ( z > ), h = h ( x, y ) the liquid depth and u the vector field of depth-averagedvelocity. In what follows, we will assume that z is a Lipschitz continuous function. Assuming that the z -dependency of the flow can be reduced to the knowledge of the water height, and supposing no grounddrag forces, the mass and momentum balance equations read ∂ t h + ∇ · ( h u ) = 0 , (1) ∂ t ( h u ) + ∇ · ( h u ⊗ u ) + ∇ p = − gh ∇ z, (2)where p = p ( h ) = g h (3)Here the notation “ p ” is used to emphasize the analogy with the usual Euler equations of compressiblegas dynamics, even it is not a pressure from a thermodynamical point of view.It is well-known that this system of nonlinear partial differential equations is hyperbolic. The quantity c = (cid:114) dpdh = (cid:112) gh (4)plays the role of a “speed of sound” in the system. For smooth solutions, the energy quantity E definedby E = h | u | g h . (5)is conserved and satisfies the additional balance equation ∂ t E + ∇ · ( E u ) = − gh ∇ z · u . As the energy is a convex function of the conservative variables h and q = h u , it can also be consideredas an entropy of the system. For general, possibly nonsmooth discontinuous solutions, we look for thephysical (entropy) weak solution that fulfills the partial differential inequality ∂ t E + ∇ · ( E u + p u ) + gh ∇ z · u ≤ (6)in the sense of distributions.For numerical discretization, stable conservative entropy schemes are searched in order to ensureconvergence toward the entropy weak solution. In this work, we propose to use the recent family ofLagrange-flux finite volume schemes [6] that provide both numerical stability and efficiency. The con-struction of the Lagrange-flux is based on a Lagrangian remapping process with a particular treatment oftime discretization. This is detailed in the next section.4 Lagrange-flux schemes for compressible Euler-type equations
For complex compressible flows involving multiphysics phenomenons like e.g. high-speed elastoplastic-ity, multimaterial interaction, plasma, gas-particles, multiphase flows etc., a Lagrangian description ofthe flow is generally preferred for easier physical coupling. To ensure robustness, some spatial remap-ping on a regular mesh may be added. A particular case is the family of the so-called Lagrange+remapschemes [7], also referred to as Lagrangian remapping that apply a remap step on a reference (sayEulerian) mesh after each Lagrangian time advance step. Acknowledged legacy codes implementingremapped Lagrange solvers usually define thermodynamical variables at cell centers and velocity vari-ables at mesh nodes (see figure 2).Figure 2: “Legacy” staggered Lagrange-remap scheme: thermodynamical variables are located at cellcenters (circles) whereas velocity variables are located at cell nodes (squares).
In Poncet et al. [8], we have achieved a node-based performance analysis of a reference legacy Lagrange-remap hydrodynamics solver used in industry. By analyzing each kernel of the whole algorithm, usingroofline-type models [9] on one side and refined Execution Cache Memory (ECM) models [10], [11]on the other side, we have been able not only to quantitatively predict the performance of the wholealgorithm — with relative errors in the single digit range — but also to identify a set of features that limitthe whole global performance. This can be roughly summarized into three features:1. Staggered velocity variables involve a rather big amount of communication to/from CPU cachesand memory with low arithmetic intensity, thus lowering the whole performance;2. Alternating direction (AD) strategies (see the appendix in [12]) or more specifically AD remappingprocedures also generate too much communication with a loss of CPU occupancy and a rather poormulticore scalability.3. For multimaterial hydrodynamics using VOF-based interface reconstruction methods, there is astrong loss of performance due to some array indirections and noncoalescent data in memory.Vectorization of such algorithms is also not trivial.From these observations and as a result of the analysis, we decided to “rethink” Lagrange-remap schemes,with possibly modifying some aspects of the solver in order to improve node-based performance of the5ydrocode solver. We have searched for alternative formulations that lower communication and improveboth arithmetic intensity and SIMD property of the algorithm. This redesign methodology has given usideas of innovative Eulerian solvers. The so-called
Lagrange-flux schemes appear to be very promisingin the Computational Fluid Dynamics (CFD) extended community, including geophysical problems.Starting from a “legacy” staggered Lagrange-remap solver and related observed performance mea-surements, we want to improve the performance by modifying the computational approach under thefollowing constraints and requirements:1. A Lagrangian solver must be used (for multiphysics coupling issue).2. To reduce communication, we prefer to use collocated cell-centered variables rather than a stag-gered scheme.3. To reduce communication, we prefer use a direct multidimensional remap solver rather than split-ted alternating direction AD projections.4. The method should be simply extended to second-order accuracy (in both space and time).Before going further, let us first comment the above requirements. The second requirement should implythe use of a cell-centered Lagrange solver. Fairly recently, Després and Mazeran in [14] and Maire and etal. [15] have proposed pure cell-centered Lagrangian solvers based on the reconstruction of nodal veloc-ities. In our study, we will examine if it is possible to use approximate and simpler Lagrangian solvers inthe Lagrange+remap context, in particular for the sake of performance. The fourth assertion requires afull multidimensional remapping step, probably taking into account geometric elements (deformation ofcells and edges) if we want to ensure high-order accuracy remapping. We have to find a good trade-offbetween simplifications-approximations and accuracy (or properties) of the numerical solver.
As example, let us consider the compressible Euler equations for two-dimensional problems. Denoting ρ, u = ( u i ) i , i ∈ { , } , p and E the density, velocity, pressure and specific total energy respectively,the mass, momentum and energy conservation equations are ∂ t U (cid:96) + ∇ · ( u U (cid:96) ) + ∇ · π (cid:96) = 0 , (cid:96) = 1 , . . . , , (7)where U = ( ρ, ( ρu i ) i , ρE ) , π = (cid:126) , π = ( p, T , π = (0 , p ) T and π = p u . For the sake ofsimplicity, we will use a perfect gas equation of state p = ( γ − ρ ( E − | u | ) , γ ∈ (1 , . The speedof sound c is given by c = (cid:112) γp/ρ .For any volume V t that is advected by the fluid, from the Reynolds transport theorem we have ddt (cid:90) V t U (cid:96) d x = (cid:90) ∂ V t { ∂ t U (cid:96) + ∇ · ( u U (cid:96) ) } d x = − (cid:90) ∂ V t π (cid:96) · ν dσ where ν is the normal unit vector exterior to V t . This leads to a natural explicit finite volume scheme inthe form | K n +1 ,L | ( U (cid:96) ) n +1 ,LK = | K | ( U (cid:96) ) nK − ∆ t n (cid:88) A n + 12 ,L ⊂ ∂K n + 12 ,L | A n + ,L | π n + ,LA · ν n + ,LA . (8)6n expression (2), the superscript “L” indicates the Lagrange evolution of the quantity. Any Euleriancell K is deformed into the Lagrangian volume K n + ,L at time t n + , and into the Lagrangian vol-ume K n +1 ,L at time t n +1 . The pressure flux terms through the edges A n +1 / ,L are evaluated at time t n + in order to get second-order accuracy in time. Of course, that means that we need a predictor step for thevelocity field u n + ,L at time t n + (not written here for simplicity).From now on, we will use the simplified notation v n + = u n + ,L . The remapping step consists in projecting the fields U (cid:96) defined at cell centers K n +1 ,L onto the initial (ref-erence) Eulerian mesh with cells K . Starting from an interpolated vector-valued field I n +1 ,L U n +1 ,L ,we project the field on piecewise-constant function on the Eulerian mesh, according to the integral for-mula U n +1 K = 1 | K | (cid:90) K I n +1 ,L U n +1 ,L ( x ) d x . (9)Practically, they are many ways to deal with the projection operation (9). One can assemble elementaryprojection contributions by computing the volume intersections between the reference mesh and the de-formed mesh. But this procedure requires the computation of all the geometrical elements. Moreover,the projection needs local tests of projection with conditional branching (think about the very differentcases of compression, expansion, pure translation, etc). Thus the procedure is not SIMD and with po-tentially a loss of performance. The incremental remapping can also interpreted as a transport/advectionprocess, as already emphasized by Dukowicz and Baumgardner [16]. Let us now write a different original formulation of the remapping process that does not explicitelyrequires the use of geometrical elements. In this step, there is no time evolution of any quantity, and insome sense we have ∂ t U = 0 , that we rewrite ∂ t U = ∂ t U + ∇ · ( − v n + U ) + ∇ · ( v n + U ) = 0 . We decide to split up this equation into two substeps, a backward convection and a forward one:i) Backward convection: ∂ t U + ∇ · ( − v n + U ) = 0 . (10)ii) Forward convection: ∂ t U + ∇ · ( v n + U ) = 0 . (11)Each convection problem is well-posed on the time interval [0 , ∆ t n ] under a standard CFL condition.Let us now focus into these two steps and the way to solve them.7 .4.1 Backward convection in Lagrangian description After the Lagrange step, if we solve the backward convection problem (4) over a time interval ∆ t n usinga Lagrangian description, we have | K | ( U (cid:96) ) n,(cid:63)K = | K n +1 ,L | ( U (cid:96) ) n +1 ,LK . (12)Actually, from the cell K n +1 ,L we go back to the original cell K with conservation of the conservativequantities. For (cid:96) = 1 (conservation of mass), we have | K | ρ n,(cid:63)K = | K n +1 ,L | ρ n +1 ,LK showing the variation of density by volume variation. For (cid:96) = 2 , , , it is easy to see that both velocityand specific total energy are kept unchanged is this step: u n,(cid:63) = u n +1 ,L , E n,(cid:63) = E n +1 ,L . Thus, this step is clearly computationally inexpensive.
From the discrete field ( U n,(cid:63)K ) K defined on the Eulerian cells K , we then solve the forward convectionproblem over a time step ∆ t n under an Eulerian description. A standard Finite Volume discretization ofthe problem will lead to the classical time advance scheme U n +1 K = U n,(cid:63)K − ∆ t n | K | (cid:88) A ⊂ ∂ K | A | U n + ,(cid:63)A ( v n + A · ν A ) (13)for some interface values U n + ,(cid:63)A defined from the local neighbor values U n,(cid:63)K . We finally get the ex-pected Eulerian values U n +1 K at time t n +1 .Notice that from (6) and (7) we have also | K | U n +1 K = | K n +1 ,L | U n +1 ,LK − ∆ t n (cid:88) A ⊂ ∂K | A | U n + ,(cid:63)A ( v n + A · ν A ) (14)thus completely defining the remap step under the finite volume scheme form (14). Let us emphasizethat we do not need any mesh intersection or geometric consideration to achieve the remapping process.The finite volume form (14) is now suitable for a straightforward vectorized SIMD treatment. From (14)it is easy to achieve second-order accuracy for the remapping step by usual finite volume tools (MUSCLreconstruction + second-order accurate time advance scheme for example). Let us note that the Lagrange+remap scheme is actually a conservative finite volume scheme: putting (8)into (14) gives for all (cid:96) : ( U (cid:96) ) n +1 K = ( U (cid:96) ) nK − ∆ t n | K | (cid:88) A n + 12 ,L ⊂ ∂K n + 12 ,L | A n + ,L | ( π (cid:96) ) n + ,LA · ν n + ,LA − ∆ t n | K | (cid:88) A ⊂ ∂K | A | ( U (cid:96) ) n + ,(cid:63)A ( v n + A · ν A ) (15) ( U (cid:96) ) n +1 K = ( U (cid:96) ) nK − ∆ t n | K | (cid:88) A ⊂ ∂K | A | (cid:32) | A n + ,L || A | ( π (cid:96) ) n + ,LA · ν n + ,LA (cid:33) − ∆ t n | K | (cid:88) A ⊂ ∂K | A | (cid:16) ( U (cid:96) ) n + ,(cid:63)A ( v n + A · ν A ) (cid:17) . (16) We recognize into (16) pressure-related fluxes and convective numerical fluxes.
From conclusions of the discussion above, we would like to be free from any “complex” collocatedLagrangian solver involving complex geometric elements. Another difficult point is to define the defor-mation velocity field v n + at time t n + , in an accurate manner.In what follows, we are trying to deal with time accuracy in a different manner. Let us come backto the Lagrange+remap formula (16). Let us consider a “small” time step t > that fulfills the usualstability CFL condition for explicit schemes. We have ( U (cid:96) ) K ( t ) = ( U (cid:96) ) nK − t | K | (cid:88) A ⊂ ∂K | A | (cid:18) | A L ( t/ || A | ( π (cid:96) ) LA ( t/ · ν LA ( t/ (cid:19) − t | K | (cid:88) A ⊂ ∂K | A | ( U (cid:96) ) (cid:63)A ( t/ v A ( t/ · ν A . By making t tend to zero, ( t > ), we have A L ( t/ → A , ( π (cid:96) ) L ( t/ → π (cid:96) , v ( t/ → u , ( U (cid:96) ) (cid:63) → U (cid:96) ,then we get a semi-discretization in space of the conservation laws. That can be seen as a method-of-linesdiscretization (see [17]): d ( U (cid:96) ) K dt = − | K | (cid:88) A ⊂ ∂K | A | (( π (cid:96) ) A · ν A ) − | K | (cid:88) A ⊂ ∂K | A | ( U (cid:96) ) A ( u A · ν A ) . (17) We get a classical finite volume method in the form dU K dt = − | K | (cid:88) A ⊂ ∂K | A | Φ A with a numerical flux Φ A whose components are (Φ (cid:96) ) A = ( U (cid:96) ) A ( u A · ν A ) + ( π (cid:96) ) A · ν A . (18)In (17), pressure fluxes ( π (cid:96) ) A and interface normal velocities ( u A · ν A ) can be computed from an approx-imate Riemann solver in Lagrangian coordinates (for example the Lagrangian HLL solver, see Toro [18]for example). Then, the interface states ( U (cid:96) ) A can be computed from a upwind process according to thesign of the normal velocity ( u A · ν A ) . To get higher-order accuracy in space, one can use a standardMUSCL reconstruction + slope limiting process. At this stage, because there is no time discretization,all acts on the Eulerian mesh and fluxes are defined at the edges the the Eulerian cells.To get high-order accuracy in time, one can then apply a standard high-order time advance scheme(Runge-Kutta 2, etc.). For the second-order Heun scheme for example, we have the following algorithm:1. Compute the time step ∆ t n subject to some CFL condition;9. Predictor step . MUSCL reconstruction + slope limitation: from the discrete values U nK , computea discrete gradient for each cell K .3. Use a Lagrangian approximate Riemann solver to compute pressure fluxes π nA and interface ve-locities u nA
4. Compute the upwind edge values ( U (cid:96) ) nA according to the sign of ( u nA · ν A ) ;5. Compute the numerical flux Φ nA as defined in (12);6. Compute the first order predicted states U (cid:63),n +1 K : U (cid:63),n +1 K = U nK − ∆ t n | K | (cid:88) A ⊂ ∂K | A | Φ nA Corrector step . MUSCL reconstruction + slope limitation: from the discrete values U (cid:63),n +1 K ,compute a discrete gradient for each cell K .8. Use a Lagrangian approximate Riemann solver to compute pressure fluxes π (cid:63),n +1 A and interfacevelocities u (cid:63),n +1 A
9. Compute the upwind edge values ( U (cid:96) ) (cid:63),n +1 A according to the sign of ( u (cid:63),n +1 A · ν A ) ;10. Compute the numerical flux Φ (cid:63),n +1 A as defined in (12);11. Compute the second-order accurate states U n +1 K at time t n +1 : U n +1 K = U nK − ∆ t n | K | (cid:88) A ⊂ ∂K | A | Φ nA + Φ (cid:63),n +1 A . One can appreciate the simplicity of the numerical solver.
A HLL approximate Riemann solver [18] in Lagrangian coordinates can be used to easily computeinterface pressure and velocity. For a local Riemann problem made of a left state U L and a right state U R , the contact pressure p (cid:63) is given by the formula p (cid:63) = ρ R p L + ρ L p R ρ L + ρ R − ρ L ρ R ρ L + ρ R max( c L , c R ) ( u R − u L ) , (19)and the normal contact velocity u (cid:63) by u (cid:63) = ρ L u L + ρ R u R ρ L + ρ R − ρ L + ρ R p R − p L max( c L , c R ) (20)leading to very simple operations. 10igure 3: Second-order Lagrange-flux scheme on reference Sod’s 1D shock tube problem. Time is T = 0 . , 384 mesh points. Here Sweby’s slope limiter with coefficient 1.5 is used. Shock tube problems.
An as example, we test the Lagrange-flux scheme presented in section 4.6on the reference one-dimensional Sod shock tube problem [19]. We use a Runge-Kutta 2 (RK2) timeintegrator and a MUSCL reconstruction using the second-order Sweby slope limiter [20] φ ( a, b ) = ( ab > sign ( a ) max (cid:0) min( | a | , β | b | ) , min( β | a | , | b | ) (cid:1) with coefficient β = 1 . . We use a uniform grid made of 384 points. The final time is T = 0 . and aCFL number equal to 0.25 . On figure 3, one can observe a good behaviour of the Eulerian solver, withrather sharp discontinuities and low numerical diffusion into rarefaction fans. This section is dedicated to the derivation of suitable Lagrange-flux schemes for the Saint-Venant equa-tions with gravity source terms.
For smooth solutions the momentum balance equation can be written ρD t u + ∇ ( g ( h + z )) = 0 . A particular solution of the Saint Venant equations is the so-called “lake-at-rest” state, meaning that afluid at rest ( u = 0 everywhere) with topography z satisfies the condition h + z = C (21)for each connected part of the fluid domain. It is expected that numerical solver also satisfies the fluid-at-rest condition at the discrete level, meaning that ( h + z ) K = C K belonging to the same connected part of the discrete domain Ω h . This conditionis name the “well-balanced” property. We are looking for Lagrange-flux schemes that satisfy the well-balance property. For debris dynamics and flooding applications, it is particularly important to have a robust numericalscheme able to deal with small water depth and/or to handle wet-dry state transitions. Rather thanconsidering wet-dry interface reconstruction techniques, for performance purpose we have decided usean interface capturing scheme that computes the variables in all the wet+dry domain. That means that thenumerical solver has to deal with vanishing water depth, see h = 0 exactly. There are three computationalissues: first we have to ensure positivity of the water depth; Secondly, when h vanishes, we have also todeal with vanishing propagation speed c = √ gh . Finally, because conservative variables are h and h u respectively, the fluid velocity in ordinarily computed as u := h u h . Of course this leads to ill-posed division operations for vanishing depths. Moreover, the division has nosense when h = 0 exactly and something different has to be done to compute the velocity. For the sake of simplicity, only first-order accurate schemes are considered in this document. We usestandard finite volume notations with K for a generic finite volume, A for a generic edge, ν A a normalunit vector at edge A , ∆ t n the time step at time t n . For Saint Venant equations, we have looking forfinite volumes explicit schemes in the form h n +1 K = h nK − ∆ t n | K | (cid:88) A ⊂ ∂K | A | Φ h ( U nK , U nK A , ν A ) , (22) ( h u ) n +1 K = ( h u ) nK − ∆ t n | K | (cid:88) A ⊂ ∂K | A | Φ hu ( U nK , U nK A , ν A ) − ∆ t n | K | (cid:88) A ⊂ ∂K | A | p nA ν A − ∆ t n | K | g ¯ h nK (cid:88) A ⊂ ∂K | A | z A ν A , (23)where Φ h represents the numerical mass flux, Φ hu represents the convective flux related to the momen-tum variable, p A (resp. z A ) is the pressure flux (resp. topography value) at the edge A . In (23), the lastterm is a finite volume discretization of the source term − ∆ t n gh ∇ z into cell K .To entirely define the numerical scheme, we have to characterize both convective fluxes, pressureflux as well as the values for z A and ¯ h nK .The Lagrange-flux scheme methodology exposed above explains how to discretize convective fluxesand pressure fluxes as soon as an approximate Riemann solver is determined for the system. The choicefor ¯ h nK and z A is guided by the well-balanced property. This is developed in the next section.12 .3 Lagrangian approximate Riemann solver and well-balanced property For simplicity purposes, we consider in this section the one-dimensional Saint-Venant equations. Wehave first to achieve an approximate Riemann solver that takes into account the gravity source termeffect. Let us rewrite the momentum equation: ∂ t ( hu ) + ∂ x ( hu ) + ∂ x p = − gh∂ x z with p = p ( h ) = g h as usual. For a lake-at-rest solution, we have ∂ x p + hg∂ x z = h ( h + z ) ∂ x ( h + z ) − gz∂ x ( h + z ) = ∂ x (cid:20) g ( h + z ) (cid:21) − g z ∂ x ( h + z ) = 0 . Let us now consider a Riemann problem made of two constant states U L = ( h L , ( hu ) L ) and U R =( h R , ( hu ) R ) with a discontinuous topography function with left and right values z L and z R . In thissection, we do not consider dry conditions, i.e. we assume that h L , h R > . Let us denote z (cid:63) a localmean topography value, function of both z L and z R . The replace the initial momentum balance equationof the local conservation law ∂ t ( hu ) + ∂ x ( hu ) + ∂ x Π = 0 (24)with the pseudo pressure Π defined by Π = Π( h, z ) = g ( h + z ) − gz (cid:63) ( h + z ) . (25)Let us emphasize that the modified system still has “lake-at-rest” solutions with constant Π giving theproperty h + z = C .We retrieve a system of conservation laws that is similar to the isentropic Euler equations with apseudo pressure defined by (25), and a propagation speed ˜ c such that ˜ c = ∂ Π ∂h = gh + g ( z − z (cid:63) ) . (26)It is important to notice that g ( z − z (cid:63) ) has no sign, so this local approximate problem has only a sense ifinvolved depths are such that h ≥ z (cid:63) − z . The Lagrangian form of the equation writes hD t τ − ∂ x u = 0 ,hD t u + ∂ x Π = 0 with τ = h − . Introducing the mass variable m such that dm = hdx , we have the conservative script D t τ − ∂ m u = 0 ,D t u + ∂ m Π = 0 . We now introduce an approximate Riemann solver made of two waves of respective speeds a L and a R and intermediate state with depth h (cid:63) and velocity u (cid:63) . By writing the Rankine-Hugoniot jump conditionson the second equation for the two waves, we have the two relations a L ( u (cid:63) − u L ) + Π (cid:63) − Π L = 0 , (27) a R ( u R − u (cid:63) ) + Π R − Π (cid:63) = 0 . (28)13he “acoustic solver” approximation consider a L = − σh L and a R = σh R for a propagation speed σ .The sub-characteristic conditions requires that σ ≥ max( c L , c R ) . From (27) and (28) we get the meanvelocity u (cid:63) = h L u R + h R u R h L + h R − σ h L + h R (Π R − Π L ) . (29)There is a convenient choice for z (cid:63) that gives an interesting and simple value for the difference of pseudo-pressure Π R − Π L : Proposition 1.
For z (cid:63) = ( z L + z R ) , we have Π R − Π L = g h L + h R h + z ) R − ( h + z ) L ] . Proof.
It is easy to check that Π R − Π L = p R − p R + gh R ( z R − z (cid:63) ) − gh L ( z L + z (cid:63) ) + g (cid:18) z L + z R − z (cid:63) (cid:19) ( z R − z L )= p R − p L − gz (cid:63) ( h R − h L ) + g (( hz ) R − ( hz ) L ) + g (cid:18) z L + z R − z (cid:63) (cid:19) ( z R − z L ) and, because p R − p L = g h L + h R ( h R − h L ) , one finds the result.Following the choice from proposition 1, we have the formula for u (cid:63) : u (cid:63) = h L u R + h R u R h L + h R − σ g [( h + z ) R − ( h + z ) L ] . (30)We then have the trivial result: Proposition 2.
For lake-at-rest conditions, i.e. u L = u R = 0 and ( h + z ) L = ( h + z ) R , we have u (cid:63) = 0 . As a direct consequence of Proposition 2, for lake-at-rest conditions, convective fluxes will be zero.Let us go now to the computation of the mean depth h (cid:63) . By writing the Rankine-Hugoniot jumpconditions on the first equation for the two waves, we have the relations h L σ ( τ (cid:63) − τ L ) = u (cid:63) − u L , − h R σ ( τ R − τ (cid:63) ) = u R − u (cid:63) , thus giving τ (cid:63) = 2 h L + h R + 12 σ h L + h R ( u R − u L ) or equivalently h (cid:63) = h L + h R σ ( u R − u L ) . (31)Then the mean pressure p (cid:63) is simply computed as p (cid:63) = g ( h (cid:63) ) . We have the following well-balancedproperty : Proposition 3.
Choosing z (cid:63) = ( z L + z R ) , for lake-at-rest conditions u L = u R = 0 and ( h + z ) L =( h + z ) R , we have h (cid:63) = h L + h R so that h (cid:63) + z (cid:63) = ( h + z ) L = ( h + z ) R . emark 1. In (31) , one can notice a “compressibility” term κ = 1 + 12 σ ( u R − u L ) . For u R − u L < which expresses local compression conditions, in order to keep κ positive, we need tochoose σ such that σ > | u R − u L | . On can choose for example σ ← max( σ, − min(0 , u R − u L )) . Now we discuss the discretization of the gravity source term. For one-dimensional problems and auniform mesh grid, the numerical scheme writes ( hu ) n +1 j = ( hu ) nj − ∆ t ∆ x (cid:16) Φ nhu,j +1 / − Φ nhu,j − / (cid:17) − ∆ t ∆ x (cid:16) p nj +1 / − p nj − / (cid:17) − ∆ t ∆ x g ¯ h nj (cid:0) z j +1 / − z j − / (cid:1) , (32)where the Φ nhu,j +1 / are the momentum convective fluxes and the p nj +1 / are the pressure at cell inter-faces, prescribed by the approximate Riemann solver. In (32), we have still to find a convenient choicefor both quantities ¯ h nj and interface topography z j +1 / . This choice is guided by the following result: Proposition 4.
Let us consider approximate Riemann values u (cid:63) and h (cid:63) given by the formulas (30) and (31) respectively. Then, the initial lake-in-rest conditions, i.e., u j = 0 and ( h + z ) j = c for all j ,then the choice z j +1 / = z j + z j +1 , ¯ h nj = h nj − / + h nj +1 / (33) provides a well-balanced numerical scheme, i.e. ( h + z ) nj = C for all j , for all n ∈ N .Proof. According to Proposition 2, the approximate Riemann solver returns a null velocity value u (cid:63) under lake-in-rest conditions. Thus, convective fluxes are zero. Let us assume that we have lake-in-restconditions at time t n . Then we have ( hu ) n +1 j = − ∆ t ∆ x (cid:16) p nj +1 / − p nj − / − g ¯ h nj (cid:0) z j +1 / − z j − / (cid:1)(cid:17) = − ∆ t ∆ x (cid:32) g (cid:32) ( h nj +1 / ) − ( h nj − / ) (cid:33) − g ¯ h nj (cid:0) z j +1 / − z j − / (cid:1)(cid:33) = − ∆ t ∆ x (cid:32) g h nj − / + h nj +1 / (cid:16) h nj +1 / − h nj − / (cid:17) − g ¯ h nj (cid:0) z j +1 / − z j − / (cid:1)(cid:33) Then for ¯ h nj defined in (33), we have ( hu ) n +1 j = − ∆ t ∆ x g h nj − / + h nj +1 / (cid:16) ( h nj +1 / + z j +1 / ) − ( h nj − / − z j − / ) (cid:17) . According the Proposition 3, the right hand side is 0.15 .5 Dealing with dry-wet phase transitions
As mentioned above, when h vanishes, we have to deal with vanishing propagation speed c = √ gh andby the fact that approximate Riemann problem quantities may leads to some divisions by zero, makingthe computational method unstable. Moreover, the division u := h u h (34)required to compute the velocity may be arbitrary ill-conditioned for h close to zero. This generallyproduces spuriously large velocities at wet-dry transition fronts, making the time step collapse or simplyleading to a code crash. This section addresses these issues.Let us first clarify some behavior of the solutions at vanishing water depth h → ( h (cid:54) = 0 ). Forsmooth solutions, the momentum equation can be written hD t u + gh ∇ ( h + z ) = 0 . Dividing by h , we find the Burgers-like equation ∂ t u + u · ∇ u + g ∇ ( h + z ) = 0 . Neglecting h in the equation gives ∂ t u + u · ∇ u = − g ∇ z. (35)We find out an autonomous equation in u . We want to insist on the fact that, at vanishing h , the valueof velocity is not arbitrary, but solution of the above equation which is subject to gravity acceleration.In particular, u has not to be zero as often considered and encountered in the literature. This is becausefriction terms are not taken into account so that fluid is sliding on the ground and accelerating according togravity forces. Let us finally remark that equation (35) can be written in conservative form by consideringany function η > solution of the conservation law ∂ t η + ∇ · ( η u ) = 0 . (36)Then we can rewrite (35) as ∂ t ( η u ) + ∇ · ( η u ⊗ u ) = − gη ∇ z. (37) Let us recall the interface velocity and water depth obtained with our Lagrangian approximate Riemannsolver, for “wet-wet” conditions: u (cid:63) = h L u R + h R u R h L + h R − σ g [( h + z ) R − ( h + z ) L ] ,h (cid:63) = h L + h R σ ( u R − u L ) σ chosen for example as σ = max ( c L , c R , − ( u R − u L ) − ) .There are two difficulties: either ( h L + h R ) is small or σ is small, what may lead to some ill-conditioneddivisions. A way to fix the problem of “speed of sound” is to use a larger non-vanishing propagationspeed σ . This can be theoretically justified by a relaxation technique [22, 23].Let us remark that the expression for u (cid:63) still holds for h L h R = 0 but h L + h R > . The remainingdifficulty is the “dry-dry” case or vanishing h on both two sides. In this case we have to use somethingdifferent. As mentioned above, at least for smooth solutions the momentum equation can be rewritten as D t u + ∂ x ( g ( h + z )) = 0 . Looking for an approximate Riemann solver considering this equation leads to the mean velocity u (cid:63) = u L + u R − σ g [( h + z ) R − ( h + z ) L ] (38)or u (cid:63) = u L + u R − σ g ( z R − z L ) (39)if h L and h R are neglected. From the practical computational point of view, we have to define a thresholdon h L + h R to use either (30) or (39) for computing u (cid:63) . Remark that expression (39) can be seen as aparticular case of (30) with h L = h R = 1 . The division operation u := ( h u ) /h cannot be used for pure dry conditions ( h = 0 ) or vanishing waterdepth ( h (cid:28) . We have to do something different. In this work, we propose to use an additional“dry velocity” equation. Let us denote u dry the “dry velocity”. It is solution of the transport-relaxationequation ∂ t u dry + u dry · ∇ u dry = − g ∇ z + u − u dry µ (40)with µ > , µ (cid:28) as relaxation coefficient. There are two strategies: either we use a non-conservativescheme that discretizes the equation (40), or we use a conservative scheme on the conservative form ∂ t ( η u dry ) + ∇ · ( η u dry ⊗ u dry ) = − gη ∇ z + η u − u dry µ . (41)using the additional variable η solution of (36). Let us comment is more details the numerical procedure.First, we have to define a field η n and a field u dry at time t n . Because η n is arbitrary, it is simply chosenas η nK = 1 ∀ K ∈ T h . We project the initial dry velocity onto the actual velocity field u n : ( u dry ) nK = u nK ∀ K ∈ T h . Then we apply the Lagrange-flux explicit scheme on the system ∂ t η + ∇ · ( η u dry ) = 0 ,∂ t ( η u dry ) + ∇ · ( η u dry ⊗ u dry ) = − gη ∇ z. η n +1 j = 1 − ∆ t ∆ x (cid:16) Φ nη,j +1 / − Φ nη,j − / (cid:17) ,η n +1 j ( u dry ) n +1 j = u nj − ∆ t ∆ x (cid:16) Φ nηu,j +1 / − Φ nηu,j − / (cid:17) − ∆ t ∆ x g (cid:0) z j +1 / − z j − / (cid:1) . Remark 2.
For water depth vanishing conditions, the auxiliary variable η can be interpreted as a waterdepth renormalization, where a scaling is operated in order to return a rescaled depth of order O (1) .Because there is an uncertainty on the good scale to apply, we simply express the rescaling by η n = 1 . Once the candidate velocity u dry in case of dry or near-dry conditions is computed, we have todefine smoothed switching strategies to compute the actual velocity. We proceed as follows: for a givenmomentum value ( h u ) and a given dry velocity u dry , we search for a velocity vector u (cid:63) solution of theTykhonov-regularized mean square problem ( P ε ) min v (cid:107) ( h u ) − h w (cid:107) + ε (cid:107) w − u dry (cid:107) (42) u (cid:63) := h ( h u ) + ε u dry h + ε . (43)Of course, for non vanishing h and rather small ε , we have u (cid:63) ≈ u and for vanishing h , we get u (cid:63) ≈ u dry . This is similar to usual velocity fixes found in the literature, as in Kurganov and Petrova [21]for example, but the authors actually consider u dry = 0 . The choice u dry = 0 creates an artificialviscosity/damping into depth-vanishing regions and in particular water cannot completely drain away.For the choice of ε , it can be adapted to the cell size as suggested by Kurganov-Petrova. This section is dedicated to the modeling of debris dynamics carried away by flooding or tsunami waves,including possible coalescence and deposition effects. For foreseen risk analysis purposes, the modelhas to be able to return high-level debris features and damage functions with sufficient fidelity. One mustfind a good trade-off between computational efficiency and ability to return quantities of interest likedamage functions.
Conditions of use of the debris models for risk analysis lead to the following requirements:1. The computational debris models have to be rather efficient in terms of computational complexity,preferably at the order of the Saint-Venant numerical solvers;2. Moreover, the model complexity has to be independent of the number of debris, allowing forlarge-scale debris computations;3. The model has to be able to return quantities of interest like damage functions for risk analysis.18 .2 Derivation of a model and model analysis
For transport-dominated problems, we have the choice between Lagrangian or Eulerian models. La-grangian models track trajectories debris particles or group of particles and, because of that, are moreaccurate than Eulerian models. Unfortunately, the algorithmic complexity is proportional to the numberof debris particles and thus do not fulfill the requirement nb 2. Thus we rather move toward Eulerianmodels, expressed in terms of averaged density quantities.In the literature, one can find volume-averaged multiphase-based debris models [25, 26, 27, 28, 29].But these models are known to be highly computationally expensive and thus are irrelevant in our case.It is here proposed to emulate well-known vehicle-based traffic flow for modeling debris dynamics,up to some adjustments.
For a one-dimensional liquid+debris flow, debris are driven by the main liquid flow and are followingthemselves. Let us first consider a Lagrangian description of the flow. For a debris particle j of mass m ,located at position x j ( t ) at time t at velocity v j ( t ) and following the debris particle number ( j + 1) , wehave the motion equations dx j dt = v j , (44) m dv j dt = − D j ( d ) − F j ( f ) + I j ( t ) , (45)where the term D j ( t ) represent a drag force due to the driving shallow water flow, F j ( t ) is a groundfriction term for debris in contact with ground and I j ( t ) is a term that represents the interaction of debrisnb j with its neighboring ones. These three terms have to be modelized and closed with respect to thevariables of computations.For the drag term, one can simply use the relaxation term D j ( t ) = m u ( x j ( t )) − v j ( t ) τ D where τ D > is a characteristic relaxation time. For the ground friction term, one can also use arelaxation term F j ( t ) = − ω F m v j ( t ) , but the relaxation rate ω F is a function of the water depth h . Let us introduce h f > the characteristicdebris plunge depth. It is expected that a debris rapidly stops when h < h f by ground contact. On theother hand, when h (cid:29) h f , there is no friction with ground and the characteristic relaxation time shouldbe infinite. We then propose to use the rather simple function ω F = ω F ( h ) = mτ F max (cid:18) , h f h (cid:19) max (cid:18) , − hh f (cid:19) β , (46)where β > and τ F > is a constant friction characteristic time.Debris interaction terms are needed to take into account unresolved small scale local water flowtowards debris (rear debris recirculation, vortexes, suction, ...). For one-dimensional problems, debris19bjects cannot overtake themselves and the debris order is preserved. This can be expressed by anacceleration/slowing down term as used in vehicular traffic flow [1, 2] (referred to as an anticipationterm in car-following models). Here, we empirically define the interaction term as a velocity relaxationbetween neighboring debris: if x j +1 ( t ) − x j ( t ) > , for a flow moving to the right, one can define forexample I j ( t ) as I j ( t ) = m a v j +1 ( t ) − v j ( t ) x j +1 ( t ) − x j ( t ) (47)where a > is homogeneous to a speed. The speed quantity a can also be modeled. One could forexample take it proportional to the local debris velocity, i.e. a = λv j ( t ) , where λ > is a dimensionlessconstant. For a general debris, we then consider following interaction term: I j ( t ) = λ m v j ( t ) v j +1 ( t ) − v j ( t ) x j +1 ( t ) − x j ( t ) (48)As a summary, we consider the following discrete dynamical system: ˙ x j = v j , (49) ˙ v j = u ( x j ( t )) − v j ( t ) τ D − τ F max (cid:18) , h f h (cid:19) max (cid:18) , − hh f (cid:19) β v j + λ v j ( t ) v j +1 ( t ) − v j ( t ) x j +1 ( t ) − x j ( t ) . (50) This discrete model allow us to derive a continuous model by a scaling process in both space and time.Let us introduce the density of debris denoted ρ . The local density can be linked to inter-debris distanceaccording to the formula ( ρ d )( x j ( t )) = ρ (cid:96)x j +1 ( t ) − x j ( t ) (51)where (cid:96) is a characteristic length of debris and ρ is the maximum debris density (without overlapping).In particular, from equation (49) we have ddt ( x j +1 − x j ) = v j +1 − v j x j +1 − x j ( x j +1 − x j ) . (52)Applying the scaling to the expression (52) leads to a continuous medium equation of conservation ofthe number of debris D (1 /ρ ) Dt = 1 ρ ∂v∂x (53)where ρ and v are now functions of both space and time and D t = ∂ t + v∂ x is the Lagrangian derivative.Equation (53) can be written in conservative form as ∂ t ρ + ∂ x ( ρv ) = 0 . (54)The scaling applied to the expression (50) leads to the following partial differential momentum-likeequation ρD t v − λρv ∂v∂x = u − vτ D − τ F max (cid:18) , h f f (cid:19) max (cid:18) , − hh f (cid:19) β v, (55)20r in Eulerian description ∂ t ( ρv ) + ∂ x ( ρv ) − λρv ∂v∂x = u − vτ D − τ F max (cid:18) , h f f (cid:19) max (cid:18) , − hh f (cid:19) β v. (56)In all what follows, for simplicity we will denote by S = S ( h, u, v ) the right hand side of equation (56)and will be referred to as the source term. For smooth solutions, equation (56) can be rewritten under theconservative form ∂ t v + ∂ x (cid:18) (1 − λ ) v (cid:19) = Sρ (57)that resembles a Burgers-like equation. In this section, we discuss the important role of the interaction term I = − λρv ∂v∂x . to get a well-posed mathematical problem. For simplicity purpose let us consider the particular case S =0 (no drag and no friction forces). The system of partial differential equation reads ∂ t ρ + ∂ x ( ρv ) = 0 ,∂ t v + (1 − λ ) ∂ x ( v /
2) = 0 . As soon as λ (cid:54) = 1 and v (cid:54) = 0 , the system is clearly hyperbolic with two distinct eigenvalues λ = (1 − λ ) v and λ = v .Let us consider now the case λ = 0 (non-existent interaction term). Then the system is ∂ t ρ + ∂ x ( ρv ) = 0 ,∂ t ( ρv ) + ∂ x ( ρv ) = 0 . It is strictly equivalent to the so-called “pressureless Euler equations” (see for example Bouchut [35]).This system is not hyperbolic due to the unique eigenvalue λ = v of multiplicity 2. Actually this systemcan have weak measure solutions with the appearance of what is called delta-shocks that are nothing elsebut Dirac measures. Indeed, the second equation can be rewritten in v -variable as an inviscid Burgersequations, which is autonomous and can develop discontinuities in velocity. At discontinuity lines, themass conservation equation says that there is mass concentration. From the point of view of debrismodel, this would represents Dirac concentration of debris, what is not really realistic from a physicalpoint of view. This model would be the lowest-regularity description of debris concentration.In this sense, the interaction term acts as a regularization term that “smooths” debris concentrationwaves. Notice that in the particular case λ = 1 , the interaction term kills the convective term and we getthe interesting (conservative) simple model ∂ t ρ + ∂ x ( ρv ) = 0 ,∂ t v = Sρ . ∂ t ( ρv ) + v ∂ x ( ρv ) = S, so that we have the balance law in conservation form ∂ t ( ρ v ) + ∂ x ( ρ v ) = S. For two-dimensional problems, the models of drag and friction effects are kept unchanged, except thatthey are now vector-valued : S = u − v τ D − τ F max (cid:18) , h f f (cid:19) max (cid:18) , − hh f (cid:19) β v . For the modeling of debris interaction, the one-dimensional debris-following approach is no more avail-able. We simply empirically extend the formula found in the 1D case, and propose to use I = − λ ρ ( ∇ · v ) v . (58)This empirical extension can be justified by the fact that ∇ · v represent a compressible factor that is pos-itive for expansion conditions and negative for compressive configuration. By this way, the regularizinginteraction term acts in the opposite direction of sgn ( ∇ · v ) v .We get the following system of PDEs: ∂ t ρ + ∇ · ( ρ v ) = 0 , (59) ∂ t ( ρ v ) + ∇ · ( ρ v ⊗ v ) − λ ρ ( ∇ · v ) v = S . (60)For smooth solutions, the second equation can be written in velocity variable as ∂ t v + v · ∇ v − λ ( ∇ · v ) v = S ρ . Unfortunately, unlike the one-dimensional case, the system cannot be written in conservative form, ex-cept in the case λ = 1 . For λ = 1 , it is easy to check that the momentum variable ρ v is solution of thesimple transport-reaction equation ∂ t ( ρ v ) + v · ∇ ( ρ v ) = S . (61)Combining it with the continuity equation, we get the balance law in conservative form ∂ t ( ρ v ) + ∇ · ( ρ v ⊗ v ) = ρ S . (62)Denoting v = ( v, w ) , the homogeneous part of the system can be written in quasi-linear form ∂ t v + A x ∂ x v + A y ∂ y v = 0 (63)22ith A x = (cid:32) − w v (cid:33) , A y = (cid:32) w − v (cid:33) . For any unit vector ν = ( ν x , ν y ) , we have A ν := A x ν x + A y ν y = (cid:32) wν y − vν y − wν x vν x (cid:33) . It is clear that tr ( A ν ) = v · ν and det ( A ν ) = 0 so that the eigenvalues of A ν are and v · ν . Accumulation of debris can create debris hills and thus act on the global flooding flow. In that case,we have a two-way shallow water equations-debris dynamics coupling. One can add to the groundtopography z a rising due to the debris, which is proportional to the debris density (with factor µ < inthe next equations). Thus, the coupled system is ∂ t h + ∇ · ( h u ) = 0 , (64) ∂ t ( h u ) + ∇ · ( h u ⊗ u ) + ∇ · (cid:18) g h (cid:19) = − g h ∇ ( z + µρ ) , (65) ∂ t ρ + ∇ · ( ρ v ) = 0 , (66) ∂ t ( ρ v ) + ∇ · ( ρ v ⊗ v ) − λ ρ ( ∇ · v ) v = S . (67) Visualization of debris dynamics added to the water flooding may be a complementary tool for decisionsupport and risk analysis. One can visualize the debris density field but it is not so easy to merge bothwater surface and debris density. Another practical and humanly comprehensible way is to visualizedebris by debris particles themselves. In some sense we may go back to a Lagrangian description fromthe Eulerian computations, for visualization purposes only. While both ρ and v are computed. We definea set of debris particles can can be initially instantiated according to the initial debris density. Then wecompute the trajectories of the set of debris particles ˙ x j = v ( x j ( t ) , t ) , t > , (68) x j (0) = x j . (69)From a user perspective, this will give him an overview of debris dynamics, zones of concentration oraccumulation, particular pathways, assessment of human safety and risks for infrastructures. Damages are caused by stresses and forces exerted on structures by both water and debris. Consideringdebris, the accumulated momentum of debris at a given point point x and at time T is D ( x, T ) = (cid:90) T ρ v ( x, t ) dx. (70)23f course, debris density ρ ( ., T ) itself is another indicator of possible damage. The case for example ofdensity of debris releases into the sea is a good indicator of risk assessment for future coastal shippingactivities. As a pioneer work, we decide to use the simplest debris model, i.e. considering λ = 1 without two-waycoupling. For numerical expectations, we rather use the conservative formulation of debris model, thatwe recall here again: ∂ t ρ + ∇ · ( ρ v ) = 0 ,∂ t ( ρ v ) + ∇ · ( ρ v ⊗ v ) = ρ S . For numerical discretization, one can use a Lagrange-flux scheme to solve the convective part of thesystem. Let us now focus on the numerical treatment of the source term in the system. At each timestep ∆ t n and each grid point x , we have to solve the differential problem d v dt = u − v τ D − τ F max (cid:18) , h f h (cid:19) max (cid:18) , − hh f (cid:19) β v , (71) v ( t n ) = v n . (72)For envisaged applications, both relaxation times τ D and τ F can be rather small compared to the timestep, leading to stiff source terms. For that reason, it appears important to use at least semi-implicit timeintegration schemes.By freezing up exogenous terms with values taken at time t n , we have the following linear differentialequation d v dt = u n τ D − (cid:34) τ D + 1 τ F max (cid:18) , h f h n (cid:19) max (cid:18) , − h n h f (cid:19) β (cid:35) v . We get the solution at time t n +1 = t n + ∆ t n v n +1 = v (∆ t ) = v n e − η n ∆ t n + (cid:0) − e − η n ∆ t n (cid:1) u n τ D η n (73)with τ n = 1 τ D + 1 τ F max (cid:18) , h f h n (cid:19) max (cid:18) , − h n h f (cid:19) β . To show the effectiveness of the numerical model, we consider a two-dimensional case defined on therectangle spatial domain [0 , × [0 , . The topography/bathymetry is that one of the figure 4 below.The tidal wave is generated by an artificial discontinuous water elevation. The flow dynamics is made acoastal tidal wave that submerges all the surface, then a pull-back wave starts and goes back to the sea.The three obstacles also create backward waves. For the debris, we initially define a nonzero uniformdebris density ρ = 2 in the sub-zone [0 . , × [0 , .24igure 4: Configuration at initial time with topography/bathymetry (in brown color) and water depth(green-blue color).Numerical results showing the sequence of the flow dynamics is given in figure 5. One can observefirst the compression of the debris caused by the tidal wave, then the transportation of the debris that eitherland on the borders of the three hills or concentrate between the hills. We also compute a cumulativedamage function as the time integral of the norm of the debris momentum: D ( ., t ) = (cid:90) t ρ | v | ( ., s ) ds. One can observe that most damages are located as intuitively guessed between the obstacles where theflow rate is maximal. This promising qualitative result let us think that our debris model is appropriatefor damage assessment caused by debris. Of course, for a quantitative computation, we would need tocalibrate debris flotation heights and drag coefficients from real measurements or laboratory experiments.
The next step is to introduce this debris model into a high-performance code like GPU-VOLNA [3,4], then use a computer design of experiment (DoCE) with damage analysis, sensitivity analysis anduncertainty quantification.A way to reduced data dimensionality is to consider reduced-order models (ROM) for spatially-defined damage functions like those proposed in the document. One can notice from the numerical resultsthat that damage function defined by the time integral of debris momentum is a rather smooth functionin space so that model-order reduced can be envisaged. This approach will simplify the generation ofemulators for damage functions. 25e will also try to propose original numerical methodologies but also think about suitable parallelsoftware environments to process the data and extract knowledge in this risk assessment framework.
Acknowledgments
I would like to warmly thank Professor Serge Guillas from the Department of Statistical Science of UCLwho supported my candidature as UCL Big Data Institute invited Researcher. This work is also partlysupported and granted by the french CNRS multidisciplinary program “Défi Littoral” 2015-2016.