**About the Editors**

**Anargiros I. Delis** is a Mathematician, Associate Professor of "Computational Mathematics" at the School of Production Engineering & Management of the Technical University of Crete, Greece, where he has been working since 2003. He is a graduate of the Department of Mathematics of the University of Crete (1993); a holder of an M.Sc. in Numerical Analysis & Computing from the University of Manchester Institute of Science & Technology (UMIST), U.K.; and received his PhD from the University of the West of England, Bristol, U.K., in Applied and Computational Mathematics in 1998. His research focuses on the areas of computational hydrodynamics of free surface flows, computational fluid dynamics, and more generally in the area of scientific computing.

**Ioannis K. Nikolos** is a full professor with the School of Production Engineering & Management, Technical University of Crete, Greece, and Director of the Turbomachines & Fluid Dynamics Laboratory (TurboLab: TUC). Prof. Nikolos has obtained his Diploma (1990) from the School of Mechanical Engineering, National Technical University of Athens (N.T.U.A.), and a Ph.D. degree (1996) from the same school in Turbomachines (Lab. of Thermal Turbomachines, N.T.U.A.). He joined Technical University of Crete in 1998. Prof. Nikolos has more than 30 years of experience in R&D projects funded by the EU, industry, and the Greek State. His research work is in the fields of computational fluid dynamics (CFD), computational heat transfer, computational engineering, turbomachines, and engineering design optimization using AI (artificial intelligence) tools.

### *Editorial* **Shallow Water Equations in Hydraulics: Modeling, Numerics and Applications**

**Anargiros I. Delis \* and Ioannis K. Nikolos**

School of Production Engineering & Management, Technical University of Crete, University Campus, 73100 Chania, Greece; jnikolo@dpem.tuc.gr

**\*** Correspondence: adelis@science.tuc.gr

**Abstract:** This Special Issue aimed to provide a forum for the latest advances in hydraulic modeling based on the use of non-linear shallow water equations (NSWEs) and closely related models, as well for their novel applications in practical engineering. NSWEs play a critical role in the modeling and simulation of free surface flows in rivers and coastal areas and can predict tides, storm surge levels and coastline changes from hurricanes and ocean currents. NSWEs also arise in atmospheric flows, debris flows, internal flows and certain hydraulic structures such as open channels and reservoirs. Due to the important scientific value of NSWEs, research on effective and accurate numerical methods for their solutions has attracted great attention in the past two decades. Therefore, in this Special issue, original contributions in the following areas, though not exclusively, have been considered: new conceptual models and applications; flood inundation and routing; open channel flows; irrigation and drainage modeling; numerical simulation in hydraulics; novel numerical methods for shallow water equations and extended models; case studies; and high-performance computing.

**Keywords:** shallow water equations; free surface flows; modeling; hydraulic engineering; computational methods; simulation

#### **1. Introduction**

In hydraulic engineering, the modeling and simulation of free surface water flows plays an important role in many real-life practical applications. An important feature of free surface water flows is that they are unbound in space, the limits of the spatial domain being an unknown of the problem to be solved. Problems in which the limits of the fluid are unknown and unsteady include, among others, dam break and flooding flows, tidal flows in coastal water regions, nearshore wave propagation with complex bathymetry structure, tsunami wave propagation and ocean modeling. The full flow field in most such applications can be described by Navier–Stokes equations. However, qualitative and/or quantitative approximations of the actual solution are given by approaches based on simplified equations. This is done in a systematic effort usually to overcome the need for excessively demanding numerical techniques to resolve Navier–Stokes equations (possibly supplied with appropriate turbulence closure models). A widely used approach is that of the 2D depth-averaged models [1,2]. The 2D character of the free surface flow is usually enforced by a horizontal length scale which is much larger than the vertical one, and by a velocity field which is quasi-homogeneous over the water depth. This small ratio between the vertical and horizontal length scales characterizes many engineering applications, mainly in river and coastal engineering. As such, many geophysical flows can be modeled by the well-known non-linear shallow water equations (NSWEs), also called Saint-Venant (SVEs) equations, and closely related models. From a mathematical point of view, these models constitute a system of partial differential equations, mainly of the hyperbolic type, also referred to as hyperbolic systems of balance laws.

Despite the relative simplicity of NSWEs, the shallow water assumption is valid in many applications in hydraulics and as such has a long tradition of providing a scientific

**Citation:** Delis, A.I.; Nikolos, I.K. Shallow Water Equations in Hydraulics: Modeling, Numerics and Applications. *Water* **2021**, *13*, 3598. https://doi.org/10.3390/w13243598

Academic Editor: Helena M. Ramos

Received: 15 November 2021 Accepted: 14 December 2021 Published: 15 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

basis for engineering practice. To this end, shallow water equations arise in modeling water flows in rivers, canals, lakes, reservoirs, coastal and urban areas and many other situations in which the water depth is much smaller than the horizontal length scale of motion. Hence, shallow water and closely related models are widely used in oceanography and atmospheric sciences to model, among others, flood waves and hazardous phenomena such as hurricanes/typhoons and tsunamis. In recent years, new conceptual models and applications based on NSWEs have emerged. These include, among others, models for flood inundation and routing, sediment transport and morphodynamic models, avalanches, pollutant transport in water models and dispersive wave propagation ones. Furthermore, in recent decades, many efforts have been devoted to develop one- and two-dimensional numerical models for unsteady free surface flows, with the most widely used mathematical framework to be the one based on the 2D NSWEs [3]. The need for accurate and efficient numerical methods has resulted in an ongoing quest for such methods that are of higher-order spatial and temporal accuracy that can effectively predict the most relevant physical phenomena.

In this Special Issue, we attempted to include manuscripts related to recent advances in terms of modeling, numerical methods and applications closely related to the implementation of the NSWEs.

#### **2. Summary of This Special Issue**

In [4], the authors presented an extensive numerical investigation for studying the accuracy and efficiency of three common shallow water finite volume solvers, namely the HLLC, Roe, and Central-Upwind (CU) schemes. Four cases dealing with shock waves and the wet–dry phenomenon were selected. All schemes were provided in an in-house code NUFSAW2D, the model of which was of second-order accuracy in space wherever the regimes were smooth and robust when dealing with strong shock waves—and of fourth-order accuracy in time. To give a fair comparison, all source terms of the 2D NSWEs were treated similarly for all schemes, namely the bed-slope terms were separately computed from the convective fluxes using a Riemann-solver-free scheme—and the friction terms were semi-implicitly computed within the framework of the RKFO method. Two important findings were shown by the presented simulations. Firstly, highly efficient vectorization could be applied to the three solvers on all hardware used. This was achieved by guided vectorization, where a cell-edge reordering strategy was employed to ease the vectorization implementations and support the aligned memory access patterns. Secondly, it was shown that for the four cases simulated, strong agreements by all schemes were obtained between the numerical results and observed data, where there were no significant differences. However, in terms of efficiency, the CU scheme was able to outperform the HLLC and Roe schemes with average factors of 1.4× and 1.25×, respectively. Finally, it was demonstrated that the edge-driven level, especially the reconstruction technique and solver computations, were the most time-consuming parts, which required 65–75% of the entire simulation time. This shows that some more "aggressive" optimization techniques will still constitute a hot topic for future studies to make shallow water simulations more efficient, particularly in the edge-driven level. Since simulating shallow water flows—especially complex phenomena that require performing long real-time computations as part of disaster planning such as dam-break or tsunami cases—on modern hardware nowadays, and even in the future, is becoming increasingly common, focusing simulations only on numerical accuracy but ignoring the performance efficiency may not be an option anymore. Wasting the performance is obviously undesirable the excessive time required by such long realtime simulations. Modern hardware offers many features for gaining efficiency, one of which is vectorization, which can be regarded as the "easiest" way for benefiting from the vector-level parallelism, and is thus non-trivial.

The authors in [5] aimed to investigate the dynamic process of the transformation between different channel patterns with different control variables. They proposed an extended 2D depth-averaged numerical model which incorporates the hydrodynamic, sediment transport and a river morphological adjustment sub-model. The model was solved in an orthogonal curvilinear grid system by using the Beam and Warming alternatingdirection implicit (ADI) scheme. The sediment transport sub-model includes the influence of non-uniform sediment with bed surface armoring and a correction for the direction of bed-load transport due to secondary flow and transverse bed slope. The bank erosion sub-model incorporates a simple simulation method for updating bank geometry during either degradational or aggradational bed evolution. Then, the extended 2D model was applied to duplicate the evolution of the channel pattern with variations in flow discharge, sediment supply and bank vegetation. Complex interactions among the flow discharge, sediment supply and bank vegetation leads to a transition from the braided pattern to the meandering one. The analysis of the simulation process indicated that (1) a decrease in the flow discharge and sediment supply can lead to the transition; and (2) the riparian vegetation helps stabilize the cut bank and bar surface but it is not key to the transition. The results are in agreement with the criterion proposed in previous research, confirming the 2D numerical model's potential in predicting the transition between different channel patterns and improving our understanding of the fluvial process. It was concluded that, further studies are needed to research the fundamental equation that governs the evolution of alluvial river, which has not been fully understood to ensure the availability of the numerical model.

The research presented in [6] aimed to present improvements in the computational efficiency of a 2D hydrodynamic model that can lead to a significant bearing on its engineering applications. In the presented study, an improved 2D shallow water dynamic model applying a local time stepping (LTS) algorithm was established using a finite volume method (FVM) on a triangular mesh. Results from the simulation and analysis of three canonical test cases showed that the LTS scheme has a satisfactory ability for adapting real complex applications and long simulations while meeting the required accuracy. The following conclusions were drawn: (1) based on the FVM for unstructured grids, the implemented LTS algorithm improved the computational efficiency of the model, while satisfying water conservation conditions. In the anti-symmetric dam break case, a speedup ratio of 2.1 was achieved, which saved 53% in execution time. The speedup ratio of the non-flat bottom dam break case was 1.3, which represented a shortening of 26% in the calculation time. The numerical simulation of the navigable flow of the river reach between the Three Gorges and Gezhouba Dams achieved a speedup ratio of 1.9, which represented a saving of 49% in modeling time; (2) the proportions of coarse to refined meshes on the acceleration effect of the LTS algorithm were noticeable. It was evident that a higher speedup ratio was obtained when the proportion of the refined mesh was minimized. When the proportion of the refined mesh was high, the acceleration effect was not significant. It is clear that the LTS algorithm is best suited to situations in which refinement is only required in small regions; and (3) when using the LTS algorithm on non-uniform unstructured grids, the larger the grid scale difference, the more obvious the grid layering became. This led to increased acceleration effects. However, computational accuracy was slightly impaired by excessive differences in grid mesh size. These results can provide technical guidelines for reducing computational time for 2D hydrodynamic models on non-uniform unstructured grids.

In [7], the main goal was to use various numerical investigations to figure out the response of coastal tides on the west coast of Korea (WCK) to the open-boundary conditions and bottom roughness using an open source computational fluids dynamics tool—the TELEMAC model. Three well-known assimilated tidal models were used to obtain a detailed tidal forcing at open boundaries—the finite element solution (FES2014); the Oregon State University TOPEX/Poseidon Global Inverse Solution Tidal Model (TPXO9.1); and the National Astronomical Observatory of Japan (NAO.99Jb). It was shown that there were no significant differences between the responses in tidal amplitudes in the WCK induced by three open-boundary conditions obtained from three assimilated tidal models. In addition, the numerical simulation of the tidal flow in the WCK should not use a uniform bottom roughness coefficient. Due to the complicated bathymetry, indented coastlines

and bed variability of the WCK, it caused strong local effects on the tides in this region. Therefore, a non-uniform bottom roughness should be applied to the modeling whereby the smallest value can be applied for Incheon, a larger value for Mokpo, and the largest value for Gunsan. The largest value of the bottom roughness coefficient was applied to the Gunsan region because its bed form was more variable than in other regions. The numerical results showed that the accuracy of the modeling of the tidal elevation around the WCK was strongly dependent on the bottom roughness rather than the offshore tidal boundary conditions. Moreover, the numerical results can provide not only a better fit to the observations but also higher spatial resolutions in comparison to the results obtained from assimilated tidal models around the WCK. However, it should be noted that the numerical results obtained from this study were still limited due to the coarse resolution (30 arcs/s) of the bathymetry obtained from GEBCO2014—which was not sufficient to capture the real geometries whose sizes were less than such resolution. Therefore, a further study with a higher resolution is necessary in order to obtain a more precise prediction of the tidal current and its elevation around the WCK. Furthermore, the wind forcing on the sea surface and the tidal energy dissipation should be taken into account.

In [8], the one-dimensional Boussinesq equations were numerically solved using the MacCormack as well as the Dissipative two–four finite difference schemes for the simulation of hydraulic jump formation in a horizontal rectangular open channel and for upstream Froude numbers Fr in the range of 2.44–5.38. The governing equations were enriched with additional terms if compared to the Saint-Venant equations, to account for the non-hydrostatic pressure distribution in the regime of rapidly varied flow. Terms related to the energy loss and the gravity forces were also included. The initial condition was a steady supercritical gradually varied flow along the whole length of the modeled channel. The upstream and downstream boundary conditions regarding the flow depth remained constant during the iteration process and equal to the values measured in experiments. The method of specified intervals was used for the calculation of the velocity at the downstream end, assuming that the positive characteristic through a point does not intersect with already established grid points. The variable time step was used in every iteration according to the CFL (Courant–Friedrichs–Lewy) stability criterion, along with artificial viscosity for the smoothing of the oscillations occurring in the jump. The computational results compared well with the experiments since the specific force was computed from the depth and mean velocity at both ends of the hydraulic jump with acceptable tolerance, and the mass conservation equation was verified for all numerical schemes and all test cases. From such a model, one can determine the sequent depth ratio as well as the length of the jump, results that are useful in the design of stilling basins (geometrical properties). Given a stilling basin with a known inflow Froude number and flow depth, the engineer must decide the end sill dimensions and the basin length, so that the hydraulic jump is contained in the stilling basin. Finally, from the comparison of the numerical results and experiments, it can be concluded that the aforementioned numerical modeling schemes can predict the basic features of the classical hydraulic jump with acceptable accuracy.

In [9], the authors introduced several innovations, including the use of the non-linear Krylov accelerator in open-channel flows, an evolutionary domain algorithm and the use of CasADi to solve steady 1D flows using the Saint-Venant model equations. These improvements led to an algorithm that is able to quickly solve steady open-channel flows. Therefore, optimization problems and uncertainty analyses that require many evaluations become more tractable. An original algorithm was implemented in order to significantly improve the computation time of a steady 1D open-channel flow problem. This includes two main optimizing strategies: a non-linear Krylov accelerator and an evolutionary domain algorithm. This new algorithm was validated against the academic benchmarks of flows over a bump. The results showed a good agreement between the numerical and analytical values. The performance of the suggested algorithm was evaluated against the non-linear optimization software CasADi. It showed good scalability properties. Indeed, the execution time of the proposed algorithm linearly evolves with the number of nodes. This is not the case with other techniques when the mesh is refined and/or when the number of nodes is increased. Finally, the capabilities of the proposed algorithm were validated on a real-world case. The optimized algorithm was used in order to quickly compute the initial condition required by the operational model for the Romanche River in France. The technique was able to provide a steady-state solution to the unsteady model in a very short period of time.

In [10], the performances of several modeling approaches were compared in order to evaluate their results and computational requirements in a transient river flow event in a reach of the Ebro river (Spain) that includes a reservoir covering a large area. A 2D distributed shallow water model solved over a triangular grid and a 1D shallow water model was used to discretize the full domain. Additionally, an aggregated volume balance model was implemented to model a reservoir region in order to allow computational saving. This led to a coupled 1D–0D model. Finally, a proportional–integral–derivative (PID) control algorithm was implemented as a regulation technique at the dam location and combined with both the 1D model and the 1D–0D model. From the comparison of the performances of the 2D and 1D models, it was concluded that the results of the 1D model for the recent flooding events at the considered Ebro River reach were very similar to those provided by the 2D model. The water level and discharge data predicted by both models follow the same trend. The cross-sections used to build the 1D model computational mesh were carefully located to reproduce the river curvature in detail, which is important to obtain a realistic evolution of the hydraulic variables. This effort is justified by the immense computational saving that the use of the 1D model offers, as long as there is no interest in representing the floodplain flow that the 1D model does not take into account. The coupling of the 1D model for the river flow at the upstream reach and the 0D model for the reservoir (1D–0D model) offers results very similar to those from the full 1D model. There was some lag due to the instantaneous propagation of the hydrograph in the reservoir assumed by the 0D model but this is acceptable considering the computational savings that the use of this model implies compared to the full 1D model. The computational times observed with the 1D–0D model justifies the use of this combined approach. Therefore, the coupling of a 0D model for the reservoir with the 2D model for the upstream river reach was envisaged as future work since this will lead to high computational savings, something very positive for simulations with 2D models as well as the possibility to simulate the floodplain flow behavior. The PID control algorithm was implemented with the objective to ensure a fixed-surface water level at the dam. The results showed that this target level value was never reached despite the time variable discharge, which means that the implementation of the control algorithm is a correct security measure to avoid exceeding certain levels in the reservoir. It will be convenient in the future to implement an algorithm that takes into account more realistic and complex objectives.

**Author Contributions:** A.I.D. and I.K.N. conceived, designed, and wrote the editorial. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** Thanks to all of the contributions to the Special Issue, the time invested by each author, as well as to the anonymous reviewers and editorial managers who have contributed to the development of the articles in this Special Issue.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

