Numerical shockwave anomalies in the resolution of the Shallow Water Equations with bed variations

. The presence of numerical shockwave anomalies appearing in the resolution of hyperbolic systems of conservation laws is a well-known problem in the scientiﬁc community. The most common anomalies are the carbuncle and the slowly-moving shock anomaly. They have been studied for decades in the framework of Euler equations, but only a few authors have considered such problems for the Shallow Water Equations (SWE). In this work, the SWE are considered and the aforementioned anomalies are studied. They arise in presence of hydraulic jumps, which are transcritical shockwaves mathematically modelled as a pure discontinuity. When solving numerically such discontinuities, an unphysical intermediate state appears and gives rise to a spurious spike in the momentum. This is observed in the numerical solution as a spike in the discharge appearing in the cell containing the jump. The presence of the spurious spike in the discharge has been taken for granted by the scientiﬁc com-munity and has even become a feature of the solution. Even though it does not disturb the rest of the solution in steady cases, it produces an undesirable shed-ding of spurious oscillations downstream when considering transient events. We show how it is possible to deﬁne a coherent spike reduction technique that reduces the magnitude of this anomaly and ensures convergence to the exact solution with mesh reﬁnement. Concerning the carbuncle, which may also appear in presence of strong hydraulic jumps, a combination of Riemann solvers is proposed to circumvent it. Also, it will be shown how there is still room from improvement when treating anomalies in moving hydraulic jumps over variable topography.


Introduction
It has been widely reported in the literature that numerical anomalies arise in presence of shock waves. Examples are the Carbuncle [1,2], the slowly-moving shock [3] and the wallheating phenomenon, all of them leading to spurious numerical solutions. In this work we focus on the Carbuncle and the slowly-moving shock anomaly.
Shock-type solutions are typically found when solving nonlinear hyperbolic systems of conservation laws. Their computation is not trivial and has been a challenging task for more than 50 years, since von Neumann, Richtmyer, Lax [4] and Godunov [5], among others, established the foundations of modern CFD. At present, the numerical resolution of shocks is still an issue under investigation, not only in the framework of gas dynamics but also for geophysical problems governed by the SWE.
Physical shockwaves have a finite width which is determined by the physical dissipation processes. On the other hand, shocks are mathematically represented by pure discontinuities in hyperbolic systems. When considering numerical shockwaves, a numerical width much greater than the physical width, is enforced by the grid size [6]. This leads to intermediate states which cannot be given a direct physical interpretation. Such states cannot be removed even when refining the grid, hence numerical schemes must be designed in a particular way to overcome such flaw. Up to the present time, most studies have been carried out in the framework of Euler equations, but the growing needs for the computation of complex geophysical flows with a morphodynamical component motivate their application to the SWE. In this work, we will extend those techniques to the SWE with bed elevation in 1D and 2D.

Mathematical and numerical models
The SWE is a nonlinear hyperbolic system of conservation laws of the form U t + F x = S, with the conserved variables, the fluxes and sources, where h is the water depth, u is the depth averaged velocity, hu the discharge and g is the acceleration of gravity. The source term S z involves the variations in bed geometry S z = −gh dz dx where z stands for the bed elevation. All the numerical techniques studied and proposed here are based on the Finite Volume method (Godunov's scheme), which reads where U n i are the discrete cell-average values at time t = t n , F − i+1/2 and F + i−1/2 the numerical fluxes at cell interfaces and ∆t and ∆x the time step and cell size, respectively.

Slowly moving shockwaves: hydraulic jumps
Here, we focus on the slowly-moving shock problem, which is associated to hydraulic jumps in the SWE. The slowly-moving shock problem was first investigated by Roberts in [3], who defined it as numerical noise generated in the discrete shock transition layer which is transported downstream. Such noise will be hereafter referred to as post-shock oscillations. The slowly-moving shock problem is related to nonlinearities of the Hugoniot curves. Such nonlinearities are found in those Hugoniot curves related to hydraulic jump-type solutions in the SWE.
Hydraulic jumps occur when a supercritical flow suddenly changes to subcritical conditions, generating a steep free surface elevation where intense mixing takes place and a large amount of mechanical energy is dissipated. Mathematically, hydraulic jumps are defined as:

Definition 1 (Hydraulic jump). Let the following discontinuous solution
As outlined before, it is worth recalling that when solving a hydraulic jump, an intermediate state appears in the numerical solution provided by Godunov's scheme with independence of the solver. The presence of this intermediate state, hereafter denoted by U M , is not of any physical relevance. It provides an unrealistic estimation of the average discharge in the intermediate cell, which is even not bounded by the left and right discharges.
As a first approach, let us compare analytically the solution for an ideal steady hydraulic jump, also called 2-state jump (U L and U R ), with the solution of a 3-state jump (U L , U M and U R ), which resembles the discrete solution provided by Godunov's scheme. Both solutions are weak solutions of the equations and they are both valid, as they satisfy the Rankine-Hugoniot (RH) conditions. Such comparison is presented in Figure 1. As observed in Figure 1, the 2-state jump is feasible because the states U L and U R lay on the intersection of the Hugoniot curve with the constant discharge curve. On the other hand, the 3-state jump is also possible as the state U M also lays on the Hugoniot curve, however, does not match the constant discharge value anymore as the curve is not linear between U L and U R . This is the explanation for the spurious intermediate state when computing hydraulic jumps, also referred to as spike.

Overcoming the slowly moving shock anomaly
The idea of flux interpolation, first presented by Zaide and Roe [6], is here applied to the SWE and extended for the non-homogeneous case. Such method is based on computing the fluxes in the untrustworthy intermediate cells by extrapolation from trustworthy neighbors using novel interpolation functions, which enforce a linear shock structure. The authors claim that by doing so, numerical shockwave anomalies are dramatically reduced.
Prior to the construction of the novel numerical fluxes F ⋆ i+1/2 , physical fluxes (which are the cell centered fluxes, F i ) are used to construct a novel approximation of the fluxes in every cell. Cell-centered fluxes, F i , are re-computed by means of extrapolation from neighboring cells.
Given by the following steps: 1. Flux interpolation according to Zaide and Roe [6]: 2. Correction ofF i by computing a novel flux estimationF i , as follows: where x S is the shock position inside the cell andS i−1,i+1 andS i−1/2 are suitable approximations of the source term that, for the sake of consistence, must satisfȳ It is worth pointing out that a shock detection algorithm is required. 3. Upwind the corrected fluxesF i at the interface to construct the numerical fluxes using an augmented solver: whereγ are the components ofΓ i+1/2 = P −1 i+1/2 δF i+1/2 , the projection of the jump in the extrapolated fluxes across cell interfaces,F i+1/2 =F i+1 −F i and β are the components ofB i+1/2 = P −1 i+1/2S i−1/2 .

Test case: travelling jump over an irregular bed profile
In this test case, a traveling hydraulic jump over an irregular bed profile is computed. To construct a solution consisting of a single jump traveling across the domain, we first compute a steady transcritical solution over the bed profile by imposing a constant discharge upstream of q = 0.6 m 2 /s. When the steady regime is reached, the boundary condition upstream is redefined, imposing now q = 0.556749458405104 m 2 /s and h = 0.12 m, which generates a supercritical state that is connected with the original subcritical state by means of a traveling hydraulic jump. The computational domain is [0, 560] and the solution is computed at t = 610 s. The CFL number is set to 0.45, g = 9.8 m/s 2 and the domain is discretized in 140 computational cells.
The bed profile will be constructed as where g(x) is an additional geometric function defined as Numerical results presented in Figures 2 show the numerical solution at t = 610 s for the water surface elevation and discharge provided by the ARoe scheme and by the proposed spike-reducing method. Major differences are observed in the solution of the discharge, which is much more oscillatory when computed by the ARoe method. On the other hand, differences on the water surface elevation are less sensitive to the spike.

Extension to 2D
The 2D extension of the spike-reducing solver is carried out by applying the 1D methodology to each direction independently. When considering a Cartesian mesh, it is possible to define the x and y interpolated fluxes as done in the 1D case. In [6], the author outlines that in the stationary case, each intermediate shock state is adjacent to at least two end states of the shock, but not necessarily aligned in the x or y direction. This would require a genuinely twodimensional method, using interpolated fluxes computed from information in both directions, however the dimension-by-dimension method proposed here is powerful enough to provide the sought results.

Test case: 2D shock wave over a sinusoidal inclined plane
This numerical experiment consist of a supercritical flow that hits a circular obstacle and generates an steady bow shock around the obstacle. The computational domain is Ω = x | (x − 80) 2 + (y − 50) 2 ≤ 400 , x ∈ Ω � and the bed elevation is given by The solution is computed at t = 150 s using CFL=0.4. In Figure 3, the numerical discharges computed by the traditional ARoe solver and the spike-reducing solver, using two different grids with using ∆x = 1 and ∆x = 0.5, are presented. A strong influence of the topography is observed in the shape of the shock. The reduction of the spike along the shock profile, provided by the spike-reducing solver, is again remarkable. The spike is marked by the dashed line.

Overcoming the Carbuncle
When computing 2D transcritical shocks such as hydraulic jumps undergoing a strong regime transition, a numerical instability in the discrete shock profile may appear. This is known as the Carbuncle and was first observed in simulations of air flow around blunt bodies by Peery and Imlay [1]. This anomaly is more likely to appear when using Cartesian grids, even though we can also find it for particular unstructured triangular grids in a weaker form.
The explanation for the Carbuncle is still not clear. Let us consider that the initial shock profile is given by a vertical straight line. After some time steps, as the flow adapts to the geometry of the solid body, the shock profile begins to curve. In a Cartesian grid, this bending of the profile is represented by part of the profile moving forward and other part moving backwards. This leads to the appearance of new RPs in the transverse direction that generate a cross-flow and induces a recirculation, which eventually produces the so-called Carbuncle instability of the discrete shock profile [8].
Most strategies to suppress the carbuncle instability are based on the detection of strong shocks and addition of artificial viscosity on such regions. To this end, the flux function is substituted by a different flux approximation provided by an incomplete Riemann solver. Incomplete solvers provide a poorer resolution of shear waves, hence introduce artificial viscosity. However, when using this technique, there is a risk that physical carbuncles could also be eliminated [8].
The easiest way to implement the aforementioned solution is to use the HLLS solver for the resolution of those RPs around the hydraulic jump. Unlike the ARoe solver, which is a complete solver and considers the complete eigenstructure of the problem, the HLLS solver only considers the two waves associated to the genuinely non-linear fields plus the wave associated to the source term. The contact wave associated to the shear velocity is neglected. This helps to damp the cross-flow as the shear waves are not resolved. The shock detection algorithm is a 2D extension of the one presented in [7].

Test case: fixing the carbuncle
This numerical experiment consist of a supercritical flow that hits a circular obstacle and generates an steady bow shock around the obstacle where the Carbuncle is likely to appear. The solution for the water surface elevation computed by the traditional ARoe solver and the modified ARoe+HLLS solver is presented in Figure 4. The map of the cells containing the shock is also included.

Conclusions
This work is devoted to the generation of numerical fixes that overcome numerical shockwave anomalies in the SWE, such as the slowly-moving shock that appears in presence of hydraulic jumps. Such anomaly has been taken for granted for a long time but in some particular cases may eventually ruin the solution. A theoretical framework of study for this anomaly has been provided and a 1D/2D spike reducing solver has been presented, based on a previous work by Zaide and Roe [6]. The resulting scheme has been exercised in a variety of scenarios and outperforms, by far, the traditional ARoe scheme. Moreover, convergence to the exact solution in presence of hydraulic jumps (measured with the L ∞ error norm), is achieved for the first time to the knowledge of the author. The proposed method will be extended to other SWE-based systems that model more complex phenomena such as geomorphological changes, interaction of several layers of fluid, granular flows, etc.