Shallow Water and Sediment Transport with Kelvin–Voigt Seabed: Numerical Insights from Theoretical Case Studies

Forecasting water waves and currents in near shore and off shore regions of the seas and oceans is essential to maintain and protect our environment and man made structures. In wave hydrodynamics, waves can be classified as shallow and deep water waves based on its water depth. The mathematical models of these waves are shallow water and free surface gravity water wave equations which describe the hydrodynamics of waves and currents near shore and off shore regions of seas and oceans. The complexity in these models exist as moving boundaries whose position depends on the solution of the governing equations. For shallow water waves, it is the shore line boundary where the water depth falls dry and for deep water waves, it is the free surface which separates the sea or ocean from atmospheric air. It is often difficult to solve these wave equations analytically while solving them numerically in an efficient and accurate way is a challenging task because of the moving boundaries. The numerical challenges are two fold: one is to develop a numerical method which is accurate and efficient for deforming grids and the other is to design a numerical algorithm for the grid adaptation following the moving boundaries. In this thesis, we aimed at first developing space-time discontinuous Galerkin finite element schemes for shallow water and free surface gravity water wave equations which are accurate and efficient for deforming grids. The shallow water equations are a leading order hydrodynamic model for coastal waves and currents. This is because they can exhibit the complicated flooding and drying phenomena due to the moving shore line boundary, and the wave breaking phenomena in the form of bores. A new space-time discontinuous Galerkin (DG) discretization is first presented for the (rotating) shallow water equations over varying topography and fixed boundaries. We formulated the discretization in an efficient and conservative way with the numerical HLLC flux on the finite element boundaries. We also designed a novel way to apply numerical dissipation around discontinuities, that are present in the form of bores, with the help of Krivodonova's discontinuity indicator such that the spurious oscillations are suppressed. The non-linear algebraic system resulting from the space-time discretization is solved using a pseudo-time integration method. A thorough verification of the space-time DG finite element method is undertaken by comparing the numerical and exact solutions. We carried out a discrete Fourier analysis of the one dimensional linear rotating shallow water equations to show that the method is unconditionally stable with minimal dispersion and dissipation error. The numerical scheme is verified and validated for a number of problems arising in geophysical flows. To demonstrate that the space-time DG method is particularly suitable for problems with dynamic grid motion, we simulated nonlinear waves generated by a wave maker and verified these for low amplitude waves where linear theory is approximately valid. Free surface gravity water wave equations is widely used in marine and offshore engineering to model waves. The mathematical nature of these equations is complex because it consists of a potential flow equation which is of elliptic nature and nonlinear free surface boundary conditions which are hyperbolic in nature. Hence, a space-time discontinuous Galerkin finite element method is presented for simplified linear free surface gravity water waves. The free surface gravity water wave equations also arise from Luke's variational formulation which is associated with the conservation of energy and phase space, under suitable boundary conditions. This variational formulation also provided a basis to obtain a novel space-time variational (dis)continuous Galerkin finite element method. Both the space-time discontinuous Galerkin and the space-time variational finite element discretizations result in an algebraic linear system of equations with a very compact stencil, i.e., the algebraic equations from each element is coupled to its immediate neighboring elements only. Thus, the linear system of equations are built using an efficient block sparse matrix storage routine and solved by using iterative linear solvers using a well-tested PETSc package. Numerical schemes are verified for harmonic waves in a periodic domain and generated in a wave basin. Extension of the space-time discontinuous Galerkin method for flooding and drying in shallow water waves and nonlinear free surface evolution of deep water waves will be the topic of future research.

Global vorticity budgets in C-grid shallow water (SW) and quasi-geostrophic (QG) models of winddriven ocean circulation with free-slip boundary conditions are considered. For both models, it is pointed out that the discretized vorticity equation is defined only over a subdomain that excludes boundary grid nodes. At finite resolution, this implies an advective flux of vorticity across the perimeter of the discretized vorticity domain. For rectangular basins where grid axes are aligned with the basin walls, this flux tends to zero as resolution is increased. We also consider the case in which the grid is rotated with respect to the basin, so that a step-like coastline results. Increased resolution then leads to more steps and, because the advective flux of vorticity out of the domain is particularly large at steps, it is no longer obvious that increased resolution should reduce the advective flux. Results are found to be sensitive to numerical details. In particular, we consider different formulations for the non-linear terms (for both the SW and QG models) and two formulations of the viscous stress tensor for the SW model [the conventional five-point Laplacian and the δ—ζ stress tensor suggested by Madec et al. (J. Phys. Oceanogr. 21, 1349—1371)]. For the SW model, the overall circulation and the behavior of the flux term are dependent on both the formulation of the viscous stress tensor and the non-linear terms. The best combination is found to be the δ—ζ tensor with an enstrophy-preserving advection scheme. With this combination, the circulation of the non-rotated basin is recovered in rotated basins and the advective flux tends to converge towards zero with increasing resolution. The poorest combination is the δ—ζ tensor with the conventional advective scheme. In this case, the advective flux term diverges with increasing resolution for some rotation angles and the model crashes for some others. For the QG model, the convergence order of the advective flux term of absolute vorticity is near unity (roughly the same as with the SW model). Most of the error (especially at high resolution) is related to errors in the β term (which is hidden in the advective contribution in the SW model). However, the overall circulation is less sensitive to the rotation of the grid with respect to the basin, especially when the Jacobian proposed by Arakawa (J. Comput. Phys. 1, 119—143) is used.

Exact solutions of the nonlinear shallow water wave equations for forced flow involving linear bottom friction in a region with quadratic bathymetry have been found. These solutions also involve moving shorelines. The motion decays over time. In the solution of the three simultaneous nonlinear partial differential shallow water wave equations it is assumed that the velocity is a function of time only and along one axis. This assumption reduces the three simultaneous nonlinear partial differential equations to two simultaneous linear ordinary differential equations. The analytical model has been tested against a numerical solution with good agreement between the numerical and analytical solutions. The analytical model is useful for testing the accuracy of a moving boundary shallow water numerical model. References Balzano, A., Evaluation of methods for numerical simulation of wetting and drying in shallow water flow models, Coastal Engineering, 34, 1998, 83--107. Carrier, G. F. and Greenspan, H. P., Water waves of finite amplitudes on a sloping beach, Journal of Fluid Mechanics, 4, 1958, 97--109. Holdahl, R., Holden, H., and Lie,K-A., Unconditionally Stable Splitting Methods For the Shallow Water Equations, BIT, 39, 1998, 451--472. Johns, B., Numerical integration of the shallow water equations over a sloping shelf, International Journal for Numerical Methods in Fluids, 2, 1982, 253--261. Kawahara, M., Hirano, H., and Tsubota, K., Selective lumping finite element method for shallow water flow, {International Journal for Numerical Methods in Fluids}, 2, 1982, 89--112. Lewis, C. H. III and Adams, W. M., Development of a tsunami-flooding model having versatile formation of moving boundary conditions, The Tsunami Society Monograph Series, 1983, No. 1, 128 pp. Parker, B. B., Frictional Effects on the Tidal Dynamics of a Shallow Estuary, PhD thesis, The John Hopkins University, 1984. Peterson P., Hauser J., Thacker W. C., Eppel D., An Error-Minimizing Algorithm for the Non-Linear Shallow-Water Wave Equations with Moving Boundaries. In Numerical Methods for Non-Linear Problems, editors C. Taylor, E. Hinton, D. R. J. Owen and E. Onate, 2, Pineridge Press, 1984, 826--836, http://www.cle.de/hpcc/publications/ Sachdev, P. L., Paliannapan, D. and Sarathy, R., Regular and chaotic flows in paraboloidal basins and eddies, Chaos, Solitons and Fractals, 7, 1996, 383--408 . Sampson, J., Easton, A., and Singh, M., Moving Boundary Shallow Water Flow in Circular Paraboloidal Basins. Proceedings of the Sixth Engineering Mathematics and Applications Conference, 5th International Congress on Industrial and Applied Mathematics, at the University of Technology, Sydney, Australia, editors R. L. May and W. F. Blyth, 2003, 223--227. Sampson, J., Easton, A., and Singh, M., Moving boundary shallow water flow in parabolic bottom topography, Australian and New Zealand Industrial and Applied Mathematics Journal, 47 (EMAC2005), C373--C387, 2006, http://anziamj.austms.org.au/V47EMAC2005/Sampson Sampson, Joe, Easton, Alan and Singh, Manmohan, A New Moving Boundary Shallow Water Wave Numerical Model, Australian and New Zealand Industrial and Applied Mathematics Journal, 48 (CTAC2006), C605--C617, 2007, http://anziamj.austms.org.au/ojs/index.php/ANZIAMJ/article/view/78 Thacker, W. C., Some exact solutions to the nonlinear shallow-water wave equations, J. Fluid. Mech., 107, 1981, 499--508. Vreugdenhil, C. B., Numerical Methods for Shallow-Water Flow, Kluwer Academic Publishers, 1998. Yoon S. B., and Cho J. H., Numerical simulation of Coastal Inundation over Discontinuous Topography, Water Engineering Research, 2(2), 2001, 75--87

Current topographic survey technology provides high-resolution (HR) datasets for urban environments. Incorporating this HR information in models aiming to provide flood risk assessment is desirable because the flood wave propagation is depending on the urban topographic features, i.e. buildings, bridges and street networks. Conceptual, numerical and practical challenges arise from the application of shallow water models to HR urban flood modeling. For instance, numerical challenges are occurrence of wet-dry fronts, geometric discontinuities in the urban environment and discontinuous solutions, i.e. shock waves. These challenges can be overcome by using a Godunov-type scheme. However, the computational cost of this type of schemes is high, such that HR two-dimensional shallow water simulations with practical relevance have to be run on supercomputers. The porous shallow water model is an alternative approach that aims to reduce computational cost by using a coarse resolution and accounting for unresolved processes by means of the porosity terms. Usually, a speedup between two and three orders of magnitude in comparison to HR simulations can be obtained. This study reports preliminary results of a practical test case concerning pluvial flooding in a district of the city of Nice, France, caused by the intense rainfall event on October 3rd, 2015. HR topography data set on a 1 m resolution is available for the district, whereby street features of infra-metric dimensions have been included. A reference solution is calculated by a HR shallow water model on a 1 m by 1 m structured computational grid. The porous shallow water model is run on a 10 m by 10 m grid and the influence of the drag source term is studied. The model results show a large deviation, which is caused by the poor meshing strategy of the porous shallow water (AP) model. The study also summarizes practical challenges that arise during the application of the AP and HR models to a large urban catchment. The main difficulty is to obtain a good mesh. In smaller scale investigations, the mesh is currently constructed by hand such that the cell edges align with buildings. This approach is not feasible for large scale urban catchments with a large number of buildings. Future steps that have to be taken, such as a strategy for automatic mesh generation, are reported on.

Global vorticity budgets in C-grid shallow water (SW) and quasi-geostrophic (QG) models of winddriven ocean circulation with free-slip boundary conditions are considered. For both models, it is pointed out that the discretized vorticity equation is defined only over a subdomain that excludes boundary grid nodes. At finite resolution, this implies an advective flux of vorticity across the perimeter of the discretized vorticity domain. For rectangular basins where grid axes are aligned with the basin walls, this flux tends to zero as resolution is increased. We also consider the case in which the grid is rotated with respect to the basin, so that a step-like coastline results. Increased resolution then leads to more steps and, because the advective flux of vorticity out of the domain is particularly large at steps, it is no longer obvious that increased resolution should reduce the advective flux. Results are found to be sensitive to numerical details. In particular, we consider different formulations for the non-linear terms (for both the SW and QG models) and two formulations of the viscous stress tensor for the SW model [the conventional five-point Laplacian and the δ—ζ stress tensor suggested by Madec et al. (J. Phys. Oceanogr. 21, 1349—1371)]. For the SW model, the overall circulation and the behavior of the flux term are dependent on both the formulation of the viscous stress tensor and the non-linear terms. The best combination is found to be the δ—ζ tensor with an enstrophy-preserving advection scheme. With this combination, the circulation of the non-rotated basin is recovered in rotated basins and the advective flux tends to converge towards zero with increasing resolution. The poorest combination is the δ—ζ tensor with the conventional advective scheme. In this case, the advective flux term diverges with increasing resolution for some rotation angles and the model crashes for some others. For the QG model, the convergence order of the advective flux term of absolute vorticity is near unity (roughly the same as with the SW model). Most of the error (especially at high resolution) is related to errors in the β term (which is hidden in the advective contribution in the SW model). However, the overall circulation is less sensitive to the rotation of the grid with respect to the basin, especially when the Jacobian proposed by Arakawa (J. Comput. Phys. 1, 119—143) is used.

Zero-Inertia (ZI) models are used in overland flow simulation due to their mathematical simplicity, compared to more complex formulations such as Shallow Water (SW) models. The main hypothesis in ZI models is that the flow is driven by water surface and friction gradients, neglecting local accelerations. On the other hand, SW models are a complete dynamical formulation that provide more information at the cost of a higher level of complexity. In realistic problems, the usually huge number of cells required to ensure accurate spatial representation implies a large amount of computing effort and time. This is particularly true in 2D models. Hence, there is an interest in developing efficient numerical methods. In general terms, numerical schemes used to solve time dependent problems can be classified in two groups, attending to the time evaluation of the unknowns: explicit and implicit methods. Explicit schemes offer the possibility to update the solution at every cell from the known values but are restricted by numerical stability reasons. This can lead to very slow simulations in case of using fine meshes. Implicit schemes avoid this restriction at the cost of generating a system of as many equations as computational cells multiplied by the number of variables to solve. In this work, an implicit finite volume numerical scheme has been used to solve the 2D equations in both ZI and SW models. The scheme is formulated so that both quadrilateral and triangular meshes can be used. A conservative linearization is done for the flux terms, leading to a non-structured matrix for unstructured meshes thus requiring iterative methods for solving the system. A comparison between 2D SW and 2D ZI is done in terms of performance, efficiency and mesh requirements, in which both models benefit of an implicit temporal discretization in steady and nearly-steady situations.

Lake circulation and sediment transport in Lake Myvatn have been calculated using AQUASEA, a numerical model developed by Vatnaskil Consulting Engineers. The goal of the modelling was to calculate changes in sediment transport within the lake due to changes in lake bathymetry caused by diatomite mining. The model uses the Galerkin finite element method and consists of a hydrodynamic flow model and a transport-dispersion model. The flow model is based on the shallow water equations and the wave equation. The transport model is based on the conservation of mass for suspended sediment. The model was calibrated against measurements performed during the summer of 1992. These included measurements of water elevation, current velocity, wave height, and concentration of suspended sediment. After calibration, the model was run for different mining scenarios to determine their impact on the sediment transport in the lake.

A two‐dimensional shallow water model based on the anisotropic porosity method is developed for predictions of water flows, sediment transport and bed morphological changes through rigid vegetation. Effects of the vegetation resistance and the spatial occupation are accounted for by adding an extra vegetation drag force and introducing porosity parameters. By defining a cell‐based porosity for the volumetric occupation and an edge‐based porosity for the flux exchange, the anisotropic property of the preferential flow can be well captured. Based on the finite volume method, the model is solved explicitly with a hybrid local time step/global maximum time step method and is parallelized using the Open MP techniques. The numerical grid convergence and quantitative accuracy of the model have been tested against a series of flume experiments for various configurations of rigid emergent vegetation over fixed or mobile beds. It is shown that using shallow water equations with anisotropic porosity, a constant drag coefficient can lead to numerical solutions of comparable accuracy as those complex empirical relations. Moreover, appropriate quantification of the stem‐scale turbulence effects brings significant improvement in modeling of the vegetation‐affected sediment transport.

A sediment transport theory based on distortion-free-boundary nonhomogeneous fluid flows

Hydrodynamic coefficient estimation for ship manoeuvring in shallow water using an optimal truncated LS-SVM

In order to accurately identify the ship's roll model parameters in shallow water, and solve the problems of difficult estimating nonlinear damping coefficients by traditional methods, a novel Nonlinear Least Squares - Support Vector Machine (NLS-SVM) is introduced. To illustrate the validity and applicability of the proposed method, simulation and decay tests data are combined and utilized to estimate unknown parameters and predict the roll motions. Firstly, simulation data is applied in the NLS-SVM model to obtain estimated damping parameters, compared with pre-defined parameters to verify the validity of the proposed method. Subsequently, decay tests data are used in identifying unknown parameters by utilizing traditional models and the new NLS-SVM model, the results illustrate that the intelligent method can improve the accuracy of parametric estimation, and overcome the conventional algorithms' weakness of difficult identification of the nonlinear damping parameter in the roll model. Finally, to show the wide applicability of the proposed model in shallow water, experimental data from various speeds and Under Keel Clearances (UKCs) are applied to identify the damping coefficients. Results reveal the potential of using the NLS-SVM for the problem of the roll motion in shallow water, and the effectiveness and accuracy are verified as well.

The independent set perturbation method for efficient computation of sensitivities with applications to data assimilation and a finite element shallow water model

The shallow water model in FLOW-3D can consist, analog to the 3D model, of multiple mesh blocks. The mesh blocks must comprise of two vertical cells. It is important to arrange the mesh in a way, that all the geometry and all the water is mapped in the bottom cell and the cell is not filled more than 90% (component and water). This is easy to achieve by defining the appropriate cell height of the bottom mesh layer. Another requirement is the definition of gravity in negative z-coordinate. Only the horizontal velocities in x-and y-directions are solved. The calculation of the pressure is reduced to a hydrostatic profile. To calculate the bottom shear stress the following approach is used.

Gravity waves, or imbalanced motions, that develop during the evolution of vortical flows in numerical models of the shallow water (SW) equations are examined in detail. The focus here is on nearly‐balanced flows, with small but non‐zero gravity‐wave activity. For properly initialized flows, it is reasonable to expect small GW activity when Froude numbers Fr < 1 and Rossby numbers Ro ≲ 1.The guiding principle in the present study is that an accurate representation of potential vorticity (PV) is the pre‐requisite to a fair assessment of the generation of gravity waves. The contour‐advective semi‐Lagrangian (CASL) algorithm for the SW equations is applied, as it shows a remarkable improvement in the simulation of PV. However, it is shown that the standard CASL algorithm for SW leads to a noticeable numerical generation of gravity waves. The false generation of gravity waves can equivalently be thought of as the false, or numerical, breakdown of balance.In order to understand the maintenance of balance in the SW equations, a hierarchy of CASL algorithms is introduced. The main idea behind the new hierarchy is to implement PV inversion partially, balancing algorithms directly within the SW algorithm, while still permitting imbalanced motions. The results of the first three members of the hierarchy, CA0 (standard CASL), CA1, and CA2, are described and are compared with the results of two other SW algorithms, a pseudo‐spectral and a semi‐Lagrangian one. The main body of results is obtained for a highly ageostrophic regime of flow, with|Ro|max ∼ 1 and Frmax ∼ 0.5, where sub‐index 'max' denotes maximum over the domain. Other flow regimes in the relevant parts of the Ro‐Fr parameter space are also explored. It is found that, for a given resolution and Froude number, there is an optimal CASL algorithm, i.e. one which gives rise to the least numerical generation of gravity waves.

A technique for determining the geoacoustic models in shallow water is described. For a horizontally stratified ocean and bottom, the method consists of measuring the magnitude and phase versus range of the pressure field due to a cw point source and numerically Hankel transforming these data to obtain the depth-dependent Green's function versus horizontal wavenumber. In shallow water, the Green's function contains prominent peaks at horizontal wavenumbers corresponding to the eigenvalues for any trapped and virtual modes excited in the waveguide. From the Green's function, one can obtain the geoacoustic model via either forward modeling or perturbative inverse techniques. In the forward modeling approach, a geoacoustic model for the bottom is obtained by computing the theoretical Green's function for various values of the bottom parameters and determining the parameter set which provides the best agreement with the experimental Green's function, particularly in the positions and relative magnitudes of the modal peaks. In the perturbative inverse technique, one uses the differences between the measured modal peaks and those predicted by a background model as input data to an integral equation, which is solved for the bottom geoacoustic parameters. These techniques are demonstrated using experimental data at 140 and 220 Hz. [Work supported by ONR.]