A Numerical Simulation of the Shallow Water Flow on a Complex Topography

Abstract

In current chapter, we have thoroughly described a numerical integration scheme of nonstationary 2D equations of shallow water. The scheme combines the smoothed particle hydrodynamics (SPH) and the total variation diminishing (TVD) methods, which are sequentially used at various steps of the combined SPH-TVD algorithm. The method is conservative and well balanced. It provides stable through calculations in presence of nonstationary “water-dry bottom” boundaries on complex irregular bottom topography including the transition of such a boundary between wet and dry bottom through the computational boundary. Multifarious tests demonstrate the effectiveness of the combined SPH-TVD scheme application for a solution of diverse problems of the engineering hydrology.

1. Introduction

The usage of the shallow water models (SWMs) allows to solve a wide range of engineering tasks related to the dynamics of surface waters (a seasonal flooding [1], a drain and rain flows [2]), the emergence and expansion of the marine nonlinear waves [3] (problems of tsunami impact on a shore, nonlinear waves formation due to earthquakes, and meteorological waves generation by an open ocean resonance [4]), and the inundation in a coastal area by storm surge [5]. The SWM modifications are effective for studying various geophysical problems such as the dynamics of the pyroclastic flows [6] and the riverbed processes including the sediment dynamics and the diffusion of pollutant particles in reservoirs. The multilayer models utilization significantly expands the opportunities of the shallow water approach [7]. The tasks associated with flooding research of river valleys or interfluves [1] are stood out among the hydrological problems. In the framework of the SWM, the important results have also been obtained for the atmospheric phenomena, meteorological forecasts, and global climate models [8, 9].

It should be noted that the SWM is also actively applied and developed for theoretical research of various cosmic gas flows. The shallow water approximation allows describing a whole series of astrophysical objects: the protoplanetary and circumstellar disks [10, 11], the accretion disks around the compact relativistic objects [12], the cyclonic movements in the giant planets’ atmospheres [13], and the spiral galaxies gas disk components [14]. The gravitational fields play the topography role in such astrophysical problems.

A lot of numerical methods have been proposed for the shallow water dynamics modeling employed for diverse tasks and conditions. Despite the fact that conservative methods of finite volume poorly describe stationary states, they allow correctly calculating shock waves and contact discontinuities. The latter problem could be overcome by the so-called well-balanced (WB) circuits [15, 16, 17].

Our main aim is through calculation for a flow with various Froude number within 0 ≤ Fr < 100 (Fr = u/c_g, where $c_g=\sqrt{\mathit{gH}}$ is an analog of the sound speed, H is the depth, u is the flow velocity, g is the specific gravitational force) in order to simulate subcritical (Fr <  < 1), transcritical (Fr~1), and supercritical (Fr >  > 1) flows. Highly heterogeneous terrain topography including vertical discontinuities and small-scale inhomogeneities at the computational domain boundary makes the calculations noticeably complicated. The latter leads to special quality requirements in numerical algorithms. Hence, the modern numerical schemes should simulate fluid movements along the dry bottom and correctly describe the interfaces between wet and dry bottom [1, 18, 19].

Among the numerous numerical methods solving shallow water equations (SWEs), the following methods should be mentioned: the discontinuous Galerkin method based on triangulation [8], the weighted surface-depth gradient method for the MUSCL scheme [18], and the modified finite difference method [20]. The so-called constrained interpolation profile/multimoment finite-volume method utilizing the shallow water approximation is developed to simulate geophysical currents on a rotating planet in spherical coordinate system ([9] and see the references there). The particle-mesh method demonstrates good opportunities for calculation of rotating shallow water [21]. As a rule, numerical schemes of the second-order accuracy give satisfactory results and allow to correctly solve a wide range of tasks for diverse applications [17]. Special attention should be focused on the numerical way of a source term setting, since in the case of discontinuous topography, it may induce an error at the shock wave front [22].

We consider the original CSPH-TVD (combined smoothed particle hydrodynamics—total variation diminishing) algorithm of numerical integration of the Saint-Venant equations. It accounts for the new modifications improving the computational properties of the scheme. A detailed description of the numerical scheme is the main aim of the chapter.

2. Mathematical models

In this section, the mathematical models for terrestrial hydrology, which accounts for the maximum number of physical and meteorological factors and can be described by the shallow water equations, are considered. The SWE or the Saint-Venant equations are fairly simple model describing the free surface dynamics of incompressible fluid (see the details in the review of [23]). The model assumes that vertical equilibrium in a medium exists at every moment of time. The conservation laws of mass and momentum in the integral form for a thin layer of moving substance with additional sources and forces are:

$$\frac{d}{\mathit{dt}}\underset{S\left(t\right)}{\iint }H\left(\mathbf{r},t\right)\mathit{dS}=Q\left(\mathbf{r},t\right),\frac{d}{\mathit{dt}}\underset{S\left(t\right)}{\iint }H\mathbf{v}\mathit{dS}=-g\underset{S\left(t\right)}{\iint }H\nabla \eta \mathit{dS}+\underset{S\left(t\right)}{\iint }H\mathbf{f}\mathit{dS},$$E1

where Q is the sources function, ∇ = e_x∂/∂x + e_y∂/∂y is the nabla operator, S(t) is the cross-sectional area of the “liquid particle”, r = xe_x + ye_y is the radius vector, v = ue_x + ve_y is the velocity vector, η(r, t) = H(r, t) + b(r) is the free surface elevation, b(r) is the bottom profile, f(r, t) = f_b + f_fr + f_Cor + f_w + f_s is the sum of the external forces, accounting for the bottom friction f_b, the internal friction (viscosity) f_fr, the Coriolis force f_Cor, the effect of atmospheric wind f_w, and the force f_s determined by the liquid momentum due to the action of the sources Q. The surface density of sources is σ(r, t) = σ⁽⁺⁾ + σ^(-) = dQ/dS, where σ⁽⁺⁾ and σ⁽⁻⁾ are the liquid sources (rain, melting snow, flows through the hydro-constructions, groundwater, etc.) and sinks (infiltration, evaporation, etc.), respectively.

The digital elevation model (DEM) provides the quality of numerical simulations of real hydrological objects to a significant extent. The DEM is determined by the height matrix b_ij = b(x_i, y_j) on numerical grid {x_i, y_j} (i = 1, …, N_x, j = 1, …, N_y). The DEM elaboration utilizes diverse geoinformation methods for the processing of spatial data obtained from various sources. A matrix of heights has been built in several stages. At the beginning stage, the remote sensing data are accounted by the function b(x_i, y_j). The river sailing directions and the actual water depth measurements allow to construct the DEM for large river beds (for example, for the Volga River and the Akhtuba River). To improve our DEM, the data on various small topography objects such as small waterways, roads, small dams, etc., should be included into the consideration. The numerical simulation results of the shallow water dynamics reveal flooding areas which may be compared with real observational data (both from the remote methods and our own GPS measurements). Such approach qualitatively improves the DEM as a result of the iterative topography refinement.

3. Numerical method СSPH-TVD

In current paragraph, the CSPH-TVD numerical method is thoroughly examined. The method combines classical TVD and SPH methods at various stages of numerical integration of hyperbolic partial differential equations and uses the particular benefits of both. The advantages of graphical processing units (GPUs) with CUDA acceleration for hydrological simulations are also discussed.

The computational domain has been covered by a fixed uniform grid with a spatial step h, while mobile “liquid particles” (hereinafter particles) are placed in the centers of cells. The time layers t_n have a nonuniform step τ_n = t_n − t_n − 1, where n denotes the index of the time layer. The vector space index i = (i, j) characterizes the radius vectors of the fixed centers of the cells ${\mathbf{r}}_{\mathbf{i}}^0=\left(x_i^0,y_j^0\right)$. The total number of cells and particles is N_p = N_x ⋅ N_y. At the initial time moment, the particles are placed in the centers of cells. The main characteristics of the particles are the volume $V_{\mathbf{i}}^n=\underset{S_{\mathbf{i}}\left(t\right)}{\iint }\mathit{HdS}$, momentum ${\mathbf{P}}_{\mathbf{i}}^n=\underset{S_{\mathbf{i}}\left(t\right)}{\iint }H\mathbf{v}\mathit{dS}$, coordinates of their centers ${\mathbf{r}}_{\mathbf{i}}^n$, and linear sizes ${\ell }_{\mathbf{i}}^n=\sqrt{S_{\mathbf{i}}^n}$.

The CSPH-TVD method main stages are the following:

1. The Lagrangian’s stage. It includes an application of the modified CSPH approach [1, 24]. At this stage, the changes of the integral characteristics and particles positions ‍($V_{\mathbf{i}}^n\to V_{\mathbf{i}}^{n+1}$, ${\mathbf{P}}_{\mathbf{i}}^n\to {\mathbf{P}}_{\mathbf{i}}^{n+1}$, and ${\mathbf{r}}_{\mathbf{i}}^n\to {\mathbf{r}}_{\mathbf{i}}^{n+1}$) due to the hydrodynamic and external forces acting on them are calculated. Figure 1 demonstrates the deformations of particles occurred during their motion, while the positions of particles are shifted relatively to the cells centers (it is assumed that ${\ell }_{\mathbf{i}}^{n+1}\ne h$). The difference between the CSPH-TVD method and the traditional SPH approach [25] consists in the fact that the particles at each time interval [t_n, t_n + 1] of the numerical time integration at the moment of time t_n are found in the initial state (${\mathbf{r}}_{\mathbf{i}}^n\equiv {\mathbf{r}}_{\mathbf{i}}^0$ and ${\ell }_{\mathbf{i}}^n=h$).

2. The Euler’s stage. At the Euler’s stage, a fixed grid is used to calculate the mass and momentum fluxes through the cells boundaries at time moment t_n + 1/2 = t_n + τ_n/2. The corresponding changes of the particles integral characteristics are proportional to the difference between the inflow and outflow fluxes. If the region shaded at Figure 1 is inside of the cell, then the substance flows into the cell, otherwise it flows out. The modified TVD approach [1, 24] and approximate solution of the Riemann’s problem [26] are applied to calculate the flows. At the end of the Euler’s stage, the particles return to their initial state.

Figure 1.

The OX-projection of the main stages scheme of the CSPH-TVD method. The dash-dotted lines x_i(t) correspond to the trajectories of the particles due to their displacement $x_i^n\to x_i^{n+1}$, the dashed black lines show the boundaries of cells. The particles boundaries deformed during the motion are shown by solid lines. The shaded regions correspond to a change in the particles integral characteristics caused by a flow of matter through the cells boundaries.

It should be noted that when the dry bottom regions are considered, the particles with zero depth and momentum are placed in the corresponding cells. For such particles, the Lagrangian’s stage is skipped. At the Euler’s stage, the mass and momentum fluxes flowing into the corresponding cell with zero depth are calculated in case of the water presence in the neighboring cells. The changes in the particles integral characteristics are determined too. Such a method permits to carry out an effective through calculation in the presence of nonstationary “water-dry bottom” boundaries in the computational domain.

3.1. The Lagrangian’s stage

Let us define the following quantities

$${\phi }_{\mathbf{i}}\left(t\right)=\frac{1}{h^2}\underset{S_{\mathbf{i}}\left(t\right)}{\iint }\phi \left(t,x,y\right)\mathit{\text{dxdy}}$$E2

being the analogous of the function mean values φ = {H, Hu, Hv} in cells in the finite volume approximation on a fixed grid. Taking into the account Eq. (2), the integral characteristics of the particles can be written in the form:

$$V_{\mathbf{i}}\left(t\right)=H_{\mathbf{i}}\left(t\right)h^2,{\mathbf{P}}_{\mathbf{i}}\left(t\right)={\left(H\mathbf{v}\right)}_{\mathbf{i}}\left(t\right)h^2$$E3

After substituting relations (2) and (3) into the system of Eq. (1), we obtain

$$\frac{d{\mathbf{U}}_{\mathbf{i}}}{\mathit{dt}}={\mathbf{\Phi }}_{\mathbf{i}}$$E4

where ${\mathbf{U}}_{\mathbf{i}}=\left(\begin{array}{c}\hfill H_{\mathbf{i}}\hfill \\ \hfill {\left(H\mathbf{v}\right)}_{\mathbf{i}}\hfill \end{array}\right)$, ${\mathbf{\Phi }}_{\mathbf{i}}=\left(\begin{array}{c}\hfill {\sigma }_{\mathbf{i}}\hfill \\ \hfill -{\mathit{gH}}_{\mathbf{i}}{\text{∇}}_{\mathbf{i}}\eta +H_{\mathbf{i}}{\mathbf{f}}_{\mathbf{i}}\hfill \end{array}\right)$, ${\text{∇}}_{\mathbf{i}}\eta ={\left.\text{∇}H\right|}_{\mathbf{r}={\mathbf{r}}_{\mathbf{i}}}+{\left.\text{∇}b\right|}_{\mathbf{r}={\mathbf{r}}_{\mathbf{i}}^0}$, σ_i = σ(r_i), f_i = f(r_i). The law of motion of a particle is

$$\frac{d{\mathbf{r}}_{\mathbf{i}}}{\mathit{dt}}={\mathbf{v}}_{\mathbf{i}}$$E5

where v_i = (Hv)_i/H_i is the velocity of the i-th particle.

For the numerical integration of the system of Eq. (4), an approximation for the spatial derivatives is required. In accordance with the SPH-approach [25], any characteristic of the medium φ and its derivatives ∇φ are replaced by their smoothed values in the flow area Ω:

$$\widehat{\phi }\left(t,\mathbf{r}\right)=\sum _{\mathbf{k}=\left(1,1\right)}^{\left(N_x,N_y\right)}\psi \left(t,{\mathbf{r}}_{\mathbf{k}}\right)W\left(|\mathbf{r}-{\mathbf{r}}_{\mathbf{k}}|,h\right),\nabla \widehat{\phi }\left(t,\mathbf{r}\right)=\sum _{\mathbf{k}=\left(1,1\right)}^{\left(N_x,N_y\right)}\psi \left(t,{\mathbf{r}}_{\mathbf{k}}\right)\nabla W\left(|\mathbf{r}-{\mathbf{r}}_{\mathbf{k}}|,h\right)$$E6

where W is the smoothing kernel function, h is the smoothing length, k is the vector index similar to the spatial index i. Accounting for Eq. (2), the quantity ψ can be determined as

$$\psi \left(t,{\mathbf{r}}_{\mathbf{k}}\right)=\underset{S_{\mathbf{k}}\left(t\right)}{\iint }\phi \left(t,\mathbf{r}\right)\mathit{dS}={\phi }_{\mathbf{k}}\left(t\right)h^2$$E7

Our modification of the SPH-method consists in the generalization of Eq. (7). In case of the traditional SPH-approach for the SWE [27] expression $\psi \left(t,{\mathbf{r}}_{\mathbf{k}}\right)=V_{\mathbf{k}}\left(t\right){\widehat{\phi }}_{\mathbf{k}}\left(t\right)/{\widehat{H}}_{\mathbf{k}}\left(t\right)$ is used instead of (7). Substituting the SPH-approximation (6) and (7) in the right-hand side of the expressions for system (4), we obtain

$${\mathbf{\Phi }}_{\mathbf{i}}\approx \left(\begin{array}{c}\hfill {\sigma }_{\mathbf{i}}\hfill \\ \hfill -{\mathit{gH}}_{\mathbf{i}}\sum _{\mathbf{k}}\left(H_{\mathbf{k}}{\nabla }_{\mathbf{i}}{\overline{W}}_{\mathbf{ik}}+b_{\mathbf{k}}{\nabla }_{\mathbf{i}}{\overline{W}}_{\mathbf{ik}}^0\right)+H_{\mathbf{i}}{\mathbf{f}}_{\mathbf{i}}\hfill \end{array}\right)$$E8

where ${\overline{W}}_{\mathbf{ik}}=h^{2}W\left(|{\mathbf{r}}_{\mathbf{i}}-{\mathbf{r}}_{\mathbf{k}}|,h\right)$, ${\overline{W}}_{\mathbf{ik}}^0=h^{2}W\left(|{\mathbf{r}}_{\mathbf{i}}^0-{\mathbf{r}}_{\mathbf{k}}^0|,h\right)$ are the functions computed in the cells centers and afterward used to approximate the gradient of the fixed bottom relief b(r). The approximation (8) ensures the conservatism of the CSPH-TVD method at the Lagrangian’s stage and the fulfillment of the WB-condition on the inhomogeneous topography.

There are three conditions superimposed on the smoothing kernels W:

1. the kernel is finiteness;

2. the normalization condition is fulfilled ${\displaystyle \underset{\Omega }{\iint }}W\left(|\mathbf{r}-{\mathbf{r}}^{\mathbf{\prime }}|,h\right)d{\mathbf{r}}^{\mathbf{\prime }}=1$;

3. the zero spatial step limit ${\displaystyle \underset{h}{\text{lim}}}W\left(|\mathbf{r}-{\mathbf{r}}^{\mathbf{\prime }}|,h\right)=\delta \left(|\mathbf{r}-{\mathbf{r}}^{\mathbf{\prime }}|,h\right)$ is the Dirac delta-function.

For the smoothing kernel W, spline functions of different orders or Gaussian distribution are used in the SPH hydrodynamic simulations [25, 28, 29]. The Monaghan’s cubic spline

$$W\left(q\right)=A_w\left\{\begin{array}{c}\hfill 1-1.5q^2+0.75q^3,\hfill \\ \hfill 0.25{\left(2-q\right)}^3,\hfill \\ 0,\hfill \end{array}\begin{array}{c}\hfill 0\le q\le 1\hfill \\ \hfill 1\le q\le 2\hfill \\ \hfill q>2\hfill \end{array}\right.,W^{\prime }\left(q\right)=-\frac{A_w}{h}\left\{\begin{array}{c}\hfill 3q-2.25q^2,\hfill \\ \hfill 0.75{\left(2-q\right)}^2,\hfill \\ 0,\hfill \end{array}\begin{array}{c}\hfill 0\le q\le 1\hfill \\ \hfill 1\le q\le 2\hfill \\ \hfill q>2\hfill \end{array}\right.$$E9

is the most commonly applied approximation, where q =  ∣ r_i − r_k ∣ /h is the relative distance between the particles, $A_w=\left(\frac{2}{3h},\frac{10}{7\pi h^2},\frac{1}{\pi h^3}\right)$ is the normalization constant for the 1D, 2D and 3D cases, respectively. In the CSPH-TVD method, the smoothing kernel W approximates only the gradients of physical quantities in Eq. (8). Therefore, the value of the normalization constant may be corrected to increase the accuracy of the numerical solution. For example, it was done by the authors of [1] when for the two-dimensional case, the normalization constant was chosen to be $A_w\approx 0.987\left(\frac{10}{7\pi h^2}\right)$. The smoothing kernel (9) can lead to unphysical clustering of particles in the simulation of hydrodynamic flows by the SPH method; therefore, an approximation for the smoothing kernel of the pressure gradient is applied in form [29]:

$$W\left(q\right)=B_w\left\{\begin{array}{c}\hfill {\left(2-q\right)}^3,\hfill \\ \hfill 0,\hfill \end{array}\right.\begin{array}{c}\hfill 0\le q\le 2\hfill \\ \hfill q>2\hfill \end{array},W^{\prime }\left(q\right)=-\frac{B_w}{h}\left\{\begin{array}{c}\hfill 3{\left(2-q\right)}^2,\hfill \\ 0,\hfill \end{array}\right.\begin{array}{c}\hfill 0\le q\le 2\hfill \\ \hfill q>2\hfill \end{array},$$E10

where $B_w=\left(\frac{1}{8h},\frac{5}{16\pi h^2},\frac{15}{64\pi h^3}\right)$ is the normalization constant. Figure 2 shows the comparison between two smoothing kernels (9) and (10). As seen from Figure 2b for q < 0.7, the characteristic feature of the smoothing kernel (9) is the impetuous decrease of the derivative value to zero as the particles approach each other. At small q values, this feature leads to a significant weakening of the hydrodynamic repulsion force between the particles and appearance of nonphysical clusters of particles.

Figure 2.

The spatial distributions of smoothing kernels (a) and their derivatives (b). The smoothing kernels (9) and (10) are shown by solid and dashed lines, correspondingly.

For the numerical integration of the set of differential Eqs. (4) and (5), the method of a predictor-corrector type of second-order accuracy (so-called leapfrog method) is employed. The main steps of the leapfrog method for the Lagrange stage of CSPH-TVD method are:

1. At the predictor step, the water depth H_i and velocity v_i are calculated at time moment t_n + 1 = t_n + τ_n:

   $${\mathbf{U}}_{\mathbf{i}}^{\ast }\left(t_{n+1}\right)={\mathbf{U}}_{\mathbf{i}}\left(t_n\right)+{\tau }_n{\mathbf{\Phi }}_{\mathbf{i}}\left({\mathbf{r}}_{\mathbf{i}}\left(t_n\right),{\mathbf{U}}_{\mathbf{i}}\left(t_n\right)\right)$$E11

2. The Particles’ spatial positions r_i are determined at time t_n + 1:

   $${\mathbf{r}}_{\mathbf{i}}\left(t_{n+1}\right)={\mathbf{r}}_{\mathbf{i}}\left(t_n\right)+\frac{{\tau }_n}{2}\left({\mathbf{v}}_{\mathbf{i}}\left(t_n\right)+{\mathbf{v}}_{\mathbf{i}}^{\ast }\left(t_{n+1}\right)\right)$$E12

3. During the corrector step, the water depth H_i and velocity v_i values are recalculated (t_n + 1):

   $${\tilde{\mathbf{U}}}_{\mathbf{i}}\left(t_{n+1}\right)=\frac{{\mathbf{U}}_{\mathbf{i}}\left(t_n\right)+{\mathbf{U}}_{\mathbf{i}}^{\ast }\left(t_{n+1}\right)}{2}+\frac{{\tau }_n}{2}{\mathbf{\Phi }}_{\mathbf{i}}\left({\mathbf{r}}_{\mathbf{i}}\left(t_{n+1}\right),{\mathbf{U}}_{\mathbf{i}}^{\ast }\left(t_{n+1}\right)\right)$$E13

To reduce the rounding errors during the numerical implementation of the algorithm in recurrent relations (11)–(13), the value of the relative displacement of particle ${\mathbf{\xi }}_{\mathbf{i}}\left(t\right)={\mathbf{r}}_{\mathbf{i}}\left(t\right)-{\mathbf{r}}_{\mathbf{i}}^0$ (here ${\mathbf{r}}_{\mathbf{i}}^0$ is the initial position) should be used instead of r_i(t).

3.2. The Euler’s stage

At this stage, the relationship between the values of ${\tilde{\mathbf{U}}}_{\mathbf{i}}^{n+1}$ and ${\mathbf{U}}_{\mathbf{i}}^{n+1}$ is determined at time moment t_n + 1. According to the CSPH-TVD approach, the changes in the particles integral characteristics at their return to the initial state ${\mathbf{r}}_{\mathbf{i}}^0$ are computed. In order to do it, the following expression

$$\frac{d}{\mathit{dt}}\underset{S_{\mathbf{i}}\left(t\right)}{\iint }\mathbf{U}\mathit{\text{dxdy}}=\underset{S_{\mathbf{i}}\left(t\right)}{\iint }\left(\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}}{\partial x}+\frac{\partial \mathbf{G}}{\partial y}\right)\mathit{\text{dxdy}},\text{where}\mathbf{F}=\left(\begin{array}{c}\hfill \mathit{Hv}\hfill \\ \hfill {\mathit{Hu}}^2\hfill \\ \hfill \mathit{Hvu}\hfill \end{array}\right),\mathbf{G}=\left(\begin{array}{c}\hfill \mathit{Hv}\hfill \\ \hfill \mathit{Huv}\hfill \\ \hfill {\mathit{Hv}}^2\hfill \end{array}\right)$$

should be integrated over time from t_n to t_n + 1 . Taking into the account Eq. (2), we obtain

$${\mathbf{U}}_{\mathbf{i}}^{n+1}={\tilde{\mathbf{U}}}_{\mathbf{i}}^{n+1}-\frac{{\tau }_n}{h}\left({\mathbf{F}}_{i+1/2,j}^{n+1/2}-{\mathbf{F}}_{i-1/2,j}^{n+1/2}+{\mathbf{G}}_{i,j+1/2}^{n+1/2}-{\mathbf{G}}_{i,j-1/2}^{n+1/2}\right)$$E14

where ${\mathbf{F}}_{i\pm 1/2,j}^{n+1/2}$ and ${\mathbf{G}}_{i\pm 1/2,j}^{n+1/2}$ are the average values (in the interval [t_n, t_n + 1]) of mass and momentum fluxes through the cell boundaries in x and y directions, respectively. To determine these fluxes, we use the TVD-approach [30] and approximate methods of the Riemann’s Problem (RP) solving for an arbitrary discontinuity decay (Lax-Friedrichs (LF) method, Harten-Lax-van Leer (HLL) method [31]). The Riemann’s problem is solved separately for each boundary of the Euler’s cells.

The values of the flow parameters to the left U^L and right U^R of the considered boundary are the initial conditions for RP. These quantities define magnitude of the jump and can be obtained by a piecewise polynomial reconstruction of the function U(t, r). We limit the consideration by the piecewise linear reconstruction providing the second order of accuracy by space coordinate for a numerical scheme [32]. The modification of the TVD-approach used in the CSPH-TVD method is concerned with the reconstruction of the value of U(t, r) relatively to the position of the particles at time t_n + 1/2 at distance ξ_i(t_n + 1/2) from the centers of the cells. Thus, the piecewise linear profile of the function U(t, r) inside of the i-th cell at time t_n + 1/2 in the projection onto the Ox axis (Figure 3) has the following form

$$\mathbf{U}\left(t_{n+1/2},x\right)={\tilde{\mathbf{U}}}_{i,j}^{n+1/2}+\left(x-x_i^{n+1/2}\right){\mathrm{\Theta }}_{i,j}^{n+1/2},x\in \left[x_{i-1/2}^0,x_{i+1/2}^0\right]$$E15

Figure 3.

The projection onto the Ox axis of the piecewise linear reconstruction of the function U(t, r) with the respect to the mass centers of the particles for the CSPH-TVD scheme at time t_n + 1/2.

where ${\mathrm{\Theta }}_{i,j}^{n+1/2}$ is the vector of slopes (angular coefficients), $x_i^{n+1/2}=x_i^0+{\xi }_i^{n+1/2}$. Taking into the account Eq. (15), the expressions for U^Land U^Rat the boundary $x_{i+1/2}^0$ can be written as

$${\mathbf{U}}_{i+1/2}^L={\tilde{\mathbf{U}}}_{i,j}^{n+1/2}+\left(\frac{h}{2}-{\xi }_i^{n+1/2}\right){\mathrm{\Theta }}_{i,j}^{n+1/2},{\mathbf{U}}_{i+1/2}^R={\tilde{\mathbf{U}}}_{i+1,j}^{n+1/2}-\left(\frac{h}{2}+{\xi }_{i+1}^{n+1/2}\right){\mathrm{\Theta }}_{i+1,j}^{n+1/2}$$E16

The slopes of the piecewise linear distribution (16) should satisfy the TVD-condition [30]. To fulfill this condition, the limiter function

$${\mathrm{\Theta }}_{i,j}^{n+1/2}=\mathrm{\Lambda }\left(\frac{{\tilde{\mathbf{U}}}_{i+1}^{n+1/2}-{\tilde{\mathbf{U}}}_i^{n+1/2}}{h+{\xi }_{i+1}^{n+1/2}-{\xi }_i^{n+1/2}},,,\frac{{\tilde{\mathbf{U}}}_i^{n+1/2}-{\tilde{\mathbf{U}}}_{i-1}^{n+1/2}}{h+{\xi }_i^{n+1/2}-{\xi }_{i-1}^{n+1/2}}\right)$$

has to be used. The analysis shows that the limiters of the minmod [33], van Leer [32], and van Albada et al. [34] satisfy the TVD-condition and suppress the nonphysical oscillations of the numerical solution in the vicinity of discontinuities.

3.3. The stability condition and parallel implementation on GPUs

For the stability of the numerical scheme CSPH-TVD the following conditions have to be fulfilled during the integration time τ_n: (1) at the Lagrangian’s stage the shift of particles’ center of mass should not exceed h/2 relatively to the initial position; (2) at the Euler’s stage the perturbations should not propagate over a distance greater than the cell size h. These requirements provide the numerical stability condition for the time step τ_n of the CSPH-TVD algorithm:

$${\tau }_n=Kmin\left(\frac{h}{2\underset{\mathbf{i}}{max}\mid {\mathbf{v}}_{\mathbf{i}}^n\mid },\frac{h}{\underset{\mathbf{i}}{max}\left(|{\mathbf{v}}_{\mathbf{i}}^n|+\sqrt{{\mathit{gH}}_{\mathbf{i}}^n}\right)}\right)$$E17

where 0 < K < 1 is the Courant number.

The CSPH-TVD method demonstrates a good ability for parallelization on graphics processors [35]. Figure 4 shows the computer implementation of the CSPH-TVD method on GPUs. It indicates the execution sequence of CUDA kernels at various stages of our algorithm.

Figure 4.

The flow diagram for each of the calculation module: K1 determines the presence of water in the CUDA block; K2 calculates the forces at the time moment t_n at the Lagrangian’s stage; K3 calculates the time step t_n + 1; K4 calculates the new positions of the particles and their integral characteristics at time t_n + 1/2; K5 defines the forces on the time layer t_n + 1/2; K6 calculates the positions of the particles and their integrated characteristics for the next time layer t_n + 1; K7 calculates the flux of physical quantities through the cells boundaries at time moment t_n + 1/2; and K8 determines the final hydrodynamic parameters at time moment t_n + 1 (see details in Ref. [35]).

4. The set of test problems verifying the numerical models

In the paragraph, the basic set of various tests for one- and two-dimensional problems suitable for revelation of positive and negative sides of the numerical shallow water models is reviewed. According to Toro [26], there are four main test types appropriate to evaluate the numerical solutions (Figure 5).

1. The comparison with the exact solutions for several 1D problems.

   1. Dam-break waves propagating on the wet bed [16, 36] with the initial distribution of the water depth H(x ≤ 0) = H₀₁ and H(x > 0) = H₀₂, (H₀₁ > 0, H₀₂ > 0). This solution is an analog of the pressure jump decay in an ordinary hydrodynamics. It also contains a hydraulic shock (an analog of a shock wave) and a rarefaction wave for which the exact solutions of the RP for a self-similar type f(x, t) = f(ξ) ≡ f(x/t) exist [37]. The bora (hydraulic jump) front for the CSPH-TVD scheme is smeared into three cells and does not contain dispersion oscillations typical for sufficiently good numerical schemes.

   2. Dam-break waves propagating on the dry bed [16]: H(x ≤ 0) = H₀₁ and H(x > 0) = 0. This is an important test imposing special requirements on a numerical scheme for a correct description of the moving boundary between the water and dry bottom. It also allows comparison of the numerical result with exact solution of the RP [37].

   3. The dam-break problem in a sloping channel is more complicated test comparing the dynamics of numerical wetting and drying fronts with the analytic solutions [18].

   4. The two shock waves collision is described in the terms of the RP and contains a rarefaction wave [36].

   5. It is possible to check the coincidence of the numerical and exact solutions for the flow over the bed profiles of different shapes (Figure 5). The transcritical mode with and without a shock for the flow over bed profile in the form of a parabolic bump (see line 3, Figure 5) [16, 17, 18, 38] or triangular obstacle with a break (see line 1) [24] have been studied.

   6. The more complicated bottom profile (see line 2 in Figure 5) helps identifying the best advantages of numerical schemes [22]. We have exact stationary solutions on a stepped bottom for subcritical (Froude number Fr < 1) and supercritical regimes (Fr > 1) [24, 39].

   7. To verify the properties of well-balanced numerical model, the oscillating lake with a heterogeneous bottom [15] or the basin with a stationary fluid with an island in the middle of the numerical domain are suitable to modeling [19]. The study of a propagation of small perturbation in stationary water with a nonuniform bottom allows to determine the quality of the well-balanced scheme [16, 19] and guarantees the absence of a numerical storm produced by the scheme [40].

   8. The same problem of assessing the quality of the WB-scheme is solved by studying the oscillations of the liquid level on the 1D parabolic bottom [24], when the profile of the oscillating free surface of the liquid should remain flat at any time. The two-dimensional analog of such problem is the oscillations in the frictionless paraboloid basin b(x, y) = b_L(x² + y²)/L², for which the free surface remains flat [18, 41]. The analytical solution for the Thacker test with curved water surface in the paraboloidic bottom basin is described in [41]. The initial circular paraboloidic water volume oscillates periodically changing the shape of the free surface in a complex manner (see Figure 5, lines 5a and 5b corresponds to the bottom profile and the initial shape of the water volume) [18].

2. The second approach is based on the tests with reliable numerical solutions for 1D equations. A simple and effective two-dimensional test is circular dam-break [16, 36], since for the exact solution can be obtained from the one-dimensional radial equation in the polar coordinate system [18]. It is necessary to distinguish the circular dam break on the wet and dry bed.

3. The third direction is represented by tests using other numerical schemes.

4. Eventually, the numerical SWM solutions can be compared with the observables data or 3D numerical solutions.

   1. Dam-break flow over the initially dry bed with a bottom obstacle is a tough test for a simple SWM due to the formation of negative wave propagation [42].

   2. The experimental results of the dam-break on a dry bed channel with varying width are used to examine numerical models [16]. The laboratory measurements of flow parameters for the dam-break in the nonprismatic converging-diverging channel are described in Ref. [38].

   3. The so-called right-angled open channel junction test gives a complex two-dimensional flow structure which is highly inhomogeneous in the cross section of the main channel [43]. Such a kind of comparison with essentially 2D currents is an important stage in testing of the numerical model [16].

   4. The dam-break flow over a triangular obstacle contains a large number of characteristic features [38], including the downstream over the dry channel bed, the smooth rarefaction wave, the upstream reflected shock, and the downstream shock after passing the triangular obstacle. The formation of a reflected shock wave from the boundary wall passing back through the hump causes complex nonlinear oscillations in such a basin.

   5. An analytical solution for steady flow case with spatially varying width, bed, and friction slope of special kind is known [38, 44]. There is a singular point in the flow in the transition from the inclined profile to the shallow bottom. A fine structure of the solution in the vicinity of this point is a good test for numerical models.

   6. The use of various software, for example, FLOW-3D, FLOW-3D-SWM, MIS2D, is a fast and effective way for the verification of new numerical models of shallow water [45].

Figure 5.

The bottom profiles with known exact solutions of SWE.

Figure 6 shows the results of the numerical simulations of shear flow instability (Kelvin-Helmholtz instability) in shallow water using the CSPH-TVD method. The characteristic stages of the instability development of the tangential velocity discontinuity in shallow water (see Figure 6a) are obvious: (I) the formation of the eigenmode; (II) the linear stage of instability development, when perturbations increase according to the law ∝ exp(t/T)(T ≈ 9.6 for Fr = 1); (III) the nonlinear mode of perturbations evolution, when the growth of the perturbation amplitude ceases and the characteristic vortex structures are formed (so-called cat eyes, Figure 6b). Thus, in the CSPH-TVD method, there is a lack of important disadvantage of classical SPH algorithm, which does not allow correctly modeling the tangential discontinuity.

Figure 6.

The problem of the tangential velocity discontinuity is represented by: (a) the time dependence of the minimum (shown by solid line) and maximum (shown by dashed line) of the amplitude of the water depth disturbances H_w(t) = H(t) − H₀  (H₀ ≡ H(0)); (b) the spatial distribution of the vorticity of the velocity vector field ω = ∇ × v/2 at t = 110. The dashed line shown the position of the tangential discontinuity.

5. Engineering applications

We have developed software to solve some engineering applications utilizing the parallel implementation of the CSPH-TVD method on GPUs.

1. The hydrological regimes of the spring flooding in the Volga-Akhtuba interfluve have been studied thoroughly [1, 35]. The optimal hydrographs of the Volga Hydroelectric Power Station have been constructed. And finally a new approach optimizing the hydrotechnical projects has been developed [46, 47, 48].

2. Our numerical experiments have reproduced the dynamics of the catastrophic flooding wave in the city of Krymsk area (Russia, Caucasus Mountains) in July 2012, leading to massive losses of life. A number of features of the hydrological regime during the flash flood of 2012 [49], associated with the landscape and the distribution of rainfall has been revealed.

3. Figure 7 shows the more detailed results of the numerical simulations of the tsunami formation in the Pacific Ocean during an earthquake off the coast of Japan in 2011. A bilateral hydraulic shock (tsunami) is formed due to the displacement of tectonic plates and then the waves propagate to both sides of the discontinuity. The tsunami reaches the shores of Japan in 20–40 min after the earthquake, while the wave height H_w reaches 10–20 m at the coast line.

Figure 7.

Numerical modeling of the tsunami occurred in 2011 off the coast of Japan: (a) the digital terrain model, where the size of the simulation area is 1500 × 1500 cells; the yellow line shows the area of tsunami formation; (b) the 3D wave structure 19 min after the earthquake; frames (c) and (d) are the tsunamis at different moments of time.

6. Conclusion

We have described in detail the numerical method for solving the equations of hydrodynamics in the approximation of the SWM. Our method is a hybrid scheme successfully combining positive properties of Euler’s and Lagrange’s approaches. The CSPH-TVD method allows to model nonstationary flows on a complex heterogeneous bottom relief, containing kinks and sharp jumps of the bottom profile. The numerical scheme is conservative, well balanced and provides a stable through calculation in the presence of non-stationary “water-dry bottom” boundaries on the irregular bottom relief, including the transition through the computing boundary.

On the SPH stage, various smoothing cores can be applied, as well as various TVD-delimiters and methods for the RP solution depending on the features of the problem being solved. The scheme has the second order of convergence on smooth solutions and the first order on discontinuities that corresponds to the accuracy of the Godunov-type schemes. In the case of a non-uniform topography the CSPH-TVD method requires less computational resources than the Godunov’s type schemes, when the WB-condition is necessary to fulfill for the numerical scheme. Comparing with various SPH-method modifications, the numerical CSPH-TVD scheme has higher accuracy and computational speed for the same number of particles; it is less dissipative and better balanced.

Acknowledgments

The first author has been supported by the Ministry of Education and Science of the Russian Federation (government task No. 2.852.2017/4.6). The second author is thankful to the RFBR (grants 16-07-01037, 15-45-02655, 16-45-340152 and 16-48-340147).