Tube wave modeling by the finite-difference method with varying grid spacing

A finite-difference approach of aP-SV modeling scheme is applied to compute seismic wave propagation in heterogeneous isotropic media, including fluid-filled boreholes. The discrete formulation of the equation of motion requires the definition of the material parameters at the grid points of the numerical mesh. The grid spacing is chosen as coarse as possible with respect to the accurate representation of the shortest wavelength. If we assume frequencies lower than 250 Hz then the grid spacing is usually chosen in the range of a few meters. One encounters difficulties because of the large-scale difference between the grid spacing and the size of the borehole, usually several centimeters.These difficulties can be overcome by a grid refinement technique. This technique provides the construction of grids with varying grid spacing. The grid spacing in the vicinity of the borehole is chosen such that the borehole is properly represented. An example demonstrates the accuracy of this technique by comparisons with other methods. Unlike many analytical methods, the FD method can handle complex subsurface geometries. Further numerical examples of walk-awayVSP configurations show tube wave propagation within fluid-filled boreholes of realistic diameters.     

Simulation of dam-break flow with grid adaptation

The unsteady free surface flow caused by sudden collapse of a dam produces discontinuities in the flow variables. As the flow surges downstream, it forms a moving bore front with steep gradients of water height and velocity. In the numerical simulation of this flow, proper grid distribution can play a crucial part in the prediction and resolution of the solutions. The use of presently available numerical schemes to solve this problem on a uniform course grid system fails to resolve the characteristic flow features and hence do a poor job in simulating this flow. In this paper, an adaptive grid which adjusts itself as the solution evolves is used for a better resolution of the flow properties. Rai and Anderson's¹² method is used to determine the grid speed; however, a different partial differential equation based on the conservative principle of grid arc lengths for clustering grids in one-dimensional flow is used along with the St. Venant equations to numerically simulate the flow. Both the subcritical and the supercritical flows under extreme boundary conditions are solved using this technique. With a specified number of grid points, this provides better quality solutions as compared to those obtained with uniformly distributed grids.     

Using the pseudospectral technique on curved grids for 2D acoustic forward modelling1

The Fourier pseudospectral method has been widely accepted for seismic forward modelling because of its high accuracy compared to other numerical techniques. Conventionally, the modelling is performed on Cartesian grids. This means that curved interfaces are represented in a ‘staircase fashion‘causing spurious diffractions. It is the aim of this work to eliminate these non-physical diffractions by using curved grids that generally follow the interfaces. A further advantage of using curved grids is that the local grid density can be adjusted according to the velocity of the individual layers, i.e. the overall grid density is not restricted by the lowest velocity in the subsurface. This means that considerable savings in computer storage can be obtained and thus larger computational models can be handled. One of the major problems in using the curved grid approach has been the generation of a suitable grid that fits all the interfaces. However, as a new approach, we adopt techniques originally developed for computational fluid dynamics (CFD) applications. This allows us to put the curved grid technique into a general framework, enabling the grid to follow all interfaces. In principle, a separate grid is generated for each geological layer, patching the grid lines across the interfaces to obtain a globally continuous grid (the so-called multiblock strategy). The curved grid is taken to constitute a generalised curvilinear coordinate system, where each grid line corresponds to a constant value of one of the curvilinear coordinates. That means that the forward modelling equations have to be written in curvilinear coordinates, resulting in additional terms in the equations. However, the subsurface geometry is much simpler in the curvilinear space. The advantages of the curved grid technique are demonstrated for the 2D acoustic wave equation. This includes a verification of the method against an analytic reference solution for wedge diffraction and a comparison with the pseudospectral method on Cartesian grids. The results demonstrate that high accuracies are obtained with few grid points and without extra computational costs as compared with Cartesian methods.     

贴体网格各向异性对坐标变换法求解起伏地表下地震初至波走时的影响

笛卡尔坐标系中的经典程函方程在静校正、叠前偏移、走时反演、地震定位、层析成像等很多地球物理工作中都有应用,然而用其计算起伏地表的地震波走时却比较困难.本文通过把曲线坐标系中的矩形网格映射到笛卡尔坐标系的贴体网格,推导出曲线坐标中的程函方程,而后,用Lax-Friedrichs快速扫描算法求解曲线坐标系的程函方程.研究表明本文方法能有效处理地表起伏的情况,得到准确稳定的计算结果.由于地表起伏,导致与之拟合的贴体网格在空间上的展布呈各向异性,且这种各向异性的强弱对坐标变换法求解地震初至波的走时具有重要影响.本文研究表明,随着贴体网格的各向异性增强,用坐标变换法求解地表起伏区域的走时计算误差增大,且计算效率降低,这在实际应用具有指导意义.     

用于地震波场模拟的变网格边长声格固体模型

提出一种新的用于地震波场模拟的变网格边长声格固体模型（Phononiclatticesolidwithvariousgridlength,简称PLSVL模型）．该模型通过改变不同介质中网格的边长来体现介质速度的变化,粒子在所有网格中的运动均保持每一步运行一个网格边长．详细推导了该模型的Boltzmann方程,证明该方程及其散射项的表达式与PLS模型相同,因而从该方程出发,可寻出宏观变量所满足的波动方程．新模型在同种网格内部的传输过程不存在误差,在界面产生的误差不会随时间而积累．本文同时给出了一个理论模型的波场模拟结果．     

Fast transport simulation with an adaptive grid refinement

One of the main difficulties in transport modeling and calibration is the extraordinarily long computing times necessary for simulation runs. Improved execution time is a prerequisite for calibration in transport modeling. In this paper we investigate the problem of code acceleration using an adaptive grid refinement, neglecting subdomains, and devising a method by which the Courant condition can be ignored while maintaining accurate solutions. Grid refinement is based on dividing selected cells into regular subcells and including the balance equations of subcells in the equation system. The connection of coarse and refined cells satisfies the mass balance with an interpolation scheme that is implicitly included in the equation system. The refined subdomain can move with the average transport velocity of the subdomain. Very small time steps are required on a fine or a refined grid, because of the combined effect of the Courant and Peclet conditions. Therefore, we have developed a special upwind technique in small grid cells with high velocities (velocity suppression). We have neglected grid subdomains with very small concentration gradients (zero suppression). The resulting software, MODCALIF, is a three-dimensional, modularly constructed FORTRAN code. For convenience, the package names used by the well-known MODFLOW and MT3D computer programs are adopted, and the same input file structure and format is used, but the program presented here is separate and independent. Also, MODCALIF includes algorithms for variable density modeling and model calibration. The method is tested by comparison with an analytical solution, and illustrated by means of a two-dimensional theoretical example and three-dimensional simulations of the variable-density Cape Cod and SALTPOOL experiments. Crossing from fine to coarse grid produces numerical dispersion when the whole subdomain of interest is refined; however, we show that accurate solutions can be obtained using a fraction of the execution time required by uniformly fine-grid solutions.     

Solution of the transport equations using a moving coordinate system

A convection-diffusion equation arises from the conservation equations in miscible and immiscible flooding, thermal recovery, and water movement through desiccated soil. When the convection term dominates the diffusion term, the equations are very difficult to solve numerically. Owing to the hyperbolic character assumed for dominating convection, inaccurate, oscillating solutions result. A new solution technique minimizes the oscillations. The differential equation is transformed into a moving coordinate system which eliminates the convection term but makes the boundary location change in time. We illustrate the new method on two one-dimensional problems: the linear convection-diffusion equation and a non-linear diffusion type equation governing water movement through desiccated soil. Transforming the linear convection diffusion equation into a moving coordinate system gives a diffusion equation with time dependent boundary conditions. We apply orthogonal collocation on finite elements with a Crank-Nicholson time discretization. Comparisons are made to schemes using fixed coordinate systems. The equation describing movement of water in dry soil is a highly non-linear diffusion-type equation with coefficients varying over six orders of magnitude. We solve the equation in a coordinate system moving with a time-dependent velocity, which is determined by the location of the largest gradient of the solution. The finite difference technique with a variable grid size is applied, and a modified Crank-Nicholson technique is used for the temporal discretization. Comparisons are made to an exact solution obtained by similarity transformation, and with an ordinary finite difference scheme on a fixed coordinate system.     

Optimization of Cell Parameterizations for Tomographic Inverse Problems

—?We develop algorithms for the construction of irregular cell (block) models for parameterization of tomographic inverse problems. The forward problem is defined on a regular basic grid of non-overlapping cells. The basic cells are used as building blocks for construction of non-overlapping irregular cells. The construction algorithms are not computationally intensive and not particularly complex, and, in general, allow for grid optimization where cell size is determined from scalar functions, e.g., measures of model sampling or a priori estimates of model resolution. The link between a particular cell j in the regular basic grid and its host cell k in the irregular grid is provided by a pointer array which implicitly defines the irregular cell model. The complex geometrical aspects of irregular cell models are not needed in the forward or in the inverse problem. The matrix system of tomographic equations is computed once on the regular basic cell model. After grid construction, the basic matrix equation is mapped using the pointer array on a new matrix equation in which the model vector relates directly to cells in the irregular model. Next, the mapped system can be solved on the irregular grid. This approach avoids forward computation on the complex geometry of irregular grids. Generally, grid optimization can aim at reducing the number of model parameters in volumes poorly sampled by the data while elsewhere retaining the power to resolve the smallest scales warranted by the data. Unnecessary overparameterization of the model space can be avoided and grid construction can aim at improving the conditioning of the inverse problem. We present simple theory and optimization algorithms in the context of seismic tomography and apply the methods to Rayleigh-wave group velocity inversion and global travel-time tomography.     

Elastic modelling on a grid with vertically varying spacing1

We present a discrete modelling scheme which solves the elastic wave equation on a grid with vertically varying grid spacings. Spatial derivatives are computed by finite-difference operators on a staggered grid. The time integration is performed by the rapid expansion method. The use of variable grid spacings adds flexibility and improves the efficiency since different spatial sampling intervals can be used in regions with different material properties. In the case of large velocity contrasts, the use of a non-uniform grid avoids spatial oversampling in regions with high velocities. The modelling scheme allows accurate modelling up to a spatial sampling rate of approximately 2.5 gridpoints per shortest wavelength. However, due to the staggering of the material parameters, a smoothing of the material parameters has to be applied at internal interfaces aligned with the numerical grid to avoid amplitude errors and timing inaccuracies. The best results are obtained by smoothing based on slowness averaging. To reduce errors in the implementation of the free-surface boundary condition introduced by the staggering of the stress components, we reduce the grid spacing in the vertical direction in the vicinity of the free surface to approximately 10 gridpoints per shortest wavelength. Using this technique we obtain accurate results for surface waves in transversely isotropic media.     

复杂地下异常体的可控源电磁法积分方程正演

可控源电磁法具有分辨率高及抗干扰能力强等特点,是一种重要的地电磁勘探方法.目前,可控源电磁法的高精度正演计算一直是其核心研究问题之一.传统积分方程法一般采用近似积分公式、简单矩形网格和近似的奇异性体积分计算技术,制约了体积分方程法处理复杂地下异常体的能力,降低了计算精度.针对上述问题,本文基于完全积分公式、四面体非结构化网格和奇异体积分的精确解析解来高精度求解复杂可控源电磁模型的正演响应.首先,从电场积分公式出发,推导了可控源电磁问题满足的积分方程;其次,借助于非结构化四面体网格离散技术,实现了地下复杂异常体的有效模拟.最后,利用散度定理把强奇异值体积分转换为一系列弱奇异性的面积分公式,并通过推导获得了这些弱奇异性的面积分公式的解析解,从而最终实现三维可控源电磁问题的高精度积分求解.以块状低阻体地电模型为测试模型,采用本文提出的积分方程方法获得的数值解与其他公开数值算法解进行对比分析,其对比结果具有高度的吻合性,验证了算法的正确性;同时,设计了球状及复杂地电模型进行算法收敛性测试,进一步验证算法的正确性以及能够处理地下复杂模型的能力.     

A Fourth Order Accurate SH-Wave Staggered Grid Finite-difference Algorithm with Variable Grid Size and VGR-Stress Imaging Technique

This article presents a new approach for the implementation of a planar-free surface boundary condition. It is based on a vertical grid-size reduction above the free surface during the explicit computation of a free surface boundary condition. This technique is very much similar to the well-known stress imaging technique. VGR-stress imaging technique name is proposed for this new free surface boundary condition (VGR stands for ‘vertical grid-size reduction’). To study the performance of the proposed VGR-stress imaging technique, it was implemented in a newly developed second order accurate in time and fourth-order accurate in space (2, 4) staggered grid SH-wave finite-difference (FD) algorithm with variable grid size. It was confirmed that the effective thickness (ETH) of first soil layer becomes less by one-half of vertical grid size than the assigned thickness (ATH), if stress imaging technique is used as a free surface boundary condition. The qualitative and quantitative results of various numerical experiments revealed that the proposed VGR-stress imaging technique is better than the stress imaging technique since it is free from the thickness discrepancy arising due to the use of images of stress components across the free surface. On the basis of iterative numerical experiments, it was confirmed that the stability condition for this FD scheme with variable grid size is It was also inferred that at least five to six grid points per shortest wavelength are required to avoid the grid dispersion. The maximum grid-spacing ratio up to 12.5 or even more did not affect the accuracy of (2,4) SH-wave algorithm. The obtained reduction of 10.46 and 5.38 folds in the requirement of computational memory and time for a particular basin-edge model, as compared with the homogeneous grid size, reflects the efficacy of the new FD algorithm.     

Parameterization of Heat and Moisture Exchange on Land Surface in Areas with Moderate Continental Climate

A technique is proposed for parameterization of the process of heat and moisture exchange between the underlying surface and the atmosphere at a regional scale under moderate continental climate, which implies a long cold period with stable snow cover and seasonally frozen soil. The technique allows the heterogeneity of the area under consideration to be explicitly accounted for in the model through properly selecting the model grid. The details at a scale below the grid size, such as the presence of field or forest ecosystems in any grid cell, also can be incorporated in the model. The proposed technique has been tested against unique observational data along with published scientific generalization for the Polomet River basin (Valdai Hills), which is 1215 km² in area.     

ACCURACY AND STABILITY OF FINITE DIFFERENCE EQUATIONS FOR SHALLOW WATER FLOWS IN ORTHOGOPNAL COORDINATES

IINTRODUCTIONNumericalmethodsasatooltosimulateflowsandpollutanttransportareincreasinglyimportantinhydraulicandenvironmentalengineering.AveryusefulapplicationofthenumericalmethodologyinengineeringproblemswouldbetosolvethesystemofZDdepth-integratedshallowwaterequations.ManysolutionsofthegoverningequationsarederivedusingtraditionalfinitedifferencemethodonCartesianregulargrids.ThedisadvantageofthismethodseemstobetheinflexibilityofCartesiangridstocomplywithirregularorcurvedperimeterswhichsur…     

三维泊松方程数值模拟的多重网格方法

本文简要介绍了多重网格方法的基本思想和原理,然后应用多重网格(MG)方法求解三维泊松方程,网格尺度从17×17×17逐次增加至257×257×257,并与不完全Chelesky共轭梯度法(ICCG)、Gauss直接解法进行比较,结果表明,MG方法计算速度明显优于ICCG、Gauss方法,对于129×129×129网格的三维数值模拟费时43s,比ICCG法快7倍,而对于257×257×257超大型网格的三维数值模拟也仅需412s.     

用于声波方程数值模拟的时间-空间域有限差分系数确定新方法

声波方程数值模拟已广泛应用于理论地震计算,同时构成了地震逆时偏移成像技术的基础.对于有限差分法而言,在满足一定的稳定性条件时,普遍存在着因网格化而形成的数值频散效应.如何有效地缓解或压制数值频散是有限差分方法研究的关键所在.为精确求解空间偏导数,相继发展了高阶差分格式优化方法和伪谱方法.近期,为更好地缓解数值频散,提出了时间-空间域有限差分方法,该方法采用了泰勒展开近似方法来确定有限差分格式系数,因而只能保证在一定的小范围内很好的拟合波场传播规律.为进一步压制数值频散效应,本文引入了时间-空间域特定波数点满足频散关系的方法,根据震源、波速和网格间距确定波数范围,同时考虑了多个传播角度,然后建立方程确定了相应的有限差分格式系数,使得差分系数能在更大范围符合波场传播规律.通过频散分析和正演模拟,验证了本文方法的有效性.     

The unsteady plane flow of ice-sheets: A parabolic problem with two moving boundaries

Abstract

Finite difference algorithms have been developed to solve a one-dimensional non-linear parabolic equation with one or two moving boundaries and to analyse the unsteady plane flow of ice-sheets. They are designed to investigate the response of an ice-sheet to changes in climate, and to reconstruct climatic changes implied by past ice-sheet variations inferred from glacial geological data. Two algorithms are presented and compared. The first, a fixed domain method, replaces time as an independent variable with span. The grid interval in real space is kept constant, and thus the number of grid points changes with span. The second, a moving mesh method, retains time as one of the independent variables, but normalises the spatial variable relative to the span, which now enters the diffusion and advection coeficients in the parabolic equation for the surface profile.

Crank-Nicholson schemes for the solution of the equations are constructed, and iterative schemes for the solution of the resulting non-linear equations are considered.

Boundary (margin) motion is governed by the surface slope at the margin. Differentiation of the evolution equations results in an evolution equation for the margin slopes. It is shown that incorporation of this evolution equation, while not formally increasing the accuracy of the finite difference schemes, in practice increases accuracy of the solution.     

Treatment of time derivative and calculation of flow when solving groundwater flow problems by Galerkin finite element methods

Two techniques connected with the use of the finite element Galerkin method for solving the linear parabolic differential equation describing unsteady groundwater flow in an anisotropic non-homogeneous aquifer are introduced. The first is a mode superposition technique for dealing with the time derivative which involves computing the smallest eigenvalues and associated eigenvectors of the matrices arising from the Galerkin method. It is shown how such a technique allows us to interpret the response of the groundwater level to input in terms of parallel linear reservoirs. It is further argued that if properly implemented, the technique will have computational advantages over standard finite difference methods, e.g. in the case when the input function is constant over relatively large time subintervals. The second is a technique based on so-called generalized flow formulae for calculating flow values across external or internal boundaries, posterior to obtaining the groundwater level values. The implementation of the technique in the case of linear triangular elements on an irregular grid is discussed. It is finally argued from simplified cases that, apart from guaranteeing a match with prescribed input, the technique may often be expected to give more accurate flow values than those obtained directly from the groundwater gradients.     

波场模拟中的数值频散分析与校正策略

波动方程有限差分法正演模拟,对认识地震波传播规律、进行地震属性研究、地震资料地质解释、储层评价等,均具有重要的理论和实际意义.但有限差分法本身固有存在着数值频散问题,数值频散在正演模拟中是一种严重的干扰,会降低波场模拟的精度与分辨率.针对TI介质波场模拟的交错网格有限差分方法,本文从空间网格离散、时间网格离散和算子近似等三个方面对其产生的数值频散进行了分析,并结合其他学者的研究成果给出了TI介质波场模拟中压制数值频散的方法与策略：在已知介质频散关系时,对差分算子可实施算子校正;通过提高差分方程的阶数来提高波场模拟精度;采用流体力学中守恒式方程的通量校正传输方法来压制波场模拟中的数值频散;在实际正演模拟时,采用交错网格高阶有限差分方程,不仅在空间上采用高阶差分,而且在时间上也要采用高阶差分,否则只在单一方向上（空间或时间）提高方程的阶数对压制数值频散也不会取得理想的效果.     

Finite difference calculation of traveltime on non-orthogonal grid

Finite difference methods have been widely employed in solving the eikonal equation so as to calculate traveltime of seismic phase. Most previous studies used regular orthogonal grid. However, much denser grid is required to sample the interfaces that are undulating in depth direction, such as the Moho and the 660 km discontinuity.Here we propose a new finite difference algorithm to solve the eikonal equation on non-orthogonal grid(irregular grid).To demonstrate its efficiency and accuracy, a test was conducted with a two-layer model. The test result suggests that the similar accuracy of a regular grid with ten times grids could achieve with our new algorithm, but the time cost is only about 0.1 times. A spherical earth model with an undulant660 km discontinuity was constructed to demonstrate the potential application of our new method. In that case, the traveltime curve fluctuation corresponds to topography. Our new algorithm is efficient in solving the first arrival times of waves associated with undulant interfaces.