Modelling estuarine and coastal flows using an unstructured triangular finite volume algorithm

Details are given of the development and application of a numerical model for predicting free-surface flows in estuarine and coastal basins using the finite volume method. Both second- and third-order accurate and oscillation free explicit numerical schemes have been used to solve the shallow water equations. The model deploys an unstructured triangular mesh and incorporates two types of mesh layouts, namely the ‘cell centred’ and ‘mesh vertex’ layouts, and provides a powerful mesh generator in which a user can adjust the mesh-size distribution interactively to create a desirable mesh. The quality of mesh has been shown to have a major impact on the overall performance of the numerical model.The model has been applied to simulate two-dimensional dam break flows for which transient water level distributions measured within a laboratory flume were available. In total 12 model runs were undertaken to test the model for various flow conditions. These conditions include: (1) different bed slopes (ranging from zero to 0.8%), (2) different upstream and downstream water level conditions, and (3) initially wet and dry bed conditions, downstream of the dam. Detailed comparisons have been made between model predicted and measured water levels and good agreement achieved between both sets of results. The model was then used to predict water level and velocity distributions in a real estuary, i.e. the Ribble Estuary, where the bed level varies rapidly at certain locations. In order to model the whole estuary, a 1-D numerical model has also been used to model the upper part of the estuary and this model was linked dynamically to the 2-D model. Findings from this application are given in detail.     

A 3-D phase-averaged model for shallow-water flow with waves in vegetated water

A 3-D shallow-water flow model has been developed to simulate the flow in coastal vegetated waters with short waves. The model adopts the 3-D phase-averaged shallow-water flow equations with radiation stresses induced by short waves. It solves the governing equations using an implicit finite volume method based on quadtree rectangular mesh in the horizontal plane and stretching mesh in the vertical direction. The flow model is coupled with a spectral wave deformation model called CMS-Wave. The wave model solves the spectral wave-action balance equation and provides wave characteristics to the flow model. The model considers the effects of vegetation on currents and waves by including the drag and inertia forces of vegetation in the momentum equations and the wave energy loss due to vegetation resistance in the wave-action balance equation. The model has been tested using several sets of laboratory experiments, including steady flows in a straight channel with submerged vegetation and in a compound channel with vegetated floodplain and random waves through a vegetated channel and on a vegetated beach slope. The calculated water levels, current velocities, and wave heights are in general good agreement with the measured data.     

基于局部加密非结构网格的海洋可控源电磁法三维有限元正演

本文基于非结构网格实现了海洋可控源电磁法三维有限元正演模拟.该算法采用完全非结构网格剖分,可以模拟任意起伏地形和复杂地电模型.为了避免场源的奇异性,采用一次场/二次场分解算法,一次场由基于Schelkunoff势函数的一维解析公式得到.为了提高算法的精度和效率,采用对测点附近单元和异常体区域进行体积约束加密的方法,实现了非结构网格的局部加密.一、二维模型计算和分析表明,本文采用的局部加密方法能够明显地改善算法的精度,最大相对误差基本在1%以内.对三维模型计算及对比分析,说明了该算法对三维可控源电磁正演的实用性.复杂海底地形模型的正演模拟表明,海底地形对电磁场的影响很大,在进行海洋可控源电磁资料解释时,地形的影响有必要考虑在内.     

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.     

Numerical modelling of bifurcation and localisation in cohesive-frictional materials

Methods are reviewed for analysing highly localised failure and bifurcation modes in discretised mechanical systems as typically arise in numerical simulations of failure in soils, rocks, metals and concrete. By the example of a plane-strain biaxial test it is shown that strain softening and lack of normality in elasto-plastic constitutive equations and the ensuing loss of ellipticity of the governing field equations cause a pathological mesh dependence of numerical solutions for such problems, thus rendering the results effectively meaningless. The need for introduction of higher-order continuum models is emphasised to remedy this shortcoming of the conventional approach. For one such a continuum model, namely the unconstrained Cosserat continuum, it is demonstrated that meaningful and convergent solutions (in the sense that a finite width of the localisation zone is computed upon mesh refinement) can be obtained.     

MAST-2D diffusive model for flood prediction on domains with triangular Delaunay unstructured meshes

A new methodology for the solution of the 2D diffusive shallow water equations over Delaunay unstructured triangular meshes is presented. Before developing the new algorithm, the following question is addressed: it is worth developing and using a simplified shallow water model, when well established algorithms for the solution of the complete one do exist?The governing Partial Differential Equations are discretized using a procedure similar to the linear conforming Finite Element Galerkin scheme, with a different flux formulation and a special flux treatment that requires Delaunay triangulation but entire solution monotonicity. A simple mesh adjustment is suggested, that attains the Delaunay condition for all the triangle sides without changing the original nodes location and also maintains the internal boundaries. The original governing system is solved applying a fractional time step procedure, that solves consecutively a convective prediction system and a diffusive correction system. The non linear components of the problem are concentrated in the prediction step, while the correction step leads to the solution of a linear system of the order of the number of computational cells. A semi-analytical procedure is applied for the solution of the prediction step. The discretized formulation of the governing equations allows to handle also wetting and drying processes without any additional specific treatment. Local energy dissipations, mainly the effect of vertical walls and hydraulic jumps, can be easily included in the model.Several numerical experiments have been carried out in order to test (1) the stability of the proposed model with regard to the size of the Courant number and to the mesh irregularity, (2) its computational performance, (3) the convergence order by means of mesh refinement. The model results are also compared with the results obtained by a fully dynamic model. Finally, the application to a real field case with a Venturi channel is presented.     

基于矢量位和标量位的空间波数混合域电磁三维正演模拟

高效、高精度电磁三维数值模拟是制约大规模电磁数据精细化三维反演成像、人机交互定量解释的核心问题.针对一问题,本文提出一种基于矢量位和标量位的空间波数混合域电磁场三维数值模拟方法.该方法利用沿水平方向的二维傅里叶变换将空间域矢量位和标量位耦合偏微分方程组转换为波数之间相互独立的常微分耦合方程组,将一个大规模三维问题分解为多个一维小问题,具有高度并行性,由此大大减少了计算量和存储量;保留垂向为空间域,浅层网格剖分适当加密,深层网格剖分适当稀疏,有效兼顾了计算精度与计算效率;采用有限单元法求解不同波数的常微分方程,充分利用追赶法求解定带宽线性方程组的高效性进一步提高数值模拟效率.在模型算例中,设计棱柱体模型验证了本文方法的正确性、计算精度和计算效率.数值试验结果表明本文方法具有数值精度高、并行度高、占用内存小、计算效率高的特性,比传统有限单元法三维数值模拟方法计算效率高1~2个数量级,且网格剖分规模越大,该方法计算效率优势越明显.

     

A discontinuous finite element baroclinic marine model on unstructured prismatic meshes

We describe the space discretization of a three-dimensional baroclinic finite element model, based upon a discontinuous Galerkin method, while the companion paper (Comblen et al. 2010a) describes the discretization in time. We solve the hydrostatic Boussinesq equations governing marine flows on a mesh made up of triangles extruded from the surface toward the seabed to obtain prismatic three-dimensional elements. Diffusion is implemented using the symmetric interior penalty method. The tracer equation is consistent with the continuity equation. A Lax–Friedrichs flux is used to take into account internal wave propagation. By way of illustration, a flow exhibiting internal waves in the lee of an isolated seamount on the sphere is simulated. This enables us to show the advantages of using an unstructured mesh, where the resolution is higher in areas where the flow varies rapidly in space, the mesh being coarser far from the region of interest. The solution exhibits the expected wave structure. Linear and quadratic shape functions are used, and the extension to higher-order discretization is straightforward.     

Evaluating Titan2D mass-flow model using the 1963 Little Tahoma Peak avalanches,Mount Rainier,Washington

The Titan2D geophysical mass-flow model is evaluated by comparing its simulation results and those obtained from another flow model, FLOW3D, with published data on the 1963 Little Tahoma Peak avalanches on Mount Rainier, Washington. The avalanches, totaling approximately 10×10⁶ m³ of broken lava blocks and other debris, traveled 6.8 km horizontally and fell 1.8 km vertically (H/L=0.246). Velocities calculated from runup range from 24 to 42 m/s and may have been as high as 130 m/s while the avalanches passed over Emmons Glacier.Titan2D is a code for an incompressible Coulomb continuum; it is a depth-averaged, ‘shallow-water’, granular-flow model. The conservation equations for mass and momentum are solved with a Coulomb-type friction term at the basal interface. The governing equations are solved on multiple processors using a parallel, adaptive mesh, Godunov scheme. Adaptive gridding dynamically concentrates computing power in regions of special interest; mesh refinement and coarsening key on the perimeter of the moving avalanche. The model flow initiates as a pile defined as an ellipsoid by a height (z) and an elliptical base defined by radii in the x and y planes. Flow parameters are the internal friction angle and bed friction angle. Results from the model are similar in terms of velocity history, lateral spreading, location of runup areas, and final distribution of the Little Tahoma Peak deposit. The avalanches passed over the Emmons Glacier along their upper flow paths, but lower in the valley they traversed stream gravels and glacial outwash deposits. This presents difficulty in assigning an appropriate bed friction angle for the entire deposit. Incorporation of variable bed friction angles into the model using GIS will help to resolve this issue.     

基于二次场的海洋可控源电磁法三维有限元正演

根据库伦规范势的定义,推导出关于磁矢量势和电标量势的偏微分方程,为了克服由电流源引起的奇异性和数值模拟计算困难,将电磁总场分解为一次场和二次场,一次场由基于Schelkunoff势函数的一维正演算法得到,二次场由有限元法计算得到,实现了海洋可控源电磁法三维有限元正演算法。通过一维数值模拟实例,验证该算法的计算精度。然后,利用该算法对带海底地形的三层储层模型进行正演,分析了海底地形对海洋控源电磁场各分量产生的影响。     

Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction

A finite-element method for computing the electric field in a 3-D conductivity model of the Earth for plane wave sources, thus enabling magnetotelluric responses to be calculated, is presented. The method incorporates in the iterative solution of the electric-field system of equations the divergence correction technique introduced for finite-difference solutions by Smith (1996). The correction technique accelerates the development of the discontinuity of the normal component of the approximate electric field across conductivity discontinuities. The convergence rate of the iterative solution is improved significantly, especially for low frequencies. The correction technique involves computing the divergence of the current density for the approximate electric field, computing the static potential whose source is this divergence of the current density, and ‘correcting’ the approximate electric field by subtracting from it the gradient of the potential. This is repeated at regular intervals during the iterative solution of the electric-field system of equations. For the method presented here, the Earth model is discretised using a rectilinear mesh comprising uniform cells. Edge-element basis functions are used to approximate the electric field and nodal basis functions are used to approximate the correction potential. The Galerkin method is used to derive the systems of equations for the approximate electric field and correction potential from the respective differential equations. A bi-conjugate gradient solver was found to be adequate for the system of equations for the correction potential; a generalised minimum residual solver was found to be better for the electric-field system of equations. The method is illustrated using the COMMEMI 3D-1A and 3D-2A models.     

The GeoClaw software for depth-averaged flows with adaptive refinement

Many geophysical flow or wave propagation problems can be modeled with two-dimensional depth-averaged equations, of which the shallow water equations are the simplest example. We describe the GeoClaw software that has been designed to solve problems of this nature, consisting of open source Fortran programs together with Python tools for the user interface and flow visualization. This software uses high-resolution shock-capturing finite volume methods on logically rectangular grids, including latitude-longitude grids on the sphere. Dry states are handled automatically to model inundation. The code incorporates adaptive mesh refinement to allow the efficient solution of large-scale geophysical problems. Examples are given illustrating its use for modeling tsunamis and dam-break flooding problems. Documentation and download information is available at www.clawpack.org/geoclaw.     

Probabilistic modeling of aquifer heterogeneity using reliability methods

A probabilistic model of groundwater contaminant transport is presented. The model is based on coupling first- and second-order reliability methods (FORM and SORM) with a two-dimensional finite element solution of groundwater transport equations. Uncertainty in aquifer media is considered by modeling hydraulic conductivity as a spatial random field with a prescribed correlation structure. FORM and SORM provide the probability that a contaminant exceeds a target level at a well, termed the probability of failure. Sensitivity of the probability of failure to basic variabilities in grid block conductivity is also obtained, at no additional computational effort. The effect of the choice of the predetermined target level at the observation well is provided, along with its effect on the relevant sensitivity information. Considerable saving in computational time was achieved by superimposing a coarse random variables mesh on a finer numerical mesh. The presence of regions of lower conductivity on the probabilistic event is analyzed, and the regions in which conductivity most affects the results are identified.     

有限元模拟弹性位错的等效体力方法

计算地震位错造成的位移场和应力场,对于估计大地震引起的后续地震活动发展趋势十分重要.本文提出在弹性位错问题的有限元模拟中,用等效体力代替位错源,从而在构建几何模型时不用包含断层,却可以处理包含任意复杂断层的问题,极大降低建模的难度.使用此方法,本文计算并讨论了在球形地球模型下2011年日本Tohoku-Oki特大地震对华北地区断层的影响,结果表明此次地震使华北主要断层趋于稳定.     

Sensitivities of an adjoint, unstructured mesh, tidal model on the European Continental Shelf

Unstructured mesh models can resolve the model domain with a variable and very fine mesh resolution. Nevertheless, tuning the model setup is still required (for example because of parametrized sub-grid processes). Adjoint models are commonly used to calculate sensitivities of ocean models and optimize their parameters so that better agreement is achieved between model simulations and observations. One major obstacle in developing an adjoint model is the need to update the reverse code after each modification of the forward code, which is not always straightforward. Automatic differentiation is a tool to generate the adjoint model code without user input. So far this method has mainly been used for structured mesh ocean models. We present here an unstructured mesh, adjoint, tidal model using this technique, and discuss the sensitivities of the misfit between simulated and observed elevations with respect to open boundary values, the bottom friction coefficient and the bottom topography. The forward model simulates tides on the European Continental Shelf and we show that the tidal model dynamics in the adjoint simulations can be used to define regions, where parameters or mesh has to be optimized. We analyze the dependence of the sensitivities on the wave type and mesh resolution to specify whether the model misfit originates from physical or numerical model deficiencies. In the sensitivity patterns, it is possible to identify islands not resolved in the mesh. We propose to refine the mesh prior to the parameter optimization.     

Combination of FEM and CMA‐ES algorithm for transmissivity identification in aquifer systems

In this paper, we propose a coupling of a finite element model with a metaheuristic optimization algorithm for solving the inverse problem in groundwater flow (Darcy's equations). This coupling performed in 2 phases is based on the combination of 2 codes: This is the HySubF‐FEM code (hydrodynamic of subsurface flow by finite element method) used for the first phase allowing the calculation of the flow and the CMA‐ES code (covariance matrix adaptation evolution strategy) adopted in the second phase for the optimization process. The combination of these 2 codes was implemented to identify the transmissivity field of groundwater by knowing the hydraulic head in some point of the studied domain. The integrated optimization algorithm HySubF‐FEM/CMA‐ES has been validated successfully on a schematic case offering an analytical solution. As realistic application, the integrated optimization algorithm HySubF‐FEM/CMA‐ES was applied to a complex groundwater in the north of France to identify the transmissivity field. This application does not use zonation techniques but solves an optimization problem at each internal node of the mesh. The obtained results are considered excellent with high accuracy and fully consistent with the hydrogeological characteristics of the studied aquifer.However, the various numerical simulations performed in this paper have shown that the CMA‐ES algorithm is time‐consuming. Finally, the paper concludes that the proposed algorithm can be considered as an efficient tool for solving inverse problems in groundwater flow.     

Three‐dimensional modelling of controlled‐source electromagnetic surveys using an edge finite‐element method with a direct solver

A fully three‐dimensional finite‐element algorithm has been developed for simulating controlled‐source electromagnetic surveys. To exploit the advantages of geometric flexibility, frequency‐domain Maxwell's equations of the secondary electric field were discretised using edge‐based finite elements while the primary field was calculated analytically for a horizontally layered‐earth model. The resulting system of equations for the secondary field was solved using a parallel version of direct solvers. The accuracy of the algorithm was successfully verified by comparisons with integral‐equations and iterative solutions, and the applicability to models containing large conductivity contrasts was verified against published data. The advantages of geometry‐conforming meshes have been demonstrated by comparing different mesh systems to simulate an inclined sheet model. A comparison of the performance between direct and iterative solvers demonstrated the superior efficiency of direct solvers, particularly for multisource problems.     

Modelling tide–surge interaction effects using finite volume and finite element models of the Irish Sea

An unstructured mesh model of the west coast of Britain, covering the same domain and using topography and open boundary forcing that are identical to a previous validated uniform grid finite difference model of the region, is used to compare the performance of a finite volume (FV) and a finite element (FE) model of the area in determining tide–surge interaction in the region. Initial calculations show that although qualitatively both models give comparable tidal solutions in the region, comparison with observations shows that the FV model tends to under-estimate tidal amplitudes and hence background tidal friction in the eastern Irish Sea. Storm surge elevations in the eastern Irish Sea due to westerly, northerly and southerly uniform wind stresses computed with the FV model tend to be slightly higher than those computed with the FE model, due to differences in background tidal friction. However, both models showed comparable non-linear tide–surge interaction effects for all wind directions, suggesting that they can reproduce the extensive tide–surge interaction processes that occur in the eastern Irish Sea. Following on from this model comparison study, the physical processes contributing to surge generation and tide–surge interaction in the region are examined. Calculations are performed with uniform wind stresses from a range of directions, and the balance of various terms in the hydrodynamic equations is examined. A detailed comparison of the spatial variability of time series of non-linear bottom friction and non-linear momentum advection terms at six adjacent nodes at two locations in water depths of 20 and 6 m showed some spatial variability from one node to another. This suggests that even in the near coastal region, where water depths are of the order of 6 m and the mesh is fine (of order 0.5 km), there is significant spatial variability in the non-linear terms. In addition, distributions of maximum bed stress due to tides and wind forcing in nearshore regions show appreciable spatial variability. This suggests that intensive measurement campaigns and very high-resolution mesh models are required to validate and reproduce the non-linear processes that occur in these regions and to predict extreme bed stresses that can give rise to sediment movement. High-resolution meshes will also be required in pollution transport problems.     

Dry Granular Flow Modelling Including Erosion and Deposition

In this paper, we considered dense avalanches made of dry and cohesion-less snow as gravitational granular flows. We used an effective friction law derived from the recent progress on granular behaviour research to represent momentum loss. We formulated a simple erosion and deposition model where the eroded and deposited mass fluxes are estimated from the knowledge of the starting and stopping conditions on a given slope angle. We adopted the shallow water hypothesis and introduced the new formulation of erosion and deposition processes. We solved the equations thanks to an accurate finite volumes solver scheme, using a non-structured mesh. This allows an easy integration of defence structures and singularities. Each component of this model was tested using analytical solutions and experimental laboratory data. Finally, we used the February 1999 Taconnaz avalanche event to test this model on a real snow avalanche case. These tests gave promising results.     

复杂固-液介质的速度应力方程三角形单元间断伽辽金地震波模拟算法

间断伽辽金法可用于存在复杂边界条件的模型,同时具有高阶精度和易于并行计算的优点,因此近年来在地震波传播模拟研究中得到了快速发展.数值通量是间断伽辽金法的关键组成部分之一.相比于其他通量,基于Rankine-Hugoniot跳跃条件的通量（RH-condition通量）在固体和液体介质边界具有更宽的稳定性,可以使用更大的时间步长,尤其是液体和固体之间的波阻抗差异较大时.它已被用于基于四边形网格和速度-应变方程的间断伽辽金地震波固液介质模拟.本文为了模拟陆地自然水体等存在复杂固-液界面形状的地震波传播,并在未来与其他数值方法耦合,发展了基于三角形网格的RH-condition通量间断伽辽金方法,使用一阶速度-压力声波和一阶速度-应力弹性波方程模拟复杂固-液介质中地震波的传播.通过水平层状和sin型起伏固-液模型,验证了该方法模拟结果的准确性.通过数值模拟展示了所提出的间断伽辽金法在不同网格大小和阶数下的准确性和效率.最后通过一个复杂模型的例子表明存在复杂固-液界面时该方法可以准确施加固-液边界条件.