Correction to: Inertial instability and phase error in Euler forward predictorcorrector time integration schemes: Improvement of modeling Great Lakes thermal structure and circulation using FVCOM

This study investigates the inertial stability properties and phase error of numerical time integration schemes in several widely-used ocean and atmospheric models. These schemes include the most widely used centered differencing (i.e., leapfrog scheme or the 3-time step scheme at n-1, n, n+1) and 2-time step (n, n+1) 1^st-order Euler forward schemes, as well as 2^nd-stage and 3^rd- and 4^th-stage Euler predictor-corrector (PC) schemes. Previous work has proved that the leapfrog scheme is neutrally stable with respect to the Coriolis force, with perfect inertial motion preservation, an amplification factor (AF) equal to unity, and a minor overestimation of the phase speed. The 1^st-order Euler forward scheme, on the other hand, is known to be unconditionally inertially unstable since its AF is always greater than unity. In this study, it is shown that 3^rd- and 4^th-order predictor-corrector schemes 1) are inertially stable with weak damping if the Coriolis terms are equally split to n+1 (new value) and n (old value); and 2) introduce an artificial computational mode. The inevitable phase error associated with the Coriolis parameter is analyzed in depth for all numerical schemes. Some schemes (leapfrog and 2^nd-stage PC schemes) overestimate the phase speed, while the others (1^st-order Euler forward, 3^rd- and 4^th-stage PC schemes) underestimate it. To preserve phase speed as best as possible in a numerical model, alternating a scheme that overestimates the phase speed with a scheme that underestimates the phase speed is recommended. Considering all properties investigated, the leapfrog scheme is still highly recommended for a time integration scheme. As an example, a comparison between a leapfrog scheme and a 1^st-order Euler forward scheme is presented to show that the leapfrog scheme reproduces much better vertical thermal stratification and circulation in the weakly-stratified Great Lakes.

     

On a general fluid dynamical theory of discrete unstable spiral modes in disk-shaped galaxies

Abstract

A general fluid dynamical theory of discrete unstable spiral modes in disk-shaped galaxies is described. This formulation of modes includes a radiation boundary condition and an exact numerical treatment of the Poisson equation. Thus, the modes are maintained by an outward transport of angular momentum, but they may be composed of both leading and trailing waves. A numerical scheme based on this formulation is described, and examples of modes obtained with this scheme are presented. These examples compare favorably with calculations based on the original asymptotic theory of Bertin, Lau, Lin and Mark. The implications of the present formulation of modes in galactic models support the hypothesis of a quasi-stationary spiral structure.     

Improvements and Corrections to AT123D Code

Although based on exact analytical solutions, semi‐analytical solute transport models can have significant numerical error in applications with high frequency oscillatory source terms and when parameter value combinations cause series solution approximations to converge slowly. Methods for correcting these numerical errors are presented and implemented in the AT123D code, which employs Green's functions to represent point, linear, and rectangular prismatic source zones. In order to increase its computational accuracy, a Romberg numerical integration scheme was added to AT123D with prespecified error criteria, variable time stepping, and partitioning of the integral to handle rapidly changing source terms. More rapidly converging series solution approximations for the Green's functions were also incorporated to improve both accuracy and computational efficiency for finite‐depth aquifers. AT123D also has been modified to eliminate redundant calculations at points where approximate steady‐state conditions have been reached to improve computational efficiency during numerical integration. These modifications help to decrease computer run times that can be excessive for three‐dimensional problems with large numbers of computational points, small time steps, and/or long simulation time periods. Errors in the original AT123D code also were corrected in this modified version, AT123D‐AT, in order to accurately simulate finite‐duration (pulse) source releases.     

Validation and Application of the Second Generation Three Dimensional Hydrodynamic Model of Chesapeake Bay

The validation and subsequent application of the current three-dimensional numerical hydrodynamic model of Chesapeake Bay is presented. The numerical model solves conservation equations for water mass, momentum, salinity, and heat on a boundary-fitted grid in the horizontal plane and a Cartesian z-grid in the vertical. A generalized ADI finite difference scheme is employed in conjunction with mode splitting technique, solving external and the internal modes. The 10-year boundary conditions including tide, slinity, temperature, wind, heat exchange coefficient, river and non-point source flows were constructed. Model validation was accomplished by demonstrating the model's ability to reproduce observed data over time scales ranging from tidal to seasonal periods. The major parameters compared include tidal elevation, intra-tidal and residual velocities, salinity, temperature, stratification, and flux calculated through the Bay mouth.After validation, the model was applied to simulate bay hydrodynamics for the 10 years of 1985–94. These results were used to drive the three-dimensional water quality model of Chesapeake Bay, which is discussed in a companion paper.     

On the computation of the barotropic mode of a free-surface world ocean model

The free-surface formulation of the equations of our world ocean model is briefly described. The barotropic mode equations are solved according to the split-explicit method, using different time steps for the external and internal modes. Because the numerical algorithm is implemented on the B-grid, a spurious, free-surface, two-grid interval mode may develop. This mode must be filtered out. The properties of two filters are theoretically investigated and their actual performance is tested in a series of numerical experiments. It is seen that one of these filters may severely perturb the local mass conservation, rendering it impossible to enforce the impermeability of the surface or the bottom of the ocean. The dynamics of the external mode is also examined, by studying the depth-integrated momentum equations. The depth-integral of the pressure force due to the slope of the ocean surface is approximately balanced by the depth-integral of the force ensuing from the horizontal variations of the density. The depth-integral of the Coriolis force is an order of magnitude smaller, except in the Southern Ocean. Two variational principles are resorted to for computing the fictitious ocean surface elevation corresponding to the approximate equilibrium between the dominant forces of the barotropic momentum equations.     

A split operator approach to reactive transport with the forward particle tracking Eulerian Lagrangian localized adjoint method

A forward particle tracking Eulerian Lagrangian localized adjoint method (ELLAM) is applied to the multicomponent reactive transport problem using a split operator approach. Two split operator algorithms are compared, the Strang algorithm and the sequential non-iterative algorithm (SNIA). The reaction equations are integrated using a coupled predictor corrector algorithm with adaptive time stepping. Reaction time steps are adjusted at the inflow boundary to reflect the actual time of transport inside the solution domain.Results show that split operator ELLAM formulations are competitive with direct or fully coupled ELLAM solutions for reactive transport problems. The SNIA algorithm is more accurate than the Strang splitting algorithm when large time steps are used. The reaction algorithm employed dominates computational effort in runs with large time step sizes. To illustrate the use of the method in practical problems, the model is fitted to aerobic aniline degradation data from laboratory scale column experiments. Model inversion is achieved using non-linear regression with a shuffled complex evolution optimization algorithm and parameter uncertainty is assessed using a Bayesian uncertainty analysis procedure.     

地壳剪切波分裂研究及地震定位误差对剪切波分裂分析的影响

地壳介质的各向异性主要是由大量充满液体的微裂隙的定向排列引起的, 剪切波穿过这种含有微裂隙的介质时会产生分裂。 剪切波分裂参数与地壳介质的物理性质有关, 与区域应力场有关。 在地壳剪切波分裂研究中, 地震的定位误差对剪切波分裂参数有影响。 本文简单介绍了几种影响定位精度的因素; 并对定位误差特别是深度方向的误差所引起的时间延迟变化进行了讨论, 比较了不同大小的深度误差对归一化后的时间延迟的影响, 定量分析结果显示, 深度定位误差对慢剪切波时间延迟影响很大。 剪切波分裂分析方法有多种, 本文还对不同的剪切波分裂分析方法做了简介。     

Adaptive time stepping and error control in a mass conservative numerical solution of the mixed form of Richards equation

Adaptive time stepping with embedded error control is applied to the mixed form of Richards equation. It is the first mathematically based adaptive scheme applied to this form of Richards equation. The key to the method is the approximation of the local truncation error of the scheme in terms of the pressure head, although, to enforce mass conservation, the principal time approximation is based on the moisture content. The time stepping scheme is closely related to an implicit Thomas–Gladwell approximation and is unconditionally stable and second-order accurate. Numerical trials demonstrate that the new algorithm fully automates stepsize selection and robustly constrains temporal discretisation errors given a user tolerance. The adaptive mechanism is shown to improve the performance of the non-linear solver, providing accurate initial solution estimates for the iterative process. Furthermore, the stepsize variation patterns reflect the adequacy of the spatial discretisation, here accomplished by linear finite elements. When sufficiently dense spatial grids are used, the time step varies smoothly, while excessively coarse grids induce stepsize oscillations.     

Comparative error analysis in finite element formulations of the advection-dispersion equation

The behaviour of numerical solutions of the one-dimensional advection-dispersion equation is investigated and comparisons between the consistent and the lumped formulations of Galerkin finite element schemes are made. Well-known criteria for the control of accuracy in the lumped (finite difference) formulation are reviewed. It is found that, because the numerical error produced by the consistent formulation is generally less than that produced by the lumped formulation, these criteria can also be used for the control of numerical dispersion in the consistent formulation. However, because the error in both types of solutions decreases in time when the discretization is invariant, the criteria can be relaxed with advancing simulation time. For the consistent formulation it is found that beyond some initial time period, the numerical error depends only on the temporal discretization. This suggests that constant accuracy can be maintained throughout the simulation period while allowing the time step length to grow.     

数字化水位仪观测的苏门答腊大地震激发的地球球型自由振荡

利用LN-3型数字化水位仪的观测资料,采用功率谱密度估计方法,在没有对资料进行去固体潮处理的情况下,获得了2004年12月26日苏门答腊大地震激发的0S3～0S30基频球型自由振荡,并检测到2S1、1S2、1S3、2S5和2S6等5个谐频球型振荡,同时记录到一些明显的振型的谱线分裂现象。用水位观测资料提取了0S3、0S4和2S1等3个振型的谱线分裂周期,虽然3个振型均只得到2个单峰,但地球自由振荡谱线分裂现象的检测为研究地球的深内部结构提供了新的、不同类型的观测资料。研究表明,用数字化水位仪的观测资料可以提取到苏门答腊大地震激发的地球球型振荡,为地球自由振荡的研究提供了一种新的观测手段。     

Alternative split-operator approach for solving chemical reaction/groundwater transport models

Various schemes are available to solve coupled transport/reaction mathematical models, one of the most efficient and easy to apply being the two-step split-operator method in which the transport and reaction steps are performed separately. Operator splitting, however, does not solve exactly the fully coupled numerical model derived from the governing partial differential and algebraic equations describing the transport and reaction processes. An error, proportional to Δt (the time step used in the numerical solution) is introduced. Thus, small time steps must be used to ensure that accurate solutions result. An alternative scheme is presented, which iterates to the exact solution of the fully coupled numerical model. The new scheme enables accurate solutions to be calculated more efficiently than the two-step method, while maintaining separation of the transport and reaction steps in the calculations. As in the two-step method, the reaction calculations are performed node-wise throughout the computation grid. However, because the scheme relies on LU factorisation of the coefficient matrix in the transport equation solution, the reaction calculations must be performed in sequence, the sequence order being determined by the ordering of the nodes in the grid. Also, because LU factorisation is used, the scheme is limited to solute transport problems for which LU factorisation is a practical solution method.     

Time domain dynamic analysis of piles under impact loading

In this paper, time domain dynamic analysis of piles under impact loading is presented. For this purpose a hybrid boundary element technique is implemented. Linear beam column finite elements are used to model the piles and resulting governing equations are solved using an implicit integration scheme. The continuum is assumed to be elastic and an efficient step-by-step time integration scheme is implemented by using an approximate half space integral formulation. By enforcing displacement equilibrium conditions at each time step, a system of equations is generated which yields the solution. Results of this time domain formulation under linear material behavior are compared with Laplace domain results to validate the methods.     

Mesh type tradeoffs in 2D hydrodynamic modeling of flooding with a Godunov-based flow solver

The effect of mesh type on the accuracy and computational demands of a two-dimensional Godunov-type flood inundation model is critically examined. Cartesian grids, constrained and unconstrained triangular grids, constrained quadrilateral grids, and mixed meshes are considered, with and without local time stepping (LTS), to determine the approach that maximizes computational efficiency defined as accuracy relative to computational effort. A mixed-mesh numerical scheme is introduced so all grids are processed by the same solver. Analysis focuses on a wide range of dam-break type test cases, where Godunov-type flood models have proven very successful. Results show that different mesh types excel under different circumstances. Cartesian grids are 2–3 times more efficient with relatively simple terrain features such as rectilinear channels that call for a uniform grid resolution, while unstructured grids are about twice as efficient in complex domains with irregular terrain features that call for localized refinements. The superior efficiency of locally refined, unstructured grids in complex terrain is attributable to LTS; the locally refined unstructured grid becomes less efficient using global time stepping. These results point to mesh-type tradeoffs that should be considered in flood modeling applications. A mixed mesh model formulation with LTS is recommended as a general purpose solver because the mesh type can be adapted to maximize computational efficiency.     

二维各向异性介质中地震波场的高阶同位网格有限差分模拟

本文将DRP/opt MacCormack有限差分格式用于模拟二维各向异性介质中的地震波传播.DRP/opt MacCormack是一种同位网格下的差分格式,避免了传统的交错网格在计算各向异性问题时由于变量插值而导致的误差.而且相对于低阶同位网格差分格式,它具有低色散、低耗散的优点.此格式将中心差分算子分成前向和后向两个空间单边差分,然后在4-6步Runge-Kutta时间积分中使用单边差分组合.在具有垂直对称轴的横向各向同性（VTI）模型下,通过对比DRP/opt MacCormack有限差分和谱元方法的模拟结果,验证了前者具有很高的精度和稳定性.由于实际地质条件下TI介质的对称轴通常是倾斜的(TTI),本文在二维三分量框架下模拟TTI介质中的地震波场.结果显示横波分裂和切平面/反平面运动耦合的特征.数值实验表明DRP/opt MacCormack是一种有效的研究各向异性介质中地震波传播规律的差分格式.     

Secondary flows induced by wind forcing in the Rhône region of freshwater influence

Secondary flows induced by the blocking effect of a river plume on a transverse upwelling are investigated in a microtidal region of freshwater influence (ROFI). A nested version of the SYMPHONIE primitive-equation free-surface model for 3-D baroclinic coastal flows has been developed for the Rhône ROFI. The main characteristics of the model are a generalized sigma coordinate system in finite differences, using a time splitting for external and internal modes and high-order numerical advection schemes for density fields in combination with an modified turbulence closure scheme. The nesting system consists of two grids forced by the high-resolution ALADIN model atmospheric data. The coarse grid of 3 km resolution for the whole Gulf of Lions allows the forcing of the Liguro-Provençal large-scale current when the fine mesh of 1-km resolution is centred on the river mouth of the Grand Rhône. Documented field experiments from the Biodypar 3 field campaign performed during March 1999 are used for validation. Numerical results, CTD profiles and a SPOT TSM visible image are in good agreement concerning the shape and structure of the river plume. Other coastal flow features can be observed from satellite imagery. Computations of realistic situations recover these main secondary structures. Complementary process-oriented runs give an explanation of how the coastal upwelling induced by an inhomogeneous offshore wind is destabilized by the combination of the river plume and along-shelf current-blocking effects. In the end, a factor-separation analysis provides evidence that the locally non-linear effects in momentum contribute to the occurrence of secondary vortices.Responsible Editor: Phil Dyke     

三角网格有限元法声波与弹性波模拟频散分析

本文对声波与弹性波方程进行有限元法离散,构造有限元法频散关系的一般特征值问题,分析了时间离散格式为中心差分的三角网格有限元法声波与弹性波模拟的频散特性. 比较了三种质量矩阵即分布式质量矩阵、集中质量矩阵和混合质量矩阵对有限元法频散的影响;选取四种典型三角网格,分析了混合质量矩阵有限元（MFEM）频散的方向各向异性;数值频散、方向各向异性随插值阶数的增加逐渐减弱,当空间为三阶插值时,频散主要表现为随采样率的变化而几乎无明显方向各向异性, 其频散幅值也较小. 控制其他影响因素不变的情况下,研究了不同波速比介质中弹性波的数值频散. 最后给出了三角网格MFEM的数值耗散性.     

Adaptive split-operator methods for modeling transport phenomena in porous medium systems

Split-operator methods are commonly used to approximate environmental models. These methods facilitate the tailoring of different approximation approaches to different portions of the differential operator and provide a means to split large coupled problems into pieces that are more amenable to parallel computation than the original fully-coupled problem. However, split-operator methods introduce an additional source of approximation error into the solution, which is typically either ignored or controlled heuristically. In this work, we develop two methods to estimate and control the error in split-operator methods, which lead to a dynamic adjustment of the temporal splitting step based upon the error estimators. The proposed methods are shown to yield robust solutions that provide the desired control of error. In addition, for a typical nonlinear reaction problem, the new methods are shown to reduce the solution error by more than two orders of magnitude compared to standard methods for an identical level of computational effort. The algorithms introduced and evaluated have widespread applicability in environmental modeling.     

Numerical modelling of stratified flow: a spectral approach

A numerical formulation is developed to solve the three-dimensional hydrodynamic equations which describe flow in a stratified sea.Arbitrary continuous physically realistic variations of density and eddy viscosity can be included in the model, which is sufficiently flexible to be applicable to sea areas of any horizontal extent and depth. A continuous current profile from sea surface to sea bed, is computed with the model. A method for expanding computed current profiles in terms of vertical modes is proposed and the contribution of these modes to the current profiles is considered.The time variation of the wind-induced circulation of a stratified lake in response to a suddenly imposed and maintained wind stress is examined. Calculations show that the wind-driven surface current is modulated by the internal seiche motion of the lake.     

Shallow water SPH for flooding with dynamic particle coalescing and splitting

In this paper an adaptive algorithm for Smoothed Particle Hydrodynamics (SPH) for the Shallow Water Equations (SWEs) is presented. The area of a particle is inversely proportional to depth giving poor resolution in small depths without particle refinement. This is a particular limitation for flooding problems of interest here. Higher resolution is created by splitting the particles, while particle coalescing (or merging) improves efficiency by reducing the number of the particles when acceptable. The new particle coalescing procedure merges two particles together if their area becomes less than a predefined threshold value. Both particle splitting and coalescing procedures conserve mass and momentum and the smoothing length of new particles is calculated by minimizing the density error of the SPH summation. The new dynamic particle refinement procedure is assessed by testing the numerical scheme against analytical, experimental and benchmark test cases. The analytical cases show that with particle splitting and coalescing typical convergence rates remain faster than linear. For the practical test case, in comparison to using particle splitting alone, the particle coalescing procedure leads to a significant reduction of computational time, by a factor of 15. This makes the computational time of the same order as mesh-based methods with the advantage of not having to specify a mesh over a flood domain of unknown extent a priori.     

Detection of spheriodal free oscillation excited by Peru 8.2 M_s earthquake with five international superconducting gravimeter data

In this paper, authors obtain the spectral peaks of the earth free oscillation and check all normal modes from 0S0 to 0S48 accurately, with the Fourier analysis and the maximum entropy spectrum method dealing jointly with six groups of the observational residual data from five international superconducting gravimeter stations. By comparing the observational results in this paper with three former groups of observations or models, authors notice an extra discrepancy between two observational 0S2 modes excited separately by Peru earthquake and Alaska earthquake, which probably mirrors the anisotropy of the Earth's inner core. The analysis on the splitting 1S2 mode shows that the asymmetric factor of rotationwise spectral splitting is possible to be different from that of anti-rotationwise spectral splitting.