An accurate time integration method for simplified overland flow models

An accurate time integration method for the diffusion-wave and kinematic-wave approximated models for the overland flow is proposed. The discretization of the first- and second-order spatial derivatives in the basic equation is obtained by using the second-order Lax–Wendroff and the three-point centred finite difference schemes, respectively. For the solution in time, the system of ordinary differential equations, obtained by the linearization of the celerity and of the hydraulic diffusivity by Taylor series expansions, is integrated analytically. The stability and the numerical dissipation and dispersion are investigated by the Fourier analysis. A proper Courant number, and the corresponding time step for the numerical simulations can be established. In addition, the proposed diffusion-wave and kinematic-wave models are straightforwardly extended to the two-dimensional flow. Test cases for both one- and two-dimensional problems, compare the solutions of the diffusion-wave and kinematic-wave models with analytical solutions, with experimental results and with numerical solutions obtained by the Saint–Venant equations. These simulations show that the proposed numerical–analytical models accurately predict the overland flow for several situations, in particular for unsteady rainfall rate and for spatial variations of the surface roughness.     

Impact of partial steps and momentum advection schemes in a global ocean circulation model at eddy-permitting resolution

Series of sensitivity tests were performed with a z-coordinate, global eddy-permitting (1/4°) ocean/sea-ice model (the ORCA-R025 model configuration developed for the DRAKKAR project) to carefully evaluate the impact of recent state-of-the-art numerical schemes on model solutions. The combination of an energy–enstrophy conserving (EEN) scheme for momentum advection with a partial step (PS) representation of the bottom topography yields significant improvements in the mean circulation. Well known biases in the representation of western boundary currents, such as in the Atlantic the detachment of the Gulf Stream, the path of the North Atlantic Current, the location of the Confluence, and the strength of the Zapiola Eddy in the south Atlantic, are partly corrected. Similar improvements are found in the Pacific, Indian, and Southern Oceans, and characteristics of the mean flow are generally much closer to observations. Comparisons with other state-of-the-art models show that the ORCA-R025 configuration generally performs better at similar resolution. In addition, the model solution is often comparable to solutions obtained at 1/6 or 1/10° resolution in some aspects concerning mean flow patterns and distribution of eddy kinetic energy. Although the reasons for these improvements are not analyzed in detail in this paper, evidence is shown that the combination of EEN with PS reduces numerical noise near the bottom, which is likely to affect current–topography interactions in a systematic way. We conclude that significant corrections of the mean biases presently seen in general circulation model solutions at eddy-permitting resolution can still be expected from the development of numerical methods, which represent an alternative to increasing resolution.     

Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation

We extend lattice Boltzmann (LB) methods to advection and anisotropic-dispersion equations (AADE). LB methods are advocated for the exactness of their conservation laws, the handling of different length and time scales for flow/transport problems, their locality and extreme simplicity. Their extension to anisotropic collision operators (L-model) and anisotropic equilibrium distributions (E-model) allows to apply them to generic diffusion forms. The AADE in a conventional form can be solved by the L-model. Based on a link-type collision operator, the L-model specifies the coefficients of the symmetric diffusion tensor as linear combination of its eigenvalue functions. For any type of collision operator, the E-model constructs the coefficients of the transformed diffusion tensors from linear combinations of the relevant equilibrium projections. The model is able to eliminate the second order tensor of its numerical diffusion. Both models rely on mass conserving equilibrium functions and may enhance the accuracy and stability of the isotropic convection–diffusion LB models.The link basis is introduced as an alternative to a polynomial collision basis. They coincide for one particular eigenvalue configuration, the two-relaxation-time (TRT) collision operator, suitable for both mass and momentum conservation laws. TRT operator is equivalent to the BGK collision in simplicity but the additional collision freedom relates it to multiple-relaxation-times (MRT) models. “Optimal convection” and “optimal diffusion” eigenvalue solutions for the TRT E-model allow to remove next order corrections to AADE. Numerical results confirm the Chapman–Enskog and dispersion analysis.     

A radiation boundary condition for finite-element models: 1-D linear advection test

Summary A linear 1-D advection equation is used to study the utilization of a finite-element method and open lateral boundary conditions. Two possible implementations of a radiation (Sommerfeld) boundary condition are tested for the case of a solitary wave passing through a computational domain — the case corresponding to zero lateral boundary values — and for the case of a simple sinusoidal wave which corresponds to a non-zero boundary forcing.     

An adaptive hybrid EnKF-OI scheme for efficient state-parameter estimation of reactive contaminant transport models

Reactive contaminant transport models are used by hydrologists to simulate and study the migration and fate of industrial waste in subsurface aquifers. Accurate transport modeling of such waste requires clear understanding of the system’s parameters, such as sorption and biodegradation. In this study, we present an efficient sequential data assimilation scheme that computes accurate estimates of aquifer contamination and spatially variable sorption coefficients. This assimilation scheme is based on a hybrid formulation of the ensemble Kalman filter (EnKF) and optimal interpolation (OI) in which solute concentration measurements are assimilated via a recursive dual estimation of sorption coefficients and contaminant state variables. This hybrid EnKF-OI scheme is used to mitigate background covariance limitations due to ensemble under-sampling and neglected model errors. Numerical experiments are conducted with a two-dimensional synthetic aquifer in which cobalt-60, a radioactive contaminant, is leached in a saturated heterogeneous clayey sandstone zone. Assimilation experiments are investigated under different settings and sources of model and observational errors. Simulation results demonstrate that the proposed hybrid EnKF-OI scheme successfully recovers both the contaminant and the sorption rate and reduces their uncertainties. Sensitivity analyses also suggest that the adaptive hybrid scheme remains effective with small ensembles, allowing to reduce the ensemble size by up to 80% with respect to the standard EnKF scheme.     

A modified TVD scheme for the advection of two or more variables with consideration for their sum

Total variation diminishing (TVD) advection schemes are known to produce results that are free from some of the numerical artifacts (no overshooting, no spurious oscillation, small diffusion) that can spoil the physical significance of the results. When two or more tracers are advected separately using a TVD scheme, the sum of these variables can however exhibit some inappropriate behaviors. The total variation of the sum will not necessarily be non- increasing and local artificial oscillations and extrema can appear. We show that these can be avoided with only minor perturbations of the original solution by adjusting the slope limiters used for the different variables. If the sum of these variables has some physical significance, for instance as refinement of a larger model compartment, the correction procedure introduced in this paper should be used to ensure a physically meaningful solution.     

Resolving frontal structures: on the payoff using a less diffusive but computationally more expensive advection scheme

This article presents some advantages using a shape-preserving total variation diminishing (TVD) advection scheme in an ecosystem model. The superbee flux-limiter has been used to the second-order Lax–Wendroff advection scheme to make it TVD. We performed simulations for three shelf sea regions with different characteristic time scales, namely, the North Sea, the Barents Sea, and the Baltic Sea. To explore the advantages, we also performed reference runs with the much simpler and computationally cheaper upwind advection scheme. Frontal structures are much better resolved with the TVD scheme. The bottom salinity in the Baltic Sea stays at realistic values throughout the 10-year simulation with the TVD scheme, while with the upwind scheme, it drifts towards lower values and the permanent haline stratification in the Baltic is almost completely eroded within one seasonal cycle. Only when applying TVD for both the vertical and horizontal advections the model succeeded to preserve haline stratification in the decadal simulation. Lower trophic level patterns are far better reproduced with the TVD scheme, and for the estimated cod larval survival, the advantages seem to be even stronger. Simulations using the TVD-derived prey fields identified distinct regions such as Dogger Bank to favor potential larvae survival (PLS), which did not appear as particularly favorable in the upstream simulations. The TVD scheme needs about 25 % more time on the central processing unit (CPU) in case of a pure hydrodynamic setup with only two scalar state variables (Barents Sea application). The additional CPU time cost increases for a coupled physical–biological model application with a large number of state variables. However, this is more than compensated by all the advantages found, and, hence, we conclude that it is worthwhile to use the TVD scheme in our ecosystem model.

     

Resolving frontal structures: on the payoff using a less diffusive but computationally more expensive advection scheme

This article presents some advantages using a shape-preserving total variation diminishing (TVD) advection scheme in an ecosystem model. The superbee flux-limiter has been used to the second-order Lax–Wendroff advection scheme to make it TVD. We performed simulations for three shelf sea regions with different characteristic time scales, namely, the North Sea, the Barents Sea, and the Baltic Sea. To explore the advantages, we also performed reference runs with the much simpler and computationally cheaper upwind advection scheme. Frontal structures are much better resolved with the TVD scheme. The bottom salinity in the Baltic Sea stays at realistic values throughout the 10-year simulation with the TVD scheme, while with the upwind scheme, it drifts towards lower values and the permanent haline stratification in the Baltic is almost completely eroded within one seasonal cycle. Only when applying TVD for both the vertical and horizontal advections the model succeeded to preserve haline stratification in the decadal simulation. Lower trophic level patterns are far better reproduced with the TVD scheme, and for the estimated cod larval survival, the advantages seem to be even stronger. Simulations using the TVD-derived prey fields identified distinct regions such as Dogger Bank to favor potential larvae survival (PLS), which did not appear as particularly favorable in the upstream simulations. The TVD scheme needs about 25?% more time on the central processing unit (CPU) in case of a pure hydrodynamic setup with only two scalar state variables (Barents Sea application). The additional CPU time cost increases for a coupled physical–biological model application with a large number of state variables. However, this is more than compensated by all the advantages found, and, hence, we conclude that it is worthwhile to use the TVD scheme in our ecosystem model.     

曲线坐标系下的块体运动与应变模型研究

本文建立了顾及地球扁率和局部切标架随点变化特性的椭球坐标系下的刚体运动模型和块体运动与应变模型,以及球坐标系下顾及局部切标架随点变化特性的严密的块体运动与应变模型,分析了球坐标系下块体运动与应变模型及椭球坐标系下的块体运动与应变模型间的差异;通过计算具体讨论了地球扁率和曲线坐标系的局部切标架随点变化特性对欧拉矢量与应变张量的影响.结果表明：地球扁率对刚体欧拉矢量和应变参数的影响甚小,具体计算时可以不予考虑,但曲线坐标系的局部切标架随点变化特性对两者的影响较大,在建模过程中需要顾及,常用的Savage模型需要修正.     

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.     

An optimized finite difference method based on a polar coordinate system for regional-scale irregular topography

The finite difference method (FDM) is an important numerical approach for simulating the propagation of seismic waves, and some FDMs can be used to study the impact of the Earth’s curvature and topography over large distances. To efficiently model the effects of the Earth’s irregular topography on the propagation of seismic waves, here we optimize a previously proposed grid mesh method and develop a novel two-dimensional boundary-conforming FDM based on a curvilinear polar coordinate system. This method efficiently simulates the propagation of seismic waves in an arc-shaped model with large variations in surface topography. Our method was benchmarked against other reported methods using several global-scale models. The consistency of the results confirms the validity of our proposed optimization strategy. Furthermore, our findings indicate that the proposed optimization strategy improves computational efficiency.     

A finite volume upwind scheme for the solution of the linear advection–diffusion equation with sharp gradients in multiple dimensions

A finite volume upwind numerical scheme for the solution of the linear advection equation in multiple dimensions on Cartesian grids is presented. The small-stencil, Modified Discontinuous Profile Method (MDPM) uses a sub-cell piecewise constant reconstruction and additional information at the cell interfaces, rather than a spatial extension of the stencil as in usual methods. This paper presents the MDPM profile reconstruction method in one dimension and its generalization and algorithm to two- and three-dimensional problems. The method is extended to the advection–diffusion equation in multiple dimensions. The MDPM is tested against the MUSCL scheme on two- and three-dimensional test cases. It is shown to give high-quality results for sharp gradients problems, although some scattering appears. For smooth gradients, extreme values are best preserved with the MDPM than with the MUSCL scheme, while the MDPM does not maintain the smoothness of the original shape as well as the MUSCL scheme. However the MDPM is proved to be more efficient on coarse grids in terms of error and CPU time, while on fine grids the MUSCL scheme provides a better accuracy at a lower CPU.     

An upwind fast sweeping scheme for calculating seismic wave first‐arrival travel times for models with an irregular free surface

The topography‐dependent eikonal equation formulated in a curvilinear coordinate system has recently been established and revealed as being effective in calculating first‐arrival travel times of seismic waves in an Earth model with an irregular free surface. The Lax–Friedrichs sweeping scheme, widely used in previous studies as for approximating the topography‐dependent eikonal equation viscosity solutions, is more dissipative and needs a much higher number of iterations to converge. Furthermore, the required number of iterations grows with the grid refinement and results in heavy computation in dense grids, which hampers the application of the Lax–Friedrichs sweeping scheme to seismic wave travel‐time calculation and high‐resolution imaging. In this paper, we introduce a new upwind fast sweeping solver by discretising the Legendre transform of the numerical Hamiltonian of the topography‐dependent eikonal equation using an explicit formula. The minimisation related to the Legendre transform in the sweeping scheme is solved analytically, which proved to be much more efficient than the Lax–Friedrichs algorithm in solving the topography‐dependent eikonal equation. Several numerical experiments demonstrate that the new upwind fast sweeping method converges and achieves much better accuracy after a finite number of iterations, independently of the mesh size, which makes it an efficient and robust tool for calculating travel times in the presence of a non‐flat free surface.     

An accurate and robust multigrid algorithm for 2D forward resistivity modelling

We present an adaptation of the full multigrid algorithm in DC resistivity modelling in an effort to increase its accuracy. There is a great difficulty with conventional multigrid solvers in representing the physics of an arbitrary distribution of electrical conductivity on a very coarse grid. In general, conventional rectangular finite‐difference or 5‐point approximations of Poisson's equation cannot represent, at a coarse grid level, the effective anisotropy on a coarse scale which results from fine structure in the model. An exception to this generalization occurs when the principal axes of structural anisotropy are aligned with the coordinate axis. Additional and similarly generated problems arise when a coarse cell is obliged to represent fine structure containing very high conductivity contrasts. We have developed an adaptation of the usual resistive‐network representation of the continuum, which avoids some of these problems, and have compared it with the traditional resistive network currently used. The network adaptation consists of replacing the usual 5‐point Laplacian operator stencil used on the finite‐difference grid with a 9‐point stencil, and the conductivity scalar with a 6‐parameter conductivity parametrization. This parametrization permits representation of arbitrarily orientated anisotropy as well as more complex behaviour related to high conductivity contrasts. The importance of multigrid solvers does not lie in their speed at forward modelling (which is comparable with other methods), but rather in their potential for inverse modelling. Inverse solvers which proceed by refinement of an initially very coarse solution can, in principle, take time only linearly proportional to the number of gridpoints in the final desired model.     

Analytical solutions for benchmarking cold regions subsurface water flow and energy transport models: One-dimensional soil thaw with conduction and advection

Numerous cold regions water flow and energy transport models have emerged in recent years. Dissimilarities often exist in their mathematical formulations and/or numerical solution techniques, but few analytical solutions exist for benchmarking flow and energy transport models that include pore water phase change. This paper presents a detailed derivation of the Lunardini solution, an approximate analytical solution for predicting soil thawing subject to conduction, advection, and phase change. Fifteen thawing scenarios are examined by considering differences in porosity, surface temperature, Darcy velocity, and initial temperature. The accuracy of the Lunardini solution is shown to be proportional to the Stefan number. The analytical solution results obtained for soil thawing scenarios with water flow and advection are compared to those obtained from the finite element model SUTRA. Three problems, two involving the Lunardini solution and one involving the classic Neumann solution, are recommended as standard benchmarks for future model development and testing.     

A 4D variational assimilation scheme with partition method for nearshore wave models

This paper summarizes the development steps of a 4D variational assimilation scheme for nearshore wave models. A partition method is applied for adjusting both wave boundary conditions and wind fields. Nonstationary conditions are assimilated by providing defined correlations of model inputs in time. The scheme is implemented into the SWAN model. Twin experiments covering both stationary and nonstationary wave conditions are carried out to assess the adequacy of the proposed scheme. Stationary experiments are carried out considering separately windsea, swells, and mixed sea. Cost functions decline to less than 5% and RMS spectrum errors are reduced to less than 10%. The nonstationary experiment covers 1 day simulation under mixed wave conditions with assimilation windows of 3 h. RMS spectrum errors are reduced to less than 10% after 30 iterations in most assimilation windows. The results show that for spacially uniform model inputs, model accuracy is improved notably by the assimilation scheme throughout the computational domain. It is found that under wave conditions in which observed spectra can be well classified, the assimilation scheme is able to improve model results significantly.     

水平界面上P-SV转换波转换点的精确解

转换波共转换点的叠加和道集选取都需要准确地计算转换点的位置.Tessmer和Behle、Taylor分别给出了水平单层介质中P-SV转换波转换点坐标的解析解,由于其表达式的复杂性,在应用中几乎不被采用.在本文中,运用Snell定律重新建立了在水平反射界面上反射的P-SV转换波的转换点坐标的四次方程,并严格地推导出与纵波速度、横波速度、炮检距和反射深度有关的转换点坐标的解析解,确定了惟一的解析表达式.将这一结果应用于P-SV转换波的速度分析和叠加处理中.简化的公式有较好的应用价值.     

两斜交柱坐标系间谐波函数坐标变换的辅助平面方法.

针对利用波函数展开法进行三维地震响应研究中的坐标变换问题,提出了在两斜交柱坐标系间谐波函数空间坐标变换的辅助平面方法．通过建立与空间点对应的一系列辅助平面,将地震波三维散射研究中的谐波函数表达式从柱坐标系变换至辅助平面内的极坐标系,然后在辅助平面内运用Graf加法公式,将柱坐标系下的波函数表达式变换至与该柱坐标系斜交的另一柱坐标系下,从而将三维问题转换成二维问题进行处理,最终得到了斜交柱坐标系下内域问题和外域问题的坐标变换公式．