Lattice Boltzmann simulations of the transient shallow water flows

A two-dimensional lattice Boltzmann model (LBM) is presented for transient shallow water flows. The model is based on the shallow water equations coupled with the large eddy simulation model. In order to obtain accurate results efficiently, a multi-block lattice scheme is applied at the area where a local finer grid is needed for strong change in physical variables. The model is verified by applying to five cases with transient processes: (a) a tidal wave over steps; (b) a perturbation over a submerged hump; (c) partial dam break flow; (d) circular dam break flow; (e) interaction between a dam break surge and four square cylinders. The objectives of this study are to validate the two-dimensional LBM in transient flow simulation and provide the detailed transient processes in shallow water flows.     

Multi-block lattice Boltzmann simulations of solute transport in shallow water flows

We report a two-dimensional multi-block lattice Boltzmann model for solute transport in shallow water flows, which is developed based on the advection–diffusion equation for mass transport and the shallow water equations for the flows. A weighting factor is included in the centered scheme for improved accuracy. The model is firstly verified by simulating three benchmark tests: wind-driven circulation in a dish-shaped lake, jet-forced flow in a circular basin, and flow formed by two parallel streams containing different uniform concentrations at the same constant velocity; and then it is applied to a practical wind-induced flow, Baiyangdian Lake, which is characterized by irregular geometries and complex bathymetries. The numerical results have shown that the model is able to produce accurate and detailed results for both water flows and solute transport, which is attractive, especially for flows in narrow zones of practical terrains and certain areas with largely varying pollutant concentrations.     

Multilayer shallow water flow using lattice Boltzmann method with high performance computing

A multilayer lattice Boltzmann (LB) model is introduced to solve three-dimensional wind-driven shallow water flow problems. The multilayer LB model avoids the expensive Navier–Stokes equations and obtains stratified horizontal flow velocities as vertical velocities are relatively small and the flow is still within the shallow water regime. A single relaxation time BGK method is used to solve each layer coupled by the vertical viscosity forcing term. To increase solution stability, an implicit step is suggested to obtain flow velocities. The main advantage of using the LBM is that after selecting appropriate equilibrium distribution functions, the LB algorithm is only slightly modified for each layer and retains all the simplicities of the LBM within the high performance computing (HPC) environment. The performance of the parallel LB model for the multilayer shallow water equations is investigated on CPU-based HPC environments using OpenMP. We found that the explicit loop control with cache optimization in LBM gives better performance on execution time, speedup and efficiency than the implicit loop control as the number of processors increases. Numerical examples are presented to verify the multilayer LB model against analytical solutions. We demonstrate the model’s capability of calculating lateral and vertical distributions of velocities for wind-driven circulation over non-uniform bathymetry.     

A kinetic flux vector splitting scheme for shallow water flows

A kinetic flux vector splitting (KFVS) scheme for shallow water flows based on the collisionless Boltzmann equation is formulated and applied. The scheme is explicit and first order in space and time with stability governed by the Courant condition. The consistency of the KFVS scheme with the shallow water equations is proven using the equivalent differential equations approach. The accuracy and efficiency of the KFVS scheme in modeling complex flow features are compared to those of the Boltzmann Bhatnagar–Gross–Krook (BGK) scheme as well as a Riemann-based scheme. In particular, all schemes are applied to (i) strong shock waves, (ii) extreme expansion waves, (iii) a combination of strong shock waves and extreme expansion waves, and (iv) a one-dimensional dam break problem. Additionally, the KFVS, BGK and Riemann schemes are applied to a one-dimensional dam break problem for which laboratory data is available. These test cases reveal that all three schemes provide solutions of comparable accuracy, but the KFVS model is 1.5–2 times faster to execute than the BGK scheme and 2–3 times faster than the Riemann-based scheme. The absence of the collision term from the Boltzmann equation not only makes the mathematical formulation of KFVS easy but also helps elucidate this approach to the novice. The accuracy, efficiency, and simplicity of the KFVS scheme indicate its potential in modeling an array of water resources problems. Due to the scalar nature of the Boltzmann equation, the extension of the KFVS scheme to 2-D surface water flows is straightforward.     

Modelling dam-break flows over mobile beds using a 2D coupled approach

Dam-break flows usually propagate along rivers and floodplains, where the processes of fluid flow, sediment transport and bed evolution are closely linked. However, the majority of existing two-dimensional (2D) models used to simulate dam-break flows are only applicable to fixed beds. Details are given in this paper of the development of a 2D morphodynamic model for predicting dam-break flows over mobile beds. In this model, the common 2D shallow water equations are modified, so that the effects of sediment concentrations and bed evolution on the flood wave propagation can be considered. These equations are used together with the non-equilibrium transport equations for graded sediments and the equation of bed evolution. The governing equations are solved using a matrix method, thus the hydrodynamic, sediment transport and morphological processes can be jointly solved. The model employs an unstructured finite volume algorithm, with an approximate Riemann solver, based on the Roe-MUSCL scheme. A predictor–corrector scheme is used in time stepping, leading to a second-order accurate solution in both time and space. In addition, the model considers the adjustment process of bed material composition during the morphological evolution process. The model was first verified against results from existing numerical models and laboratory experiments. It was then used to simulate dam-break flows over a fixed bed and a mobile bed to examine the differences in the predicted flood wave speed and depth. The effects of bed material size distributions on the flood flow and bed evolution were also investigated. The results indicate that there is a great difference between the dam-break flow predictions made over a fixed bed and a mobile bed. At the initial stage of a dam-break flow, the rate of bed evolution could be comparable to that of water depth change. Therefore, it is often necessary to employ the turbid water governing equations using a coupled approach for simulating dam-break flows.     

Numerical calculation for bed variation in compound-meandering channel using depth integrated model without assumption of shallow water flow

In a compound meandering channel, patterns of flow structures and bed variations change with increasing water depth owing to complex momentum exchange between high-velocity flow in a main channel and low-velocity flows in flood plains. We have developed a new quasi-three-dimensional model without the shallow water assumption, i.e., hydrostatic pressure distribution; our method is known as the general bottom velocity computation (BVC) method. In this method, a set of depth-integrated equations, including depth-integrated momentum and vorticity equations, are prepared for evaluating bottom velocity and vertical velocity distributions. The objective of this study is to develop a bed variation calculation method for both single and compound meandering channels by using the BVC method coupled with a sediment transport model. This paper shows that the BVC method can reproduce the pattern change of bed variation in a compound meandering channel flow with increasing relative depth. The variation in sediment transport rate due to overbank flow is explained by experimental and computational results.     

Generic boundary conditions for lattice Boltzmann models and their application to advection and anisotropic dispersion equations

We address a “multi-reflection” approach to model Dirichlet and Neumann time-dependent boundary conditions in lattice Boltzmann methods for arbitrarily shaped surfaces. The multi-reflection condition for an incoming population represents a linear combination of the known population solutions. The closure relations are first established for symmetric and anti-symmetric parts of the equilibrium functions, independently of the nature of the problem. The symmetric part is tuned to build second- and third-order accurate Dirichlet boundary conditions for the scalar function specified by the equilibrium distribution. The focus is on two approaches to advection and anisotropic-dispersion equations (AADE): the equilibrium technique when the coefficients of the expanded equilibrium functions match the coefficients of the transformed dispersion tensor, and the eigenvalue technique when the coefficients of the dispersion tensor are built as linear combinations of the eigenvalue functions associated with the link-type collision operator. As a particular local boundary technique, the “anti-bounce-back” condition is analyzed. The anti-symmetric part of the generic closure relation allows to specify normal flux conditions without inversion of the diffusion tensor. Normal and tangential constraints are derived for bounce-back and specular reflections. The bounce-back closure relation is released from the non-physical tangential flux restriction at leading orders. Solutions for the Poisson equation and for convection–diffusion equations are presented for isotropic/anisotropic configurations with specified Dirichlet and Neumann boundary conditions.     

Use of Large-Eddy Simulation for the bed shear stress estimation over a dune

Environmental flows are generally characterized by complex bed morphology and high current speeds. Such configurations favor the formation of vortex structures that strongly affect hydrody-namics and sediment transport. Large-Eddy Simulation (LES) enables investigation of the dynam-ics of the largest turbulence scales and, thanks to enhanced calculation resources, has now become applicable for simulating environmental flows. In this paper, a LES approach is developed in a CFD code (TELEMAC-3D), which was originally developed to simulate free surface flows using RANS methods. The present developments involve implementing subgrid models, boundary con-ditions and numerical schemes suitable for LES. The LES version of TELEMAC-3D was validated by comparing results on the model with experimental data for flow past a cylinder. Then, the model was applied to a test case representing flow over dunes. After validating the hydrodynamics, the model was used to assess the bottom shear stress, using both a RANS and a LES approach. Com-parison highlighted the potential contribution of LES to investigating the hydrodynamic forces acting on the bottom.     

Two-dimensional numerical model of two-layer shallow water equations for confluence simulation

This study presents a finite-volume explicit method to solve 2D two-layer shallow water equations. This numerical model is intended to describe two-layer shallow flows in which the superposed layers differ in velocity, density and rheology in a two-dimensional domain. The rheological behavior of mudflow or debris flow is called the Bingham fluid. Thus, the shear stress on rigid bed can be derived from the constitutive equation. The computational approach adopts the HLL scheme, a novel approach for the purpose of computing a Godunov flux and solving the Riemann problem approximately proposed by Harten, Lax and van Leer, as a basic building block, treats the bottom slope by lateralizing the momentum flux, and refines the scheme using the Strang splitting to manage the frictional source term. This study successfully performed 2D two-layer shallow water computations on a rigid bed. The proposed numerical model can describe the variety of depths and velocities of substances including water and mud, when the hyperconcentrated tributary flows into the main river. The analytical results in this study will be valuable for further advanced research and for designing or planning hydraulic engineering structures.     

A reduced model for nonlinear dispersive waves in a rotating environment

Abstract

The simplest model for geophysical flows is one layer of a constant density fluid with a free surface, where the fluid motions occur on a scale in which the Coriolis force is significant. In the linear shallow water limit, there are non-dispersive Kelvin waves, localized near a boundary or near the equator, and a large family of dispersive waves. We study weakly nonlinear and finite depth corrections to these waves, and derive a reduced system of equations governing the flow. For this system we find approximate solitary Kelvin waves, both for waves traveling along a boundary and along the equator. These waves induce jets perpendicular to their direction of propagation, which may have a role in mixing. We also derive an equivalent reduced system for the evolution of perturbations to a mean geostrophic flow.     

Bubble coalescence in basalt flows: comparison of a numerical model with natural examples

Gas accumulation in magma may be aided by coalescence of bubbles because large coalesced bubbles rise faster than small bubbles. The observed size distribution of gas bubbles (vesicles) in lava flows supports the concept of post-eruptive coalescence. A numerical model predicts the effects of rise and coalescence consistent with observed features. The model uses given values for flow thickness, viscosity, volume percentage of gas bubbles, and an initial size distribution of bubbles together with a gravitational collection kernel to numerically integrate the stochastic collection equation and thereby compute a new size spectrum of bubbles after each time increment of conductive cooling of the flow. Bubbles rise and coalesce within a fluid interior sandwiched between fronts of solidification that advance inward with time from top and bottom. Bubbles that are overtaken by the solidification fronts cease to migrate. The model predicts the formation of upper and lower vesicle-rich zones separated by a vesicle-poor interior. The upper zone is broader, more vesicular, and has larger bubbles than the lower zone. Basaltic lava flows in northern California exhibit the predicted zonation of vesicularity and size distribution of vesicles as determined by an impregnation technique. In particular, the size distribution at the tops and bottoms of flows is essentially the same as the initial distribution, reflecting the rapid initial solidification at the bases and tops of the flows. Many large vesicles are present in the upper vesicular zones, consistent with expected formation as a result of bubble coalescence during solidification of the lava flows. Both the rocks and model show a bimodal or trimodal size distribution for the upper vesicular zone. This polymodality is explained by preferential coalescence of larger bubbles with subequal sizes. Vesicularity and vesicle size distribution are sensitive to atmospheric pressure because bubbles expand as they decompress during rise through the flow. The ratio of vesicularity in the upper to that in the lower part of a flow therefore depends not only on bubble rise and coalescence, but also on flow thickness and atmospheric pressure. Application of simple theory to the natural basalts suggests solidification of the basalts at 1.0±0.2 atm, consistent with the present atmospheric pressure. Paleobathymetry and paleoaltimetry are possible in view of the sensitivity of vesicle size distributions to atmospheric pressure. Thus, vesicular lava flows can be used to crudely estimate ancient elevations and/or sea level air pressure.     

A robust well-balanced finite volume model for shallow water flows with wetting and drying over irregular terrain

An unstructured Godunov-type finite volume model is developed for the numerical simulation of geometrically challenging two-dimensional shallow water flows with wetting and drying over convoluted topography. In the framework of sloping bottom model, a modified formulation of shallow water equations is used to preserve mass conservation during flooding and recession. The key ingredient of the model is the use of this combination of the sloping bottom model and the modified shallow water equations to provide a robust technique for wet/dry fronts tracking and, together with centered discretization of the bed slope source term, to exactly preserve the static flow on irregular topographies. The variable reconstruction technique ensures nonnegative reconstructed water depth and reasonable reconstructed velocity, and the friction terms are solved by semi-implicit scheme that does not invert the direction of velocity components. The robustness and accuracy of the proposed model are assessed by comparing numerical and reference results of extensive test cases. Moreover, the results of a dam-break flooding over real topography are presented to show the capability of the model on field-scale application.     

High-resolution numerical simulation of turbulence in natural waterways

We develop an efficient and versatile numerical model for carrying out high-resolution simulations of turbulent flows in natural meandering streams with arbitrarily complex bathymetry. The numerical model solves the 3D, unsteady, incompressible Navier-Stokes and continuity equations in generalized curvilinear coordinates. The method can handle the arbitrary geometrical complexity of natural streams using the sharp-interface curvilinear immersed boundary (CURVIB) method of Ge and Sotiropoulos (2007) [1]. The governing equations are discretized with three-point, central, second-order accurate finite-difference formulas and integrated in time using an efficient, second-order accurate fractional step method. To enable efficient simulations on grids with tens of millions of grid nodes in long and shallow domains typical of natural streams, the algebraic multigrid (AMG) method is used to solve the Poisson equation for the pressure coupled with a matrix-free Krylov solver for the momentum equations. Depending on the desired level of resolution and available computational resources, the numerical model can either simulate, via direct numerical simulation (DNS), large-eddy simulation (LES), or unsteady Reynolds-averaged Navier-Stokes (URANS) modeling. The potential of the model as a powerful tool for simulating energetic coherent structures in turbulent flows in natural river reaches is demonstrated by applying it to carry out LES and URANS in a 50-m long natural meandering stream at resolution sufficiently fine to capture vortex shedding from centimeter-scale roughness elements on the bed. The accuracy of the simulations is demonstrated by comparisons with experimental data and the relative performance of the LES and URANS models is also discussed.     

A finite element model for simulating surface run‐off and unsaturated seepage flow in the shallow subsurface

This work introduces water–air two‐phase flow into integrated surface–subsurface flow by simulating rainfall infiltration and run‐off production on a soil slope with the finite element method. The numerical model is formulated by partial differential equations for hydrostatic shallow flow and water–air two‐phase flow in the shallow subsurface. Finite element computing formats and solution strategies are presented to obtain a numerical solution for the coupled model. An unsaturated seepage flow process is first simulated by water–air two‐phase flow under the atmospheric pressure boundary condition to obtain the rainfall infiltration rate. Then, the rainfall infiltration rate is used as an input parameter to solve the surface run‐off equations and determine the value of the surface run‐off depth. In the next iteration, the pressure boundary condition of unsaturated seepage flow is adjusted by the surface run‐off depth. The coupling process is achieved by updating the rainfall infiltration rate and surface run‐off depth sequentially until the convergence criteria are reached in a time step. A well‐conducted surface run‐off experiment and traditional surface–subsurface model are used to validate the new model. Comparisons with the traditional surface–subsurface model show that the initiation time of surface run‐off calculated by the proposed model is earlier and that the water depth is larger, thus providing values that are closer to the experimental results.     

Coupling water flow and solute transport into a physically-based surface-subsurface hydrological model

A distributed-parameter physically-based solute transport model using a novel approach to describe surface-subsurface interactions is coupled to an existing flow model. In the integrated model the same surface routing and mass transport equations are used for both hillslope and channel processes, but with different parametrizations for these two cases. For the subsurface an advanced time-splitting procedure is used to solve the advection-dispersion equation for transport and a standard finite element scheme is used to solve Richards equation for flow. The surface-subsurface interactions are resolved using a mass balance-based surface boundary condition switching algorithm that partitions water and solute into actual fluxes across the land surface and changes in water and mass storage. The time stepping strategy allows the different time scales that characterize surface and subsurface water and solute dynamics to be efficiently and accurately captured. The model features and performance are demonstrated in a series of numerical experiments of hillslope drainage and runoff generation.     

A reduced 3D hydrodynamic model of a shallow,long, and weakly curved stream

A reduced three-dimensional mathematical model of a free-flow stream in nondeformable channels (a model of a long shallow flow), proposed earlier, has been studied analytically and numerically. The reduced model has been verified by direct numerical simulation of the flow by full hydrodynamic models in COMSOL finite-element software complex for laminar and turbulent flows of a viscous fluid. The obtained results show that the proposed reduced model of a shallow weakly curved stream flow adequately describes its hydrodynamics, so it can be used in systems of complex simulation of the ecology of water objects and the use of water resources.     

Balancing the source terms in a SPH model for solving the shallow water equations

A shallow flow generally features complex hydrodynamics induced by complicated domain topography and geometry. A numerical scheme with well-balanced flux and source term gradients is therefore essential before a shallow flow model can be applied to simulate real-world problems. The issue of source term balancing has been exhaustively investigated in grid-based numerical approaches, e.g. discontinuous Galerkin finite element methods and finite volume Godunov-type methods. In recent years, a relatively new computational method, smooth particle hydrodynamics (SPH), has started to gain popularity in solving the shallow water equations (SWEs). However, the well-balanced problem has not been fully investigated and resolved in the context of SPH. This work aims to discuss the well-balanced problem caused by a standard SPH discretization to the SWEs with slope source terms and derive a corrected SPH algorithm that is able to preserve the solution of lake at rest. In order to enhance the shock capturing capability of the resulting SPH model, the Monotone Upwind-centered Scheme for Conservation Laws (MUSCL) is also explored and applied to enable Riemann solver based artificial viscosity. The new SPH model is validated against several idealized benchmark tests and a real-world dam-break case and promising results are obtained.     

A Cartesian cut cell method for shallow water flows with moving boundaries

A new computational method for the calculation of shallow water flows with moving physical boundaries is presented. The procedure can cope with shallow water problems having arbitrarily complex geometries and moving boundary elements. Although the method provides a fully boundary-fitted capability, no mesh generation is required in the conventional sense. Solid regions are simply cut out of a background Cartesian mesh with their boundaries represented by different types of cut cell. Moving boundaries are accommodated by up-dating the local cut cell information on a stationary background mesh as the boundaries move. No large-scale re-meshing is required. For the flow calculations, a multi-dimensional high resolution upwind finite volume scheme is used in conjunction with an efficient approximate Riemann solver at flow interfaces, and an exact Riemann solution for a moving piston at moving boundary elements. The method is validated for test problems that include a ship's hull moving at supercritical velocity and two hypothetical landslide events where material plunges laterally into a quiescent shallow lake and a fiord.     

A new analytical solution for assessing climate change impacts on subsurface temperature

Groundwater temperature is an important water quality parameter that affects species distributions in subsurface and surface environments. To investigate the response of subsurface temperature to atmospheric climate change, an analytical solution is derived for a one‐dimensional, transient conduction–advection equation and verified with numerical methods using the finite element code SUTRA. The solution can be directly applied to forward model the impact of future climate change on subsurface temperature profiles or inversely applied to produce a surface temperature history from measured borehole profiles. The initial conditions are represented using superimposed linear and exponential functions, and the boundary condition is expressed as an exponential function. This solution expands on a classic solution in which the initial and boundary conditions were restricted to linear functions. The exponential functions allow more flexibility in matching climate model projections (boundary conditions) and measured temperature–depth profiles (initial conditions). For example, measured borehole temperature data from the Sendai Plain and Tokyo, Japan, were used to demonstrate the improved accuracy of the exponential function for replicating temperature–depth profiles. Also, the improved accuracy of the exponential boundary condition was demonstrated using air temperature anomaly data from the Intergovernmental Panel on Climate Change. These air temperature anomalies were then used to forward model the effect of surficial thermal perturbations in subsurface environments with significant groundwater flow. The simulation results indicate that recharge can accelerate shallow subsurface warming, whereas upward groundwater discharge can enhance deeper subsurface warming. Additionally, the simulation results demonstrate that future groundwater temperatures obtained from the proposed analytical solution can deviate significantly from those produced with the classic solution. Copyright © 2013 John Wiley & Sons, Ltd.     

The groundwater–land-surface–atmosphere connection: Soil moisture effects on the atmospheric boundary layer in fully-coupled simulations

This study combines a variably-saturated groundwater flow model and a mesoscale atmospheric model to examine the effects of soil moisture heterogeneity on atmospheric boundary layer processes. This parallel, integrated model can simulate spatial variations in land-surface forcing driven by three-dimensional (3D) atmospheric and subsurface components. The development of atmospheric flow is studied in a series of idealized test cases with different initial soil moisture distributions generated by an offline spin-up procedure or interpolated from a coarse-resolution dataset. These test cases are performed with both the fully-coupled model (which includes 3D groundwater flow and surface water routing) and the uncoupled atmospheric model. The effects of the different soil moisture initializations and lateral subsurface and surface water flow are seen in the differences in atmospheric evolution over a 36-h period. The fully-coupled model maintains a realistic topographically-driven soil moisture distribution, while the uncoupled atmospheric model does not. Furthermore, the coupled model shows spatial and temporal correlations between surface and lower atmospheric variables and water table depth. These correlations are particularly strong during times when the land-surface temperatures trigger shifts in wind behavior, such as during early morning surface heating.