ADMESH:?An advanced, automatic unstructured mesh generator for shallow water models

In this paper, we present the development and application of a two-dimensional, automatic unstructured mesh generator for shallow water models called Admesh. Starting with only target minimum and maximum element sizes and points defining the boundary and bathymetry/ topography of the domain, the goal of the mesh generator is to automatically produce a high-quality mesh from this minimal set of input. From the geometry provided, properties such as local features, curvature of the boundary, bathymetric/topographic gradients, and approximate flow characteristics can be extracted, which are then used to determine local element sizes. The result is a high-quality mesh, with the correct amount of refinement where it is needed to resolve all the geometry and flow characteristics of the domain. Techniques incorporated include the use of the so-called signed distance function, which is used to determine critical geometric properties, the approximation of piecewise linear coastline data by smooth cubic splines, a so-called mesh function used to determine element sizes and control the size ratio of neighboring elements, and a spring-based force equilibrium approach used to improve the element quality of an initial mesh obtained from a simple Delaunay triangulation. Several meshes of shallow water domains created by the new mesh generator are presented.     

Simulating the hydraulic characteristics of the lower Yellow River by the finite‐volume technique

The finite‐volume technique is used to solve the two‐dimensional shallow‐water equations on unstructured mesh consisting of quadrilateral elements. In this paper the algorithm of the finite‐volume method is discussed in detail and particular attention is paid to accurately representing the complex irregular computational domain. The lower Yellow River reach from Huayuankou to Jiahetan is a typical meandering river. The generation of the computational mesh, which is used to simulate the flood, is affected by the distribution of water works in the river channel. The spatial information about the two Yellow River levee, the protecting dykes, and those roads that are obviously higher than the ground, need to be used to generate the computational mesh. As a result these dykes and roads locate the element interfaces of the computational mesh. In the model the finite‐volume method is used to solve the shallow‐wave equations, and the Osher scheme of the empirical function is used to calculate the flux through the interface between the neighbouring elements. The finite‐volume method has the advantage of using computational domain with complex geometry, and the Osher scheme is a method based on characteristic theory and is a monotone upwind numerical scheme with high resolution. The flood event with peak discharge of 15 300 m³/s, occurring in the period from 30 July to 10 August 1982, is simulated. The estimated result indicates that the simulation method is good for routing the flood in a region with complex geometry. Copyright © 2002 John Wiley & Sons, Ltd.     

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.     

Modelling the influence of spatially varying hydrodynamics on the cross-sectional stability of double inlet systems

The cross-sectional stability of double inlet systems is investigated using an exploratory model that combines Escoffier’s stability concept for the evolution of the inlet’s cross-sectional area with a two-dimensional, depth-averaged (2DH) hydrodynamic model for tidal flow. The model geometry consists of four rectangular compartments, each with a uniform depth, associated with the ocean, tidal inlets and basin. The water motion, forced by an incoming Kelvin wave at the ocean’s open boundary and satisfying the linear shallow water equations on the f -plane with linearised bottom friction, is in each compartment written as a superposition of eigenmodes, i.e. Kelvin and Poincaré waves. A collocation method is employed to satisfy boundary and matching conditions. The analysis of resulting equilibrium configurations is done using flow diagrams.

Model results show that internally generated spatial variations in the water motion are essential for the existence of stable equilibria with two inlets open. In the hydrodynamic model used in the paper, both radiation damping into the ocean and basin depth effects result in these necessary spatial variations. Coriolis effects trigger an asymmetry in the stable equilibrium cross-sectional areas of the inlets. Furthermore, square basin geometries generally correspond to significantly larger equilibrium values of the inlet cross-sections. These model outcomes result from a competition between a destabilising (caused by inlet bottom friction) and a stabilising mechanism (caused by spatially varying local pressure gradients over the inlets).

     

A wetting and drying algorithm with a combined pressure/free-surface formulation for non-hydrostatic models

A wetting and drying method for free-surface problems for the three-dimensional, non-hydrostatic Navier–Stokes equations is proposed. The key idea is to use a horizontally fixed mesh and to apply different boundary conditions on the free-surface in wet and dry zones. In wet areas a combined pressure/free-surface kinematic boundary condition is applied, while in dry areas a positive water level and a no-normal flow boundary condition are enforced. In addition, vertical mesh movement is performed to accurately represent the free-surface motion. Non-physical flow in the remaining thin layer in dry areas is naturally prevented if a Manning–Strickler bottom drag is used. The treatment of the wetting and drying processes applied through the boundary condition yields great flexibility to the discretisation used. Specifically, a fully unstructured mesh with any finite element choice and implicit time discretisation method can be applied. The resulting method is mass conservative, stable and accurate. It is implemented within Fluidity-ICOM [1] and verified against several idealized test cases and a laboratory experiment of the Okushiri tsunami.     

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.     

A lattice Boltzmann model coupled with a Large Eddy Simulation model for flows around a groyne

This present paper proposes a two-dimensional lattice Boltzmann model coupled with a Large Eddy Simulation (LES) model and applies it to flows around a non-submerged groyne in a channel. The LES of shallow water equations is efficiently performed using the Lattice Boltzmann Method (LBM) and the turbulence can be taken into account in conjunction with the Smagorinsky Sub-Grid Stress (SGS) model. The bounce-back scheme of the non-equilibrium part of the distribution function is used to determine the unknown distribution functions at inflow boundary, the zero gradient of the distribution function is set normal to outflow boundary to obtain the unknown distribution functions here and the bounce-back scheme, which states that an incoming particle towards the boundary is bounced back into fluid, is applied to the solid wall to ensure non-slip boundary conditions. The initial flow field is defined firstly and then is used to calculate the local equilibrium distributions as initial conditions of the distribution functions. These coupled models successfully predict the flow characteristics, such as circulating flow, velocity and water depth distributions. The comparisons between the simulated results and the experimental data show that the model scheme has the capacity to solve the complex flows in shallow water with reasonable accuracy and reliability.     

Analysis of high-frequency local coupling instability induced by multi-transmitting formula—P-SV wave simulation in a 2D waveguide

Local coupling instability will occur when the numerical scheme of absorbing boundary condition and that of the field wave equation allow energies to spontaneously enter into the computational domain. That is, the two schemes support common wave solutions with group velocity pointed into the computation domain. The key to eliminate local coupling instability is to avoid such wave solutions. For lumped-mass finite element simulation of P-SV wave motion in a 2D waveguide, an approach for stable implementation of high order multi-transmitting formula is provided. With a uniform rectangular mesh, it is proven and validated that high-frequency local coupling instability can be eliminated by setting the ratio of the element size equal to or greater than \(\sqrt 2 \) times the ratio of the P wave velocity to the S wave velocity. These results can be valuable for dealing instability problems induced by other absorbing boundary conditions.     

Application of a finite element model (TELEMAC) to computing the wind induced response of the Irish Sea

Initially a brief overview of the problem of computing the wind-induced circulation on the west coast of Britain is reviewed together with storm surge modelling. To date this work has primarily been performed with finite difference models. However, here new work is presented using a finite element model with a range of mesh refinements in shallow water regions to examine the influence of mesh resolution upon the wind-induced circulation off the west coast of Britain. Steady state current fields are computed for uniform westerly and southerly winds and compared with a uniform grid (of order 7 km) finite difference model solution. Calculations show that in deep water regions away from the coastal influence, the large-scale circulation features in the finite element solution are in good agreement with those found in the finite difference model. This suggests that they can be adequately resolved on a 7 km mesh. In the nearshore region and within estuaries a significantly finer mesh is required, with the variable mesh finite element model showing significant small scale variability in the nearshore area. Refining the mesh in the Mersey and using an accurate topographic data set, shows that although the larger scale features in the estuary can be resolved in the coarser mesh model, accurate topography is required to model their exact location. In addition smaller scale features are found that were not resolved in the coarser mesh models. Due to the effects of “wetting and drying” and the importance of non-linear processes in shallow regions difficulties occurred in de-tiding the full solution in order to determine the wind forced residual. Determining the wind forced solution in shallow water from a calculation in which wind and tidal forcing are included poses problems as to how to “de-tide” the solution in such a highly non-linear region. An approach based upon the harmonic analysis of the total solution, rather than subtracting a “tide only” solution is shown to be most effective and has implications for storm surge prediction.General and specific conclusions on the importance of highly accurate bathymetry, good mesh resolution and de-tiding method upon the accuracy of the wind forced solution in nearshore regions are summarized in the final part of the paper. The implications for storm surge prediction together with suggestions for future research to enhance the accuracy of storm surge prediction, namely “the way forward” are given at the end of the paper.     

Conservative identification of nodal transmissivities: theory,computational experiments and application to the Rharb coastal aquifer (Morocco)

Abstract

In order to calculate the transmissivity from the inverse problem corresponding to the groundwater flow in an isotropic horizontal aquifer, a numerical conservative approach is tested. The method deals with triangulation of the domain and applies the conservation of mass to elements of the mesh using the harmonic mean for internodal transmissivities. An optimal sweeping algorithm is used to evaluate nodal transmissivities from one element to another with a minimal relative error accumulation. The practical importance of the method is demonstrated through two synthetic examples representing those experienced in the field, then through application to a Moroccan aquifer. The computed hydraulic head is well fitted to the reference one, which confirms the validity of the identified transmissivity model.     

On the well posedness of some wave formulations of the shallow water equations

According to the theory of characteristics, the number of boundary conditions required for the adequate definition of a PDE problem is equal to the number of characteristic half planes entering the domain associated with the PDE problem.This theory was applied to the primitive form of the shallow water equations two decades ago to determine the number of initial and boundary conditions required by these equations. The results of this early study are remembered here. Subsequently, the same theory is applied to some wave formulations of the shallow water equations (the wave continuity and primitive momentum equation model, and the wave continuity and wave momentum equation model).Circumstances under which the number of boundary conditions required by the mathematical model can be reduced to the number of available boundary conditions are discussed for both the primitive and wave formulations of the shallow water equations.     

Numerical simulation of groundwater flow patterns using flux as upper boundary

Since the 1960s, most of the studies on groundwater flow systems by analytical and numerical modelling have been based on given‐head upper boundaries. The disadvantage of the given‐head approach is that the recharge into and discharge from a basin vary with changes in hydraulic conductivity and/or basin geometry. Consequently, flow patterns simulated with given‐head boundaries but with different hydraulic conductivities and/or basin geometry may not reflect the effects of these variables. We conducted, therefore, numerical simulations of groundwater flow in theoretical drainage basins using flux as the upper boundary and realistically positioned fluid‐potential sinks while changing the infiltration intensity, hydraulic conductivities, and geometric configuration of the basin. The simulated results demonstrate that these variables are dominant factors controlling the flow pattern in a laterally closed drainage basin. The ratio of infiltration intensity to hydraulic conductivity (R_ic) has been shown to be an integrated pattern‐parameter in a basin with a given geometric configuration and possible fluid‐potential‐sink distribution. Successively, the changes in flow patterns induced by stepwise reductions in R_ic are identical, regardless of whether the reductions are due to a decrease in infiltration intensity or an increase in hydraulic conductivity. The calculated examples show five sequential flow patterns containing (i) only local, (ii) local–intermediate, (iii) local–intermediate–regional, (iv) local–regional, and (v) just regional flow systems. The R_ic was found to determine also whether a particular sink is active or not as a site of discharge. Flux upper boundary is preferable for numerical simulation when discussing the flow patterns affected by a change of infiltration, the hydraulic conductivity, or the geometry of a basin. Copyright © 2012 John Wiley & Sons, Ltd.     

Seismic analysis of semi-sine shaped alluvial hills above subsurface circular cavity

In this study, a seismic analysis of semi-sine shaped alluvial hills above a circular underground cavity subjected to propagating oblique SH-waves using the half-plane time domain boundary element method (BEM) was carried out. By dividing the problem into a pitted half-plane and an upper closed domain as an alluvial hill and applying continuity/boundary conditions at the interface, coupled equations were constructed and ultimately, the problem was solved step-by-step in the time domain to obtain the boundary values. After solving some verification examples, a semi-sine shaped alluvial hill located on an underground circular cavity was successfully analyzed to determine the amplification ratio of the hill surface. For sensitivity analysis, the effects of the impedance factor and shape ratio of the hill were also considered. The ground surface responses are illustrated as three-dimensional graphs in the time and frequency domains. The results show that the material properties of the hill and their heterogeneity with the underlying half-space had a significant effect on the surface response.

     

A low‐dimensional hillslope‐based catchment model for layered groundwater flow

Despite the strong interaction between surface and subsurface waters, groundwater flow representation is often oversimplified in hydrological models. For instance, the interplay between local or shallow aquifers and deeper regional‐scale aquifers is typically neglected. In this work, a novel hillslope‐based catchment model for the simulation of combined shallow and deep groundwater flow is presented. The model consists of the hillslope‐storage Boussinesq (hsB) model representing shallow groundwater flow and an analytic element (AE) model representing deep regional groundwater flow. The component models are iteratively coupled via a leakage term based on Darcy's law, representing delayed recharge to the regional aquifer through a low conductivity layer. Simulations on synthetic single hillslopes and on a two‐hillslope open‐book catchment are presented, and the results are compared against a benchmark three‐dimensional Richards equation model. The impact of hydraulic conductivity, hillslope plan geometry (uniform, convergent, divergent), and hillslope inclination (0.2%, 5%, and 30%) under drainage and recharge conditions are examined. On the single hillslopes, good matches for heads, hydrographs, and exchange fluxes are generally obtained, with the most significant differences in outflows and heads observed for the 30% slope and for hillslopes with convergent geometry. On the open‐book catchment, cumulative outflows are overestimated by 1–4%. Heads in the confined and unconfined aquifers are adequately reproduced throughout the catchment, whereas exchange fluxes are found to be very sensitive to the hillslope drainable porosity. The new model is highly efficient computationally compared to the benchmark model. The coupled hsB/AE model represents an alternative to commonly used groundwater flow representations in hydrological models, of particular appeal when surface–subsurface exchanges, local aquifer–regional aquifer interactions, and low flows play a key role in a watershed's dynamics. Copyright © 2011 John Wiley & Sons, Ltd.     

Capturing the residence time boundary layer—application to the Scheldt Estuary

At high Peclet number, the residence time exhibits a boundary layer adjacent to incoming open boundaries. In a Eulerian model, not resolving this boundary layer can generate spurious oscillations that can propagate into the area of interest. However, resolving this boundary layer would require an unacceptably high spatial resolution. Therefore, alternative methods are needed in which no grid refinement is required to capture the key aspects of the physics of the residence time boundary layer. An extended finite element method representation and a boundary layer parameterisation are presented and tested herein. It is also explained how to preserve local consistency in reversed time simulations so as to avoid the generation of spurious residence time extrema. Finally, the boundary layer parameterisation is applied to the computation of the residence time in the Scheldt Estuary (Belgium/The Netherlands). This timescale is simulated by means of a depth-integrated, finite element, unstructured mesh model, with a high space–time resolution. It is seen that the residence time temporal variations are mainly affected by the semi-diurnal tides. However, the spring–neap variability also impacts the residence time, particularly in the sandbank and shallow areas. Seasonal variability is also observed, which is induced by the fluctuations over the year of the upstream flows. In general, the residence time is an increasing function of the distance to the mouth of the estuary. However, smaller-scale fluctuations are also present: they are caused by local bathymetric features and their impact on the hydrodynamics.     

Local transmitting boundaries for transient elastic analysis

The aim of this paper is to investigate and develop alternative methods of analyzing problems in dynamic soil–structure-interaction (SSI). The interaction means that the amplitude of structural response is effected by additional energy dissipation through radiation and material damping in the soil. The surrounding soft soil behaves as a natural damper for a massive and stiff structure supported or embedded in it. The main focus is the major difficulty posed by such an analysis — the phenomena of waves that radiate outward from the excited structures towards infinity. In numerical calculations only a finite region of the foundation medium is analyzed and something is done to prevent the outgoing radiation waves from reflecting at the boundary region.Development of a simple and efficient finite element (FE) procedure for the solution directly in the time domain of transient SSI problems is the main concern. The central feature of the procedure is local absorbing boundaries used to render the computational domain finite. These boundaries are local in both time and space and are completely defined by a pair of symmetric stiffness and damping matrices. As the effort for implementing them is the same as for the impedance boundary condition (BC) considering the angle of incidence, standard assembly procedure can be used. Due to the local nature they also preserve the overall structure of the global equations of motion. Even though the focus is in the time domain the same equations of motions can be used to determine the solution under time-harmonic excitation directly in the frequency domain. Explicit formulae for the element matrices are included in the paper and numerical examples for transient radiation model problems to illustrate the validity and accuracy of the new procedures, are given.     

基于卷积完全匹配层的非规则网格时域有限元探地雷达数值模拟

复频移完全匹配层(Complex Frequency-Shifted PML,CFS-PML)在长时间时域计算中对凋落波、倏失波具有好的吸收效果,并被广泛应用于时域有限差分模拟中.而本文采用卷积方法将CFS-PML应用于时域有限元求解GPR波动方程的数值模拟中.论文以TM波为例,推导了基于CPML(Convolutional PML)边界的时域有限元GPR波动方程求解公式,采用Newmark-β方法对时间导数进行离散,有效改善了时域有限元GPR数值计算程序的稳定性.并以狭长模型为例,开展了CPML边界中关键参数m、R和κ的选取实验,通过对比反射误差大小确定了综合最优参数组合.相同时刻UPML与CPML波场快照、3个检测点的反射误差比较,说明CPML较UPML具有更好的吸收效果.最后,采用非规则四边形网格对1个复杂GPR模型进行剖分,应用加载CPML边界条件的FETD程序对该模型进行了正演,得到了二维剖面法、宽角法正演GPR剖面图,说明非规则四边形对复杂模型的良好适应性,基于CPML边界条件的FETD可有效减少边界反射误差,能实现对任意复杂不规则模型的正演模拟.     

Application of an unstructured mesh model to storm surge propagation in the Mersey estuary region of the Irish Sea

The storm surge period of 13–16 November 1977 when there was a major positive surge followed by a negative surge in the Irish Sea is investigated using a two-dimensional unstructured mesh model of the west coast of Britain. The model accounts for tidal and external surge forcing across its open boundaries which are situated in the Celtic Sea and off the west coast of Scotland. Although this period has been examined previously using a uniform finite-difference model, and a finite element model, neither of these could resolve the Mersey estuary which is the focus of the present study. By using a finite element model with very high mesh resolution within the Mersey, the spatial variability of surge elevations and currents within the Mersey to rapidly changing surge dynamics can be examined. The mesh in the model varies from about 7 km in deep water, to the order of 100 m in the Mersey, with the largest mesh length reaching 17 km in deep offshore regions, and smallest of order 26 m occurring in shallow coastal regions of the Mersey estuary. The model accounts for wetting/drying which occurs in shallow water coastal areas. Calculations showed that during the positive surge period, the amplitude and speed of propagation of the surge was largest in the deep water channels. This gave rise to significant spatial variability of surge elevations and currents within the estuary. As wind stresses decreased over the Irish Sea, a negative surge occurred over Liverpool Bay and at the entrance to the Mersey. However, within the Mersey there was a local positive surge which continued to propagate down the estuary. This clearly showed that although the large scale response of the Irish Sea to changing wind fields occurred rapidly, the response in the Mersey was much slower. These calculations with a west coast variable mesh model that included a high-resolution representation of the Mersey revealed for the first time how elevations and currents within the Mersey responded to Irish Sea surges that rapidly changed from positive to negative.     

Influence of advection on scales of ecological studies in a coastal equilibrium flow

Spatial patterns are generated as a result of the coupling between biogeochemical and physical processes and the ability to capture and reproduce patchiness is crucial for the better comprehension of an ecosystem and its response to external perturbations. A 1D reaction–diffusion–advection equation is used to investigate the formation of patterns and relevant time and spatial scales and thus define an approach for the determination of a critical domain size that allows differentiation of the role of local and internal cycling from advective fluxes across the open boundaries in a shallow coastal ecosystem. By using a 3D numerical model, in conjunction with an extensive field data set, it is shown that domain sizes must be larger than this critical value in order to capture the patterns generated within the system. For smaller domains, the evolution of the system is controlled by transport processes across the boundaries misleading the interpretation of the internal ecological dynamics. The study of the influence of boundary fluxes on ecological patchiness was motivated by the need to define the size of the domain necessary for the assessment of the impact of a sewage outflow on a coastal regime.