FMRIJ
Fluid Mechanics Research International Journal
Review Article
Volume 1 Issue 3  2017
Simulation of the Propagation of Tsunamis in Coastal Regions by a TwoDimensional NonHydrostatic Shallow Water Solver
Costanza Arico*
Department of Civil Engineering, University of Palermo, Italy
Received: November 14, 2017  Published: November 27, 2017
*Corresponding author:
Costanza Arico, Department of Civil Engineering, Environmental, Aerospace, Materials, University of Palermo, Avenue of Sciences, 90128, Palermo, Italy, Email:
Citation:
Arico C (2017) Simulation of the Propagation of Tsunamis in Coastal Regions by a TwoDimensional NonHydrostatic Shallow Water Solver.
Fluid Mech Res Int 1(3): 00011. DOI:
10.15406/fmrij.2017.01.00011
Abstract
Due to the enormous damages and losses of human lives in the inundated regions, the simulation of the propagation of tsunamis in coastal areas has received an increasing interest of the researchers. We present a 2D depthintegrated, nonhydrostatic shallow waters solver to simulate the propagation of tsunamis, solitary waves and surges in coastal regions. We write the governing continuity and momentum equations in conservative form and discretize the domain with unstructured triangular Generalized Delaunay meshes. We apply a fractionaltimestep procedure, where two problems (steps) are consecutively solved. In the first and in the second step, we hypothesize a hydrostatic and a nonhydrostatic distribution of the pressure, respectively. Several literature models, which solve the same set of equations, are based on a fractionaltimestep procedure. In the hydrostatic step of these literature solvers, the flow field does not satisfy the depthintegrated continuity equation, since the momentum equations are solved independently of the continuity equation. In both steps of the proposed numerical scheme, the computed flow field satisfies the depthintegrated continuity equation discretized in each computational cell. This is obtained by solving together the governing equations, according to the different numerical procedures adopted in the steps of the algorithm. Wet/dry problems are implicitly embedded in the proposed numerical solver. We present several model applications where analytical or measured reference solutions are available. The computational effort required by the proposed procedure is investigated.
Keywords: Shallow waters; Nonhydrostatic pressure; Unstructured mesh; Wetting/drying; Tsunami propagation; Long waves; Voronoi cells; RungeKutta method; Galerkin scheme; Manning equation; Dirichlet condition; OpenFOAM
Abbrevations
BTMs: BoussinesqType Models; NLSWEs: Nonlinear Shallow Water Equations; FV: Finite Volumes; WD: Wetting/Drying; PDEs: Partial Differential Equations; DC: Diffusive Correction; MAST: Marching in Space and Time; CC: Convective Correction; GD: Generalized Delaunay; ODEs: Ordinary Differential Equations; CP: Convective Prediction
Introduction
During the last years, great effort has been addressed by several authors to simulate the propagation of solitary waves/tsunamis, tides or surges, due to the tremendous damages and losses of human lives in the inundated rural and residential areas. Tsunamis are sea waves usually generated by undersea landslides and earthquakes. They can be regarded as long/solitary waves with small amplitude and long wavelength, travelling with high speed over long distances. Approaching the coast, their amplitude increases, becoming potentially destructive. The propagation of tsunami in coastal regions can be studied by investigating the shoaling and breaking of solitary waves over inclined bottoms [1]. Chanson [2] asserted that the front of tsunamis over dry plains becomes a shock wave, and presented a similarity between the propagation of the tsunamis over dry coastal areas and the classical dambreak problems. Generally, 3D simulations of the processes described as above, e.g., by RANS models [3], or Smoothed Particle Hydrodynamics methods [4], or Volume of Fluids methods [5], require very high computational costs. For these reasons, a significant amount of literature based on depthintegrated equations has been published for simulations of long waves/tsunamis propagation. Generally, the two modeling approaches proposed in literature are the BoussinesqType Models (BTMs) and the NonHydrostatic Nonlinear Shallow Water Equations (NLSWEs) models.
The classical formulation of the BTMs [6] involves weak nonlinearity and dispersion. Highorder BTMs have been proposed to overcome these difficulties, but suffer from complex numerical discretization, which leads to high computational efforts [7,8]. Generally, the BTMs also suffer from the inclusion, in the governing equations, of extra terms, based on empirical considerations, introduced for the simulation of wave breaking with the associated energy dissipation, and wetting/drying transition [8].
Hereafter, we call "hydrostatic NLSWEs" and "nonhydrostatic NLSWEs" (or "dynamic NLSWEs") respectively the NLSWEs with hydrostatic and nonhydrostatic (or dynamic) pressure distribution. Due to simplicity and accuracy, the hydrostatic NLSWEs, are largely applied for simulating wave processes in rivers, lakes, estuaries. The hydrostatic NLSWEs lack frequency dispersion and are not able in simulating some dispersive features of long waves/tsunamis propagation [813] as well as freesurface undulations at front/tail of bores and shock waves of tsunamis and dambreaks [14,15]. Essentially two aspects allow reproducing the dispersive effects in the NLSWEs, 1) the decomposition of the pressure into its hydrostatic and dynamic components and 2) the inclusion of the vertical momentum equation in the governing equations [16]. Two fractionalstep methodologies are usually applied to perform the decomposition of the pressure, the pressure projection and the pressure correction method, respectively [17].
One of the most demanding topics of the NLSWEs models is the simulation of the wave breaking [18]. As widely recognized, an accurate simulation of discontinuous flows (e.g., bores or jumps), is obtained if the numerical model solves the weak form of the hydrostatic NLSWEs written in conservative form [1921]. Numerous shockcapturing procedures, in the framework of the Finite Volumes (FV), Finite Elements (FE) and Finite Differences (FD) schemes, have been proposed for the solution of the hydrostatic NLSWEs [22,23]. Only a few nonhydrostatic models are equipped with this shockcapturing feature [18,24,25].
In the NLSWEs applications involving fast transitions between wet and dry surfaces, e.g., coastal flooding and runup of tsunamis, an accurate modeling of wetting/drying (WD) processes plays an important role. Several WD numerical techniques have been proposed [2628], but generally these are affected by several drawbacks, e.g., the mass imbalance and the high computational costs due to additional relationships to include in the numerical procedure at the transition wet/dry surface [27].
In the present paper, we propose a nonhydrostatic NLSWEs solver to simulate the propagation of long waves, tsunamis and surges processes. We write the governing equations in conservative form and the total pressure is decomposed into its hydrostatic and dynamic components. A hydrostatic and a nonhydrostatic (or dynamic) problem (or step), consecutively solved, are embedded in a general fractionaltimestep procedure. The domain is discretized with unstructured triangular Generalized Delaunay meshes [4]. Generally, the triangular meshes are suitable to fit irregular natural boundaries and arbitrary geometries. The hydrostatic step, where the dynamic pressure terms are neglected, is solved by a predictioncorrection procedure, recently proposed in [3033]. Unlike studies in [3033], where a cellcentered formulation is applied, in the present solver we adopt a nodecentered formulation.
We organize the paper as follows. We present the governing equations and the general formulation of the numerical solver in sections 2 and 3, respectively. The hydrostatic and dynamic steps are disentangled in sections 4 and 5, respectively, and the boundary conditions, in section 6. The properties of the proposed model (e.g. the proof of the local and global mass balance, the WD properties, the study of the phase speed, ... ) are outlined in section 7. Finally, in section 8, we present several model applications where reference analytical solutions or measured data are available. The proposed test cases highlight the skill of the model in simulating challenging features of the propagation of long, solitary waves and tsunami, as the wave breaking, runup, and wetting/drying problems. In the same section we also investigate the computational burden.
The Governing Equations
We derive the governing equations from the continuity and Reynolds equations of an incompressible fluid (water),
$\nabla \cdot u=0$
(1),
$\frac{\partial u}{\partial t}+\nabla \cdot \left(uu\right)+\frac{1}{\rho}\nabla p+g\cdot \nabla z\upsilon {\nabla}^{2}u=0$
(2),
where t is the time, u is the velocity vector (u, v, and w are its components along x, y and z directions, respectively), g is the gravitational acceleration vector,
$g=\text{}{\left(0\text{}0g\right)}^{T}$
, with g the norm of the gravitational acceleration (notation (*)^{T} indicates the transposed of vector (*)),
$\rho $
is the density of the water, p is the pressure and
$\upsilon $
is the kinematic viscosity coefficient.
We assume a physical domain bounded along the vertical directions by the free surface H(x, y, t) and the bed surface (or bottom) z_{b}(x, y, t), and h is the water depth (h = H  z_{b}) (Figure 1). We operate under the hypothesis of fixed bed condition, i.e.,
$\partial {z}_{b}/\partial t=0$
, which implies
$\partial H/\partial t=\partial h/\partial t$
. The total derivatives of the bottom and free surfaces lead to the kinematic boundary conditions [34]
${w}_{b}=\frac{D{z}_{b}}{Dt}={u}_{b}\frac{\partial {z}_{b}}{\partial x}+{v}_{b}\frac{\partial {z}_{b}}{\partial y}$
and
${w}_{s}=\frac{DH}{Dt}=\frac{\partial H}{\partial t}+{u}_{s}\frac{\partial H}{\partial x}+{v}_{s}\frac{\partial H}{\partial y}$
(3),
and the subindices b and s mark the values at the bottom and at the free surface, respectively. The (total) pressure p is split into its hydrostatic and dynamic components [16],
$p\left(x,t\right)=\underset{\text{hydrostaticcomponent}}{\underbrace{\rho g\left(H\left(x,t\right)z\right)}}+\underset{\text{dynamiccomponent}}{\underbrace{\widehat{q}\left(x,t\right)}}\text{}x={\left(\begin{array}{ccc}x& y& z\end{array}\right)}^{T}$
(4).
Figure 1: Definition sketch.
At z = H we set
$p=\widehat{q}=\text{}0$
, i.e., zero value of the atmospheric pressure.
If we integrate Eqs. (1)(2) over the water depth h, from z_{b} to H, we derive the depthintegrated NLSWEs [35]. If we merge Eqn (3) in the continuity equation, integrated along the water depth, it results
$\frac{\partial h}{\partial t}+\frac{\partial \left(Uh\right)}{\partial x}+\frac{\partial \left(Vh\right)}{\partial y}=0$
(5).
If we hypothesize that
$\widehat{q}$
changes linearly along z, from value 0 at z = H to value
${\widehat{q}}_{b}$
at z = z_{b}, dividing Eq. (2) by
$\rho $
, the depthintegrated momentum equations are
$\begin{array}{c}\underset{\begin{array}{c}\text{localderivative}\\ \text{term}\end{array}}{\underbrace{\frac{\partial \text{\hspace{0.05em}}\left(Uh\right)}{\partial t}}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\underset{\text{convectiveterms}}{\underbrace{\frac{\partial \text{\hspace{0.05em}}\left({U}^{2}h\right)}{\partial x}\text{\hspace{0.05em}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial \text{\hspace{0.05em}}\left(UVh\right)}{\partial y}}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}+\underset{\begin{array}{c}x\text{component}\\ \text{ofdepthintegrated}\\ \text{hydrostaticpressure}\end{array}}{\underbrace{gh\frac{\partial h}{\partial x}}}+\underset{\begin{array}{c}x\text{component}\\ \text{ofhydrostatic}\\ \text{pressureatbottm}\end{array}}{\underbrace{gh\frac{\partial {z}_{b}}{\partial x}}}\text{\hspace{0.17em}}+\underset{\begin{array}{c}x\text{component}\\ \text{ofdepthintegrated}\\ \text{dynamicpressure}\end{array}}{\underbrace{\frac{1}{2}\text{\hspace{0.17em}}\frac{\partial \text{}\left({q}_{b}h\right)}{\partial x}}}\text{\hspace{0.17em}}+\underset{\begin{array}{c}x\text{component}\\ \text{ofdynamic}\\ \text{pressureatbottm}\end{array}}{\underbrace{{q}_{b}\text{\hspace{0.05em}}\frac{\partial {z}_{b}}{\partial x}}}\text{\hspace{0.17em}}+\\ \text{\hspace{0.17em}}\underset{\begin{array}{c}x\text{component}\\ \text{ofbottomstress}\end{array}}{\underbrace{\frac{g{n}^{2}\text{}\left(Uh\right)\text{\hspace{0.05em}}\sqrt{{\left(Uh\right)}^{2}\text{\hspace{0.05em}}+{\left(Vh\right)}^{2}}}{{h}^{7/3}}}}=0\end{array}$
(6.a),
$\frac{\partial \text{}\left(Vh\right)}{\partial t}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial \text{}\left(UVh\right)}{\partial x}\text{\hspace{0.17em}}+\frac{\partial \text{}\left({V}^{2}h\right)}{\partial y}\text{\hspace{0.17em}}+gh\frac{\partial h}{\partial y}\text{\hspace{0.17em}}+gh\frac{\partial {z}_{b}}{\partial y}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{1}{2}\frac{\partial \text{}\left({q}_{b}h\right)}{\partial y}\text{\hspace{0.17em}}+{q}_{b}\text{\hspace{0.17em}}\frac{\partial {z}_{b}}{\partial y}+\text{\hspace{0.17em}}\frac{g{n}^{2}\text{\hspace{0.05em}}\left(Vh\right)\text{}\sqrt{{\left(Uh\right)}^{2}+{\left(Vh\right)}^{2}}}{{h}^{7/3}}=0$
(6.b),
$\frac{\partial \left(Wh\right)}{\partial t}={q}_{b}$
(6,c),
and the components of the depthintegrated velocity vector along x, y and z directionsare U, V and W, respectively, and the corresponding components of the specific flow rate vector (i.e. the flow rate per unitary width) are Uh, Vh and Wh,
${q}_{b}={\widehat{q}}_{b}/\rho $
and n is the Manning bed roughness coefficient. The Leibniz rule is applied to perform the integration along the depth of the terms proportional to q_{b} [11,18]. We refer the reader to [36] for further details of the depthintegration of the governing equations. For brevity, we only explicit the terms in the xmomentum equation (6.a), and the terms in the other two momentum equations have similar meanings along the y and z directions, respectively. In x and y momentum equations (Eqs. (6.a)(6.b)), we neglect the stress at the free surface (e.g., due to the action of the wind), and the depthaveraged turbulent/viscous stresses. Similarly to other studies [8,11,13,17,37,38], in Eq. (6.c), we neglect the vertical advection and dissipation terms compared with the dynamic pressure term, that is, we assume that the dynamic pressure is essentially due to the local time derivative of the vertical velocity of the water column.
The linear vertical variation of q_{b} has been adopted by other authors [8,11,13,17,18,37,38] and it is consistent with the hypothesis of zero value of the atmospheric pressure. Following other authors, 8,11,13,17,18,37,38] we also hypothesize that the vertical component of the velocity changes linearly, from wb to w_{s}, and its depthintegrated W value is
$W=\left({w}_{s}+{w}_{b}\right)/2$
(7),
Both assumptions of the vertical distribution of q_{b} and w require further investigation, but beyond the purpose of this paper. The unknown variables of the governing Partial Differential Equations (PDEs) (5)(6) are h, Uh, Vh, Wh and q_{b}.
General Features of the Numerical Algorithm
The fractionaltimestep procedure
A fractionaltimestep scheme, in the general framework of the pressure projection methods [17], is applied to solve the governing equations (5)(6). The hydrostatic problem, in its turn, is solved applying a predictorcorrector scheme. If we neglect in Eqn. (6.a)(6.b) the terms with the dynamic pressure, the hydrostatic NLSWEs become:
$\frac{\partial h}{\partial t}+\frac{\partial \left(Uh\right)}{\partial x}+\frac{\partial \left(Vh\right)}{\partial y}=0$
(8),
$\frac{\partial \left(Uh\right)}{\partial t}+\frac{\partial}{\partial x}\left({U}^{2}h\right)+\frac{\partial}{\partial y}\left(UVh\right)+gh\frac{\partial H}{\partial x}+g\frac{{n}^{2}\left(Uh\right)\sqrt{{\left(Uh\right)}^{2}+{\left(Vh\right)}^{2}}}{{h}^{7/3}}=0$
(9.a),
$\frac{\partial \left(Vh\right)}{\partial t}+\frac{\partial}{\partial x}\left(UVh\right)+\frac{\partial}{\partial y}\left({V}^{2}h\right)+gh\frac{\partial H}{\partial y}+g\frac{{n}^{2}\left(Vh\right)\sqrt{{\left(Uh\right)}^{2}+{\left(Vh\right)}^{2}}}{{h}^{7/3}}=0$
(9.b),
and the vertical momentum equation becomes trivial. The unknown variables of system (8)(9) are h, Uh and Vh. In Appendix A we give more details of the predictorcorrector procedure applied for the solution of the hydrostatic problem (8)(9). The functional characteristics of the prediction and correction steps of problem (8)(9) are the ones of a convective and a diffusive process, respectively [31,32] (see also Appendix A). According to this, the prediction and correction systems are named "convective prediction" (CP_{0}) and "diffusive correction" (DC) system, respectively. The original hyperbolic shallow waters problem is split into a convective plus a diffusive one.
We solve the hydrostatic CP_{0} problem using the A Marching in Space and Time (MAST) procedure [30,31,33], where the cells are ordered at the beginning of each time iteration, according to rules indicated as follows, and consecutively solved. Originally, MAST has been used for scalar potential flow fields (where the flow vector and the spatial gradient of the potential have opposite signs). To adopt this procedure in cases where a scalar potential of the flow field does not exist, as the present one (where the velocity vector is an unknown of the problem), a further convective correction step has been added [31,32]. The original hydrostatic CP_{0} step is split into a "convective prediction" (CP) and a "convective correction" (CC) step [31,32]. We use an auxiliary scalar function, the "approximated potential" for cells ordering (details in [31,32]).
In the hydrostatic DC step, a large linear system, whose dimension is equal to the number of the computational cells, is solved for the H unknowns, and the horizontal components of the flow rate are subsequently updated. In the dynamic problem of the general fractionaltimestep scheme, we merge the terms of the momentum equations proportional to the nonhydrostatic pressure, as well as the boundary conditions at the free and bottom surfaces (Eqs. (3)), in the continuity equation (Eq. (1)), integrated over a computational cell. The resulting system is solved for the q_{b} unknowns. Hence, the components of the flow rate are corrected, while the water level (depth) is not updated.
Others numerical schemes have been recently proposed for the solution of the nonhydrostatic NLSWEs, where a fractionaltimestep procedure is applied [13,18,37,39]. In the hydrostatic step of these models, the flow field does not satisfy the depthintegrated continuity equation. This is why the momentum and the continuity equations are separately solved. Unlike these literature solvers, in both the hydrostatic and dynamic steps of the present scheme, the computed flow field satisfies the depthintegrated continuity equation discretized in each cell. As shown in sections 4 and 5, this is due 1) to the simultaneous solution of the continuity and momentum equations in the CP and CC steps, and 2) to the embedding of the momentum equations in the continuity equation in the hydrostatic DC and dynamic steps. As a result, the model preserves in all its steps the local and global mass balance, as shown in the following sections.
Computational mesh properties
To facilitate the reader, we briefly remind some basic concepts and nomenclatures of the unstructured triangular Generalized Delaunay (GD) meshes [29] adopted in the present study to discretize the domain.
$\Omega $
is a 2D bounded domain and T_{h} an unstructured triangulation of
$\Omega $
. N_{T}and N represent the total number of triangles and nodes (or vertices) of T_{h}, respectively, and T_{e}, e = 1, …, N_{T}, is a generic triangle of T_{h}, with area
$\left{T}_{e}\right.\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{r}_{i}{}_{,ip}$
is the edge vector connecting nodes i and i_{p}, oriented from i to i_{p}. Triangle T_{e} shares side r_{i},_{ip} with T_{ep},
${c}_{i,ip}^{{T}_{q}}$
is the distance with sign between the T_{q} circumcenter
${c}_{{T}_{q}}$
(q = e, ep) and the midpoint of ri,_{ip}, P_{i},_{ip} (Figure 2a) (see also Eqs. (19.b)(19.b) in [29] and Eq. (B.1) in Appendix B), given by the cross product
${c}_{i,ip}^{{T}_{q}}=\left({x}_{{c}_{q}}{x}_{i,ip}\right)\times \left({x}_{i}{x}_{ip}\right)/\left{r}_{i,ip}\right{\delta}_{q}$
(10),
where x_{i}(_{ip}), x_{i},_{ip} and
${x}_{{c}_{q}}$
are the vectors of the coordinates of i (i_{p}), P_{i},_{ip} and
${c}_{{T}_{q}}$
, respectively,
${\delta}_{q}=\text{}1$
if the direction of r_{i},_{i}_{p} is counterclockwise in triangle T_{q},
${\delta}_{q}=\text{}1$
otherwise, and "x" is the symbol of the cross product. Eq. (B.2) in Appendix B represents the condition to be satisfied to guarantee the GD condition for side r_{i},_{i}_{p}.
Figure 2: a) Notations. b) The computational cell (dual finite volume, or Voronoi polygon). In bold black the boundary of the cell.
We construct a dual mesh over T_{h}, where a dual finite volume ei, i =1, ..., N, is associated with node i. It is the Voronoi polygon, defined as in [40], shown in Figure 2b. We call the dual finite volumes, (computational) cells. A nodecentered formulation is adopted, with the unknown variables stored in each node i. As motivated at the beginning of section 4.2 and in Appendix C, we concentrate the storage capacity in the nodes in the measure of 1/3 of the area of all the triangles sharing the node.
We will prove that the GD mesh property ensures that 1) the sign of the flux among contiguous cells only depends on the flow rate vector, in the CP and CC hydrostatic steps (in section 4.1), 2) the Mproperty of the matrix of the systems in the DC hydrostatic step and in the dynamic problem, respectively (in sections 4.2 and 5), and, 3) the sign of the flux among contiguous cells are consistent with the differences of the values of H and q_{b} in the same cells, respectively, in the DC hydrostatic step and in the nonhydrostatic problem (in sections 4.2 and 5). See for example in [41] the definition of a Mmatrix.
Let Supp_{i} be the support of cell i, i.e., the set of the N_{T}_{,i} triangles sharing node i. A_{i} and Nl_{i} denote the area and the number of the sides of the boundary of the cell associated to node i, respectively. A_{i} is computed as
${A}_{i}=\frac{1}{3}{\displaystyle \sum _{{T}_{r}\in Sup{p}_{i}}\left{T}_{r}\right}$
r = 1, ..., N_{T},_{i} (11).
Hereafter we assume that cell i shares its l^{th} side (l = 1, ..., Nl_{i}) with cell ip (Figure 2b). Length with sign d_{l}_{,i} of side l is given by (Figure 2a)
${d}_{l,i}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\displaystyle \sum _{q=1,2}{c}_{i,ip}^{{T}_{{l}_{q}}}}$
with
${T}_{{l}_{q}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\in \text{\hspace{0.05em}}\text{\hspace{0.05em}}Sup{p}_{i},\text{}q\text{\hspace{0.05em}}=\{\begin{array}{l}1\\ 2\end{array}\text{}\Rightarrow \text{}{T}_{{l}_{q}}\text{\hspace{0.05em}}=\{\begin{array}{l}{T}_{e}\\ {T}_{ep}\end{array},\text{}l=1,\mathrm{...},N{l}_{i}$
(12),
with
${c}_{i,ip}^{{T}_{{l}_{q}}}$
defined in Eq. (10) and Eq. (B.1) in Appendix B. According to the definition of
${c}_{i,ip}^{{T}_{{l}_{q}}}$
, side l is orthogonal to r_{i},_{ip}. From Eqs. (10), (12) and (B.1) in Appendix B, d_{l,i} could be either positive, either negative. In a GD mesh, d_{l,i} is always positive (see Eq. (10) and Appendix B). The boundary of the cell i is defined as
${\Gamma}_{i}={\displaystyle {\sum}_{l=1,N{l}_{i}}{d}_{l,i}}$
, always positive in a GD mesh.
The Hydrostatic Problem
Eqs. (8)(9) are integrated in space (see also Eqs. (A.3)(A.5) in Appendix A), and, using the Green’s theorem, we get the integral formulation
$\frac{\partial {h}_{i}}{\partial t}{A}_{i}+{\displaystyle \sum _{l=1,N{l}_{i}}{F}_{l,i}}=0$
i = 1, …, N (13),
$\frac{\partial {q}_{s,i}}{\partial t}{A}_{i}+{\displaystyle \sum _{l=1,N{l}_{i}}\text{\hspace{0.17em}}{M}_{l,i}^{s}}+{R}_{i}^{s}=0$
s = x, y (14.a),
$s\text{\hspace{0.05em}}=\{\begin{array}{c}x\\ y\end{array}\text{}\Rightarrow \text{}{q}_{s}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}\{\begin{array}{c}Uh\\ Vh\end{array},\text{}{M}_{l,i}^{s}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}\{\begin{array}{c}{M}_{l,i}^{x}\\ {M}_{l,i}^{y}\end{array},\text{}{R}_{i}^{s}\text{\hspace{0.05em}}=\{\begin{array}{c}{R}_{i}^{x}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}{A}_{i}g\left({h}_{i}\text{\hspace{0.05em}}\frac{\partial {H}_{i}^{k}}{\partial x}+{\left(Uh\right)}_{i}\text{\hspace{0.05em}}\text{\hspace{0.05em}}st\right)\\ {R}_{i}^{y}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}{A}_{i}g\left({h}_{i}\text{\hspace{0.05em}}\frac{\partial {H}_{i}^{k}}{\partial y}+{\left(Vh\right)}_{i}st\right)\end{array}\text{}st=\frac{{n}_{i}^{2}\text{\hspace{0.05em}}\sqrt{{\left(Uh\right)}_{i}^{2}\text{\hspace{0.05em}}+{\left(Vh\right)}_{i}^{2}}}{{h}_{i}^{7/3}}$
(14.b),
where sub index i marks the variables in cell i, F_{l},_{i} and
${M}_{l,i}^{x\left(y\right)}$
, further specified, are the flux and the x(y) component of the momentum flux across side l (l = 1, ..., Nl_{i}) of cell i and
${R}_{i}^{x\left(y\right)}$
is the source term. The Manning coefficient in cell i, n_{i} is computed as the mean value of the friction coefficients in the triangles sharing node i. The spatial gradients of the water level are kept constant in the prediction step of the hydrostatic problem (see Eq (14,b)), as motivated in [31,32].
In accordance to the formulation in Eqs. (A.3)(A.5) in Appendix A, the linearized differential formulation of the DC problem is
$\frac{\partial h}{\partial t}+\frac{\partial Uh}{\partial x}+\frac{\partial Vh}{\partial y}=\frac{\partial \left(\overline{Uh}\right)}{\partial x}+\frac{\partial \left(\overline{Vh}\right)}{\partial y}$
(15),
$\frac{\partial {q}_{s}}{\partial t}+g\overline{h}\frac{\partial H}{\partial s}+g\text{\hspace{0.17em}}{n}^{2}\left({q}_{s}\right)\overline{\left(\frac{\sqrt{{\left(Uh\right)}^{2}+{\left(Vh\right)}^{2}}}{{h}^{7/3}}\right)}=g\overline{h}\frac{\partial {H}^{k}}{\partial s}+g\text{\hspace{0.17em}}{n}^{2}\left({\overline{q}}_{s}\right)\overline{\left(\frac{\sqrt{{\left(Uh\right)}^{2}+{\left(Vh\right)}^{2}}}{{h}^{7/3}}\right)}$
(16),
and the symbol (
$\overline{*}$
) marks the meanintime value of the variable (*). The initial states of the CC and DC steps are the solutions of the CP and CC steps, respectively. In Eqs. (16) we assume negligible the difference between the sum of the convective terms and the meanintime values of the same quantities obtained after the (CP + CC) steps [31,32]. Symbols (*)^{k}, (*)^{k}^{+1/3}, (*)k^{+2/3}, (
${(\tilde{*})}^{k+1}$
mark the values of the variable (*) at the beginning of the time iteration, at the end of the convective prediction, convective correction and diffusive correction steps, respectively.
In the prediction steps, we assume, inside each computational cell, piecewise constant values of h, Uh, Vh. The correction system (15)(16) is solved applying a technique similar to the linear (P1) conforming FE Galerkin scheme (see section 4.2). The procedure for the solution of the hydrostatic problem has been proposed in [32], where we adopted a cellcentered formulation with storage capacity concentrated in the mesh triangles. The nodecentered formulation with respect to the Voronoi cells, presented in this paper, leads to a different spatial discretization of the governing equations. In the following sections 4.1 and 4.2, we only report the steps of the numerical formulation of the hydrostatic problem involved in the new spatial discretization. In test 3 (section 8.3) we compare the efficiency and accuracy of the two formulations using an analytical reference solution.
The prediction problem
All the computational cells, at the beginning of each time iteration, are ordered according to the value of the approximated potential
${\varphi}_{i}^{k}$
, known inside each cell, and computed according to the procedure proposed in [31,32]. This procedure, at the beginning of the time step, minimizes the differences [31,32]. From Eq. (12), we set the flux FLl,i through side l as
$F{L}_{l.i}=\left({\widehat{n}}_{i,ip}\cdot {q}_{i}\right){d}_{l,i}$
(17),
where
${\widehat{n}}_{i,ip}$
is the unit vector parallel to side r_{i,ip} (i.e., orthogonal to side l) positive outwards from i, i.e., oriented from i to i_{p}, qi is the specific horizontal flow rate vector in cell i, qi = ((Uh)i, (Vh)i)T. and symbol
$"\cdot "$
marks the dot product. Eq. (17) becomes
$F{L}_{l,i}=\left({\left(Uh\right)}_{i}\left({x}_{ip}{x}_{i}\right)+{\left(Vh\right)}_{i}\left({y}_{ip}{y}_{i}\right)\right){d}_{l,i}/\left{r}_{i,ip}\right$
(18).
If the GD property holds for side r_{i,ip}, d_{l,i} is positive (see Eq. (10) and Appendix B), and the sign of FLl,i only depends on the dot product
$\left({\widehat{n}}_{i,ip}\cdot {q}_{i}\right)$
, and a leaving (incoming) flux from (in) cell i is positive (negative). This condition is not fulfilled if the GD property is not satisfied, since d_{l,i} is negative. The fluxes and momentum fluxes between cells i and i_{p} in Eqs. (13)(14) are defined as
${F}_{l,i}=\text{\hspace{0.17em}}F{L}_{l,i}$
if (
$F{L}_{l,i}>0$
and
$F{L}_{l,i}>F{L}_{m,ip}$
)
${F}_{l,i}=F{L}_{m,ip}$
otherwise (19),
${M}_{l,i}^{s}\text{}=\text{}{F}_{l,i}{q}_{s,i}/{h}_{i}$
if
${F}_{l,i}\text{}\text{\hspace{0.17em}}=\text{}\text{\hspace{0.17em}}F{L}_{l,i}$
,
${M}_{l,i}^{s}\text{}=\text{}\text{\hspace{0.17em}}{F}_{l,i}{q}_{s,ip}/{h}_{ip}$
otherwise (20),
where m is the side of cell ip shared with cell i. Eqs. (19)(20) guarantee the continuity of the flux and the momentum flux at each internal cell side, so that
${F}_{l,i}={F}_{m,ip}$
and
${M}_{l,i}^{x\left(y\right)}=\text{\hspace{0.17em}}{M}_{m,ip}^{x\left(y\right)}$
. For an external boundary side, the condition Fl,i = Fll,i holds for an outward oriented (positive) flux.
Thanks to the consecutive solution of the cells performed in the MAST procedure, the system of PDEs (13)(14) for each cell i are written as a system of Ordinary Differential Equations (ODEs) (more details in [3133]). After the cells ordering, in the CP step we proceed from the cell with the highest to the cell with the lowest value of
${\varphi}^{k}$
, and solve a ODEs system for each cell in the interval
$\left[0,\text{}\Delta t\right]$
(from time level tk to time level tk+1/3and
$\Delta t$
is the size of the time step). After that, in the CC step, we proceed from the cell with the lowest to the cell with the highest value of
${\varphi}^{k}$
, solving an ODEs system in each cell in the interval
$\left[0,\text{}\Delta t\right]$
(from time level tk+1/3 to time level tk+2/3).
A 5th order RungeKutta method with adaptive stepsize control [42] is applied for the solution of the ODEs system. The ODEs system for the CP step is
$\frac{d{h}_{i}}{dt}{A}_{i}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{\displaystyle \sum _{l=1,N{l}_{i}}{\delta}_{l,i}{F}_{l,i}^{p,out}\left(t\right)dt}}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{F}}_{l,i}^{p,in}}$
(21),
$\frac{d{q}_{s,i}\left(t\right)}{dt}{A}_{i}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{R}_{i}^{p,s}\left(t\right)dt}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{\displaystyle \sum _{l=1,N{l}_{i}}{\delta}_{l,i}{M}_{l,i}^{p,s,out}\left(t\right)dt}}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{M}}_{l,i}^{p,s,in}\left(t\right)}$
s = x, y (22),
with
${R}_{i}^{x\left(y\right)}$
defined in Eqs. (14,b),
${\delta}_{li}=1$
if flux across side l is oriented outward i , 0 otherwise and
${F}_{l,i}^{p,out}=\mathrm{max}\left(0,{F}_{l,i}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{M}_{l,i}^{p,s,out}={F}_{l,i}^{p,out}{q}_{s,i}/{h}_{i},\text{\hspace{1em}}s=x,\text{}y$
, if
${\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}$
(23.a),
${F}_{l,i}^{p,in}=\mathrm{min}\left(0,{F}_{l,i}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{M}_{l,i}^{p,x,in}={F}_{l,i}^{p,in}{q}_{s,ip}/{h}_{ip}$
, if
${\varphi}_{i}^{k}<{\varphi}_{ip}^{k}$
(23.b).
The right hand side (r.h.s.) of Eqs. (21)(22) are the meanintime values of the incoming flux and x/y momentum fluxes, respectively, known after the solution of the linked i_{p} cells with higher value of
${\varphi}^{k}$
, solved before i [31,32]. The gradients of H in the Rx(y) term are computed as proposed in 4.3.
Once the ODEs (21)(22) are solved for cell i, the meanintime value of the total flux
${\overline{F}}_{i}^{out}$
leaving from i, is given by the local mass balance [32],
${\overline{F}}_{i}^{out}={\overline{F}}_{i}^{in}\left({h}_{i}^{k+1/3}{h}_{i}^{k}\right)/\Delta t{A}_{i}$
(24),
where
${\overline{F}}_{i}^{in}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{F}}_{l,i}^{p,in}}$
(i.e., the r.h.s. of Eq. (21)). Following [32], we compute the meanintime value of the flux
${\overline{F}}_{l,i}^{out}$
leaving from side l of i to the neighboring cell ip with
${\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}$
, by partitioning
${\overline{F}}_{i}^{out}$
on the base of the ratio between the flux
${F}_{l,i}$
and the sum of all the leaving fluxes to the neighboring cells at the end of the time step,
${\overline{F}}_{l,i}^{out}={\overline{F}}_{i}^{out}{\left({F}_{l,i}^{p,\text{}out}\right)}^{k+1/3}{\scriptscriptstyle \raisebox{1ex}{$$}\!\left/ \!\raisebox{1ex}{$$}\right.}{\sum}_{l=1,N{l}_{i}}{\delta}_{l,i}{\left({F}_{l,i}^{p,\text{}out}\right)}^{k+1/3}$
(25.a),
and we estimate in a similar way the meanintime values of the leaving momentum flux
${\overline{M}}_{l,i}^{x\left(y\right),out}$
. Finally, we set:
${\overline{F}}_{m,ip}^{in}={\overline{F}}_{l,i}^{out}$
${\overline{M}}_{m,ip}^{s,in}={\overline{M}}_{l,i}^{s,out}\text{}s=x,\text{}y$
(25.b),
for all the neighboring ip cells with
${\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}$
and we proceed to solve system (21)(22) for the i_{p} cells among the unsolved ones. A similar algorithm is applied to solve the convective correction step, whose ODEs system for the generic cell i is
$\frac{d{h}_{i}}{dt}{A}_{i}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{\displaystyle \sum _{l=1,N{l}_{i}}{\delta}_{l,i}{F}_{l,i}^{c,out}\left(t\right)dt}}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{F}}_{l,i}^{c,in}}$
(26),
$\frac{d{q}_{s,i}\left(t\right)}{dt}{A}_{i}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{\displaystyle \sum _{l=1,N{l}_{i}}{\delta}_{l,i}{M}_{l,i}^{c,s,out}\left(t\right)dt}}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{M}}_{l,i}^{c,s,in}}$
s = x, y (27),
with
${F}_{l,i}^{c,out}=\mathrm{max}\left(0,{F}_{l,i}\right)$
${M}_{l,i}^{c,s,out}={F}_{l,i}^{c,out}{q}_{s,i}/{h}_{i}$
s = x, y if
${\varphi}_{i}^{k}<{\varphi}_{ip}^{k}$
(28,a),
and
${F}_{l,i}^{c,in}=\mathrm{min}\left(0,{F}_{l,i}\right)$
${M}_{l,i}^{c,s,in}={F}_{l,i}^{c,in}{q}_{s,ip}/{h}_{ip}$
if
${\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}$
(28,b).
The source terms is allocated in the convective prediction step [32]. The initial state in each cell of the convective correction step is the one computed in the same cell at the end of the convective prediction step. After the ODEs system (26)(27) is solved in cell i, the meanintime value of the flux/momentum flux leaving from side l of i to the neighboring cell ip with
${\varphi}_{i}^{k}<{\varphi}_{ip}^{k}$
, from tk+1/3 to tk+2/3, is computed similarly to the previous CP step. Further comments concerning the MAST procedure can be found in [31,32].
The meanintime values
$\overline{h}$
and
$\overline{\sqrt{{\left(Uh\right)}^{2}+{\left(Vh\right)}^{2}\raisebox{1ex}{$$}\!\left/ \!\raisebox{1ex}{$$}\right.{h}^{7/3}}}$
, required for the DC step (see Eqs. (16)), are obtained by numerical integration during the solution of each ODEs system during the convective prediction and correction steps [31,32].
The diffusive correction problem
After the prediction components of the water depths (levels) are computed with piecewise constant approximation (P0) inside each cell, a continuous piecewise linear (P1) shape is assigned to the water depths (levels) computed at the end of the (CP + CC) steps, obtained from the nodal values. If we move from the P0 to the P1 approximation, at a given time level, we do not alter the estimation of the global mass, if the polygons of the P0 approximation (Figure 2b) are the same as the ones of the P1 approximation. If we use, for the solution of the DC problem the conforming P1 Galerkin FE scheme, where the polygons are delimited by the centers of mass (see Appendix C), we introduce a small mass balance error. We avoid this error, using the same polygons for the P0 and P1 approximations. Other details in Appendix C.
We use a fully implicit time discretization of the DC problem (15)(16) [32]. Under the hypothesis of fixed bed condition, we get the following differential equation for the generic cell i
$\frac{{\tilde{H}}_{i}^{k+1}{H}_{i}^{k+2/3}}{\Delta t}+\frac{\partial {\left(\tilde{U}\tilde{h}\right)}^{k+1}}{\partial x}+\frac{\partial {\left(\tilde{V}\tilde{h}\right)}^{k+1}}{\partial y}=\frac{\partial \overline{\left(Uh\right)}}{\partial x}+\frac{\partial \overline{\left(Vh\right)}}{\partial y}$
(29),
$\frac{{\tilde{q}}_{s,i}^{k+1}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\tilde{q}}_{s,i}^{k+2/3}}{\Delta t}+g\text{\hspace{0.17em}}{\overline{h}}_{i}\frac{\partial {\tilde{H}}^{k+1}}{\partial s}\text{\hspace{0.05em}}+\text{\hspace{0.17em}}\text{\hspace{0.05em}}g\text{\hspace{0.17em}}{\tilde{q}}_{s,i}^{k+1}so{r}_{i}=g\text{\hspace{0.17em}}{\overline{h}}_{i}\frac{\partial {H}^{k}}{\partial s}\text{\hspace{0.05em}}+g\text{\hspace{0.17em}}{\overline{q}}_{s,i}\text{\hspace{0.17em}}so{r}_{i}$
s = x, y (30),
where
$so{r}_{i}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}{n}_{i}^{2}\overline{\left(\sqrt{{\left(Uh\right)}_{i}^{2}\text{\hspace{0.05em}}+\text{\hspace{0.05em}}{\left(Vh\right)}_{i}^{2}}/{h}_{i}^{7/3}\right)}$
. From Eq. (30), the flow rate components at the end of the hydrostatic step are
${\tilde{q}}_{s,i}^{k+1}={\tilde{D}}_{i}\frac{\partial {\tilde{H}}^{k+1}}{\partial s}+{k}_{i}^{s}+{q}_{s,i}^{k+2/3}$
s = x, y (31),
with coefficients
${\tilde{D}}_{i}$
and
${k}_{i}^{x\left(y\right)}$
defined as
${\tilde{D}}_{i}=\{\begin{array}{l}g{\overline{h}}_{i}\Delta t/\left[1+\Delta t\text{\hspace{0.17em}}g\text{\hspace{0.17em}}\text{\hspace{0.17em}}so{r}_{i}\right]\text{if}{\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}\\ g{\overline{h}}_{ip}\Delta t/\left[1+\Delta t\text{\hspace{0.17em}}g\text{\hspace{0.17em}}\text{\hspace{0.17em}}so{r}_{ip}\right]\text{if}{\varphi}_{ip}^{k}{\varphi}_{i}^{k}\end{array}\text{\hspace{0.17em}}\text{}{k}_{i}^{s}={\tilde{D}}_{i}\frac{\partial {H}^{k}}{\partial s}$
s = x, y (32).
We merge Eqs. (31) in Eq. (29) and, after space integration and application of the Green's theorem, we get the following balance law for cell i
$\underset{{A}_{i}}{\int}\frac{\partial {H}_{i}}{\partial t}d{A}_{i}}{\displaystyle \sum _{l=1,N{l}_{i}}\left({\tilde{D}}_{i}\frac{\partial {\tilde{H}}^{k+1}}{\partial {n}_{i,ip}}{d}_{l,i}\right)}={\displaystyle \sum _{l=1,N{l}_{i}}\left({\tilde{D}}_{i}\frac{\partial {\tilde{H}}^{k}}{\partial {n}_{i,ip}}{d}_{l,i}\right)}+{\displaystyle \sum _{l=1,N{l}_{i}}\left({\overline{q}}_{i}{q}_{i}^{k+2/3}\right)\cdot {\widehat{n}}_{i,ip}{d}_{l,i}$
(33),
where
$\partial H/\partial {n}_{i,ip}$
is the water level gradient along the
${\widehat{n}}_{i,ip}$
direction and the other symbols are specified as above. Concerning the second term on the r.h.s. of Eq. (33), we have
$\sum _{l=1,N{l}_{i}}\left(\text{\hspace{0.17em}}{\overline{q}}_{i}\cdot {\widehat{n}}_{i,ip}{d}_{l,i}\right)}={\overline{F}}_{i}^{in}{\overline{F}}_{i}^{out}=\frac{{H}_{i}^{k+2/3}{H}_{i}^{k}}{\Delta t}{A}_{i$
(34),
while
${q}_{i}^{k+2/3}$
is computed after the (CP + CC) steps. After discretization in time, Eq. (33) becomes
$\frac{{\tilde{H}}_{i}^{k+1}{H}_{i}^{k+2/3}}{\Delta t}{A}_{i}+{\displaystyle \sum _{l=1,N{l}_{i}}{\tilde{F}}_{l,i}}={\displaystyle \sum _{l=1,N{l}_{i}}{\tilde{b}}_{l,i}}$
i = 1, ..., N (35),
where
${\tilde{F}}_{l,i}$
is the corrective flux across side l of cell i
${\tilde{F}}_{l,i}={\tilde{D}}_{i}\frac{\partial {\tilde{H}}^{k+1}}{\partial {n}_{i,ip}}{d}_{l,i}$
(36.a),
and the term
${\tilde{b}}_{l,i}$
on the r.h.s. of Eq. (35) is
${\tilde{b}}_{l,i}={\tilde{D}}_{i}\frac{\partial {\tilde{H}}^{k}}{\partial {n}_{i,ip}}{d}_{l,i}+\left({\overline{q}}_{i}{q}_{i}^{k+2/3}\right)\cdot {\widehat{n}}_{i,ip}{}_{i,ip}{d}_{l,i}$
(36.b).
Due to the P1 spatial approximation,
$\partial {\tilde{H}}^{k+1}/\partial {n}_{i,ip}$
only depends on the values at nodes i and i_{p} (since
${\widehat{n}}_{i,ip}\text{\hspace{0.05em}}$
is parallel to r_{i,ip} and orthogonal to side l, see sections 3.2 and 4.1), and
${\tilde{F}}_{l,i}$
is written
${\tilde{F}}_{l,i}={\tilde{D}}_{i}\left({\tilde{H}}_{ip}^{k+1}{\tilde{H}}_{i}^{k+1}\right)\raisebox{1ex}{$$}\!\left/ \!\raisebox{1ex}{$$}\right.\text{\hspace{0.17em}}\left{r}_{i,ip}\right{d}_{l,i}$
(37).
Eqs. (35) constitute a linear system of order N in the Hi (i = 1, …, N) unknowns, whose matrix is sparse and symmetric. Generally the nonzero entries are about 5% or less of the total number of matrix coefficients. The diagonal and offdiagonal matrix coefficients are, respectively
${s}_{i,i}=\frac{{A}_{i}}{\Delta t}+{\displaystyle \sum _{l=1,N{l}_{i}}\frac{{\tilde{D}}_{i}}{\left{r}_{i,ip}\right}{d}_{l,i}}\text{\hspace{0.17em}}$
${s}_{i,ip}=\frac{{\tilde{D}}_{i}}{\left{r}_{i,ip}\right}\text{\hspace{0.17em}}\text{\hspace{0.05em}}{d}_{l,i}\text{\hspace{0.17em}}{\delta}_{i,ip}$
i, i_{p} = 1, ..., N (38),
where
${\delta}_{i}{}_{,ip}=\text{}1$
if nodes i and i_{p} are linked, otherwise it is equal to zero.
According to Eqs. (37)(38), the offdiagonal coefficients are nonnegative if the GD mesh property is not satisfied for side r_{i,ip} and the Mmatrix property is lost. In this case, d_{l,i} is negative (see Eq. (10) and Appendix B) and the sign of
${\tilde{F}}_{l,i}$
is not consistent with the difference
$\left({\tilde{H}}_{ip}^{k+1}{\tilde{H}}_{i}^{k+1}\right)$
in Eq. (37). On the opposite, if side r_{i,ip} satisfies the GD property, d_{l,i} is positive, and the sign of the flux in Eq. (37) only depends on the difference of the water levels in cells i and i_{p}. If the GD property is satisfied for all the sides of the mesh, the matrix of system (35) is a Mmatrix, positivedefinite and, according to [41], the system is wellconditioned.
The numerical scheme adopted for the solution of the DC problem has a structure similar to the P1 conforming FE Galerkin scheme, and we present and discuss analogies and differences between the two numerical formulations in Appendix C.
We solve system (35) by a preconditioned conjugate gradient method (with incomplete Cholesky preconditioning). Due to the symmetry of the matrix of system (35), we only save the nonzero entries of the lower triangular matrix, using a Symmetric Coordinate Storage algorithm [43 ].
After the solution of system (35) for the H unknowns, we compute the gradients of the water level
${\nabla}_{x\left(y\right)}{\tilde{H}}_{i}^{k+1}$
at the end of the DC step as shown in the next section. As specified in section 4.1, the values of the gradients of the water level are assumed constant during the solution of the two convective steps of the hydrostatic problem in the next time iteration [31,32]. Finally, the components of the flow rate
${\tilde{q}}_{s,i}^{k+1}$
(s = x, y) are updated according to Eqs. (31)(32).
According to Eq. (37) and the definition of the coefficient
${\tilde{D}}_{i}$
in (Eq. (32)), the flux
${\tilde{F}}_{l,i}$
leaving from cell i (i_{p}) to i_{p} (i) is equal to the flux entering cell i_{p} (i) from cell i (i_{p}). If
${\tilde{F}}_{l,i}$
is leaving from i to i_{p} or vice versa depends on the difference
$\left({\tilde{H}}_{i}^{k+1}{\tilde{H}}_{ip}^{k+1}\right)$
(see Eq. (37).
Computation of the spatial gradients of the water level
We compute the gradient of H in cell i as follows. Call
${\widehat{n}}_{q}$
the unit vector parallel to the flow vector in cell i. Let
$\nabla {H}_{i,ip}$
be the vector computed as (i_{p} is one of the cells linked with i),
$\nabla {H}_{i,ip}=\left({H}_{i}{H}_{ip}\right)/{r}_{i,ip}$
(39).
We assume
$\nabla {H}_{i}=\frac{{\displaystyle {\sum}_{ip=1,N{l}_{i}}\mathrm{max}\left(0,\left({\widehat{n}}_{i,ip}\cdot {\widehat{n}}_{q}\right)\nabla {H}_{i,ip}\right)}}{{\displaystyle {\sum}_{ip=1,N{l}_{i}}\mathrm{max}\left(0,\left({\widehat{n}}_{i,ip}\cdot {\widehat{n}}_{q}\right)\right)}}$
(40),
where
${\widehat{n}}_{i,ip}\cdot {\widehat{n}}_{q}$
is the dot product between the direction of side r_{i,ip} and the unitary flow rate vector. If cell i lies on an impervious boundary, the flow rate vector is parallel to the boundary and
$\nabla {H}_{i}$
, computed as in Eq. (40), is projected along the boundary direction.
The expression in Eq. (40) can be regarded as an upwind formulation, where
$\nabla {H}_{i}$
depends on the differences between Hi and the values of H in the downstream (along the flow direction) i_{p} cells. This is consistent with the procedure adopted in the (CP + CC) steps, where the unknown fluxes/momentum fluxes are the ones leaving from i to the downstream i_{p} cells.
The Hydrostatic Problem
Eqs. (8)(9) are integrated in space (see also Eqs. (A.3)(A.5) in Appendix A), and, using the Green’s theorem, we get the integral formulation
$\frac{\partial {h}_{i}}{\partial t}{A}_{i}+{\displaystyle \sum _{l=1,N{l}_{i}}{F}_{l,i}}=0$
i = 1, …, N (13),
$\frac{\partial {q}_{s,i}}{\partial t}{A}_{i}+{\displaystyle \sum _{l=1,N{l}_{i}}\text{\hspace{0.17em}}{M}_{l,i}^{s}}+{R}_{i}^{s}=0$
s = x, y (14.a),
$s\text{\hspace{0.05em}}=\{\begin{array}{c}x\\ y\end{array}\text{}\Rightarrow \text{}{q}_{s}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}\{\begin{array}{c}Uh\\ Vh\end{array},\text{}{M}_{l,i}^{s}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}\{\begin{array}{c}{M}_{l,i}^{x}\\ {M}_{l,i}^{y}\end{array},\text{}{R}_{i}^{s}\text{\hspace{0.05em}}=\{\begin{array}{c}{R}_{i}^{x}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}{A}_{i}g\left({h}_{i}\text{\hspace{0.05em}}\frac{\partial {H}_{i}^{k}}{\partial x}+{\left(Uh\right)}_{i}\text{\hspace{0.05em}}\text{\hspace{0.05em}}st\right)\\ {R}_{i}^{y}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}{A}_{i}g\left({h}_{i}\text{\hspace{0.05em}}\frac{\partial {H}_{i}^{k}}{\partial y}+{\left(Vh\right)}_{i}st\right)\end{array}\text{}st=\frac{{n}_{i}^{2}\text{\hspace{0.05em}}\sqrt{{\left(Uh\right)}_{i}^{2}\text{\hspace{0.05em}}+{\left(Vh\right)}_{i}^{2}}}{{h}_{i}^{7/3}}$
(14.b),
where sub index i marks the variables in cell i, Fl,i and
${M}_{l,i}^{x\left(y\right)}$
, further specified, are the flux and the x(y) component of the momentum flux across side l (l = 1, ..., Nli) of cell i and
${M}_{l,i}^{x\left(y\right)}$
is the source term. The Manning coefficient in cell i, ni is computed as the mean value of the friction coefficients in the triangles sharing node i. The spatial gradients of the water level are kept constant in the prediction step of the hydrostatic problem (see Eq (14,b)), as motivated in [31,32].
In accordance to the formulation in Eqs. (A.3)(A.5) in Appendix A, the linearized differential formulation of the DC problem is
$\frac{\partial h}{\partial t}+\frac{\partial Uh}{\partial x}+\frac{\partial Vh}{\partial y}=\frac{\partial \left(\overline{Uh}\right)}{\partial x}+\frac{\partial \left(\overline{Vh}\right)}{\partial y}$
(15),
$\frac{\partial {q}_{s}}{\partial t}+g\overline{h}\frac{\partial H}{\partial s}+g\text{\hspace{0.17em}}{n}^{2}\left({q}_{s}\right)\overline{\left(\frac{\sqrt{{\left(Uh\right)}^{2}+{\left(Vh\right)}^{2}}}{{h}^{7/3}}\right)}=g\overline{h}\frac{\partial {H}^{k}}{\partial s}+g\text{\hspace{0.17em}}{n}^{2}\left({\overline{q}}_{s}\right)\overline{\left(\frac{\sqrt{{\left(Uh\right)}^{2}+{\left(Vh\right)}^{2}}}{{h}^{7/3}}\right)}$
(16),
and the symbol (
$\overline{*}$
) marks the meanintime value of the variable (*). The initial states of the CC and DC steps are the solutions of the CP and CC steps, respectively. In Eqs. (16) we assume negligible the difference between the sum of the convective terms and the meanintime values of the same quantities obtained after the (CP + CC) steps [31,32]. Symbols (*)k, (*)k+1/3, (*)k+2/3, (
${(\tilde{*})}^{k+1}$
mark the values of the variable (*) at the beginning of the time iteration, at the end of the convective prediction, convective correction and diffusive correction steps, respectively.
In the prediction steps, we assume, inside each computational cell, piecewise constant values of h, Uh, Vh. The correction system (15)(16) is solved applying a technique similar to the linear (P1) conforming FE Galerkin scheme (see section 4.2). The procedure for the solution of the hydrostatic problem has been proposed in [32], where we adopted a cellcentered formulation with storage capacity concentrated in the mesh triangles. The nodecentered formulation with respect to the Voronoi cells, presented in this paper, leads to a different spatial discretization of the governing equations. In the following sections 4.1 and 4.2, we only report the steps of the numerical formulation of the hydrostatic problem involved in the new spatial discretization. In test 3 (section 8.3) we compare the efficiency and accuracy of the two formulations using an analytical reference solution.
The prediction problem
All the computational cells, at the beginning of each time iteration, are ordered according to the value of the approximated potential
${\varphi}_{i}^{k}$
, known inside each cell, and computed according to the procedure proposed in [31,32]. This procedure, at the beginning of the time step, minimizes the differences [31,32]. From Eq. (12), we set the flux FL_{l},_{i} through side l as
$F{L}_{l.i}=\left({\widehat{n}}_{i,ip}\cdot {q}_{i}\right){d}_{l,i}$
(17),
where
${\widehat{n}}_{i,ip}$
is the unit vector parallel to side r_{i},_{i}_{p} (i.e., orthogonal to side l) positive outwards from i, i.e., oriented from i to i_{p}, q_{i} is the specific horizontal flow rate vector in cell i, q_{i} = ((Uh)_{i}, (Vh)_{i})^{T}. and symbol
$"\cdot "$
marks the dot product. Eq. (17) becomes
$F{L}_{l,i}=\left({\left(Uh\right)}_{i}\left({x}_{ip}{x}_{i}\right)+{\left(Vh\right)}_{i}\left({y}_{ip}{y}_{i}\right)\right){d}_{l,i}/\left{r}_{i,ip}\right$
(18).
If the GD property holds for side r_{i},i_{p}, d_{l,i} is positive (see Eq. (10) and Appendix B), and the sign of FL_{l},_{i} only depends on the dot product
$\left({\widehat{n}}_{i,ip}\cdot {q}_{i}\right)$
, and a leaving (incoming) flux from (in) cell i is positive (negative). This condition is not fulfilled if the GD property is not satisfied, since d_{l,i} is negative. The fluxes and momentum fluxes between cells i and i_{p} in Eqs. (13)(14) are defined as
${F}_{l,i}=\text{\hspace{0.17em}}F{L}_{l,i}$
if (
$F{L}_{l,i}>0$
and
$F{L}_{l,i}>F{L}_{m,ip}$
)
${F}_{l,i}=F{L}_{m,ip}$
otherwise (19),
${M}_{l,i}^{s}\text{}=\text{}{F}_{l,i}{q}_{s,i}/{h}_{i}$
if
${F}_{l,i}\text{}\text{\hspace{0.17em}}=\text{}\text{\hspace{0.17em}}F{L}_{l,i}$
,
${M}_{l,i}^{s}\text{}=\text{}\text{\hspace{0.17em}}{F}_{l,i}{q}_{s,ip}/{h}_{ip}$
otherwise (20),
where m is the side of cell ip shared with cell i. Eqs. (19)(20) guarantee the continuity of the flux and the momentum flux at each internal cell side, so that
${F}_{l,i}={F}_{m,ip}$
and
${M}_{l,i}^{x\left(y\right)}=\text{\hspace{0.17em}}{M}_{m,ip}^{x\left(y\right)}$
. For an external boundary side, the condition F_{l},_{i} = Fl_{l},_{i} holds for an outward oriented (positive) flux.
Thanks to the consecutive solution of the cells performed in the MAST procedure, the system of PDEs (13)(14) for each cell i are written as a system of Ordinary Differential Equations (ODEs) (more details in [3133]). After the cells ordering, in the CP step we proceed from the cell with the highest to the cell with the lowest value of
${\varphi}^{k}$
, and solve a ODEs system for each cell in the interval
$\left[0,\text{}\Delta t\right]$
(from time level t^{k} to time level t^{k}^{+1/3}and
$\Delta t$
is the size of the time step). After that, in the CC step, we proceed from the cell with the lowest to the cell with the highest value of
${\varphi}^{k}$
, solving an ODEs system in each cell in the interval
$\left[0,\text{}\Delta t\right]$
(from time level t^{k}^{+1/3} to time level t^{k}^{+2/3}).
A 5^{th} order RungeKutta method with adaptive stepsize control [42] is applied for the solution of the ODEs system. The ODEs system for the CP step is
$\frac{d{h}_{i}}{dt}{A}_{i}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{\displaystyle \sum _{l=1,N{l}_{i}}{\delta}_{l,i}{F}_{l,i}^{p,out}\left(t\right)dt}}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{F}}_{l,i}^{p,in}}$
(21),
$\frac{d{q}_{s,i}\left(t\right)}{dt}{A}_{i}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{R}_{i}^{p,s}\left(t\right)dt}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{\displaystyle \sum _{l=1,N{l}_{i}}{\delta}_{l,i}{M}_{l,i}^{p,s,out}\left(t\right)dt}}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{M}}_{l,i}^{p,s,in}\left(t\right)}$
s = x, y (22),
with
${R}_{i}^{x\left(y\right)}$
defined in Eqs. (14,b),
${\delta}_{li}=1$
if flux across side l is oriented outward i , 0 otherwise and
${F}_{l,i}^{p,out}=\mathrm{max}\left(0,{F}_{l,i}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{M}_{l,i}^{p,s,out}={F}_{l,i}^{p,out}{q}_{s,i}/{h}_{i},\text{\hspace{1em}}s=x,\text{}y$
, if
${\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}$
(23.a),
${F}_{l,i}^{p,in}=\mathrm{min}\left(0,{F}_{l,i}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{M}_{l,i}^{p,x,in}={F}_{l,i}^{p,in}{q}_{s,ip}/{h}_{ip}$
, if
${\varphi}_{i}^{k}<{\varphi}_{ip}^{k}$
(23.b).
The right hand side (r.h.s.) of Eqs. (21)(22) are the meanintime values of the incoming flux and x/y momentum fluxes, respectively, known after the solution of the linked i_{p} cells with higher value of
${\varphi}^{k}$
, solved before i [31,32]. The gradients of H in the R^{x}(^{y}) term are computed as proposed in 4.3.
Once the ODEs (21)(22) are solved for cell i, the meanintime value of the total flux
${\overline{F}}_{i}^{out}$
leaving from i, is given by the local mass balance [32],
${\overline{F}}_{i}^{out}={\overline{F}}_{i}^{in}\left({h}_{i}^{k+1/3}{h}_{i}^{k}\right)/\Delta t{A}_{i}$
(24),
where
${\overline{F}}_{i}^{in}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{F}}_{l,i}^{p,in}}$
(i.e., the r.h.s. of Eq. (21)). Following [32], we compute the meanintime value of the flux
${\overline{F}}_{l,i}^{out}$
leaving from side l of i to the neighboring cell ip with
${\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}$
, by partitioning
${\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}$
on the base of the ratio between the flux
${\overline{F}}_{i}^{out}$
and the sum of all the leaving fluxes to the neighboring cells at the end of the time step,
${\overline{F}}_{l,i}^{out}={\overline{F}}_{i}^{out}{\left({F}_{l,i}^{p,\text{}out}\right)}^{k+1/3}{\scriptscriptstyle \raisebox{1ex}{$$}\!\left/ \!\raisebox{1ex}{$$}\right.}{\sum}_{l=1,N{l}_{i}}{\delta}_{l,i}{\left({F}_{l,i}^{p,\text{}out}\right)}^{k+1/3}$
(25.a),
and we estimate in a similar way the meanintime values of the leaving momentum flux
${\overline{M}}_{l,i}^{x\left(y\right),out}$
. Finally, we set:
${\overline{F}}_{m,ip}^{in}={\overline{F}}_{l,i}^{out}$
${\overline{M}}_{m,ip}^{s,in}={\overline{M}}_{l,i}^{s,out}\text{}s=x,\text{}y$
(25.b),
for all the neighboring ip cells with
${\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}$
and we proceed to solve system (21)(22) for the i_{p} cells among the unsolved ones. A similar algorithm is applied to solve the convective correction step, whose ODEs system for the generic cell i is
$\frac{d{h}_{i}}{dt}{A}_{i}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{\displaystyle \sum _{l=1,N{l}_{i}}{\delta}_{l,i}{F}_{l,i}^{c,out}\left(t\right)dt}}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{F}}_{l,i}^{c,in}}$
(26),
$\frac{d{q}_{s,i}\left(t\right)}{dt}{A}_{i}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{\displaystyle \sum _{l=1,N{l}_{i}}{\delta}_{l,i}{M}_{l,i}^{c,s,out}\left(t\right)dt}}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{M}}_{l,i}^{c,s,in}}$
s = x, y (27),
with
${F}_{l,i}^{c,out}=\mathrm{max}\left(0,{F}_{l,i}\right)$
${M}_{l,i}^{c,s,out}={F}_{l,i}^{c,out}{q}_{s,i}/{h}_{i}$
s = x, y if
${\varphi}_{i}^{k}<{\varphi}_{ip}^{k}$
(28,a),
and
${F}_{l,i}^{c,in}=\mathrm{min}\left(0,{F}_{l,i}\right)$
${M}_{l,i}^{c,s,in}={F}_{l,i}^{c,in}{q}_{s,ip}/{h}_{ip}$
if
${\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}$
(28,b).
The source terms is allocated in the convective prediction step [32]. The initial state in each cell of the convective correction step is the one computed in the same cell at the end of the convective prediction step. After the ODEs system (26)(27) is solved in cell i, the meanintime value of the flux/momentum flux leaving from side l of i to the neighboring cell ip with
${\varphi}_{i}^{k}<{\varphi}_{ip}^{k}$
, from t^{k}^{+1/3} to t^{k}^{+2/3}, is computed similarly to the previous CP step. Further comments concerning the MAST procedure can be found in [31,32].
The meanintime values
$\overline{h}$
and
$\overline{\sqrt{{\left(Uh\right)}^{2}+{\left(Vh\right)}^{2}\raisebox{1ex}{$$}\!\left/ \!\raisebox{1ex}{$$}\right.{h}^{7/3}}}$
, required for the DC step (see Eqs. (16)), are obtained by numerical integration during the solution of each ODEs system during the convective prediction and correction steps [31,32].
The diffusive correction problem
After the prediction components of the water depths (levels) are computed with piecewise constant approximation (P0) inside each cell, a continuous piecewise linear (P1) shape is assigned to the water depths (levels) computed at the end of the (CP + CC) steps, obtained from the nodal values. If we move from the P0 to the P1 approximation, at a given time level, we do not alter the estimation of the global mass, if the polygons of the P0 approximation (Figure 2b) are the same as the ones of the P1 approximation. If we use, for the solution of the DC problem the conforming P1 Galerkin FE scheme, where the polygons are delimited by the centers of mass (see Appendix C), we introduce a small mass balance error. We avoid this error, using the same polygons for the P0 and P1 approximations. Other details in Appendix C.
We use a fully implicit time discretization of the DC problem (15)(16) [32]. Under the hypothesis of fixed bed condition, we get the following differential equation for the generic cell i
$\frac{{\tilde{H}}_{i}^{k+1}{H}_{i}^{k+2/3}}{\Delta t}+\frac{\partial {\left(\tilde{U}\tilde{h}\right)}^{k+1}}{\partial x}+\frac{\partial {\left(\tilde{V}\tilde{h}\right)}^{k+1}}{\partial y}=\frac{\partial \overline{\left(Uh\right)}}{\partial x}+\frac{\partial \overline{\left(Vh\right)}}{\partial y}$
(29),
$\frac{{\tilde{q}}_{s,i}^{k+1}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\tilde{q}}_{s,i}^{k+2/3}}{\Delta t}+g\text{\hspace{0.17em}}{\overline{h}}_{i}\frac{\partial {\tilde{H}}^{k+1}}{\partial s}\text{\hspace{0.05em}}+\text{\hspace{0.17em}}\text{\hspace{0.05em}}g\text{\hspace{0.17em}}{\tilde{q}}_{s,i}^{k+1}so{r}_{i}=g\text{\hspace{0.17em}}{\overline{h}}_{i}\frac{\partial {H}^{k}}{\partial s}\text{\hspace{0.05em}}+g\text{\hspace{0.17em}}{\overline{q}}_{s,i}\text{\hspace{0.17em}}so{r}_{i}$
s = x, y (30),
where
$so{r}_{i}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}{n}_{i}^{2}\overline{\left(\sqrt{{\left(Uh\right)}_{i}^{2}\text{\hspace{0.05em}}+\text{\hspace{0.05em}}{\left(Vh\right)}_{i}^{2}}/{h}_{i}^{7/3}\right)}$
. From Eq. (30), the flow rate components at the end of the hydrostatic step are
${\tilde{q}}_{s,i}^{k+1}={\tilde{D}}_{i}\frac{\partial {\tilde{H}}^{k+1}}{\partial s}+{k}_{i}^{s}+{q}_{s,i}^{k+2/3}$
s = x, y (31),
with coefficients
${\tilde{D}}_{i}$
and
${k}_{i}^{x\left(y\right)}$
defined as
${\tilde{D}}_{i}=\{\begin{array}{l}g{\overline{h}}_{i}\Delta t/\left[1+\Delta t\text{\hspace{0.17em}}g\text{\hspace{0.17em}}\text{\hspace{0.17em}}so{r}_{i}\right]\text{if}{\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}\\ g{\overline{h}}_{ip}\Delta t/\left[1+\Delta t\text{\hspace{0.17em}}g\text{\hspace{0.17em}}\text{\hspace{0.17em}}so{r}_{ip}\right]\text{if}{\varphi}_{ip}^{k}{\varphi}_{i}^{k}\end{array}\text{\hspace{0.17em}}\text{}{k}_{i}^{s}={\tilde{D}}_{i}\frac{\partial {H}^{k}}{\partial s}$
s = x, y (32).
We merge Eqs. (31) in Eq. (29) and, after space integration and application of the Green's theorem, we get the following balance law for cell i
$\underset{{A}_{i}}{\int}\frac{\partial {H}_{i}}{\partial t}d{A}_{i}}{\displaystyle \sum _{l=1,N{l}_{i}}\left({\tilde{D}}_{i}\frac{\partial {\tilde{H}}^{k+1}}{\partial {n}_{i,ip}}{d}_{l,i}\right)}={\displaystyle \sum _{l=1,N{l}_{i}}\left({\tilde{D}}_{i}\frac{\partial {\tilde{H}}^{k}}{\partial {n}_{i,ip}}{d}_{l,i}\right)}+{\displaystyle \sum _{l=1,N{l}_{i}}\left({\overline{q}}_{i}{q}_{i}^{k+2/3}\right)\cdot {\widehat{n}}_{i,ip}{d}_{l,i}$
(33),
where
$\partial H/\partial {n}_{i,ip}$
is the water level gradient along the
${\widehat{n}}_{i,ip}$
direction and the other symbols are specified as above. Concerning the second term on the r.h.s. of Eq. (33), we have
$\sum _{l=1,N{l}_{i}}\left(\text{\hspace{0.17em}}{\overline{q}}_{i}\cdot {\widehat{n}}_{i,ip}{d}_{l,i}\right)}={\overline{F}}_{i}^{in}{\overline{F}}_{i}^{out}=\frac{{H}_{i}^{k+2/3}{H}_{i}^{k}}{\Delta t}{A}_{i$
(34),
while
${q}_{i}^{k+2/3}$
is computed after the (CP + CC) steps. After discretization in time, Eq. (33) becomes
$\frac{{\tilde{H}}_{i}^{k+1}{H}_{i}^{k+2/3}}{\Delta t}{A}_{i}+{\displaystyle \sum _{l=1,N{l}_{i}}{\tilde{F}}_{l,i}}={\displaystyle \sum _{l=1,N{l}_{i}}{\tilde{b}}_{l,i}}$
i = 1, ..., N (35),
where
${\tilde{F}}_{l,i}$
is the corrective flux across side l of cell i
${\tilde{F}}_{l,i}={\tilde{D}}_{i}\frac{\partial {\tilde{H}}^{k+1}}{\partial {n}_{i,ip}}{d}_{l,i}$
(36.a),
and the term
${\tilde{b}}_{l,i}$
on the r.h.s. of Eq. (35) is
${\tilde{b}}_{l,i}={\tilde{D}}_{i}\frac{\partial {\tilde{H}}^{k}}{\partial {n}_{i,ip}}{d}_{l,i}+\left({\overline{q}}_{i}{q}_{i}^{k+2/3}\right)\cdot {\widehat{n}}_{i,ip}{}_{i,ip}{d}_{l,i}$
(36.b).
Due to the P1 spatial approximation,
$\partial {\tilde{H}}^{k+1}/\partial {n}_{i,ip}$
only depends on the values at nodes i and i_{p} (since
${\widehat{n}}_{i,ip}\text{\hspace{0.05em}}$
is parallel to r_{i}_{,}_{ip} and orthogonal to side l, see sections 3.2 and 4.1), and
${\tilde{F}}_{l,i}$
is written
${\tilde{F}}_{l,i}={\tilde{D}}_{i}\left({\tilde{H}}_{ip}^{k+1}{\tilde{H}}_{i}^{k+1}\right)\raisebox{1ex}{$$}\!\left/ \!\raisebox{1ex}{$$}\right.\text{\hspace{0.17em}}\left{r}_{i,ip}\right{d}_{l,i}$
(37).
Eqs. (35) constitute a linear system of order N in the Hi (i = 1, …, N) unknowns, whose matrix is sparse and symmetric. Generally the nonzero entries are about 5% or less of the total number of matrix coefficients. The diagonal and offdiagonal matrix coefficients are, respectively
${s}_{i,i}=\frac{{A}_{i}}{\Delta t}+{\displaystyle \sum _{l=1,N{l}_{i}}\frac{{\tilde{D}}_{i}}{\left{r}_{i,ip}\right}{d}_{l,i}}\text{\hspace{0.17em}}$
${s}_{i,ip}=\frac{{\tilde{D}}_{i}}{\left{r}_{i,ip}\right}\text{\hspace{0.17em}}\text{\hspace{0.05em}}{d}_{l,i}\text{\hspace{0.17em}}{\delta}_{i,ip}$
i, i_{p} = 1, ..., N (38),
where
${\delta}_{i}{}_{,ip}=\text{}1$
if nodes i and i_{p} are linked, otherwise it is equal to zero.
According to Eqs. (37)(38), the offdiagonal coefficients are nonnegative if the GD mesh property is not satisfied for side r_{i,ip} and the Mmatrix property is lost. In this case, d_{l,i} is negative (see Eq. (10) and Appendix B) and the sign of
${\tilde{F}}_{l,i}$
is not consistent with the difference
$\left({\tilde{H}}_{ip}^{k+1}{\tilde{H}}_{i}^{k+1}\right)$
in Eq. (37). On the opposite, if side r_{i,ip} satisfies the GD property, d_{l,i} is positive, and the sign of the flux in Eq. (37) only depends on the difference of the water levels in cells i and i_{p}. If the GD property is satisfied for all the sides of the mesh, the matrix of system (35) is a Mmatrix, positivedefinite and, according to [41], the system is wellconditioned.
The numerical scheme adopted for the solution of the DC problem has a structure similar to the P1 conforming FE Galerkin scheme, and we present and discuss analogies and differences between the two numerical formulations in Appendix C.
We solve system (35) by a preconditioned conjugate gradient method (with incomplete Cholesky preconditioning). Due to the symmetry of the matrix of system (35), we only save the nonzero entries of the lower triangular matrix, using a Symmetric Coordinate Storage algorithm [43 ].
After the solution of system (35) for the H unknowns, we compute the gradients of the water level
${\nabla}_{x\left(y\right)}{\tilde{H}}_{i}^{k+1}$
at the end of the DC step as shown in the next section. As specified in section 4.1, the values of the gradients of the water level are assumed constant during the solution of the two convective steps of the hydrostatic problem in the next time iteration [31,32]. Finally, the components of the flow rate
${\tilde{q}}_{s,i}^{k+1}$
(s = x, y) are updated according to Eqs. (31)(32).
According to Eq. (37) and the definition of the coefficient
${\tilde{D}}_{i}$
in (Eq. (32)), the flux
${\tilde{F}}_{l,i}$
leaving from cell i (i_{p}) to i_{p} (i) is equal to the flux entering cell i_{p} (i) from cell i (i_{p}). If
${\tilde{F}}_{l,i}$
is leaving from i to i_{p} or vice versa depends on the difference
$\left({\tilde{H}}_{i}^{k+1}{\tilde{H}}_{ip}^{k+1}\right)$
(see Eq. (37).
Computation of the spatial gradients of the water level
We compute the gradient of H in cell i as follows. Call
${\widehat{n}}_{q}$
the unit vector parallel to the flow vector in cell i. Let
$\nabla {H}_{i,ip}$
be the vector computed as (i_{p} is one of the cells linked with i),
$\nabla {H}_{i,ip}=\left({H}_{i}{H}_{ip}\right)/{r}_{i,ip}$
(39).
We assume
$\nabla {H}_{i}=\frac{{\displaystyle {\sum}_{ip=1,N{l}_{i}}\mathrm{max}\left(0,\left({\widehat{n}}_{i,ip}\cdot {\widehat{n}}_{q}\right)\nabla {H}_{i,ip}\right)}}{{\displaystyle {\sum}_{ip=1,N{l}_{i}}\mathrm{max}\left(0,\left({\widehat{n}}_{i,ip}\cdot {\widehat{n}}_{q}\right)\right)}}$
(40),
where
${\widehat{n}}_{i,ip}\cdot {\widehat{n}}_{q}$
is the dot product between the direction of side r_{i,ip} and the unitary flow rate vector. If cell i lies on an impervious boundary, the flow rate vector is parallel to the boundary and
$\nabla {H}_{i}$
, computed as in Eq. (40), is projected along the boundary direction.
The expression in Eq. (40) can be regarded as an upwind formulation, where
$\nabla {H}_{i}$
depends on the differences between Hi and the values of H in the downstream (along the flow direction) i_{p} cells. This is consistent with the procedure adopted in the (CP + CC) steps, where the unknown fluxes/momentum fluxes are the ones leaving from i to the downstream i_{p} cells.
The NonHydrostatic Problem
The divergencefree continuity Eq. (1) is integrated to the control volume
${\tau}_{i}$
corresponding to the computational cell i,
${\tau}_{i}=({A}_{i}\times {\tilde{h}}_{i}^{k+1})$
and, due to the divergence theorem, we obtain the following balance relation between the horizontal and the vertical fluxes,
${\iint}_{{\Sigma}_{i}}\left({U}_{h}\cdot {\widehat{n}}_{\Sigma}\right)\text{\hspace{0.05em}}\text{\hspace{0.05em}}d\sigma}+{A}_{i}\left({w}_{s,i}{w}_{b,i}\right)=0$
(41.a),
where
${\Sigma}_{i}$
is the vertical surface of
${\tau}_{i}$
(
${\sum}_{i}{\tilde{h}}_{i}^{k+1}\times {\Gamma}_{i}$
with
${\Gamma}_{i}$
defined in section 3.2), U_{h} = (U, V)^{T} and
${\widehat{n}}_{\Sigma}$
is the unitary vector perpendicular to
${\Sigma}_{i}$
, positive outward. Eq. (41.a) can be written as
${\int}_{{\Gamma}_{i}}\left(\left({U}_{h}\text{\hspace{0.05em}}{\tilde{h}}_{i}^{k+1}\right)\cdot {\widehat{n}}_{\Sigma}\right)\text{\hspace{0.05em}}\text{\hspace{0.05em}}dl}+{A}_{i}\left({w}_{s,i}{w}_{b,i}\right)=0$
(41.b).
Let
$\sum}_{l=1,N{l}_{i}}{\tilde{F}}_{l,i$
be the first term on the left hand side (l.h.s.) of Eq. (41.b), with
${\tilde{F}}_{l,i}$
the flux across the lth side of
${\Gamma}_{i}$
, computed by Eq. (18) at the end of the hydrostatic DC step (with
${\tilde{q}}_{s}^{k+1}$
, s = x, y).
According to Eq. (7), the z momentum equation (Eq. (6.c)) becomes
$\frac{\partial \left(Wh\right)}{\partial t}=\frac{1}{2}\left(\frac{\partial \left(h{w}_{s}\right)}{\partial t}+\frac{\partial \left(h{w}_{b}\right)}{\partial t}\right)={q}_{b}$
(42).
After discretization in time, we approximate Eq. (42) as
$\frac{{\tilde{h}}^{k+1}}{2\Delta t}\left({w}_{s}^{k+1}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{w}_{s}^{k}\text{\hspace{0.05em}}+\text{\hspace{0.05em}}{w}_{b}^{k+1}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{w}_{b}^{k}\right)={q}_{b}^{k+1}$
(43),
the bottom vertical velocity as (see the first of Eqs. (3))
${w}_{b}^{k+1}\simeq {\tilde{w}}_{b}^{k+1}=\left(\frac{{\left(\tilde{U}\tilde{h}\right)}^{k+1}}{{\tilde{h}}^{k+1}}\right)\frac{\partial {z}_{b}}{\partial x}+\left(\frac{{\left(\tilde{V}\tilde{h}\right)}^{k+1}}{{\tilde{h}}^{k+1}}\right)\frac{\partial {z}_{b}}{\partial y}$
(44),
and the vertical velocity at the free surface is computed by Eq. (43), as
${w}_{s}^{k+1}={w}_{b}^{k}{\tilde{w}}_{b}^{k+1}+{w}_{s}^{k}+\frac{\text{\hspace{0.17em}}2\Delta t\text{\hspace{0.17em}}{q}_{b}^{k+1}}{{\tilde{h}}^{h+1}}$
(45).
The horizontal, dynamic components of the momentum equations are
$\frac{\partial {q}_{s}}{\partial t}+\frac{1}{2}\frac{\partial \left({q}_{b}h\right)}{\partial s}+{q}_{b}\frac{\partial {z}_{b}}{\partial s}=0$
s = x, y (46).
After discretization in time, Eqs. (46) become
${q}_{s}^{k+1}\text{\hspace{0.05em}}=\text{\hspace{0.05em}}{\tilde{q}}_{s}^{k+1}\frac{\Delta t}{2}\left({\tilde{h}}^{k+1}\frac{\partial {q}_{b}^{k+1}}{\partial s}+{q}_{b}^{k+1}\left(\frac{\partial {\tilde{H}}^{k+1}}{\partial s}+\frac{\partial {z}_{b}}{\partial s}\right)\right)=0\text{}s=\{\begin{array}{c}x\\ y\end{array}\text{}\Rightarrow \text{}{q}_{s}^{k+1}=\{\begin{array}{c}{\left(Uh\right)}^{k+1}\\ {\left(Vh\right)}^{k+1}\end{array}$
(47).
Symbol (
$\cdot $
)k+1 marks the values of the variables at the end of the nonhydrostatic step.
Substituting Eqs. (44)(45) and Eqs. (47) in Eq. (41,b), using the previous symbols, we get the following equation in the q_{b} unknown
$\begin{array}{c}{\stackrel{\u2323}{D}}_{i,1}{\displaystyle \sum _{l=1,N{l}_{i}}\frac{\partial {q}_{b}^{k+1}}{\partial {n}_{i,ip}}{d}_{l,i}}\frac{\Delta t}{2}\text{\hspace{0.17em}}{q}_{b,i}^{k+1}\text{\hspace{0.05em}}{\displaystyle \sum _{l=1,N{l}_{i}}\left(\frac{\partial {\tilde{H}}^{k+1}}{\partial {n}_{i,ip}}+\frac{\partial {z}_{b}}{\partial {n}_{i,ip}}\right){d}_{l,i}}+\frac{2\Delta t}{{\tilde{h}}_{i}^{k+1}}{q}_{b,i}^{k+1}{A}_{i}=\\ {\displaystyle \sum _{l=1,N{l}_{i}}{\tilde{F}}_{l,i}}{A}_{i}\left({w}_{s,i}^{k}+{w}_{b,i}^{k}2{\tilde{w}}_{b,i}^{k+1}\right)\end{array}$
i = 1, ..., N (48),
where
${\stackrel{\u2323}{D}}_{i,1}=\frac{\Delta t}{2}{\tilde{h}}_{i}^{k+1}\text{if}{\varphi}_{i}^{k}\ge {\varphi}_{ip}^{k}\text{}{\stackrel{\u2323}{D}}_{i,1}=\frac{\Delta t}{2}{\tilde{h}}_{ip}^{k+1}\text{if}{\varphi}_{ip}^{k}{\varphi}_{i}^{k}$
(49),
and
$\partial {q}_{b}/\partial {n}_{i,ip}$
is the q_{b} gradient along the
${\widehat{n}}_{i,ip}$
direction. We set
${\stackrel{\u2323}{F}}_{l,i}={\stackrel{\u2323}{D}}_{i,1}\frac{\partial {q}_{b}^{k+1}}{\partial {n}_{i,ip}}{d}_{l,i}$
(50),
and it can be considered a (corrective) flux due to the difference of q_{b} in cell i and in the neighboring i_{p} one (its units are cubic meters per seconds, i.e., the ones of a volumetric flux).
As for the hydrostatic DC step, we adopt a linear P1 spatial approximation of the q_{b} unknowns according to the nodal values, and the computational cell is the same as in the hydrostatic problem.
Similar to Eq. (37), we rewrite Eq. (50) as
${\stackrel{\u2323}{F}}_{l,i}={\stackrel{\u2323}{D}}_{i,1}\left({q}_{b,ip}^{k+1}{q}_{b,i}^{k+1}\right)\raisebox{1ex}{$$}\!\left/ \!\raisebox{1ex}{$$}\right.\text{\hspace{0.17em}}\left{r}_{i,ip}\right{d}_{l,i}$
(51).
As in the DC problem, if the GD property is satisfied for side r_{i,ip}, d_{l,i} is positive (see Eq. (10) and Appendix B), and the sign of
${\stackrel{\u2323}{F}}_{l,i}$
only depends on the difference of q_{b} in the two cells i and i_{p}, while the opposite occurs if the GD property is not satisfied.
Eqs. (48) form a linear system of order N in the q_{b,i} (i = 1, …, N) unknowns. The diagonal and offdiagonal entries of the matrix are, respectively
${s}_{i,i}=\frac{2\Delta t\text{\hspace{0.17em}}{A}_{i}}{{\tilde{h}}_{i}^{k+1}}+{\displaystyle \sum _{l=1,N{l}_{i}}\frac{{\stackrel{\u2323}{D}}_{i,1}}{\left{r}_{i,ip}\right}{d}_{l,i}\text{\hspace{0.17em}}}{\displaystyle \sum _{l=1,N{l}_{i}}\frac{\Delta t}{2}\left({\nabla}_{n}{\tilde{H}}^{k+1}+{\nabla}_{n}{z}_{b}\right)\text{\hspace{0.05em}}\text{\hspace{0.17em}}{d}_{l,i}}$
${s}_{i,ip}=\frac{{\stackrel{\u2323}{D}}_{i,1}}{\left{r}_{i,ip}\right}{d}_{l,i}{\delta}_{i,ip}$
i, i_{p} = 1, ..., N (52),
where
${\nabla}_{n}{\tilde{H}}^{k+1}$
and
${\nabla}_{n}{z}_{b}$
are the spatial gradients of H and z_{b} computed along the direction
${\widehat{n}}_{i,ip}$
. The gradient of z_{b} in cell i are computed as the gradient of H (see section 4.3). The other symbols are specified as above. The source term of system (48) is
${b}_{i}={\displaystyle {\sum}_{l=1,N{l}_{i}}{\tilde{F}}_{l,i}}{A}_{i}\left({w}_{s,i}^{k}+{w}_{b,i}^{k}2{\tilde{w}}_{b,i}^{k+1}\right)$
(53).
Once q_{b} is computed in each cell i, the horizontal flow rate components are updated according to Eqs. (47) and the velocity of the free surface according to Eq. (45).
Similarly to the DC step (section 4.2), the matrix of system (48) is sparse, symmetric and, for a GD mesh, positivedefinite, and system (48) is solved as system (35) (see section 4.2).
As for the corrective flux
${\tilde{F}}_{l,i}$
in the hydrostatic DC step, according to Eq. (51) and the definition of the coefficient
${\stackrel{\u2323}{D}}_{i,1}$
in (Eq. (49)), the flux
${\stackrel{\u2323}{F}}_{l,i}$
leaving from cell i (i_{p}) to i_{p} (i) is equal to the flux entering cell i_{p} (i) from cell i (i_{p}). If
${\stackrel{\u2323}{F}}_{l,i}$
is leaving from i to i_{p} or vice versa depends on the difference
$\left({q}_{b,i}^{k+1}{q}_{b,ip}^{k+1}\right)$
(see Eq. (51)).
Boundary Conditions (B.C.)
As specified in section 3, at the free surface, the stress of the wind is neglected and the (atmospheric) pressure is set to zero. The friction at the bottom is approximated by the Manning equation, (see Eqs. (6.a)(6.b)). We assign free slip condition at the impervious walls.
During the CP hydrostatic problem, we deal with the B.C. for those cells with incoming and leaving flux/momentum fluxes as in section 4.4 of [33] (for brevity we refer the reader to the mentioned paper). The external assigned incoming fluxes/momentum fluxes are considered in the CP step [31,33]. During the CC step, in the inlet/outlet cells, we set the incoming fluxes and momentum fluxes from the neighboring cells equal to the leaving ones. This means that, in the boundary cells, the solution of the convection prediction problem is not altered by the convection correction problem. The B.C. in the DC step is handled as in [33]. We set q_{b} = 0 in the inlet and outlet cells and this represents a Dirichlet condition for the solution of system (48).
Model Properties
Preservation of the local and global mass balance
The hydrostatic and dynamic problems involved in the proposed procedure rely on the local mass balance. This is imposed, for each cell i, after the solution of the ODEs system of the governing equations in the CP and CC steps (see Eq. (24)). The systems solved in the hydrostatic DC and dynamic steps are derived by a local mass balance, obtained by embedding the momentum equations in the continuity equation (see Eqs. (33)(34) for the DC step and Eqs. (41) and (48) for the dynamic problem). We show that the proposed solver also preserves the global mass balance in all its steps.
Summing all Eqs. (24) at the end of the CP problem, we get
${\sum}_{i=1,N}{A}_{i}\left({h}_{i}^{k+1/3}{h}_{i}^{k}\right)}=\Delta t{\displaystyle {\sum}_{i=1,N}\left({\overline{F}}_{i}^{in}{\overline{F}}_{i}^{out}\right)$
(54),
and, due to Eq. (25,b), for each internal couple of linked cells i and i_{p}, the r.h.s. of Eq. (54) is equal to the difference between the boundary incoming and leaving fluxes. i.e.,
${\sum}_{i=1,N}{A}_{i}\left({h}_{i}^{k+1/3}{h}_{i}^{k}\right)}=\Delta t\left({\overline{F}}_{bou}^{in}{\overline{F}}_{bou}^{out}\right)$
(55,a),
with
${\overline{F}}_{bou}^{in}$
and
${\overline{F}}_{bou}^{out}$
the sum of the incoming (assigned) and leaving boundary fluxes. In a similar way, at the end of the CC step, due to the continuity of the meanintime flux between i and i_{p} (similarly to Eq. (25.b)) and since we handle the assigned boundary incoming fluxes
${\overline{F}}_{bou}^{in}$
in the CP step (see section 6), after simple manipulations, we can easily show that
${\sum}_{i=1,N}{A}_{i}\left({h}_{i}^{k+2/3}{h}_{i}^{k}\right)}=\Delta t\left({\overline{F}}_{bou}^{in}{\overline{F}}_{bou}^{out}\right)$
(55.b),
where
${\overline{F}}_{bou}^{out}$
is the sum of the leaving boundary flux computed at the end of the CC problem. In the hydrostatic DC problem, as specified at the end of section 4.2, the flux
${\tilde{F}}_{l,i}$
leaving from cell i (i_{p}) to i_{p} (i) is equal to the flux entering cell i_{p} (i) from cell i (i_{p}). This implies that the sum of all these fluxes is zero and, after simple manipulations, it can be show that
${\sum}_{i=1,N}{A}_{i}\left({\tilde{h}}_{i}^{k+1}{h}_{i}^{k+2/3}\right)}=\Delta t\text{\hspace{0.17em}}{F}_{dir$
(55.c),
with F_{dir} the sum of the boundary fluxes computed at the cells with imposed Dirichlet conditions.
If we add Eq. (55,c) to Eq. (55,b), we obtain
${\sum}_{i=1,N}{A}_{i}\left({\tilde{h}}_{i}^{k+1}{h}_{i}^{k}\right)}=\Delta t\left({\overline{F}}_{bou}^{in}{\overline{F}}_{bou}^{out}+{F}_{dir}\right)$
(55.d),
which can be assumed as a relationships representing the global mass balance at the end of the hydrostatic problem. Because for each couple of linked internal cells i and i_{p}, the corrective dynamic flux
${\stackrel{\u2323}{F}}_{l,i}$
leaving from i (i_{p}) to i_{p} (i) (in Eq. (51)) is equal to the flux entering i_{p} (i) from i (i_{p}) (see also the consideration at the end of section 5), the sum of all these internal fluxes is zero. On the other hand, the corrective boundary flux of the cells at the boundaries is zero, since we set q_{b} = 0 as boundary Dirichlet assigned value (see section 6). This means that the nonhydrostatic step does not affect the global mass balances.
Preservation of a general steady gtate (GSS) condition
We assume that the GSS condition is achieved if the two relationships in Eq. (56) hold in each cell:
$\frac{\partial {h}_{i}}{\partial t}=0$
and
$\frac{\partial {q}_{s,i}}{\partial t}=0$
s = x, y i = 1, ..., N (56).
The numerical procedure proposed for the solution of the (CP + CC) steps for a generic cell relies on (see sections 4 and 5): 1) the balance of the incoming and leaving fluxes with the time change of the storage term, for the continuity equation and, 2) the balance among the incoming and leaving momentum fluxes with the time change of the flow rates, with the resultant of the hydrostatic pressure distribution, along cell interfaces and bottom surfaces, and the resultant of the bottom friction term, for the momentum equations.
Given an initial condition in the domain, h_{i}, and q_{s}_{,}_{i}, i = 1, ..., N and s = x, y, the above considerations imply the two following conditions,
a) If the incoming fluxes in any cell i are equal to the leaving ones, then h does not change in time, i.e., from Eqs. (21) and (26),
$\text{if}\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{\displaystyle \sum _{l=1,N{l}_{i}}{\delta}_{l,i}{F}_{l,i}^{r,out}\left(t\right)dt}}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{F}}_{l,i}^{r,in}}\text{}\Rightarrow \text{}\frac{d{h}_{i}}{dt}{A}_{i}=0$
with symbol r = p, c (57),
b) If the incoming and leaving momentum fluxes in any cell i are balanced by the resultant of the hydrostatic pressure distribution along the cell interfaces and bottom surfaces, as well as by the resultant of the bottom stress term, then the specific flow rates in cell i do not change in time, i.e., from Eqs. (22) and (27),
$\text{if}\left(\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{R}_{i}^{r,s}\left(t\right)dt}\right){\delta}_{r}+\frac{1}{\Delta t}{\displaystyle \underset{\Delta t}{\int}{\displaystyle \sum _{l=1,N{l}_{i}}{\delta}_{l,i}{M}_{l,i}^{r,s,out}\left(t\right)dt}}={\displaystyle \sum _{l=1,N{l}_{i}}\left(1{\delta}_{l,i}\right){\overline{M}}_{l,i}^{r,s,in}\left(t\right)}\text{}\Rightarrow \text{}\frac{d{q}_{s,i}\left(t\right)}{dt}{A}_{i}=0$
(58),
with
${\delta}_{r}=\text{}1$
in the CP step,
${\delta}_{r}=\text{}0$
in the CC step. Apex r in Eqs. (57)(58) is equal to "p" (in the convective prediction step) or "c" (in the convective correction step) and in this second case the source term R_{i} is equal to zero (see section 4.1).
Moreover, due to the fixed bed assumption, if the water depths in all the i cells do not change in time, neither the water levels nor their spatial gradients will do. We will show that, if conditions a) and b) of the (CP + CC) steps hold, the hydrostatic DC step and the dynamic problem do not modify this steady state condition for each cell. If the flow rate components do not change in time during the (CP + CC) steps (condition b) holds), the second sum on the r.h.s. of Eq. (33) in the DC step is zero, and, after simple manipulations, the same equation can be written as:
$\underset{{A}_{i}}{\int}\frac{\partial {\eta}_{i}}{\partial t}d{A}_{i}}{\displaystyle \sum _{l=1,N{l}_{i}}{\tilde{D}}_{i}\frac{\partial \left(\eta \vartheta \right)}{\partial {n}_{l,i}}{d}_{l,i}}=0$
(59),
with
$\eta ={\tilde{H}}^{k+1}{H}^{k}{}^{+2/3}$
and
$\theta ={H}^{k}{H}^{k}{}^{+2/3}$
. If h (H) in all the i cells does not change in time during the (CP + CC) steps (condition a) holds),
$\theta $
is zero, and the solution of Eq. (59) is identically
$\eta =0$
, i.e., zero water level correction. From Eq. (31) we get that the correction of the specific flow rate components after the DC step is zero. In brief, if conditions a) and b) hold, we get:
c)
${\tilde{H}}^{k+1}={H}^{k}{}^{+2/3}={H}^{k}$
for each cell i,
d)
${\tilde{q}}_{s}^{k+1}={q}_{s}^{k+2/3}={q}_{s}^{k}$
, s = x, y for each cell i.
We prove that the correction due to the nonhydrostatic problem is zero as well. We focus our attention on the r.h.s. of Eq. (48). We find the two following conditions e) and f),
e) If conditions b) and d) hold, the sum of the horizontal fluxes is zero for any i cell, i.e.,
${\sum}_{l=1,N{l}_{i}}{\tilde{F}}_{l,i}}=0$
,
f) If conditions a) to d) hold,
${\tilde{q}}_{s}^{k+1}/{\tilde{h}}^{k+1}={q}_{s}^{k}/{h}^{k}$
(s = x, y).
From condition(f) and Eqs. (3), we obtain that
${\tilde{w}}_{s\left(b\right)}^{k+1}={w}_{s\left(b\right)}^{k}$
. Integrating the divergencefree continuity equation (Eq. (1)) to a control volume corresponding to the generic computational cell (see section 5), From condition e) we get
${\tilde{w}}_{s}^{k+1}={\tilde{w}}_{b}^{k+1}$
.
According to these considerations, the r.h.s. of Eqs. (48) is zero and this leads to q_{b} = 0. This implies that the dynamic correction of the flow rate components, with respect to the initial condition, is zero and it is consistent with the second condition in Eq. (54).
The "water at rest" condition [44] is a particular case of the general situation described as above, where the flow rate components q_{s}_{,}_{i} (s = x, y) are zero for all the i cells. In this case, according to the argumentations as above, it easy to prove the Cproperty of the present model [44].
Handling of wetting/drying processes
In [32] the authors show 1) how negative water depths do not arise during the solution of the (CP + CC) steps and 2) how to handle small negative values of water depth calculated after the solution of the linear system of the DC step. If, after the solution of the hydrostatic step, we have, in cell i,
${h}_{i}\le 0$
, we set q_{b} = 0, as Dirichlet value in system (48).
Study of the frequency dispersion
We compare the phase speed of the proposed model with the one given by the linear wave theory. The phase speed from the linear wave theory is [34]:
${c}_{1}=\sqrt{gh}\sqrt{\left(\mathrm{tanh}\left({k}_{w}h\right)/{k}_{w}h\right)}$
(60.a)
where k_{w} is the wave number (
${k}_{w}=\text{}2\pi /{L}_{w},{L}_{w}$
is the wave length). The wave period T_{w} is equal to L_{w}/c1. We assume a standing regular wave with L_{w} = 20 m in a closed basin, 40 m long (= 2 L_{w}) and 1 m wide, with frictionless and horizontal bottom. The wave profile is:
$\eta \left(x,0\right)={a}_{w}\mathrm{cos}\left(2\pi x/{L}_{w}\right)$
(60,b),
where
$\eta $
is the distance between the still water level and the free surface and a_{w} is the wave amplitude. We keep constant the value for a_{w}, equal to 0.1 m and h ranges from 1 m to 20 m. This experiment is similar to the one performed in [11,16,18], where the authors assumed a channel length equal to L_{w} /2. We discretize the domain with equilateral triangles with side 0.1 m. The wave period is divided in 160 parts and the total simulation time is 30 Tw. We registered the computed
$\eta $
elevation at two gauges, placed along the centerline of the channel, far L_{w} /2 and 3/2 L_{w} from the upstream end, so that the distance between the two gauges is equal to the wave length. We compute the numerical phase speed c_{n} = L_{w}/T_{w}_{,}_{n}, where T_{w}_{,}_{n} is the time measured between two consecutive wave peaks of the same gauge, the first in the exact
$\eta \left(t\right)$
profile, the second in the computed
$\eta \left(t\right)$
profile (Figure 3.a). In Figure 3.b we compare the numerical and the linear wave theory phase speeds versus kh, and we also plot the phase celerity analytically computed for the shallow waters limit (e.g., c = (gh)^{1/2}. The numerical phase speed plotted in Figure 3.b is obtained as the mean value between the phases speed computed for the two gauges, which have very similar values. The numerical phase celerity begins to deviate from the linear theory at about kh = 2.74, i.e., for h/L_{w} = 0.436, corresponding to the "intermediate waters" [34]), similarly to the weakly dispersive BTMs [11,45]. Setting T_{w}_{,}_{n} as the time interval between two consecutive wave peaks computed at gauges G1 and G2, we obtained results very similar to the ones in Figure 3.b. A small increment of c_{n} has been computed between kh = 4  5 (Figure 3.b). According to the numerical results, in this interval, the numerical solution progressively deviates from the linear theory, and a slight increment of the phase speed is balanced by a reduction of the computed wave peaks. For kh > 5, c_{n} monotonically decreases.
Figure 3: a) Computation of the numerical phase celerity. b) Phase speed vs. water depth.
Model Applications
Test 1: Propagation of a solitary wave in a frictionless and horizontal channel
We deal with a solitary wave, travelling in a horizontal and frictionless channel, 850 m long and 10 m wide. The initial water level h_{0} is 10 m, the amplitude wave a_{w} is 2 m and the wave is initially centered at x = 200 m. The initial conditions are given by the analytical solution [11,13,37].
$\eta \left(x,t\right)=\frac{{a}_{w}}{{\mathrm{cosh}}^{2}\left({k}_{w}\left(x{x}_{0}ct\right)\right)}$
$U\left(x,t\right)=\frac{\eta c}{\eta +{h}_{0}}$
${w}_{s}\left(x,t\right)={z}_{z={h}_{0}+\eta}\frac{\partial U}{\partial x}$
(61),
where
$\eta $
is the elevation of the free surface measured from h_{0}, c is the wave celerity
$c=\sqrt{g\left({h}_{0}+{a}_{w}\right)}$
, k_{w} is the wave number
${k}_{w}=\sqrt{\left(3{a}_{w}/4{h}_{0}^{3}\right)}$
and the other symbols are specified as above. Left and right sides of the channel are open boundaries; the other two sides are impervious walls. The GD mesh used to discretize the channel has 160925 triangles and 84409 nodes. The mesh generator Netgen [46] has been used to create the meshes in all the applications presented in this paper. The size of the time step
$\Delta t$
is 0.025 s. The maximum value of the CFL number we computed is 1.925, where the CFL number for cell i is obtained as
$CFL=\left(\sqrt{{\left(Uh\right)}_{i}^{2}+{\left(Vh\right)}_{i}^{2}}/{h}_{i}+\sqrt{g{h}_{i}}\right)\Delta t/\sqrt{{A}_{i}}$
(62).
In Figure 4 we show the initial condition for
$\eta $
, U and w_{s} and the solutions computed after 10, 20 and 50 s. The computed absolute values of the peaks slightly decrease at the very beginning, since the initial condition is given by an analytical solution and other authors observed similar trends [8,18]. The computed peak values of
$\eta $
and U after 10 s are 1.926 m and 1.677 m/s, respectively. The maximum and minimum computed values of w_{s} at the same time are 0.546 and 0.552 m/s, respectively. These values and shape of the computed profiles are almost saved at 50 s.
Figure 4: Test 1. H (top), U (middle) and w_{s} (bottom). Initial solution and computed solutions after 10, 20 and 50 s.
Test 2: Moving shoreline inside a parabolic bowl. Investigation of the computational costs and convergence order
Sampson et al. [47] provided for this case the analytical solution for H, U and V (see Eqs. (39)(40) in [48]) The domain is [8000]^{2} m^{2} with symmetric bottom with respect to the centre of the squared domain (see Eq. (38) in [48]). The GD mesh we use for spatial discretization has 272 elements and 149 nodes, and
$\Delta t$
is 40 s (we adopt the same mesh and time discretizations as in a previous paper [33], where we present a hydrostatic NLSEWs solver). The maximum CFL value obtained during the simulation is 2.82. We refine the computational mesh three times. Each refinement is performed by halving each sides of the previous coarser mesh (i.e., each triangle is divided in four equal ones). We limit the growth of the CFL number by halving the time step at each refinement. We compute the L_{2} norms of the relative errors of h, Uh and Vh with respect to the exact values (Table 1). The relative error for the mesh level l, err_{l} and the rate of convergence r_{c}, are computed as in Eqs. (70,a) and (70,b) in [33], respectively. The values of the rc are close to 2 (Table 1). The reason is that inside each triangle of the mesh, h (H), Uh and Vh change linearly, according to the three nodal values (2_{nd} spatial approximation order). Refining the mesh, r_{c} does not change significantly and this implies that the numerical model computes accurate results also if the domain in discretized using a coarse mesh instead of a very fine one.
Refinement
Level 
N_{T} 
N 
L_{2,h} 
r_{c,h} 
L_{2,Uh} 
r_{c,Uh} 
L_{2,Vh} 
r_{c,Vh} 
0 
272 
159 
3.71E02 

5.16E+00 

3.31E+00 

1 
1088 
589 
1.10E02 
1.75E+00 
1.60E+00 
1.69E+00 
1.00E+00 
1.73E+00 
2 
4352 
2265 
2.90E03 
1.92E+00 
4.60E01 
1.80E+00 
2.70E01 
1.89E+00 
3 
17408 
8881 
7.46E04 
1.96E+00 
1.20E01 
1.94E+00 
7.00E02 
1.95E+00 
Table 1: Test 2. L_{2} norms of relative errors and rate of convergence r_{c}.
We use this test to investigate the computational costs (CPU times). Table 2 shows the specific mean values of the CPU times, CPU_{s}, required for cells ordering and for the solution of the (CP+CC) steps, DC step and nonhydrostatic problem. CPU_{s} are computed by dividing the total CPU time required for the solution of the different steps of the algorithm, by the number of the cells and the time iterations. We use a single processor Intel CORE i74770 at 3.40 GHz. The (CP+CC) steps are the most demanding ones, and the small decrement of the CPU_{s} can be motivated as in section 5.8 in [33]. The DC step, the nonhydrostatic problem and the cells ordering step require CPU_{s} approximately one magnitude order less than the one required for (CP+CC) steps.
We also measure the mean CPU times per iteration, obtained by dividing the total CPU times of the different steps of the algorithm, by the number of the time iterations. We compute the growth rate
$\beta $
of the mean CPU time as in Eq. (71) in [33]. Results are shown in Figure 5 and the growth exponents
$\beta $
are close to 1. The value of the (CP + CC) steps, slightly less than 1, can be due to the small decrement of the computational effort observed with increasing the cells number, discussed in [32]. In the same figure we also show the total CPU times, obtained by adding the mean CPU times of the four algorithm steps, and the value of its
$\beta $
is very close to 1.
Refinement Level 
N_{T} 
Cells Ordering 
CP + CC 
DC 
NonHydrostatic 
0 
272 
2.81E07 
1.19E05 
7.64E07 
5.36E07 
1 
1088 
3.15E07 
1.21E05 
8.63E07 
6.15E07 
2 
4352 
3.33E07 
1.18E05 
8.44E07 
6.72E07 
3 
17408 
3.40E07 
1.17E05 
9.05E07 
6.76E07 
Table 2: CPUs times (seconds).
Figure 5: Mean CPU times and exponents
$\beta $
.
In Table D.1 (in Appendix D), we show the L_{2} norms of the relative errors of h, Uh and Vh, as well as the corresponding convergence orders, computed by the present solver without the nonhydrostatic solver. The norms are slightly greater than the ones obtained with the present nonhydrostatic solver, probably because of the lack of the correction of the flow field due to the components of the dynamic pressure. It is also important to compare the values in Table D.1 with the ones computed by the hydrostatic solver in [33] (Table 1 in [33].). In the solver of the referred paper, a cellcentered formulation is adopted, with a constant (1^{st} order) spatial approximation order of the unknown variables h, Uh and Vh inside the computational cells (the triangles of the mesh). Changing from a nodecentered to a cellcentered formulation, the number of the computational cells and unknowns increases. So we would expect a reduction of the errors in the solution of the cellcentered formulation. The reason of the smaller values of the errors of the present algorithm, compared to the ones in [33], could be due to the difference of the spatial approximation order of the unknown variables in the two adopted formulations.
Another benefit of the present solver with respect to the one in [33], is the abatement of the mean CPU times. In Table D.2 (in Appendix D), we list the CPU_{s} of the solver in [33], a bit smaller than the ones in Table 1 of the present model. The reason could be the higher number of neighboring computational cells in the nodecentered formulation, with respect to the three (or less, for the boundary triangles) in the cellcentered formulation. Due to the greater number of the computational cells, the total CPU times required by the cellcentered formulation (in Figure D.1, Appendix D), are greater than the ones of the present nodecentered formulation (Figure 5).
Test 3: Propagation and breaking of a solitary wave over a sloping bottom
A solitary wave with ratio a_{w}/h_{0} = 0.3, travels, breaks and runs up on a sloping plane beach. The measures provided in [49] constitute the reference solution for the one provided by the present model. The channel consists in a first horizontal stretch, followed by a beach with slope 1:19.85. The solitary wave in the numerical experiment is initially centered at half wavelength L_{w} from the beach toe [49],
${L}_{w}=2/{k}_{w}\text{arccosh}\left(\sqrt{20}\right)$
(63),
For the numerical experiments, we use a 120 m long and 1.2 m wide channel, discretized with a GD mesh with 24234 triangles and 13282 nodes. The toe of the beach is 60 m far from the left channel end. We assign zero flux at the right side of the domain, and the boundary conditions at the other sides of the domain, as well as the initial conditions, are as in test 1. The size of the time step is 0.016 s and the maximum computed value of CFL is 1.17. In Figure 6 we compare the simulated and the measured free surface elevations. The length and time scales are normalized by the initial still water level h_{0} and (h_{0}/g)^{1/2}, respectively. When the wave propagates over the inclined beach, its front skews and, according to the experimental data, it breaks around t(g/h_{0})^{1/2} = 20. This is accurately reproduced by the proposed model. Several numerical models equipped with special treatment for wave breaking, not often perform very well and fail in reproducing wave breaking in time and/or amplitude [13]. After t(g/h_{0})^{1/2} = 20, the front of the wave moves towards the shoreline and the height of the wave dramatically reduces. The wave progressively runs up to the beach and reaches the maximum height around t(g/h_{0})^{1/2} = 45. The present model simulates very well the runup process. During the drawdown, a hydraulic jump is formed approximately at t(g/h_{0})^{1/2} = 50. This is accurately simulated by the present model and, generally, good agreement is obtained till the end of the wave retreat process. During the drawdown at t(g/h_{0})^{1/2} = 55, several authors [8,13,50] compute numerical results that not perfectly match the experimental data. On the opposite, the proposed model matches very well the measures at t(g/h_{0})^{1/2} = 55. Dimakopoulos et al. [51] simulated the same test case using OpenFOAM© [3]. Their results are consistent with the experimental ones. Dimakopoulos et al. [51] discretize the domain with 180000 3D cells and their total simulation time is 10s. They use 12 CPU cores with 2.6 GHz, with a total CPU time of 2 hrs. For the simulations of the present model we use a single processor Intel Q 6600 with 2.4 GHz. The total simulation time is 25 s, requiring 425 s of computational time [49].
Figure 6: Test 3. Measured (circles) and simulated (blue solid line) surface profiles. Measures.
Test 4: Runup of a solitary wave on a conical island
We deal with the runup of a tsunami on a conical island [52]. The lab flume (30 m wide and 25 m long, in Figure 7) presents an island, i.e., a truncated circular cone, 0.625 m high, with diameters 7.2 m and 2.2 m at the toe and crest, respectively, and slope of the face 1:4. The surface of the physical model is made of smooth concrete and we assume zero Manning friction coefficient. Other details of the experiments can be found in [52].
Figure 7: Test 4. Definition sketch.
Figure 8: Test 4. Snapshot of the water surface for a_{w}/h_{0} = 0.181 at different times. The black arrow indicates the direction of the incident wave.
The initial conditions are as in tests 1 and 3, and all the boundaries of the domain are open. The domain is discretized with a GD mesh with 143622 triangles and 72467 nodes. The initial water depth value h0 in our numerical run is 0.32 m and we studied three cases with relative wave heights a_{w}/h_{0} = 0.045, 0.096 and 0.181, respectively. These experiments have been widely tested in the literature [8,13,17,39,45,53,54].
For brevity, we shown in Figure 8 some snapshots of the wave evolution only for a_{w}/h_{0} = 0.181, where the dispersive waves are more pronounced, compared with the other two cases. In Figure 8a the wave reaches the maximum height of runup over the front of the island. After that, the wave withdraws below the initial water level. The refracted waves partially propagate around the island towards the downstream side and two trapped waves are generated at both sides of the cone. These two waves interact at the back face of the island and generate a second runup. Then, water waves withdraw from the back face of the island and pass through each other, moving around the cone. Figure 8e shows the first group of the dispersive waves wrapping around the island and colliding on the lee side. Figure 8f shows the second group of the dispersive waves after the collision on the back side. Similar observations have been provided in [8,54].
We compare the computed and measured water levels at gauges WG 6, 9, 16 and 22 in Figure 9. The data measured at WG 3 are used for time adjustment of the computed solutions [13]. Generally, the numerical results computed by the present solver match satisfactorily the measured time series. The rundown at the front face of the island, following the first runup, is properly simulated. The differences between numerical and experimental data become more evident at later times, especially for a_{w}/h_{0} = 0.181. Fuhrman & Madsen [53] assert that during the experiments, the front of the wave has been generated more accurately than the rear. This could motivate the above discrepancies after the initial runup. The present model accurately describes the phase and the amplitude of the peak. The peak is slightly delayed at WG 22 for a_{w}/h_{0} = 0.096 and 0.181 and at WG 9 and 16 in case a_{w}/h_{0} = 0.181. The amplitude of the peak is reasonably predicted by the present model, which slightly overestimates the amplitude of the leading wave at WG 9 for a_{w}/h_{0} = 0.096 and 0.181. The depression after peak at WG 9 is slightly underestimated for a_{w}/h_{0} = 0.096 and overestimated for a_{w}/h_{0} = 0.181. Other literature models [53,54], over predict the amplitude of the peak at WG 9, 16 and 22, expecially for a_{w}/h_{0} = 0.18.
Figure 9: Test 4. Measured (circles) and simulated (blue solid line) time evolution of water level at the gauges. Measures from [8].
${a}_{w}/{h}_{0}=\text{}0.045$
(Left panel),
${a}_{w}/{h}_{0}=\text{}0.096$
(Middle panel),
${a}_{w}/{h}_{0}=\text{}0.181$
(Right panel).
Dimakopoulos et al. [51] applied the open source code OpenFOAM© [3] to reproduce the same test. They simulated half domain, discretized with 6000000 3D cells for a total simulation time 17 sec [51]. Their results are consistent with the experimental ones and very similar to the one provided by the present model. Dimakopoulos et al. [51] used 12 CPU cores with 2.6 GHz, with a total CPU time 24 hrs. The total simulation time in our run is 20 s and, over one processor Intel Q 6600 with 2.4 GHz, the total CPU time is 1795 s.
Test 5: Tsunami in the Monai valley
This benchmark concerns a laboratory scale model of the real 1993 Hokkaido NanseiOki Tsunami even, around the Monai Valley. The experimental runs have been performed in the 1:400 scale lab flume of the Central Research Institute for Electric Power Industry (CRIEPI) [55]. Field observations of the real event refer to an initial leading depression wave, usually reproduced in the laboratory by an Nwave [9]. Figure 10, 10.a & 10.b show the bathymetry of the model and the incoming Nwave, respectively [56]. This is assumed as Dirichlet condition assigned at the offshore (western) boundary of the computational domain free slip condition is imposed at the two lateral sides (impervious walls). The GD mesh has 33004 triangles and 16730 nodes and
$\Delta t$
is 0.01 s. The Manning's coefficient n is 0.012 s/m1/3 [9]. The maximum value of the computed CFL number in our simulation is 4.93. Three gauges, located in front of the coast, behind the Muen Island (coordinates (4.521 m, 1.196 m), (4.521 m, 1.696 m) and (4.521 m, 2.196 m), respectively for gauges 1, 2 and 3), have been used for the time registration of the free surface elevation. In Figures 11, 11.a to 11.c shows some snapshots. A series of waves with short period are generated, during the leading depression, by the dispersion of the incoming Nwave over the coastal relief (Figure 11a). The incident wave shoals over the plane slope of the relief and refracts and diffracts around the Muen Island (Figure 11b). After reflection at the coast, some smallamplitude dispersive waves are generated (Figure 11c). At the same time, the waves reflected at the coast interact and generate some breaking waves (Figure 11c). In Figure 12 we compare the measured and the computed elevation of the free surface. The model matches very well the registrations. The initial drawdown, the arrival time and the amplitude of the leading wave are properly reproduced by the model, as well as the reflection at the coast. The discrepancy, approximately for the first 10 s, especially at gauge 1, are due to the initial condition inside the tank of the flume, where
$\eta $
was not uniformly equal to zero [57].
Figure 10: Test 5. a) Bathymetry of the lab model. b) The Nwave assigned at the western side [56].
Figure 11: Test 5. Snapshots of the water levels at different simulation times.
Figure 12: Test 5. Measured (circles) and simulated (red solid line) water levels at the gauges. Measures from [56].
Kesservani & Liang [58] simulated the same benchmark with different numerical solvers: a Discontinuous Galerkin (RKDG2) scheme, an adaptivegrid RKDG2 and a MUSCL2 scheme with 2^{nd} spatial approximation order. They used different resolutions of the spatial discretizations with quadrilateral cells for the RKDG2 scheme [58]. The number of the elements ranged approximately from 7000 to 42000 for the simulation of the adaptivegrid RKDG2, and they used a uniform mesh with 95648 elements for the MUSCL2 scheme [58]. Kesservani & Liang [58] performed their simulations, reproducing a total duration of 25 s, using one processor Core Duo T9400 at 2.53 GHz. The total CPU times required by the models in [58] are listed in Table 4 of the same paper, and range from 0.25 hrs to 8.8 hrs for the RKDG2, 0.75 hrs for the MUSCL2 solver and 1.52 hrs for the adaptivegrid RKDG2. We used one processor Intel Q 6600 at 2.4 GHz and the total computational time was 0.425 hrs for a total simulation time of 200 s. The outputs of the present solver are similar to the results of the adaptivegrid RKDG2 and MUSCL2 schemes in Figure 14 [58], which required higher computational burden. The simulation with the RKDG2 with (98 x 61) elements in [58] obtained CPU times comparable to ours, but returned very bad solutions (Figures 14a14c).
Conclusions
A nonhydrostatic NLSWEs model has been proposed for the simulations of long waves/tsunamis processes in shallow waters regions. It is our opinion that one of the principal merit of the present solver with respect to others literature schemes proposed for the solution of the same equations [13,19,39,49,] is that in both the hydrostatic and dynamic steps, the computed flow field satisfies the depthintegrated continuity equation discretized in each computational cell. As addressed in the previous sections, this has important implications in the local and global mass balance of the numerical solver. We provide several applications involving breaking/nonbreaking waves and runup, also over irregular topography. Generally, the computed results are in good agreement with the analytical solutions or lab data and with the results provided by other literature algorithms. Another important merit of the proposed model concerns the required computational effort. The CPU times are several magnitude order less than the ones of 3D models and, generally, much less that the ones of other 2D depthintegrated NLSWEs models. All the steps of the proposed algorithm can be parallelizable, as an important future advancement for this study, since it could reduce the computational burden. The parallelization of the hydrostatic CP and CC steps could be obtained by simultaneously solving the ODEs systems in all those cells next to already solved cells with higher and lower values of the approximated potential, respectively. The solution of the systems by the conjugate gradient method, in the hydrostatic DC and dynamic steps, is naturally parallelizable [59]. As a consequence of these considerations, the proposed solver could be embedded in a more general forecasting and prevention tool for the inundation generated by tsunamis or solitary waves.
Disclosure Statement
Some details of the proposed numerical procedure are provided in the file "supplementary material.docx" (Appendices A, B and C), as well as some of the figures and tables concerning the presented test cases (Appendix D).
Acknowledgment
This research did not receive any specific grant from funding agencies in the public, commercial, or notforprofit sectors.
Conflict of Interest
The Author of the paper, Costanza Aricò, discloses any actual or potential conflict of interest including any financial, personal or other relationships with other people or organizations within three (3) years of beginning the submitted work.
References
 Goring DG (1979) TsunamisThe propagation of long waves onto a shelf. Ph.D. Thesis, CalTech, California Institute of Technology, California, USA, pp 356.
 Chanson H (2006) Tsunami surges on dry coastal plains: application of a dam break wave equation. Coast Eng Journal 48(4): 355370.
 (2013) OpenFOAM User Guide. The Architects of OpenFOAM, OpenFOAM Foundation.
 Dalrymple RA, Rogers BD (2006) Numerical modeling of water waves with the SPH method. Coastal Engineering 53(23): 141147.
 Hirt CW, Nichols BD (1981) Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39(1): 201225.
 Peregrine DH (1967) Long waves on a beach. J Fluid Mech 27(4): 815827.
 Brocchini M (2013) A reasoned overview on Boussinesqtype models: the interplay between physics, mathematics and numerics. Proc Math Phys Eng Sci 469(2160): 20130496.
 Yamazaki Y, Kowalik Z, Cheung KF (2008) Depthintegrated, nonhydrostatic model for wave breaking and runup. Int J Num Meth Fluids 61(5): 473497.
 Yamazaki Y, Cheung KF, Kowalik Z (2011) Depth‐integrated, non‐hydrostatic model with grid nesting for tsunami generation, propagation, and run‐up. Int J Num Meth Fluids 67(12): 20812107.
 Zelt JA (1991) The runup of nonbreaking and breaking solitary waves. Coastal Eng 15(3): 205246.
 Walters RA (2005) A semiimplicit finite element model for nonhydrostatic (dispersive) surface waves. Int J Num Meth Fluids 49(7): 721737.
 Wei G, Kirby JT, Grilli ST, Subramanya R (1995) A fully nonlinear Boussinesq model for surface waves. Part 1: Highly nonlinear unsteady waves. J Fluid Mech 294: 7192.
 Wei Z, Jia Y (2014) Simulation of nearshore wave processes by a depthintegrated nonhydrostatic finite element model. Coast Eng 83: 93107.
 Kim DH, Lynett PL (2011) Dispersive and Nonhydrostatic Pressure Effects at the Front of Surge. ASCE ASCE J Hydraul Eng 137(7): 75465.
 SoaresFrazao S, Zech Y (2002) Undular bores and secondary waves experiments and hybrid finitevolume modelling. J Hydraul Res 40(1): 3343.
 Casulli V, Stelling GS (1998) Numerical simulation of 3D quasihydrostatic, freesurface flows. J Hydraul Eng 124(7): 678686.
 Cui H, Pietrzak JD, Stelling GS (2012) Improved efficiency of a nonhydrostatic, unstructured grid, finite volume model. Ocean Modelling 54(55): 5567.
 Stelling GS, Zijlema M (2003) An accurate and efficient finitedifference algorithm for nonhydrostatic freesurface flow with application to wave propagation. Int J Num Meth Fluids 43(1): 123.
 Witham GB (1974) Linear and nonlinear waves. John Wiley and Sons, UK.
 LeVeque RJ (1992) Numerical Methods for Conservation Laws. In: Birkhauser (Ed.), Lectures Notes in Mathematics (2_{nd} edn), ETH, Zurich, Springer, Switzerland, pp. 220.
 Toro EF (2009) Riemann solvers and numerical methods for fluid dynamics: A practical introduction. (3^{rd} edn), Springer, Dordrecht, UK, USA, pp. 724.
 Toro EF (2001) ShockCapturing Methods for FreeSurface Shallow Flows. Wiley and Sons Ltd, UK, pp. 326.
 Toro EF, GarciaNavarro P (2007) Godunovtype methods for freesurface shallow flows: A review J Hydr Res 45(6): 736751.
 Fang KZ, Liu ZB, Zou ZL (2014) Modelling coastal water waves using a depth integrated, nonhydrostatic model with shockcapturing ability. J Hydraul Res 53(1): 119113.
 Zijlema M, Stelling GS (2008) Efficient computation of surf zone waves using the nonlinear shallow water equations with nonhydrostatic pressure. Coastal Engineering 55(10): 780790.
 Balzano A (1998) Evaluation of methods for numerical simulation of wetting and drying in shallow water flow models. Coastal Engineering 34(12): 83107.
 Brocchini M, Svendsen IA, Prasad RS, Bellotti G (2002) A comparison of two different types of shoreline boundary conditions. Comput Meth Appl Mech Eng 191(3940): 44754496.
 Gourgue O, Comblen R, Lambrechts J, Kärnä T, Legat V, et al. (2009) A flux limiting wettingdrying method for finiteelement shallowwater models, with application to the Scheldt Estuary. Advances in Water Resources 32(12): 17261739.
 Aricò C, Sinagra M, Begnudelli L, Tucciarelli T (2011) MAST2D diffusive model for flood prediction on domains with triangular Delaunay unstructured meshes. Advances in Water Resources 34(11): 14271449.
 Aricò C, Tucciarelli T (2007) MAST solution of advection problems in irrotational flow fields. Advances in Water Resources 30(3): 665685.
 Aricò C, Tucciarelli T (2007) A Marching in space and time (MAST) solver of the shallow water equations. Part I: The 1D case. Advances in Water Resources 30(5): 12361252.
 Aricò C, Nasello C, Tucciarelli T (2007) A Marching in space and time (MAST) solver of the shallow water equations. Part II: The 2D case. Advances in Water Resources 30(5):12531271.
 Aricò C, Sinagra M, Tucciarelli T (2013) MAST solution of “anisotropic” flow problems: application to the shallow water equations. Advances in Water Resources 62(A): 1336.
 Dean RG, Dalrymple RA (1991) Water Wave Mechanics for Engineers and Scientists. (2_{nd} edn), Advanced Series on Ocean Engineering: World Scientific Publishing Co Pte Ltd, USA, pp. 353.
 Wu W (2007) Computational River Dynamics. Taylor & Francis Group, London, UK, pp. 494.
 Smit PB, Stelling GS, Roelvink D, J van Thiel de Vries, R McCall et al. (2010) X Beach: Nonhydrostatic model. Validation, verification and model description. Delft University of Technology and Deltares, Netherlands, p. 169.
 Wei Z, Jia Y (2013) A depthintegrated nonhydrostatic finite element model for wave propagation. Int J Num Meth Fluids 73(11): 9761000.
 Lu X, Dong B, Mao B, Zhang X (2015) A two dimensional depthintegrated nonhydrostatic numerical model for near shore wave propagation. Ocean Modelling 96(2): 187202.
 Zijlema M, Stelling GS, Smit P (2011) SWASH: an operational public domain code for simulating wave fields and rapidly varied flows in coastal waters. Coastal Engineering 58(10): 9921012 .
 Putti M, Cordes C (1988) Finite element approximation of the diffusion operator on tetrahedral. SIAM J Sci Comput 19(4): 11541168.
 Li X, Huang W (2010) An anisotropic mesh adaptation method for the finite element solution of heterogeneous anisotropic diffusion problems. J Comp Phys 229(21): 80728094.
 Nag Library Manual (2005) The Numerical Algorithms Group Ltd, Oxford, UK.
 Saad Y (2000) Iterative methods for sparse linear systems. (2_{nd} edn ) Society for Industrial and Applied Mathematics.
 Bermudez A, Vazquez ME (1994) Upwind methods for hyperbolic conservation laws with source terms, Comput. Fluids 23(8): 10491071.
 Lynett PJ, Wu TR, Liu PLF (2002) Modeling wave run up with depthintegrated equations. Coastal Engineering 46(2): 89107.
 NETGEN.
 Sampson J, Easton A, Singh M (2003) Moving boundary shallow water flow in circular paraboloidal basins. Proc 6^{th} Eng Math and Appl Conf, Sydney Australia, pp. 223227.
 Hou J, Liang Q, Simons F, Hinkelmann R (2013) A 2D wellbalanced shallow flow model for unstructured grids with novel slope source term treatment. Advances in Water Resources 52: 107131.
 Synolakis CE (1987) The run up of solitary waves. J Fluid Mech 185: 523545.
 Roeber V, Cheung KF, Kobayashi MH (2010) Shockcapturing Boussinesqtype model for nearshore wave processes. Coast Eng 57(4): 407423.
 Dimakopoulos AS, Guercio A, Cuomo G (2014) Advanced numerical modeling of tsunami wave propagation, transformation and runup. Proc. of the Institution of Civil Engineers Engineering and Computational Mechanics 167(3): 139151.
 Briggs MJ, Synolakis CE, Harkins GS, Green DR (1995) Laboratory experiments of tsunami run up on a circular island. Pure and applied geophysics 144(34): 569593.
 Fuhrman DR, Madsen PA (2008) Simulation of nonlinear wave runup with a high order Boussinesq model. Coastal Engineering 55(2): 139154.
 Kazolea M, Delis AI, Nikolos IK, Synolakis EC (2012) An unstructured finite volume numerical scheme for extended 2D Boussinesqtype equations. Coastal Engineering 69: 4266.
 Matsuyama M, Tanaka H (2001) An experimental study of the highest runup height in the 1993 Hokkaido NanseiOki earthquake tsunami. Proceedings of the International Tsunami Symposium, Australia, pp. 879889.
 Tsunami Runup onto a Complex Threedimensional Beach; Monai Valley. Noaa Center for Tsunami Research.
 Nikolos IK, Delis AI (2009) An unstructured nodecentered finite volume scheme for shallow water flows with wet/dry fronts over complex topography. Comput Methods Appl Mech Engg 198(4748): 37233750.
 Kesserwani G, Liang Q (2012) Dynamically adaptive grid based discontinuous Galerkin shallow water model. Advances in Water Resources 37: 2339.
 Paglieri L, Ambrosi L, Formaggia L, Quarteroni A, Scheinine AL (1997) Parallel computation for shallow water flow: a domain decomposition approach. Parallel Computing 23(9): 12611277.
Useful Links


For Authors

For Editors

For Reviewers

Downloads

MedCrave Reprints
MedCrave Group is ardent to provide article reprints at an instant affordable
Read more...

