N Chevaugeon

J Xin

P Hu

X Li

D Cler

J E Flaherty

M S Shephard

Discontinuous Galerkin Methods Applied to Shock and Blast Problems

Keywords: Discontinuous Galerkin methods, hyperbolic conservation problems, blast, shock detection, anisotropic mesh refinement, adaptive methods

We describe procedures to model transient shock interaction problems using discontinuous Galerkin methods to solve the compressible Euler equations. The problems are motivated by blast flows surrounding cannons with perforated muzzle brakes. The goal is to predict shock strengths and blast over pressure. This application illustrates several computational difficulties. The software must handle complex geometries. The problems feature strong interacting shocks, with pressure ratios on the order of 1000 as well as weaker precursor shocks traveling rearward that also must be accurately captured. These aspects are addressed using anisotropic mesh adaptation. A shock detector is used to control the adaptation and limiting. We also describe procedures to track projectile motion in the flow by a level-set procedure.

INTRODUCTION

We propose to solve blast problems motivated by flows surrounding cannons with perforated muzzle brakes using a discontinuous Galerkin method (DGM), as described in Sec. 2. Detection of shock structures appears crucial to drive both adaptive mesh refinement and limiting and we describe such detection scheme (Sec. 3) and present some one-dimensional examples. We also describe an anisotropic adaptive mesh refinement procedure (Sec. 4), that allows us to solve large problems, in an efficient manner. We apply this to two examples involving blast simulation. In Sec. 5, we describe progress to include the presence of solid objects moving in the flow field by using a level-set/ghost-fluid approach, and present some applications.

DISCONTINUOUS GALERKIN METHODS FOR EULER EQUATIONS

The discontinuous Galerkin method was introduced by Reed and Hill in 1973 [START_REF] Reed | Triangular Mesh Methods for the Neutron Transport Equation[END_REF] as a technique to solve neutron transport problems. Lesaint and Raviart [START_REF] Lesaint | On a finite element method for solving the neutron transport equation[END_REF] presented the first numerical analysis of the method for a linear advection equation. However, the technique lay dormant for several years and has only recently become popular as a method for solving fluid dynamics and electromagnetic problems [START_REF] Cockburn | Discontinuous Galerkin Methods[END_REF]. As with all mesh based procedures, DGM uses a double discretization. First, the physical domain Ω is discretized into a collection of N e elements

T e = N e e=1 e
(2.1) called a mesh. Then, the function space V (Ω) containing the solution of the problem is approximated on each element e of the mesh by a finite-dimensional space V e (T e ).

The accuracy of the double discretization depends on geometrical and functional discretization. Classical Galerkin finite element approximation uses conforming meshes where elements share only complete boundary segments. The space V e is also constrained to be a subspace of a continuous function space, for example, H 1 , with a basis that is typically associated with element vertices, edges, faces, or interiors. These simplify the imposition of the C 0 continuity requirements of H 1 but limit choices. The DGM allows more general mesh configurations and discontinuous bases that simplify both h-and p-refinement. For example, non-conforming meshes and arbitrary bases for functional approximation [START_REF] Warburton | Spectral basis functions for 2d hybrid hp elements[END_REF] may be used. Herein, we use a L 2 -orthogonal basis as a product of Joacobi polynomials [START_REF] Dubiner | Spectral methods on triangles and other domains[END_REF] that yields a diagonal mass matrix [START_REF] Remacle | An adaptive discontinuous Galerkin technique with an orthogonal basis applied to compressible flow problems[END_REF].

The DGM can easily support arbitrary orders of spatial discretization accuracy without the need to construct complex stencils for high-order reconstruction [START_REF] Remacle | An adaptive discontinuous Galerkin technique with an orthogonal basis applied to compressible flow problems[END_REF]. Indeed, the DGM stencil remains invariant for all polynomial degrees. This greatly simplifies parallel implementation. The DGM has additional flexibility in specifying fluxes across elements faces and permits the use of fluxes based on exact or approximate solutions of Riemann problems [START_REF] Cockburn | The Runge-Kutta local projection discontinuous Galerkin finite element method for the conservation laws IV: the multidimensional case[END_REF].

We apply the DGM to hyperbolic conservation laws (Sec. 2.1) using a spatial discretization of p-degree orthogonal polynomials. Time discretization is performed by an explicit total variation diminishing Runge-Kutta method [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]. To improve the performance of the explicit integration on irregular and unstructured grid, we use a new local time stepping procedure [START_REF] Flaherty | Adaptive local refinement with octree load-balancing for the parallel solution of threedimensional conservation laws[END_REF][START_REF] Remacle | An efficient local time stepping-discontinuous Galerkin scheme for adaptive transient computations[END_REF].

Discontinuous Finite Element Formulation for Conservation Laws

Consider an open set Ω ⊂ R 3 with boundary ∂Ω. We seek to determine u(Ω, t): R 3 × R → L 2 (Ω, t) m = V (Ω, t) as the solution of a system of m conservation laws

∂ t u + div F(u) = r.
(2.2)

Here div = (div , . . . , div ) is the vector valued divergence operator and

F(u) = ( F 1 (u), . . . , F m (u)) (2.3)
is the flux vector with the ith component

F i (u): (H 1 (Ω)) m → H(div , Ω),
where H(div , Ω) consists of square integrable vector valued functions whose divergence is also square integrable, that is,

H(div , Ω) = v | v ∈ L 2 (Ω 3 ) , div v ∈ L 2 (Ω) . ( 2
(∂ t u, w) Ω -( F(u), grad w) Ω + F(u) • n, w ∂Ω = (r, w) Ω , ∀w ∈ V (Ω).
(2. F n (u e , u e k ), w ∂e k = (r, w) e , ∀w ∈ V e (e), (2.10) where n e is the number of faces of element e. Only the normal traces have to be defined on ∂e k and several options are possible (see, e.g., Refs. [START_REF] Van Leer | Flux Vector Splitting for the Euler Equations[END_REF][START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]). It is usual to define the trace as the solution of a Riemann problem across ∂e k . Herein, when we consider problems with strong shocks [START_REF] Colella | Efficient solution algorithms for the Riemann problem for real gases[END_REF][START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF], an exact Riemann solver is used to compute the numerical fluxes and a moment limiter [START_REF] Biswas | Parallel adaptive finite element method for conservation laws[END_REF] is used to restrict spurious oscillations when polynomial degrees p > 0 are used. The moment limiter is applied on each element after a new solution field is computed, so that the minimal and maximal value of the solution moments in one element are bounded by their mean values on neighboring elements.

Euler's Equations

Flows satisfy the compressible Euler equations, which, in two dimensions, have the form (2.2) with

u = [ρ, ρu, ρv, E] T (2.11)
where ρ is the density of the fluid, v is the velocity with Cartesian components u and v and E is the internal energy. The flux is where P is the pressure and e x and e y are unit vectors in the x and y coordinate directions. An equation of state (EOS) P = P (ρ, E) is necessary to close the system. The DGM software can be used with any physically consistent EOS and we choose the perfect gas law.

F(u) = F 1 , F 2 , F 3 , F 4 = ρ v,
P = (γ -1)ρ E - v 2 2 (2.13)
with gas constant γ = 1.4.

DISCONTINUITY DETECTION

A discontinuity detection procedure automatically detects regions of the flow field where the solution is discontinuous. For the Euler equations, discontinuities may be contact discontinuities, where the density is discontinuous; shock waves, where the density, velocity, and pressure are discontinuous; and expansion wave where flow gradients are discontinuous. A knowledge of discontinuity locations allows us to restrict the limiting procedure to the vicinity of discontinuities, which both saves time and improves solution accuracy in smooth flow regions. Discontinuity detection is also used to drive the anisotropic adaptation procedure (Sec. 4).

To define a discontinuity detection strategy, we partition an element boundary ∂e into portions ∂e -and ∂e + where, respectively, the flow is into ( v • n < 0) and out of ( v • n > 0) element e.

Let q be a flow variable or a derived quantity from the flow field, for example, density or entropy, for the Euler equations, and Q e be the DG value of q on element e. Consider the jump in Q e across the inflow edges (faces) of e and examine

I e = ∂e - (Q e -Q e k ) dσ = ∂e - (Q e -q) dσ + n e k=1 ∂e + k (q -Q e k ) dσ. (3.1)
In smooth flow regions, the first integral on the right is O(h p+2 ), while the second might be O(h 2(p+1) ) [START_REF] Adjerid | A posteriori error estimation for discontinuous Galerkin solutions of hyperbolic problems[END_REF][START_REF] Krivodonova | Error estimation for discontinuous Galerkin solutions of two-dimensional hyperbolic problems[END_REF][START_REF] Krivodonova | Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws[END_REF] due to super-convergence. Here, h is a measure of an element size, as defined below. Thus, in smooth regions,

I e = O(h p+2 ).
If q is discontinuous in the immediate vicinity of ∂e, then either or both of q -Q e and q -Q e k are O(1); hence,

I e = O(h p+2 ), if q| ∂e is smooth, O(h), if q| ∂e is discontinuous. (3.2)
We construct the detector by normalizing I e relative to an "average" O(h (p+1)/2 ) convergence rate and the solution on e to obtain

I e = | ∂e -(Q e -Q e k ) dσ | h (p+1)/2 ∂e -Q e . ( 3.3)
In examples, we choose h as the radius of the circumscribed circle in element e, and use a maximum norm based on local solution maxima at integration points in two dimensions and an element average in one dimension. Using (3.2), I e → 0 as either h → 0 or p → ∞ in smooth solution regions, whereas I e → ∞ near discontinuities. Thus, a detection scheme is

if I e > 1, q is discontinuous, if I e < 1, q is smooth. (3.4)
Example 3.1. We apply the discontinuity detection scheme [Eqs. (3.3) and (3.4)] to a blast problem [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF] for the one-dimensional Euler Eq. [[(2.2), (2.11)-(2.13)] with v = 0] subject to the initial data (ρ, u, P )(x, 0) =    (1, 0, 1000), 0 x < 0.1, (1, 0, 0.01), 0.1 x < 0.9, (1, 0, 100), 0.9 x 1.

(3.5)

The solution involves multiple interactions of strong shocks, expansions, and contact discontinuities. The most important interactions take place in a short period prior to t = 0.04. An exact solution of this problem is not known, so, following Woodward and Colella [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF], comparisons were done relative to a numerical solution computed with N e = 3072 uniform elements and piecewise polynomials of degree p = 2. We present solutions for the density (Fig. 1) with N e = 800, p = 2, using entropy as the quantity q for the shock detection scheme and moment limiting [START_REF] Biswas | Parallel adaptive finite element method for conservation laws[END_REF] in regions where discontinuities were detected.

We see that the shocks have been identified and limiting is confined to these regions. Comparing the numerical and "exact" solutions, we see that all the important features of the various interactions have been sharply captured.

MESH ADAPTATION

Our software is able to perform hierarchical mesh refinement where elements are refined using an isotropic refinement template. Since interelement continuity constraints are not needed with the DGM, element refinement need not consider neighboring cells. The level of mesh refinement of a cell is computed using an error indicator [START_REF] Remacle | An adaptive discontinuous Galerkin technique with an orthogonal basis applied to compressible flow problems[END_REF] based on the value of the jump in solutions between neighboring cells. Coarsening is equally simple, by reversing refinement after projecting the solution from a finer mesh to the coarser parent mesh. From left to right: initial mesh, non-conforming isotropic mesh at t = 0.6, enlargement of this mesh near the discontinuity, anisotropic conforming mesh at t = 0.6, and enlargement of this mesh near the discontinuity.

ρ = 1, if 0 x < 0.8, 2, if 0.8 x 1. (4.1)
This initial field is advected by a unit velocity in the -x direction. The boundary conditions are reflective at y = 0, 1, and prescribed at x = 1. With six levels of adaptive mesh refinement, the minimum edge has a size of order h 0 = 0.1/(2 6 ) = 0.0015625. The solution has the same accuracy as one on a uniform mesh of 819200 elements in two dimensions. Figure 2 displays the initial mesh and the adaptive mesh at t = 0.6 with six levels of mesh adaptation. The adapted mesh has 6097 elements, which is 134 times less than needed by a uniform mesh with the same resolution. This gain in degrees of freedom has to be weighted against the cost of adapting the mesh with the evolving solution.

Isotropic mesh refinement does not account for situations, such as the one at hand, where refinement is only needed in one direction. Thus, an anisotropic mesh refinement scheme has the potential to greatly improve performance [START_REF] Frey | Anisotropic mesh adaptation in 3D. In 12th International Meshing Roundtable[END_REF].

Our anisotropic mesh refinement procedure adapts a conforming mesh using local mesh modification operations such as element subdivision, edge collapsing, and edge and face swapping. A complete description of this mesh adaptation procedure appears in Refs. [START_REF] Li | Mesh Modification Procedures for General 3-D Non-Manifold Domains[END_REF][START_REF] Li | 3D anisotropic mesh adaptation by mesh modifications[END_REF].

The anisotropic mesh adaptation procedure requires the definition of an anisotropic mesh size field. The natural way to represent an anisotropic size field is to define a symmetric positive definite tensor field or a metric. Given any point x in the size field domain, and a direction a, the size field describes the size l of a mesh edge in the given direction. The mesh adaptation procedure uses a definition of the metric at each vertex and assumes a linear variation of the size field along element edges.

Building an anisotropic size field requires directional error information. For piecewise linear (p = 1) approximations, this may be based on the Hessian of the field. The principal directions of the Hessian are then the principal direction of the metric field. The eigenvalues of the Hessian field are used to define the eigenvalues of the metric field by using a normalization,obtained by selecting a small mesh size in directions where the Hessian is large.

We compute an approximation of the Hessian of the field using a procedure of Remacle et al. [START_REF] Remacle | Transient mesh adaptation using conforming and non-conforming mesh modifications[END_REF], where a stabilized gradient of the field is first computed by solving a DG problem on each element. A stabilized Hessian is obtained from this gradient by solving a second local DG problem. Near discontinuities, the Hessian of the exact solution does not exist. Thus, using an approximate Hessian near discontinuities will provide unreliable results. Our procedure uses the shock detector (Sec. 3) to decide if the Hessian should be used or not. If the Hessian cannot be used, we use the gradient of the field on each side of the shock to determine principal directions of the size field. A user defined minimal size field and a user defined maximum aspect ratio are used to specify mesh spacing in the principal direction of the size field.

The problem of Example 4.1 is solved using the anisotropic procedure. For consistency with the prior results we set h min = h 0 and select an aspect ratio α = 100. The mesh obtained at time t = 0.6 with this setting is displayed in Fig. 2. This mesh only has 100 elements, and the results are as accurate as the solution with isotropic refinement. Element edges are aligned with the discontinuity, which makes the computation of the numerical flux across edges more accurate. In realistic situations, discontinuity surfaces are not planar and the aspect ratio that the mesh adaptation procedure will be able to reach for a fixed minimum size depends on the curvature of the discontinuity surface. Nevertheless, the savings in the number of degrees of freedom makes simulation of problems involving complex shock structures less expensive.

Example 4.2. The mesh adaptation procedure is used to solve a simplified two-dimensional blast problem. The domain under consideration (Fig. 3) is a slice of the three-dimensional problem with a tube having a diameter of 7.65 mm, and a length of 30 mm. The left side of the domain represents an axis of symmetry. On the tube walls we prescribe slip boundary conditions. The edge inside the tube at the bottom left is the inlet, and other edges are subsonic outflow surfaces. At t = 0, the external pressure is atmospheric at 101325 Pa, the velocity is 0 m/s and temperature is 300 K. At the inflow boundary, we impose a pressure of 6 × 10 6 Pa, a temperature of 631.626 K, and a velocity of 631 m/s, which correspond to a typical precursor shock in a blast problem [START_REF] Cler | CFD application to gun muzzle blast-a validation case study[END_REF][START_REF] Klingenberg | Gun Muzzle Blast and Flash[END_REF]. Time steps are adaptively set to maintain a maximum CFL limit of 0.9, and correspond to a range of 1 × 10 -8 -5 × 10 -8 s. Mesh Adaptation is done every 1 × 10 -6 s and the problem is solved to for 1 × 10 -4 s. We set the minimum mesh size for adaptation to 0.1 mm and the maximum aspect ratio to 50. In the first phase of the simulation, the flow is mostly one-dimensional while the shock structure progresses in the tube. Thus, like a shock tube, a shock is followed by a contact discontinuity which is followed by an expansion wave traveling toward the open end of the tube. The adaptive mesh is concentrated at the contact discontinuity and the shock is aligned with the shock plane. Even when the two discontinuities are close, the mesh adaptation procedure is precise enough to separate the two structures with some coarsening between them. At time t = 6 × 10 -5 , the shock has exited the tube and forms a characteristic bubble shape. The contact is just at the edge of the tube. At time t = 8 × 10 -5 , the structure is fully developed, with a main "spherical" shock, followed by the contact, and a shear layer linking the contact with the exit of the tube. The final mesh contains 50, 000 degrees of freedom.

Example 4.3. We consider a more complex tube with a muzzle brake. The tube has a diameter of 155 mm and a set of 16 perforations of 28.6 mm diameter near the muzzle. We model a slice of it with a cut plane through the center of the holes. A shock-tube type initial condition, with a shock placed 75 mm to the left of the perforations is used to initiate the flow. The upstream pressure is 565 times the downstream pressure and the upstream temperature is 2111.5 K. Fig. 4 shows the density and meshes at time t = 1 × 10 -4 , 5 × 10 -4 and 1 × 10 -3 s. At time t = 1 × 10 -4 s, the initial shock has reached the fourth hole, and "mushroom-like" structures are forming outside the first and second holes. At time 5 × 10 -4 s, the initial shock is outside the tube, forming its own plume. Inside the tube, shear layers from each hole are propagating downstream. A main "pear-shaped" shock structure, enclosing the complex "diamond-shaped" structures resulting from the interacting shocks, is forming and moving away from the muzzle. At time t = 2 × 10 -3 s, the initial shock has exited the muzzle and formed a "spherical" structure that merges with the pear-shaped one. The mesh is closely following the evolution of the shock structures and is aligned with them. Even though the shocks are of very different strengths, they are captured by the size field computation. The mesh at t = 1 × 10 -3 s after 1001 mesh adaptations contains 900,000 degrees of freedom.

MOVING OBJECT

The next step is to include moving objects such as projectiles in the flow field. The basic approaches to handle such calculations involve (i) using deformable meshes that move with the interface [START_REF] Stein | Automatic mesh update with the solidentention mesh moving technique[END_REF][START_REF] Tezduyar | Finite element methods for flow problems with moving boundaries and interfaces[END_REF] and (ii) using a fixed mesh and a technique to capture the interface, such as a volume of fluid [START_REF] Lafauric | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF], front tracking [START_REF] Mulder | Computing interface motion in compressible gas dynamics[END_REF], or level-set [START_REF] Chang | A level set formulation of Eulerian interface capturing methods for incompressible fluid flows[END_REF][START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF]. Herein, we describe a level-set technique [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] to track the interface between the rigid body and the compressible fluid. We adapt the Ghost Fluid Method (GFM) [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] to capture the correct boundary conditions at the interface and enable us to perform flow calculations over the whole domain, without explicit regard of the object.

The level-set function φ( x, t) is defined as a signed distance function from a point x to the boundary of the moving object at time t. If we denote the boundary of the object as B, the inside of the object as I , and the fluid part of the domain as O, then φ(x, t) satisfies

φ( x, t) < 0, if x ∈ I, φ( x, t) = 0, if x ∈ B, φ( x, t) > 0, if x ∈ O.
(5.1)

The normal vector to the moving object is

n( x, t) = φ( x, t), x ∈ F, ( 5.2)
and is extended to the entire domain (I ∪ B ∪ O). This normal vector is a constant on any trajectory perpendicular to a surface with a constant value of φ. In the present work, we assume that φ( x, t) can be explicitly defined, which is reasonable for the rigid objects under consideration. Treatment of deformable objects will need an actual discretization and computation of the change of the shape of φ. The development of the following "boundary" conditions implied by the moving object would, however, be the same.

Elements inside the rigid object (φ < 0) are modeled as "ghost fluid" elements. The Euler equations are solved using the DG method just as they are in the actual fluid. The field variables in ghost elements are set to reproduce the behavior at a solid-fluid interface.

Mesh entities that lie completely in the fluid (φ > 0) require no special attention. Entities that are either in the object (φ < 0) or span the interface φ = 0 require a proper setting of fluid variables in the "ghost" (interior) part of the object to capture the correct conditions at the moving interface. The moving interface is treated as a contact discontinuity since the fluid near the interface moves with the speed of the object. Thus, the pressure and normal velocity v n = v • n are continuous at the interface and we prescribe ghost fluid values as a reflection about φ = 0 surface to enforce these conditions [START_REF] Fedkiw | Coupling an Eulerian fluid calculation to a Lagrangian solid calculation with the ghost fluid method[END_REF][START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF]. Consider any ghost point x G in the object and its reflected image x R across the interface. The density ρ G , energy E G , and velocity v G at x G , in terms of similar quantities at x R are

ρ G = ρ R , (5.3) v G = v t R -v n R + 2 U n G , ( 5.4
)

E G = E R -1 2 ρ R | v R | 2 + 1 2 ρ R | v G | 2 , ( 5.5)
where v t R and v n R are the tangential and normal velocity components and U G is the velocity of the moving rigid body at x G . For a given point x G , the closest point x I on the interface can be approximated as

x I = x G -φ( x G ) φ( x G ).
(5.6) Thus, the reflected point is

x R = x G -2φ( x G ) φ( x G ).
(5.7)

In the practice, only those ghost entities within a small distance to the interface need be considered. The bandwidth can be set to about three to four times the element cell size near the interface. Thus, the level-set function need only be defined in a small neighborhood of the fluid-solid interface. Thus, the computation is efficient even when the object has a large volume relative to the element size. Density contours are shown at t = 0.1, 0.25, and 1.6 s in Fig. 5. These results compare well with the flow about a fixed ball in a channel with the fluid having a speed corresponding to the sum of the speed of the cylinder and the speed of the fluid.

Example 5.2. We solve Example 4.2, with the flow driven by a square projectile having an initial speed of Mach 0.9. Preliminary results, presented in Fig. 6, show the pressure field at stages when the projectile is inside the tube and when it has exited. These preliminary results are encouraging. We reproduce the strong precursor shock. Had we applied a moving-mesh to explicitly track the projectile, difficulties arise when the projectile exits the tube. With the presented method, no special difficulties arise in this case.

CONCLUSION

We described procedures to model transient shock problems applied to blast flow using DGMs and mesh adaptation. The DGM enables us to sharply resolve shocks and other discontinuities. Error estimation results related to DG theory [START_REF] Adjerid | A posteriori error estimation for discontinuous Galerkin solutions of hyperbolic problems[END_REF][START_REF] Krivodonova | Error estimation for discontinuous Galerkin solutions of two-dimensional hyperbolic problems[END_REF] allow us to build robust discontinuity indicators that efficiently drive the limiting and the mesh adaptation procedures. The level-set method/GFM introduced to account for moving objects is an interesting alternative to the more complex procedure of following objects with a moving mesh. More work is necessary to assess the method's accuracy, stability, and efficiency. Future work includes comparisons with a moving-mesh approach and extensions to three-dimensional problems. Encouraging results in this direction are shown in Fig. 7. Two three-dimensional idealized blast problems were solved in a parallel computational environment. The results on the left of Fig. 7 were obtained on a four-processor system with a final mesh of about two million degrees of freedom. The idealized tube has a square section. The result on the right of Fig. 7 is an idealized tube with a hexagonal cross section and with holes in the tube wall. This computation was done on an eight-processor system with a final mesh of about eight million degrees of freedom.

. 4 ) 2 )

 42 With the aim of constructing a Galerkin form of Eq. (2.2), let (•, •) Ω and •, • ∂Ω , respectively, denote the standard L 2 (Ω) and L 2 (∂Ω) scalar products (u, w) Ω = by a test function w ∈ V (Ω), integrate over Ω and use the divergence theorem to obtain the variational formulation

Fig. 1 .

 1 Fig. 1. Density for Example 3.1 at t = 0.016 (top left), 0.026 (top right), 0.028 (bottom left), and 0.038 (bottom right) for computations performed with p = 2, N e = 800. Elements where discontinuity is detected and moment limiting is applied are marked with a + sign at unity.

Example 4 . 1 .

 41 To illustrate this process, consider a simple scalar advection equation in two dimensions on a unit square. The square is partitioned by an unstructured mesh of triangles with edge length of the order of approximately 0.1. The initial density is

Fig. 2 .

 2 Fig. 2. Meshes for Example 4.1.From left to right: initial mesh, non-conforming isotropic mesh at t = 0.6, enlargement of this mesh near the discontinuity, anisotropic conforming mesh at t = 0.6, and enlargement of this mesh near the discontinuity.

Fig. 3 .

 3 Fig. 3. Density profile (left) near the tube and adaptive mesh (right) for Example 4.2 at t = 2 × 10 -5 , 6× 10 -5 , and 8 × 10 -5 s (left to right).

Fig. 4 .

 4 Fig. 4. Density (left) and adaptive mesh (right) for Example 4.3. at t = 1 × 10 -4 , 5 × 10 -4 , and 1 × 10 -3 s (top to bottom).

Example 5 . 1 .

 51 Consider a two-dimensional Mach 3 flow in a 1 × 4 channel. Initially the channel contains a diatomic gas having a unit density and pressure. A gas having this density and pressure enters from the left boundary at Mach 3. A circular cylinder with radius 0.125 moves from right to left at Mach 3. Reflecting boundary conditions are applied to solid channel walls.

Fig. 5 .

 5 Fig. 5. Density field at t = 0.1 (upper left), 0.25 (upper right), 0.5 (lower left), and 1.6 s (lower right).

Fig. 6 .

 6 Fig. 6. Pressure field for Example 5.2 at t = 6 × 10 -5 s when the projectile is inside the tube (left) and at t = 2 × 10 -4 s when the projectile has left the tube (right).

Fig. 7 .

 7 Fig. 7. Three-dimensional results: density iso-surface. Left idealized square section cannon, right, idealized hexagonal section cannon with muzzle break.

u) = F(u) •

  • n, where n is the unit outer normal to ∂e, is not defined on ∂e. In this situation, a "numerical flux" F n (u e , u e k ) is used on each portion ∂ e k of ∂e shared by element e and neighboring element e k . The numerical flux must be continuous, so F n ∈ H(div , Ω) m , and consistent, so F n (u,

	With the DGM, V e is a "broken" function space that consists of the direct
	sum of elementary approximations u e
	V e (T e ) = u | u ∈ L 2 (Ω) m , restriction (u) | e
	= u e ∈ P p (e) m = V e (e), ∀e ∈ T e ,	(2.8)
	where P p (e) m is a polynomial of degree p on element e. Because all
	approximations are disconnected, we can solve on each element to obtain
	(∂ t u e , w) e -( F(u e ), grad w) e + F n , w ∂e = (r, w) e , ∀w ∈ V e (e). (2.9)
	A discontinuous basis implies that the normal flux F n = F(u) n. With such
	a numerical flux, Eq. (2.9) becomes
	(∂ t
		7)

u e , w) e -( F(u e ), grad w) e

	n e
	+
	k=1

  ρu v + P e x , ρv v + P e y , (ρE + P ) v(2.12)

Scientific Computation Research Center, Rensselaer Polytechnic Institute, Troy, NY 12180, USA. E-mail: flahej@rpi.edu

Benet Laboratories, TACOM/ARDEC, U. S. Army, Watervliet, NY 12189, USA.


Mohamed Nassar
email: mohamed.nassar@loria.fr

Saverio Niccolini
email: saverio.niccolini@netlab.nec.de

Radu State
email: radu.state@loria.fr

Thilo Ewald
email: thilo.ewald@netlab.nec.de

Holistic VoIP Intrusion Detection and Prevention System

VoIP security is crucial for current and future networks and services. The rapid shift from a closed and confined telephony towards an all IP network supporting end to end VoIP services provides major challenges to the security plane. Faced with multiple attack vectors, new and comprehensive defensive security solutions for VoIP must emerge from the research community.

This paper describes a multilayer intrusion detection and prevention system architecture for VoIP infrastructures. The key components of the approach are based on a VoIP-specific honeypot and on an application layer event correlation engine. While each component alone can detect only a subset of VoIP-specific attacks, the two of them together can provide an effective defense for the many class of attacks. We show in this paper, how different and complementary conceptual approaches can jointly provide an in depth defense for VoIP architectures.

INTRODUCTION

Securing VoIP infrastructures constitutes one of the major challenges for both the operational and research communities because security by design was not a key component in the early phases of VoIP research and development. VoIPspecific security solutions are currently demanded by the market while the research and standardization are still trying hard to address the issues of securing and monitoring VoIP infrastructures. Over the past few years, different approaches emerged, anyway most of them only address the defense against a subset of potential attack vectors.

We propose in this paper a holistic approach to VoIP security intrusion detection and prevention. Our approach is based on a combined VoIP-specific honeypot and application layer monitoring scheme based on SIP. Such an approach is capable of detecting multiple types of attacks: The VoIPspecific honeypot is best suited for preventing social attacks like Spam over Internet Telephony (SPIT) and VoIP Phishing (Vishing) as well as other stealthy reconnaissance actions, while the SIP correlation engine is adapted to detect Denial of service and/or fraudulent usage. Our paper is structured as follows: Section 2 starts with on overview of major threats that must be addressed by current and future VoIP infrastructures. The global architecture proposed in this paper is described in section 3, while the next two sections (4 and 5) describe in detail the two main security components of this architecture: a VoIP honeypot and respectively an anomaly based intrusion detection engine. An overview of related work is given in section 6. The paper ends with conclusions and pointers to future work in section 7.

VOIP THREATS AND PROBLEMS

VoIP security threats constitute a superset of those faced by data networks. The main source of security threats comes from the fact that both the signaling and the control plane in VoIP are carried over the IP network. Therefore, VoIP infrastructure shares the same vulnerabilities as the data networks. In addition, VoIP-specific threats (both at signaling and data layer) do also represent major causes of concern. Among the most dangerous VoIP-specific attacks, identity theft, eavesdropping, fraudlent usage and social attacks (SPIT, Vishing) are becoming reality. Denial of service (DoS) attacks can be oriented against the VoIP infrastructure (servers, proxies, agents) and lead to the crippling and total shut down of a VoIP infrastructure. Viruses, worms and backdoors can allow the remote control of IP phones and VoIP proxies allowing more than just simple malicious purposes. Password cracking bots can be launched and assure a fraudulent usage, privacy violation and SPIT/Vishing operations. In addition, with VoIP growing rapidly, hackers are becoming interested in extending and benefiting from this market. Designing a security solution for a distributed, multi-protocol and QoS-sensitive application as VoIP is a hard task. This paper addresses the issue of a global and holistic security solution capable to deal with the multiple threats faced by a VoIP infrastructure.

Interception and Modification Threats

In contrast to the difficulties encountered by the illegal interception in PSTN, a VoIP conversation reconstruction is possible using traffic captured and decoded. Free software like Vomit 1 (Voice over misconfigured internet telephones) are already existing for enterprise and convenient Cisco IP phones. Similar tools for SIP have also emerged [START_REF] Endler | Hacking VoIP Exposed[END_REF] and do not represent major technological difficulties. The attacker is left with a huge choice of actions ranging from injecting/eavesdropping and stealing sensitive financial and commercial information.

Denial of Service and Toll fraud

Flooding attacks can target the signaling plane elements (e.g. proxy, gateway, etc.) with the objective to take them down and produce havoc in the VoIP network. This is very easy to do by either flooding the signaling plane with a large quantity of messages, malformed messages or device specific vulnerabilities. Abuse of service attacks have as objective the fraudulent usage of the VoIP services. Threats that can be cited are unauthorized or unaccountable resource utilization that exploits specific vulnerabilities if the identity management of VoIP infrastructures aimed at fraudulous usage and/or spoofing call identification. For a comprehensive overview on the VoIP security threats, please check the taxonomy developed at [START_REF] Voipsa | VoIP security and privacy threat taxonomy[END_REF].

Social Threats

Social threats are attacks ranging from the generation of unsolicited communications -which are annoying and disturbing for the users-to more dangerous data stealing (Phishing) attacks. The threat is classified as social since the term "unsolicited" is strictly bound to user-specific preferences and this makes hard for system to identify this kind of attack. An example of this is a threat commonly referred to as Spam over Internet Telephony (SPIT) (similar to Spam in the email systems but delivered by mean of voice calls) which leverage on the cheap cost that VoIP has with respect of legacy phone systems (it is currently estimated that generating VoIP calls is three order of magnitude cheaper than generating PSTN calls). SPIT calls can be telemarketing calls used for guiding callees to a service deployed to sell products. A subtle variant of SPIT is the so-called Vishing (VoIP phishing) attack that aims either to make the callee dialing expensive numbers in order to get the promised prize or to collect personal data redirecting the users towards Interactive Voice Responder (IVR) pretending to be trusted. Most of these attacks are going to be generated by machines (bot-nets) programmed to do such a job. Unsolicited communications (like SPIT or Vishing) are, from a signaling point of view, technically correct transactions. It is not possible to distinguish from the INVITE (in the case of SIP) if such a transaction is SPIT or not. From a technical point of view the challenge is even more complicated since the 1 http://vomit.xtdnet.nl/ content is not available to help in the detection until the phone rings (disturbing the user) and the callee answers the call, for this reason techniques available from email spam like text filtering are hardly reusable. Even if a transaction is identified as unsolicited it depends strongly on the legal country environment how to handle it.

ARCHITECTURE FOR INTRUSION DE-TECTION AND PREVENTION

VoIP architectures are distributed (proxy servers, gateways, application servers, terminals, etc.) and thus difficult to tackle by a centralized security approach. The architecture of an intrusion detection and prevention solution should be distributed, unless the network is designed in a dedicated way. The architecture presented in this paper takes into account the possibility of using a single entry point but for the sake of generality it proposes a distributed approach combining a VoIP-specific Honeypot and application layer monitoring scheme based on SIP. Such an approach is capable of detecting multiple types of attacks basing on the combination of two solutions.

Figure 1 depicts the general architecture where the VoIPspecific Honeypot (see section 4) domain is separated by the real infrastructure domain. The distributed application layer monitoring is achieved deploying the VoIP Security Event Correlation (SEC) (see section 5) on VoIP infrastructural elements (SIP Proxy Servers, Terminals, etc.) both in the real domain and in the honeypot one.

In order to accommodate the VoIP-specific Honeypot into the intrusion detection and prevention architecture an additional domain was introduced. The task of this domain is to transparently bridge the internal domains (real and honeypot) and the external one correctly routing attacks to the honeypot domain and real requests to the real domain while hiding domain-specific information. In the architecture depicted in Figure 1 the domain is represented by an inbound / outbound proxy in a public reachable network.

The production domain (real) and the honeypot one are strictly disconnected from each other to prevent backdoors to the production environment when the VoIP-specific Honeypot is compromised. Thus a sharing information database and a black-listed users one are deployed in the inbound / outbound domain as information sharing mean between internal domains both domains. The scope of such databases is to have a sharing of information between the Honeypot domain and the real one to improve detection and prevention methods and schemes based on observations of the VoIPspecific Honeypot (see section 4.2).

VOIP-SPECIFIC HONEYPOT

In literature, a Honeypot is a trap set to detect, deflect and monitor attacks to information systems. Generally it consists of a computer, data or a network site that appears to be part of a network but which is actually isolated and monitored. A Honeypot has a specific value in attack detection and deflection. Honeypots usually are specific systems that normally should not see any legitimate traffic or activity. Whatever activity is seen on a Honeypot it can be interpreted as malicious or unauthorized.

The Honeypot concept can become very useful as specific component of a VoIP intrusion and prevention architecture. We brought the Honeypot concept into the VoIP world de- The key features of the VoIP-specific Honeypot are the mitigation of the SPIT threat with low cost infrastructure as described in 4.1 and the fact that such an infrastructure well complements all the other possible intrusion detection and prevention by information gathering useful to improve error rate of other methodologies as described in 4.2.

Architecture and Implementation

The VoIP-specific Honeypot depicted in Figure 2 is composed of standard open-source VoIP components based on the SIP protocol. The main goal of a Honeypot is to attract attacks into a secured and observed environment for analyzing their attacking schemes and deduce new kind of prevention methods against them. To reach this goal the Honeypot has to offer attractive services that are worth to be attacked. These services were simulated because of performance and resources issues. A SIP Honeypot has to simulate a whole SIP network, which can offer many fetching values: The openser advanced routing logic allows to configure complex routing schemes and its modular design can be extended with additional modules for more advanced services;

• SIP components • SIP services • SIP users
• Asterisk3 : an open source PBX. Asterisk offers a wide range of services and applications (e.g. gateway functionalities, mailboxes or Interactive Voice Response, IVR, applications) allowing a wide range of PBX scenarios emulation on a Honeypot.

The VoIP-specific Honeypot emulates a real VoIP network as follows:

• the SIP infrastructure is emulated by Openser. An instance of Openser is configured in the VoIP-specific Honeypot. Multiple instances can be used to build more complex SIP networks. Openser is configured with a random call-dispatcher for handling unknown callees. Callers initiating sessions to unknown callees gets randomly connected to the emulated services and users.

• the SIP services are emulated by Asterisk. With dedicated numbers (Asterisk extensions) the attacker gets redirected to voice-mail boxes (e.g. the SIP-PSTN gateway service is a simple prerecorded messages that signals the temporal unavailability of such a service.

• the SIP users are emulated by Asterisk mailboxes as well. Optional SIP phones can be implemented to mirror every incoming call for observation purposes.

Detection and Prevention Scheme

In order to stimulate attacks towards the Honeypot and to trap as much attackers as possible we prepared scripts that are going to publish systematically fake user URIs (in web pages, emails, newsgroups, etc.) registered to the Honeypot SIP proxy server in a way that is either not possible to be seen by humans (e.g. URIs written in white over a white background in a web-page) or easily recognizable by humans as URIs not corresponding to real users (e.g. URIs with the word do-not-call) but at the same time that can be easily harvested by bot-nets automatically looking for semantics of URIs. The objective is, with such a big number of user URIs advertised on the Internet, to make the VoIP-specific Honeypot and its users target of social attacks like SPIT and/or Vishing (since bot-nets that harvested the URIs are not going to make distinctions on them).

The rationale behind such VoIP deflectors (VoIP-specific Honeypot users) is that they should normally not experience traffic activity and therefore every activity seen by these deflectors can be interpreted as a malicious one. If any of these deflectors receives a call this is considered as an indication that the sender is an initiator of SPIT/Vishing attacks.

Information gathered by such a monitoring activity is then used as input to the prevention system reacting to attacks. A simple reaction is to add users that initiated activities to the VoIP-specific Honeypot to a black list that can be temporary or permanent. The users in such a list are not allowed to initiate communications to the real infrastructure. For this purpose we programmed the VoIP-specific Honeypot to publish its results in a database accessible for reading to the real infrastructure. VoIP-specific Honeypot federations sharing such monitoring activities among multiple site for the purpose of improved detection is currently subject of an on-going work. The idea in this case would be to have federations of domains that share the VoIP-specific Honeypot activity reports in order to build a list of untrusted users (users seen on multiple VoIP-specific Honeypots are more likely to be malicious and not to have made an error in dialing).

The big number of fake user URIs has also a second objective which is the deflection of attacks, as a matter of fact if the attacker harvests for SIP URIs (either from the web or doing a user enumeration attack) it will retrieve not only the real x ones but also the fake y ones. In this case the mitigation percentage will be proportional to the ratio between real and fake users y/(x+y)since the SPIT/Vishing calls will be distributed among all users (real and fake ones).

The VoIP-specific Honeypot has not only the objective of monitoring intrusions but has also to gather additional features of attacks characteristics for intrusion detection and prevention in the case of SPIT/Vishing attacks. The Honeypot software (mainly the terminal part) is composed of different answering machine that can:

• register the call and analyze it later: messages of the same length can be classified as SPIT/Vishing attacks, messages can be analyzed and compared to messages in the voice mailbox of real users. Such information can be used to process messages in the real voice mailbox and sort accordingly voice mailbox of users dividing messages in good ones and SPIT/Vishing ones in order to facilitate users in listening to their voice messages.

• deeper interaction with the originator in order to gather additional information on the source in order to characterize it and input such results to other identification mechanisms. An example of this are the fingerprint checks (IP addresses, software client used to initiate the call, IP path towards the source, other layer 4 ports the initiator may have open that can indicate it is initiating many calls in parallel, etc.) detailed in [START_REF] Quittek | Detecting SPIT calls by checking communication patterns[END_REF] or an extension of the audio turing test detailed in [START_REF] Yan | Incorporating Active Fingerprinting into SPIT Prevention Systems[END_REF] aimed to fingerprint the audio sent by machines. The output of such fingerprint can then be used by other detection mechanisms to further reduce the error rate (calls with characteristics observed already by the deflectors will be detected with higher probability).

• make the blocking (suppression in legal terms) of communications compliant to the suppression of communication law. The idea is to send the calls that are considered malicious by the VoIP Security Event Correlator detailed in 5 to the deflector infrastructure instead of blocking them. In this case the communication will not be suppressed and further recorded as if it were in a separate voice mailbox being available for the callee; at the same time these calls will be used to further refine the detection mechanism.

VOIP SECURITY EVENT CORRELATION (SEC)

As mentioned before, VoIP systems are distributed by nature, thus it is obvious that studying traces of just one VoIP user agent or server does not give an overall picture of what is happening to the whole system. An attack involving several domains could result in not being detected just monitoring the single traces at each domains. Aggregating traces from different sources is vital to detect useful signatures. Our approach is to deploy at every point of interest (proxy, user agent, gateways) a first layer of monitoring composed of threading and correlation. Then, events of interest are given as input from different sources to a central correlator. The event correlation hierarchy is depicted in figure 3, in order to better explain the advantage of such an approach the following example is given:

Malicious gateway.

In a normal scenario, a user sets up a call towards the PSTN by contacting a call agent. The call agent controls the gateway by using the MGCP protocol. It opens a media trunk in the gateway using a CRCX command letting the user send the RTP flow according to the specification included in the SDP protocol body. In a toll fraud attack scenario, the malicious user aims to bypass billing. The malicious user could set up the call normally, then immediately sends a BYE to the call agent. The call agent releases the call by a DLCX command, the accounting procedure stops upon receiving a 200 OK from the gateway. But if the gateway was instructed by the malicious user to not respond to the call agent and to send falsified notifications, the malicious user will continue to use the media trunk without being billed for the usage. Such intrusions could not be detected if only the MGCP call agent trace or the end point one are monitored but by a central Security Event Correlator (SEC) can detect RTP packets received by the endpoint after a DLCX command is sent by the call agent.

Event generation

Most VoIP servers and phones provide protocol logs either to standard output or as text files in order to help the debugging process. Asterisk, OpenSER and Kphone4 are some examples. However, the purpose of an intrusion detection system is not debugging VoIP agents but detecting security incidents.

The task of the event generator located at each call agent (MGCP Call Agent, SIP Proxy, etc.) is to extract predefined fields from the protocol message and to add a timestamp. SIP and MGCP signaling protocols have completely ASCII oriented headers and bodies (SDP). In case of binary based protocols as RTP for media transport and H323 for signaling, the event generator has to properly parse the relevant fields to ASCII to be coherent. Such an event generator can be configured to allow multi-layer analysis by including IP and TCP/UDP information and has all reassembly capabilities to deal with IP fragmentation, TCP segmentation and protocols decoding in order to recognize protocol messages. To better explain the event generation, the following array is given as example when a SIP INVITE message is considered: In order to detect intrusions in VoIP architectures a modelbuilding of normal behavior or malicious activities is of paramount importance. Appropriate specification of an event content and its relationship to other events in the time are the key requirements to build complex events from primitive ones. We show the importance of these requirements by mentioning some examples:

Event threading.

Besides primitive events, building patterns of activities is mandatory for security intrusion monitoring. User profiling is a possible example where user normal patterns need to be defined (e.g. user 'A' makes 5 calls on average during weekdays and 2 calls on average during weekends). To model such patterns, it is easier to deal with primitive call events than with basic protocol messages. A PBX has already support for such events in its Call Detail Records (CDRs) but a SIP proxy log needs a dedicated application matching every INVITE with the corresponding BYE to identify a session as follows: Also, when searching for some attack patterns, the signature of an attack could not be detected by matching one single event but by correlating a sequence of events appearing to be normal if analyzed separately. For example, in a flooding DoS we do like to bound a large number of similar events as just one alarming event and suppress not necessary redundancy by just adding the number of messages from which the alarming event was composed.

Temporal restrictions.

Scheduling restrictions and the inter-arrival timing of events are included in the temporal restrictions considered by VoIP SEC.

Scheduling restrictions are very useful in monitoring tasks. A MGCP call agent is supposed to send an AUdit Endpoint command (AUEP) to each endpoint under its supervision. Assuming that the normal activity is to make audits at fixed dates (e.g. every hour o'clock), any fluctuation in the call agent behavior can be detected by noticing the absence of such events at specified times.

The inter-arrival time of events is even more relevant. A short inter-arrival time of requests to a SIP target characterizes a flooding DoS, a high rate of appearance of '404 Not Found' responses characterizes a domain enumeration and a high rate of appearance of '403 Forbidden' responses characterizes a password cracking attack.

To clarify our approach regarding temporal restrictions let us take the following attack scenario:

1. One SIP caller from Internet attempts to call a PSTN number through a gateway controlled by the MGCP protocol;

2. the MGCP call agent receives the INVITE, translates it into SS7 signaling to ring the PSTN phone;

3. once the call is answered, the MGCP call agent responds to the SIP caller with a 200 OK, it sends a CRCX (create connection) command to one end point in the media gateway to link the RTP flow and a voice trunk;

4. the caller has to acknowledge such a command with the final response (200 OK). A malicious caller may do not complete the handshake to achieve one of two purposes:

• bypass billing: if the SIP call agent is compromised and it waits for the ACK to consider that the session is open and to start billing, the caller can use the media connection without being charged for it.

• resource exhaustion DoS: if the caller uses an army of bots initiating a set of messages aimed to open media connections bypassing billing without sending RTP data, the call agent will be overloaded with a high number of paralyzed state machines.

A solution addressing this attack scenario is to impose a temporal restriction on the ACK message that must be received within a few round trip times after the 200 OK.

Event Correlation using SEC

After investigating available software for event correlation, the most suitable was found to be the open source and platform independent SEC 5 . SEC fulfills VoIP events modeling requirements using static rules, accepting input from text streams named pipes and generating output events, log messages or executing shell commands. Each rule in SEC contains an event matching condition based on a Perl regular expression, an actions list, and optionally a boolean context that permits or not the application of the rule.

SEC is a lightweight online monitoring tool initially invented to fill the gap between homegrown and commercial event correlation solutions, and it already proved its efficiency in several domains as network management, intrusion detection systems, system monitoring and fraud detection. However, to the authors best knowledge, with new generation Internet applications having specific properties as VoIP, it is not been tested yet. Some VoIP threats have known patterns that can be detected by signatures matching. However, the misuse detection solution based on signatures can not accommodate previously unknown attacks and is not effective to detect service abuse. Anomaly-based intrusion detection is efficient to detect fraudulent usage, virus propagation and social threats like SPIT. In the anomaly-based approach, a statistical profile is created for normal activity of subjects 5 http://kodu.neti.ee/~risto/sec/ (e.g. user, call session, SIP server) with respect to objects (e.g. network resources, gateway trunks, server CPU and memory). This statistical profile generates anomaly records when deviations from normal behavior are occurring. SEC rules are efficient to detect attack signatures, and statistical profiles of users, group of users and traffic can be defined to monitor behavior anomalies. In table 5.2, examples of both approaches are given to show how different types of SEC rules can be designed.

Implementation

We have implemented a module in the OpenSER code that copies the SIP messages routed by the server towards a parser that is written with SEC rules. A few lines in the OpenSER configuration file have been added after installation for this scope. The SEC parser creates a line with fields Three rules of type PairWithWindow are needed to detect such an attack. PairWithWindow is a composed SEC rule that after a first defined event arrives, waits for t seconds for other defined events, and executes one list of actions if the second event arrived in time, or another list of actions if the window timeout expires. The detection scheme for this attack consists on: • the third PairWithWindow rule to match the INVITE-2OO OK-BYE event with a RTP event sent from the caller to the callee that occurs in the time window and in this case write an alarms to the GUI as shown in figure 4.

DoS using broken handshaking: this attack is based on broken SIP handshaking where the attacker sends an INVITE request and then ignores the 200 OK response refusing to send the ACK. The attacker (or an army of bots) proceeds with a large number of broken initiations in order to exhaust the target performances.

one rule of type PairWithWindow and one rule of type Single-WithThreshold are needed to detect this attack. SingleWith-Threshold is a composed SEC rule that counts matching events in a window of t seconds, and if a defined threshold is exceeded executes a list of actions, otherwise it waits the expiration of the time window to execute another list. The detection scheme for this attack consists on:

• one PairWithWindow rule to match the INVITE with the corresponding 200 OK and generate an INVITE-2OO OK event with the Call-ID, From and To tags;

• the second rule to match the INVITE-2OO OK event with the corresponding ACK event. If no ACK is received within the time window, the rule generates a broken handshaking event;

• the SingleWithThreshold rule counts the broken handshaking event in a time window and writes an alarm to the GUI in case of a defined number is exceeded as shown in figure 5.

Anomaly-based intrusion detection User profile: One method of building an account profile is to make use of histograms. A day is divided into bins of specified time (e.g. one hour). For each bin a predefined metric is calculated (e.g. number of calls, number of different recipients, average duration of a call) matching predefined events (e.g. call). In the learning phase (e.g. a month), daily statistics are built to extract a long term account profile (e.g. daily average of the number of calls for each bin). In the detecting stage (e.g. a day), a short term profile is compared to the long term one by using an appropriate distance function (e.g. Euclidean distance, quadratic distance, Mahalanobis distance). A recent profile which is quite different from the long term one indicates possible misuse. On the other hand, long term profiles can be compared to different accounts to group them into classes. This is of high importance since a class of VoIP bots send SPIT calls can be detected if a known bot profile falls into the same class. Another method is to study non stationary features of an account, for example the distribution of calls over all callees or the shape of the callees Š list size over all dialed calls. By comparing changes of a distribution over the time by using of an appropriate distance function (e.g. Hellinger distance), sudden bursts may be detected and treated as abnormalities.

In addition to appropriate scripts to build up such a profile for a given user, SEC rules can be of importance to manage the process. Calendar is a basic SEC rule that gets activated at specific times and executes a list of actions using a UNIX cronlike syntax. SingleWithScript is a composed SEC rule that matches input event and executes a list of actions depending form the return value of an external script. For a complete view of SEC design and allowed actions, please refer to [START_REF] Brown | Working with SEC -the Simple Event Correlator[END_REF]. The detection of abnormal behaviors can be implemented using SEC as follow:

• a rule of type Calendar gets activated at specified times to launch a time to compare event.

• a rule of type SingleWithScript catches the event and launches an external script that has access to short-term and long-term profiles. If the external script returns a zero value for a weak similarity, the corresponding action will notify about an abnormal behavior. The corresponding rule set is depicted in figure 6. Table 2 shows the attacks implemented so far and the number of SEC rules written for their detection. The scalability of the solution is clear when looking at the small number of rules sufficient to detect different series of attacks. The performance of the VoIP Security Event Correlator built in this way is therefore promising but was not quantitatively tested. In addition to OpenSER extensions, the current work is focused on implementing a "local" event correlator in Asterisk and in additional softphones (Kphone) and a central correlation engine as a second layer of detec-tion.

RELATED WORKS

Intrusion detection systems are a second line of defense behind intrusion prevention mechanisms as password authentication and firewalls. One of the earlier works in this domain is the model proposed in [START_REF] Denning | An Intrusion-Detection Model[END_REF]. The authors of [START_REF] Krügel | Service specific anomaly detection for network intrusion detection[END_REF] shows the necessity of domain knowledge in specific IDSs especially with web-based applications. Security is of great interest in the new VoIP generation design so several intrusion detection approaches are proposed in response to different threats. Scidive [START_REF] Wu | SCIDIVE: A Stateful and Cross Protocol Intrusion Detection Architecture for Voice-over-IP Environments[END_REF] uses signature-based statefull and cross protocol schemes. We proposed a statistical frame work based on Bayes model to recognize normal and attack SIP traffic classes in [START_REF] Nassar | Intrusion detections mechanisms for VoIP applications[END_REF]. SEC was proposed in [START_REF] Vaarandi | SEC -A Lightweight Event Correltion Tool[END_REF] as a lightweight event correlator that can serve different applications ranging from log file and system monitoring to fraud detection, network management and intrusion detection and a good overview on it can be found in [START_REF] Rouillard | Real-tile Logfile Analysis Using the Simple Event Correlator (SEC)[END_REF]. Recent papers on the detection of social attacks proposed blacklist/graylist type solutions [START_REF] Shin | Voice SPAM Control with Gray Leveling[END_REF] and centralized network entities [START_REF] Mathieu | SPIT mitigation by a network level anti SPIT entity[END_REF]. Previous work of a subset of the authors [START_REF] Niccolini | Third annual security workshop (VSW'06)[END_REF] considered network level plugins for Snort capable of detecting SPIT attacks. An interesting idea of SPIT detection based on device fingerprinting is exposed in [START_REF] Yan | Incorporating Active Fingerprinting into SPIT Prevention Systems[END_REF]. A first description of a VoIP honeypot can be found in [START_REF] Nassar | VoIP Honeypot Architecture[END_REF] where a subset of the current authors proposed network level SIP honeypot.

CONCLUSION

We have presented in this paper a holistic approach for VoIP security monitoring. The key components of our solution are a VoIP honeypot and a SIP level event correlation engine. This solution is capable to defend against both brute force denial of service attacks as well as more stealthy social type (SPIT, Vishing). We leveraged the capabilities of SEC to constitute the technical tool in building efficient VoIP IDS and showed the feasibility by developing a prototype. We extended in this paper the notion of a honeypot towards VoIP applications and showed how social attacks can be mitigated. We have implemented and tested our solution in a testbed environment, but more real life tests and performance evaluation will be done in the future on a VoIP network. Future work will address also the extension of the current event correlation techniques towards machine learning inspired paradigms.

Figure 1 :

 1 Figure 1: Network Architecture

Figure 2 :

 2 Figure 2: Honeypot Architecture

Figure 3 :

 3 Figure 3: Two layers of event correlation

Figure 4 :

 4 Figure 4: Diagram of SEC ruleset to detect BYE attack

Figure 5 :

 5 Figure 5: Diagram of SEC ruleset to detect broken handshaking flooding

•

  one PairWithWindow rule to match an INVITE with the corresponding 200 OK and generate an INVITE-2OO OK event with the Call-ID, From and To tags, the media IP and port of the caller (in the INVITE) and the media IP and port of the callee (in the 200 OK); • the second rule to match the INVITE-200 OK event with the corresponding BYE event and generate an INVITE-2OO OK-BYE event;

Figure 6 :Figure 7 :

 67 Figure 6: SEC ruleset to manage profile comparison

Table 1 :

 1 Using SEC to build efficient correlation rules Attack description Detection scheme Signature-based intrusion detection DoS using BYE attack: In a call between A and B, a malicious entity C sends a DoS on A by prematurely tearing down the call with a crafted BYE message to A. If A receives RTP packets from B after receiving the BYE, this means that B was unaware of the BYE sent. Similar attacks are call hijacking, instant messenger hijacking and RTP play out attack.

Table 2 :

 2 Size of rulesets detecting various attacks

	Name of attack	Number of rules Name of attack	Number of rules
	Request flooding	2	BYE attack	3
	Broken handshaking flooding	3	CANCEL attack	3
	REGISTER scan	4	Malicious gateway	4
	password cracking	2	SPIT	4
	type=Calendar
	time=0 0 * * *
	desc=time to compare
	action= event 0 %s user1
	type= SingleWithScript
	ptype=RegExp
	pattern=time to compare (\S+)
	script= /usr/bin/perl compare.pl $1
	desc=compare long-term and short-term profiles for user $1
	action= shellcmd notify.sh ''abnormal behavior user $1''

http://www.openser.org/

http://www.asterisk.org/

http://sourceforge.net/projects/kphone


Jie Hu
email: e0013640@u.nus.edu

A simple numerical scheme for the 2D shallow-water system

Keywords: shallow-water equations, approximate Riemann solver, Godunov-type finite volume method, well-balanced scheme, flow regime, free-surface shallow flow

, this paper extends the problem from 1D to 2D with the simplicity of application preserves. The new scheme is implemented into the code TELEMAC-2D [tel2d, 2014] and several tests are made to verify the scheme ability under an equilibrium state at rest and different types of flow regime (i.e., fluvial regime, transcritical flow from fluvial to torrential regime, transcritical flow with a hydraulic jump). The sensitivity analysis is conducted to exam the scheme convergence.

Introduction

The SWEs have been proposed by [START_REF] Saint-Venant | Theorie du mouvement non permanent des eaux, avec application aux crues des rivieres et a l'introduction de marees dans leurs lits[END_REF] to model flows in a channel.

Nowadays, they are used to model flows in a wide variety of physical phenomena, such as: overland flow, flooding, dam breaks, tsunami (e.g. [START_REF] Esteves | Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values[END_REF][START_REF] Caleffi | Finite volume method for simulating extreme flood events in natural channels[END_REF][START_REF] Valiani | Case study: Malpasset dam-break simulation using a two-dimensional finite volume method[END_REF][START_REF] Kim | Propagation and run-up of nearshore tsunamis with HLLC approximate Riemann solver[END_REF]). These equations are a time-dependent twodimensional system of non-linear partial differential equations of hyperbolic type.

In real situations (realistic geometry, sharp spatial), there is little hope to solve explicitly the SWEs, i.e. to produce analytic formula for the solutions. Therefore, it is necessary to develop specific numerical methods to compute approximate solutions of SWEs (e.g. [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF][START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. Implementation of any of such methods raises the question of the validation of the code.

Validation is a necessary step to check if a model (the numerical methods) suitably describes the considered phenomena. There exists at least three complementary types of numerical tests to ensure a numerical code is relevant for the considered systems of equations. First, one can produce convergence or stability results (e.g. by refining the mesh). This validates only the numerical method and its implementation. Second, approximate solutions can be matched to analytic solutions available for some simplified or specific cases. Finally, numerical results can be compared with experimental data, produced indoor or outdoor. This step should be done after the previous two; it is the most difficult one and must be validated by a specialist of the domain. This paper focuses on the first two steps.

A simply implementary scheme for 1D SWEs is provided in [START_REF] Audusse | A simple well-balanced and positive numerical scheme for the shallow-water system[END_REF], and this scheme is proved to be accurate and robust on several typical test cases. Enlightened by [START_REF] Audusse | A simple well-balanced and positive numerical scheme for the shallow-water system[END_REF]'s work, the present paper describes a numerical scheme for the 2D SWEs to study the free surface shallow flows.

The paper is organized as follows: the general mathematical model is described in section 2.

In section 3, the property of rotational invariance of SWEs is applied to split the governing equations into x direction, thus simplifying the problem. The classic Riemann solver is reviewed in section 4 and a Godunov-type finite volume scheme is derived for the augmented 1D SWEs. A simple treatment of the source term is adopted from [START_REF] Audusse | A simple well-balanced and positive numerical scheme for the shallow-water system[END_REF], in section 5. Several test cases are conducted to exam the ability of the new scheme under the equilibrium state and several types of flow regime, in section 6.

Mathematical model

The SWEs can be deduced from Navier-Stokes equations for an incompressible fluid by making the hypothesis of hydrostatic pressure, uniform velocities along the vertical direction.

For inviscid flow, the model can be written in its conservative form as follows

Rotational invariance

We first consider the homogeneous (no source terms) time-dependent two-dimensional SWEs.

𝑼 𝑡 + 𝐹(𝑼) 𝑥 + 𝐺(𝑼) 𝑦 = 𝟎 Eq. 2

where 𝑼 𝑡 is

The Riemann Problem and the Godunov Flux

Our concern is about solving numerically the general initial-boundary value problem (IBVP)

for the augmented one-dimension SWEs in section 3. Utilising Godunov-type methods based on the explicit conservative formula

𝑼 𝑖 𝑛+1 = 𝑼 𝑖 𝑛 - 𝛥𝑡 𝛥𝑥 [𝐹 𝒊+ 1 2 -𝐹 𝒊-1 2 ] , ∀ 𝑖 ∈ 𝕫 Eq. 5
where 4.1 The HLLC Riemann solver [START_REF] Harten | On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws[END_REF] suggested a way of solving the Riemann problem approximately by finding an approximation to the numeral flux 𝐹 𝒊+ 1 2 (flux through the interface). The mathematical bases of the approach are given in [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF]. Assuming all wave speed estimates are available, the HLLC numerical flux is shown as follows

𝑼
𝐹 𝑖+ 1 2 ℎ𝑙𝑙𝑐 = { 𝐹 𝐿 , 𝑖𝑓 0 ≤ 𝑆 𝐿 , 𝐹 𝐿 * , 𝑖𝑓 𝑆 𝐿 ≤ 0 ≤ 𝑆 * , 𝐹 𝑅 * , 𝑖𝑓 𝑆 * ≤ 0 ≤ 𝑆 𝑅 , 𝐹 𝑅 , 𝑖𝑓 0 ≥ 𝑆 𝑅 Eq. 6
where

𝐹 𝐿 * = 𝐹 𝐿 + 𝑆 𝐿 (𝑈 𝐿 * -𝑈 𝐿 ) & 𝐹 𝑅 * = 𝐹 𝑅 + 𝑆 𝑅 (𝑈 𝑅 * -𝑈 𝑅 )
Eq. 7

The states 𝑈 𝐿 * , 𝑈 𝑅 * are given by

𝑈 𝐾 * = ℎ 𝐾 ( 𝑆 𝐾 -𝑢 𝐾 𝑆 𝐾 -𝑆 * ) [ 1 𝑆 * 𝑣 𝐾
], (𝐾 = 𝐿, 𝑅) Eq. 8

For determining the numerical flux, the three approximate wave speeds 𝑆 𝐿 , 𝑆 𝑅 and 𝑆 * should a priori be known.

4.1.1 Estimation of speeds 𝑆 𝐿 , 𝑆 𝑅 and 𝑆 *

The determination of numerical flux in Eq. 6 requires the pre-known estimated wave speeds. [START_REF] Toro | Shock-Capturing Methods for Free-Surface Shallow Flows[END_REF] suggests the following choice of wave speeds that can lead to accurate and robust scheme:

𝑆 𝐿(𝑅) = 𝑢 𝐿(𝑅) -𝑎 𝐿(𝑅) 𝑞 𝐿(𝑅)
Eq. 9

where 𝑎 𝐾 = √𝑔ℎ 𝐾 and 𝑞 𝐾 (𝐾 = 𝐿, 𝑅) is given by

𝑞 𝐾 = { √ 1 2 [ (ℎ * +ℎ 𝐾 )ℎ * ℎ 𝐾 2 ] 𝑖𝑓 ℎ * > ℎ 𝐾 1 𝑖𝑓 ℎ * ≤ ℎ 𝐾
Eq. 10 here ℎ * is an estimate for the exact solution for ℎ in the intermediate region. Similar choices can also be found in [START_REF] Bello | A HLLC Riemann solver to compute shallow water equations with topography and friction[END_REF]. This new solver has a simple form. It can deal very well with situations involving very shallow water and is found to be very robust in dealing with shock waves [START_REF] Toro | Direct Riemann solvers for the time-dependent Euler equations[END_REF].

Including of the topography

The previous section (4.1) elaborates the HLLC scheme for the homogeneous SWEs. In this section, our scheme is extended to include the topography term. And Eq. 4 becomes SWEs with uphill topography. There are four approximate waves with speeds 𝑆 𝐿 , 0, 𝑆 * and 𝑆 𝑅 , while 0 is induced by the effect of topography.

𝑼 𝑡 + 𝐹(𝑼) 𝑥 = ( 0 -𝑔ℎ𝜕 𝑥 𝑏 0 ) Eq.
Integrating Eq. 12 on volume [𝑥 𝑙 , 𝑥 𝑟 ]x[0, 𝑡] and deploy the first two components (as flux is along x-axis, no impact on third component), we can have:

ℎ 𝑅 𝑢 𝑅 -ℎ 𝐿 𝑢 𝐿 = 𝑆 𝐿 (ℎ 1 -ℎ 𝐿 ) + 𝑆 * (ℎ 3 -ℎ 2 ) + 𝑆 𝑅 (ℎ 𝑅 -ℎ 3 ) Eq. 13

(ℎ 𝑅 𝑢 𝑅 2 + 𝑔ℎ 𝑅 2 2 ) -(ℎ 𝐿 𝑢 𝐿 2 + 𝑔ℎ 𝐿 2 2 ) + 𝑔∆𝑥{ℎ𝜕 𝑥 𝑏} = 𝑆 𝐿 (ℎ 1 𝑢 1 -ℎ 𝐿 𝑢 𝐿 ) + 𝑆 * (ℎ 3 𝑢 3 -ℎ 2 𝑢 2 ) + 𝑆 𝑅 (ℎ 𝑅 𝑢 𝑅 -ℎ 3 𝑢 3 )
Eq. 14

where {ℎ𝜕 𝑥 𝑏} stands for a consistent approximation of the source term ℎ𝜕 𝑥 𝑏 which will be precised later on. In order to close this system, two relations are missing and we suggest to impose:

 Hydraulic balance between states 𝑈 1 and 𝑈 2 . Eq. 20

{ ℎ 1 𝑢 1 = 𝑞 1 = ℎ 1 𝑢 1 = 𝑞 2 = 𝑞 * ℎ 1 + 𝑏 1 = ℎ 2 +
where 𝐹 𝐿 * is the HLLC flux in Eq. 6.

Hence, the numerical flux for SWEs with uphill bottom is:

𝐹 𝑖+ 1 2 ℎ𝑙𝑙𝑐 (𝑼 𝐿 , 𝑼 𝑅 ) = { 𝐹 𝑖 , 𝑖𝑓 0 ≤ 𝑆 𝐿 ,

Eq.21

It should note that similar expression can also be found in [START_REF] Audusse | A simple well-balanced and positive numerical scheme for the shallow-water system[END_REF].

An downhill bottom 𝜕 𝑥 𝑏 > 0

Analogous with 4.2.1, for the topography with the left bottom height 𝑏 1 higher than the right bottom height 𝑏 2 , it can be deduced that 𝑆 * is negative, it has been found that the numerical flux for SWEs with downhill bottom has the same expressions as in Eq. 21 except the HLLC flux 𝐹 𝐿 * is replaced by 𝐹 𝑅 * in Eq. 6.

5 Treatment of the source term {𝒉𝝏 𝒙 𝒃} An adopted method is from [START_REF] Audusse | A simple well-balanced and positive numerical scheme for the shallow-water system[END_REF] to discrete source term as follows

{ℎ𝜕 𝑥 𝑏} = { ℎ 𝐿 +ℎ 𝑅 2∆𝑥 𝑚𝑖 𝑛(ℎ 𝐿 , ∆𝑏) , 𝑖𝑓 ∆𝑏 ≥ 0, ℎ 𝐿 +ℎ 𝑅 2∆𝑥 max (-ℎ 𝑅 , ∆𝑏), 𝑖𝑓 ∆𝑏 < 0, Eq. 22
This treatment can preserve the lake at rest in the case of a wet-dry transition or the case of a dry-wet transition.

Numerical results

The convergence of the new numerical scheme is first examined on two flow regimes (a fluvial flow and a transcritial flow). This validates the numerical method and its implementation. The behaviour of our new scheme is then examined on several test cases to check if the approximate solutions are accurate. These test cases include: an equilibrium state with a bump bottom, a steady flow over this bump in two regimes (fluvial and transcritical).

Mesh sensitivity study

To ensure the convergence of our new scheme, three different meshes are built for a rectangular topography. This topography is discretized by N triangular elements, where N are 968, 1925 and 3839, respectively (see Fig 3). The free software Bluekenue is used for creating these meshes. And more information of the rectangular topography is in section 6.2.

The relative error norms are defined as following:

‖𝑒 1 (𝑁)‖= ∫ |𝐻 𝑁 (𝑥)-𝐻 𝑒𝑥𝑎𝑐𝑡 (𝑥)|𝑑𝑥 𝐿 0 𝐿 Eq. 23.1 ‖𝑒 2 (𝑁)‖= √ ∫ (𝐻 𝑁 (𝑥)-𝐻 𝑒𝑥𝑎𝑐𝑡 (𝑥)) 2 𝑑𝑥 𝐿 0 𝐿 Eq. 23.2
‖𝑒 1 (𝑁)‖=max (𝐻 𝑁 (𝑥) -𝐻 𝑒𝑥𝑎𝑐𝑡 (𝑥))

Eq. 23.3

with 𝐻 𝑁 (𝑥) the free surface (water depth h plus bottom elevation b) profile from the new scheme with N cell elements used and 𝐻 𝑒𝑥𝑎𝑐𝑡 (𝑥) is the free surface profile from exact solution as the reference result with the same cell elements. L is the total length of simulated domain and dx is the space step. The simulations are performed with different cell element numbers (i.e., 968, 1925, 3839)

Equilibrium steady state at rest

The aim of this test is to verify the scheme for equilibrium steady state. The topography is flat at both left and right sides except there is a bump in the middle. The length of the topography is 20.5 m. And the topography is discretized by 968 elements (Fig 3). This geometry is adopted for all the following tests. The free surface is initially constant with a height 1.8 m which submerges all the topography. At the right side, we impose a 1.8 m water depth which is equal to the initial condition, and at the left side, a zero discharge is imposed. In this way, we can ensure the equilibrium state. The time step is 0.01 s, and the total number of time steps is 10000 which is long enough to ensure the steady state. After the simulation, the result of free surface is shown in Fig 6. noted that the free surface remains stable since the simulation starts. We can thus believe that our scheme can assure the equilibrium steady state at rest. Next, we exam our scheme in nonequilibrium steady states which consist of three types of situation: fluvial regime, transcritical flow from fluvial regime to torrential regime, transcritical flow with a hydraulic jump. The results calculated by our new scheme are compared quantitatively with exact solution.

Fluvial regime

The initial condition is the same as that in 6.2. A 1.8 m water depth is imposed at the right side and a discharge of 8.85889 m 3 /s is imposed at the left side. After the flow is at steady state, the free surface and discharge is depicted in Fig 6.

Transcritcal flow from fluvial regime to torrential regime

The free surface is initially constant with a height 0.13 m which submerges all the topography.

Only a discharge of 0.6 m 3 /s is imposed at the left side. The time step is 0.01 s, and the total number of time steps 10000 which is long enough to ensure the steady state. The simulation result is depicted in Fig 7.

Transcritical flow with a hydraulic jump

The initial condition is the same as that in 6.3. A discharge of 0.18 m 3 /s is imposed at the left side and a water depth with the height 0.13 m imposed at the right side. The time step is 0.01 s, and the total number of time steps 30000 which is long enough to ensure the steady state.

The simulation result is depicted in Fig 8. From Fig 6-8, the numerical result of our new scheme and analytic solution are in very good agreement. By looking carefully at the differences in free surface, it appears that the difference is extremely low. While for discharge, the maximum differences are 0.27%, 0.93% and 2.22%, respectively. The fluctuation in discharge from the new scheme is due to the numerical error which can be improved by densifying the computing mesh number. It is also clear that the fluctuation is tiny and the accuracy can be guaranteed.

Conclusion

We have proposed a simple numerical scheme for the two dimensional Shallow-Water Equations (SWEs). This new scheme can ensure the well-balance condition and is proved to be accurate on several typical test cases: Fluvial regime, transcritical flow from fluvial regime to torrential regime and transcritical flow with a hydraulic jump.

  ) in Eq. 1. Following the rotational invariance property of the SWEs by Toro (2001), the two-dimensional problem can be reduced to augmented one-dimensional: 𝑼 𝑡 * + 𝐹(𝑼 * ) 𝑥 = 𝟎 Eq. 3 where 𝑼 * = [ℎ, ℎ𝑢 ̃, ℎ𝑣 ̃]𝑇 , 𝑢 ̃, 𝑣 ̃ are the velocity components after the angular rotational by The new velocities are related with the original ones by 𝑢 ̃= 𝑢𝑐𝑜𝑠𝜃 + 𝑣𝑠𝑖𝑛𝜃, 𝑣 ̃= -𝑢𝑠𝑖𝑛𝜃 + 𝑣𝑐𝑜𝑠𝜃. The tilde (. ) ̃ and (*) is dropped later for brevity.

  Fig 1 Structure of the approximate solution of the Riemann problem for the x-split homogenous twodimensional SWEs. There are three approximate wave families 𝑆 𝐿 , 𝑆 𝑅 and 𝑆 * associated with the proper eigenvalues of Eq. 4.

  uphill bottom 𝜕 𝑥 𝑏 < 0 If the topography is not flat and the left bottom height 𝑏 1 is less higher than the right bottom height 𝑏 2 , which means the left water depth ℎ 𝐿 is greater than the right water depth ℎ 𝑅 , if the free surface is at the same elevation. Then the first approximation of the height in intermediate region ℎ * is smaller than ℎ 𝐿 , but bigger than ℎ 𝑅 (ℎ 𝐿 > ℎ * > ℎ 𝑅 ). We deduce that 𝑆 * is positive. Left intermediate region 𝑈 𝐿 * in Fig. 1 is divided into two zones denoted as 𝑈 1 and 𝑈 2 due to the influence of the topography. This is equivalent to introduce a velocity whose value is zero at t-axis, shown in Fig 2.

Fig 2 .

 2 Fig 2. Structure of the approximate solution of the Riemann problem for the x-split two-dimensional

Fig 3 Fig 4

 34 Fig 3 Top view of the topography discretized by 968 triangular elements, the refined area is with more intensive elements, different colours represent different elevations of the topography.

Fig 6 A

 6 Fig 6 A front view of equilibrium steady state at rest over a bump

Fig 6

 6 Fig 6 Fluvial flow: Comparison of free surfaces h+b (left) and discharge q (right) resulting from the new scheme represented by asterisk and the exact solution represented by line for a bottom topography b.

Fig 7

 7 Fig 7 Transcritical flow: Comparison of free surfaces h+b (left) and discharge q (right) resulting from the new scheme represented by asterisk and the exact solution represented by line for a bottom topography b.

Fig 8

 8 Fig 8 Flow with hydraulic jump: Comparison of free surfaces h+b (left) and discharge q (right) resulting from the new scheme represented by asterisk and the exact solution represented by line for a bottom topography b.

  4.1.2 Estimation of water depth ℎ * and particle velocity 𝑢 * in the intermediate region The primitive variables in intermediate region are a prior needed to determine wave speeds in sub-section 4.1.1. In[START_REF] Toro | Shock-Capturing Methods for Free-Surface Shallow Flows[END_REF], a new scheme solver is given for ℎ * and 𝑢 * .

	ℎ * = 𝑢 * = 1 2 1 (ℎ 𝐿 + ℎ 𝑅 ) -2 (𝑢 𝐿 + 𝑢 𝑅 ) -(ℎ 𝑅 -ℎ 𝐿 )(𝑎 𝐿 + 𝑎 𝑅 )/(ℎ 𝐿 + ℎ 𝑅 ) 1 4 (𝑢 𝑅 -𝑢 𝐿 )(ℎ 𝐿 + ℎ 𝑅 )/(𝑎 𝐿 + 𝑎 𝑅 )	}	Eq. 11

Acknowledgements

This study is one part of my final year's internship financially supported by LNHE R&D, EDF. Thanks is given to Dr Ung for his attentive explanation of his new numerical approach solving SWEs and my advisor Dr Goutal for her kindly guidance. Thanks is also given to Prof Delestre for his sincere suggestions.