The work herein presented has been developed at the Computational Hydraulics Group at Univerisity of Zaragoza. Computational resources have been provided by LIFTECCSIC. The validation of the computational tools has been done in colaboration with the LCHEPFL.
This work is the result of my Doctoral Thesis Accurate simulation of shallow flows using arbitrary order ADER schemes and overcoming numerical shockwave anomalies, supervised by Dr. J. Murillo and presented at University of Zaragoza on the 27th April 2018. In the following sections, the most relevant aspects of this thesis are brought out and supplementary multimedia content is presented.
Generation of arbitrary order augmented schemes for hyperbolic problems with source terms: application to the SWE
We aim at the generation of fullydiscrete arbitrary order numerical schemes, based on the WENO spatial reconstruction and ADER timestepping technique, with application to hyperbolic problems with source terms. The proposed schemes are based on augmented Riemann solvers that include the contribution of source terms in the definition of the Riemann problem, allowing the presevation of equilibrium states with machine precision.
We consider the application of the aforementioned methods to the SWE with bed elevation, friction and Coriolis. In presence of such sources, the relevant equilibrium states are (see Figure 3):
 Quiescent equilibrium (wellbalanced property).
 Moving water equilibrium (energybalanced property).
 Geostrophic equilibrium (wellbalanced property in the rotating frame).
Furthermore, additional dissipation effects are accounted for by including turbulence models using the Boussinesq approximation. The RANS and URANS approaches are considered and algebraic as well as 1 and 2 equation turbulence models are used.
The WENO AR/ARLADER method in 1D: wellbalanced and energybalanced simulation of the SWE with bed variation
The first stage of the project was the development of the mathematical framework for the resolution of hyperbolic conservation laws with geometric source terms with arbitrary order of accuracy using WENOADER schemes. A new family of solvers for the Derivative Riemann problem are proposed, using the augmentedsolver methodology and based on previous work from Toro, Castro and Titarev. The novelty of the proposed methods is:

For transient cases, they converge with arbitrary order to the analytical solution. The numerical solution to a complete set of RPs, including solutions in the resonant regime, are presented in NavasMontilla, 2015 and NavasMontilla, 2016 (see Figure 4).

For steady cases, they provide the exact solution with independence of the grid thanks to the energybalanced property.
Details of the methods and more results can be found in NavasMontilla, 2015 and NavasMontilla, 2016.
The WENO ARLADER method in 2D: wellbalanced simulation of the SWE with bed variation, friction and Coriolis
The numerical schemes in the previous section were then extended to the resolution of the 2D SWE with geometric source term and their application to other shallow water models involving nongeometric sources was explored. The following issues are highlighted:
 The proposed scheme offers an arbitrary order resolution of the 2D SWE with bed elevation and friction ensuring the wellbalanced property. The application of the scheme to the resolution of a complex flow pattern around a solid obstacle is shown in Figure 5.
 The proposed scheme allows to compute the advection of solutes with arbitrary order of accuracy, a shown in Figure 6. The scheme can be coupled with diffusionreaction terms for the solutes, with application to industrial processes in water treatment plants.

The scheme is applied to the resolution of the SWE in the rotating frame (including Coriolis). In NavasMontilla, 2018, the results for the benchmark of the propagation of equatorial Rossby solitons are presented and the relative errors for the amplitude and celerity of the soliton are kept below 2% (Watch video). The performance of the model is also assessed for the resolution of an anticyclonic eddy in the beta plane (Watch video).

The proposed methods offer a remarkable gain in computational efficiency when applied to 2D shallow water scenarios with source terms. Figure 7 shows the numerical error vs. CPU time (singlethreaded/serial execution) and wall time (parallel execution in 24 threads implemented using the OpenMP paradigm) for the resolution of an smooth nonequilibrium Gaussian water surface over a smooth bottom topography (NavasMontilla, 2018). The numerical results evidence that the 3rd order scheme is able to provide the same level of accuracy than a 1st order scheme requiring a 65 times shorter computational time, for an error of around 1.e4. It is worth noting that this gain is increased when seeking lower errors. Furthermore, the plots also show that the paralellization of the code using OpenMP allows a significant speedup.
 As the order is increased, the numerical diffusion is strongly reduced and the proposed scheme is able to reproduce a larger extent of the kinetic energy cascade for 2D turbulence. Figures 8 and 9 show the numerical resolution of double shear layer, which is initially perturbed in the transverse direction. Due to the inherent instability (KelvinHelmholtz instability) of the flow across the shear layer, any initial perturbation is amplified and evolves into a hydrodynamic twodimensional turbulence. The vortices interact to each other and are combined to form bigger vortices. The mass exchange across the layers is governed by such turbulent pattern and only very high order numerical schemes, that produce low dissipative solutions, can provide trustworthy predictions of such process in an affordable time. In NavasMontilla, 2018, the performance of the 3rd order ARLADER and HLLSADER schemes for the resolution of diffusive and turbulent shear layers is assessed.
Watch the following videos for more examples:
Details of the methods and more results can be found in NavasMontilla, 2018
URANS simulation of shallow flows using the WENO ARLADER method for the SWE
The proposed scheme offers a very low numerical dissipation, allowing the computation of small turbulent structures and the reproduction of the theoretical energy cascade as the grid is refined. However, the mathematical model considers a 2D depth averaged flow, while turbulence is, in essence, three dimensional. Therefore, extra dissipation terms must be included to account for the small scale turbulent dissipation that the model cannot resolve.
In shallow flows, there is a coexistence of smallscale 3D turbulence, mainly generated by the friction on the bottom, and largescale 2D turbulence, generated by horizontal gradients. The proposed model uses a turbulence model to account for the effects of the smallscale turbulence and resolves the largescale 2D vortices. This approach is often called depth averaged URANS simulation or depth averaged LES simulation. A representation of the typical energy cascade in shallow flows is presented in Figure 10.
Turbulence modelling is of practical application when simulating the flow interaction with solid obstacles such as bridge piers and groynes, or even with strong variations in the bed elevation and bed friction.
Simulation of the flow over a submerged conical island
This test case is a benchmark for tsunami simulation models and was introduced with this purpose for the first time at the National Tsunami Hazard Mitigation Program (NTHMP) workshop (Portland, 2015). It consists of a shallow flow over a smallslope submerged island. The aim of this benchmark is to test the model’s ability to generate a separation region and the resulting oscillatory wake for an idealized and simplified case.
Figure 11 shows a snapshot of simulated tracer concentration distribution where the vonKarman vortex street can be observed. Only when including a suitable calibration of the friction coefficient and the turbulence model, the velocities in the wake are properly predicted and there is a stable periodicity in the shedding of vonKarman vortices. Figure 12 displays a comparison between numerical and experimental velocities in the wake region. The detailed description of the simulation of this experiment is described in NavasMontilla, 2019.
Application of the model to a practical problem of environmental relevance: river restoration
In the last decades, riverine and coastal habitats have degenerated because of anthropogenic activities. Nowadays, the scientific community is making a big effort to design novel approaches to recover biodiversity in such ecosystems. The utilization of fast and reliable predictive tools will suppose a breakthrough in this field as they will be able to plan efficient strategies based on the predicted quantitative variables.
Channels with lateral cavities are commonly used for river restoration purposes as the presence of cavities enhances fine sediment trapping. The flow in these channels, far from being simple, involves the presence of steady seiching waves produced by the coupling between the instability of the separated turbulent layer along the opening of the cavities and a gravity standing wave within the cavities. Such coupling is associated with largescale coherent vortical structures in the unstable shear layer and periodic oscillations of the free surface within the cavity. The complexity of such flow configuration challenges the prediction capability of simulation models.
Detailed results of the application of the proposed schemes to the simulation of channels with lateral cavities can be found in the following article NavasMontilla, 2019, done in colaboration with C. Juez (UPM), M. J. Franca (UNIHE Delft) and J. Murillo (Unizar). An animation of the results can be watched here.
Watch the following videos for more examples: