Quadrangle-grid velocity–stress finite difference method for poroelastic wave equations

A quadrangle-grid velocity–stress finite difference method, based on a first-order hyperbolic system that is equivalent to Biot's equations, is developed for the simulation of wave propagation in 2-D heterogeneous porous media. In this method the velocity components of the solid material and of the pore fluid relative to that of the solid, and the stress components of three solid stresses and one fluid pressure are defined at different nodes for a staggered non-rectangular grid. The scheme uses non-orthogonal grids, allowing surface topography and curved interfaces to be easily modelled in the numerical simulation of seismic responses of poroelastic reservoirs. The free-surface conditions of complex geometry are achieved by using integral equilibrium equations on the surface, and the source implementations are simple. The algorithm is an extension of the quadrangle-grid finite difference method used for elastic wave equations.     

The construction of effective methods for electromagnetic modelling

Summary. This paper deals with the further development of finite-difference methods for electromagnetic field modelling in two-and three-dimensional cases. The main feature of the approach suggested here is the application of generalized asymptotic boundary conditions valid with the accuracy (1/ρ^N), where ρ is the distance from the heterogeneities. The finite-difference approximation of problems under solution is made using the balance method, which results in 5-point difference schemes in the 2-D case and 7-point difference schemes in the 3-D case. To solve the linear system of difference equations the successive over-relaxation (SOR) method is used, the relaxation factor being chosen during the iteration procedure. In view of the vectorial character of the problem for the 3-D case, a successive blocked over-relaxation method (SBOR) is applied.
The model's validity is based on the comparison of the fields accounted at the ground surface with those computed by the integral transformation of excessive currents, determined in the heterogeneity region using the finite-difference scheme.     

Surface motion of a fluid planet induced by impacts

In order to approximate the free-surface motion of an Earth-sized planet subjected to a giant impact, we have described the excitation of body and surface waves in a spherical compressible fluid planet without gravity or intrinsic material attenuation for a buried explosion source. Using the mode summation method, we obtained an analytical solution for the surface motion of the fluid planet in terms of an infinite series involving the products of spherical Bessel functions and Legendre polynomials. We established a closed form expression for the mode summation excitation coefficient for a spherical buried explosion source, and then calculated the surface motion for different spherical explosion source radii a (for cases of   a / R = 0.001  to 0.035, R is the radius of the Earth) We also studied the effect of placing the explosion source at different radii r ₀ (for cases of   r ₀/ R = 0.90  to 0.96) from the centre of the planet. The amplitude of the quasi-surface waves depends substantially on a / R , and slightly on   r ₀/ R   . For example, in our base-line case,   a / R = 0.03, r ₀/ R = 0.96  , the free-surface velocity above the source is 0.26 c , whereas antipodal to the source, the peak free surface velocity is 0.19 c . Here c is the acoustic velocity of the fluid planet. These results can then be applied to studies of atmosphere erosion via blow-off caused by asteroid impacts.     

Traction image method for irregular free surface boundaries in finite difference seismic wave simulation

In this study, we propose a new numerical method, named as Traction Image method, to accurately and efficiently implement the traction-free boundary conditions in finite difference simulation in the presence of surface topography. In this algorithm, the computational domain is discretized by boundary-conforming grids, in which the irregular surface is transformed into a 'flat' surface in computational space. Thus, the artefact of staircase approximation to arbitrarily irregular surface can be avoided. Such boundary-conforming gridding is equivalent to a curvilinear coordinate system, in which the first-order partial differential velocity-stress equations are numerically updated by an optimized high-order non-staggered finite difference scheme, that is, DRP/opt MacCormack scheme. To satisfy the free surface boundary conditions, we extend the Stress Image method for planar surface to Traction Image method for arbitrarily irregular surface by antisymmetrically setting the values of normal traction on the grid points above the free surface. This Traction Image method can be efficiently implemented. To validate this new method, we perform numerical tests to several complex models by comparing our results with those computed by other independent accurate methods. Although some of the testing examples have extremely sloped topography, all tested results show an excellent agreement between our results and those from the reference solutions, confirming the validity of our method for modelling seismic waves in the heterogeneous media with arbitrary shape topography. Numerical tests also demonstrate the efficiency of this method. We find about 10 grid points per shortest wavelength is enough to maintain the global accuracy of the simulation. Although the current study is for 2-D P-SV problem, it can be easily extended to 3-D problem.     

Smooth inversion for ground surface temperature histories: estimating the optimum regularization parameter by generalized cross-validation

The inversion of recent borehole temperatures has proved to be a successful tool to determine ancient ground surface temperature histories. To take into account heterogeneity of thermal properties and their non-linear dependence on temperature itself, a versatile 1-D inversion technique based on a finite-difference approach has been developed. Regularization of the generally ill-posed problem is obtained by an appropriate version of Tikhonov regularization of variable order. In this approach, a regularization parameter has to be determined, representing a trade-off between data fit and model smoothness. We propose to select this parameter by generalized cross-validation. The resulting technique is employed in case studies from the Kola ultradeep drilling site, and another borehole from northeastern Poland. Comparing the results from both sites corroborates the hypothesis that subglacial ground surface temperatures as met in Kola often are much higher than the ones in areas exposed to atmospheric conditions (Poland).     

Modelling of scattering of seismic waves from a corrugated rough sea surface: a comparison of three methods

We compare three numerical methods to model the sea surface interaction in a marine seismic reflection experiment (the frequencies considered are in the band 10–100 Hz): the finite-difference method (FDM), the spectral element method (SEM) and the Kirchhoff method (KM). A plane wave is incident at angles of 0° and 30° with respect to the vertical on a rough Pierson–Moskowitz surface with 2 m significant wave height and the response is synthesized at 6, 10 and 50 m below the average height of the sea surface. All three methods display an excellent agreement for the main reflected arrival. The FDM and SEM also agree very well all through the scattered coda. The KM shows some discrepancies, particularly in terms of amplitudes.     

2-Dreflectivity method and synthetic seismograms for irregularly layered structures -I. SH-wave generation

Summary. The reflectivity method for complete SH seismograms has been extended to two-dimensionally layered structures. The Aki-Larner technique is generalized to solve the integral equations for 2-D boundary conditions, and propagator matrices are enlarged to express a total SH wavefield. Synthetic seismograms in a soft basin are calculated for an incident plane-wave. They compare favourably with the results of the finite-element and finite-difference methods even in the later portion where asymptotic ray and beam theories break down. Synthetic seismograms due to a line force and a point dislocation are also presented.     

Diffraction of P, S and Rayleigh waves by three-dimensional topographies

The diffraction of P, S and Rayleigh waves by 3-D topographies in an elastic half-space is studied using a simplified indirect boundary element method (IBEM). This technique is based on the integral representation of the diffracted elastic fields in terms of single-layer boundary sources. It can be seen as a numerical realization of Huygens principle because diffracted waves are constructed at the boundaries from where they are radiated by means of boundary sources. A Fredholm integral equation of the second kind for such sources is obtained from the stress-free boundary conditions. A simplified discretization scheme for the numerical and analytical integration of the exact Green's functions, which employs circles of various sizes to cover most of the boundary surface, is used.
The incidence of elastic waves on 3-D topographical profiles is studied. We analyse the displacement amplitudes in the frequency, space and time domains. The results show that the vertical walls of a cylindrical cavity are strong diffractors producing emission of energy in all directions. In the case of a mountain and incident P, SV and SH waves the results show a great variability of the surface ground motion. These spatial variations are due to the interference between locally generated diffracted waves. A polarization analysis of the surface displacement at different locations shows that the diffracted waves are mostly surface and creeping waves.     

Numerical modelling of propagation of elastic waves in anisotropic inhomogeneous media for the half-space and the sphere

Summary. A method based on a combination of partial separation of variables and finite-difference method is used for the calculation of complete theoretical seismograms for inhomogeneous anisotropic media. Examples of theoretical seismograms for several anisotropic models are presented.     

沙漠下垫面地面温度及能量收支的BATS模式参数化改进

利用陆面过程模式BATS,引入地表发射率及两种大气发射率参数化方案,同时引入不同的地表粗糙度参数化方案,对比各种参数化方案对沙漠下垫面地面温度及能量收支的模拟状况。结果表明：采用Van Bavel等发展的地表发射率及Chung等发展的大气发射率方案可以明显改进地面温度及向上长波辐射的模拟,Chung等方案在夜间与正午的模拟效果更好,减小了1 ℃左右的地面温度模拟偏差,减小了10 W·m^-2左右的向上长波辐射模拟偏差。晴天地面温度及向上长波辐射的模拟结果优于阴天。利用Zhang等发展的裸土粗糙度参数化方案也会提高模式对地表感热通量模拟的准确性。     

Split-step complex Padé-Fourier depth migration

We present a split-step complex Padé-Fourier migration method based on the one-way wave equation. The downward-continuation operator is split into two downward-continuation operators: one operator is a phase-shift operator and the other operator is a finite-difference operator. A complex treatment of the propagation operator is applied to mitigate inaccuracies and instabilities due to evanescent waves. It produces high-quality images of complex structures with fewer numerical artefacts than those obtained using a real approximation of a square-root operator in the one-way wave equation. Tests on zero-offset data from the SEG/EAGE salt data show that the method improves the image quality at the cost of an additional 10 per cent computational time compared to the conventional Fourier finite-difference method.     

Shear-wave polarizations on a curved wavefront at an isotropic free surface

Summary. We present polarization diagrams of the particle motions at the free surface of an isotropic half-space generated by incident shear waves from a local buried point source. The reflectivity technique is used to calculate synthetic seismograms from which the particle motions are plotted. The particle motions are examined over a range of epicentral distances in a uniform isotropic half-space for different source frequencies and polarization angles, and for different Poisson's ratios. The particle motions due to a curved wavefront possess different characteristics from those generated by plane wavefronts at corresponding incidence angles. A curved wavefront generates a local SP -phase: a P -headwave which propagates along the free surface, and arrives shortly before the direct S -wave. These two arrivals give rise to cruciform particle motions in the sagittal and horizontal planes, which could be misinterpreted as anisotropy-induced shear-wave splitting. An examination of the particle motion in the transverse plane, mutually orthogonal to the sagittal and horizontal planes, can be used to discriminate between isotropic and anisotropic interpretations. The amplitude of the SP -phase is enhanced when it propagates in a low-velocity surface layer overlying the source layer, and may then become the dominant phase on radial-component seismograms. The presence of even a single surface layer may introduce considerable complexity into the seismogram, and we examine the effects of layer thickness, velocity contrast, and source depth on the corresponding polarization diagrams. Reliable information on the source and propagation path characteristics of shear waves from a buried local point source can only be obtained from free-surface records if they are recorded within a very limited epicentral distance range.     

Non-linear and finite-amplitude thermal convection in a heterogeneous terrestrial planet

Summary. Equations governing non-linear and finite-amplitude convection in a heterogeneous planetary interior are developed. Using spherical harmonic expressions of variables, together with Green's function of Laplacian operator in a spherical coordinate, the equations are reduced to one-dimensional integro-differential equations and their numerical solutions are obtained by a finite-difference scheme. The theory is then applied to several lunar models and the following conclusions are obtained.
(1) The mean temperatures and velocities of convecting zones of variable viscosity models are higher than those of constant viscosity ones. This is due to the development of lithospheres with 400–500 km thicknesses in the former models, which reduce heat loss considerably.
(2) Molten regions are continuous shells in variable viscosity models whereas they become discontinuous and localized in a constant viscosity model. The continuous molten shells decrease lateral variations of temperature significantly and tend to stabilize convection.
(3) Lateral variations of viscosity have negligible effects on the thermal evolution of the models considered.     

Comparison of global synthetic seismograms calculated using the spherical 2.5-D finite-difference method with observed long-period waveforms including data from the intra-Antarctic region

We have been developing an accurate and efficient numerical scheme, which uses the finite-difference method (FDM) in spherical coordinates, for the computation of global seismic wave propagation through laterally heterogeneous realistic Earth models. In the field of global seismology, traditional axisymmetric modeling has been used widely as an efficient approach since it can solve the 3-D elastodynamic equation in spherical coordinates on a 2-D cross-section of the Earth, assuming structures to be invariant with respect to the axis through the seismic source. However, it has the severe disadvantages that asymmetric structures about the axis cannot be incorporated and the source mechanisms with arbitrary shear dislocation have not been attempted for a long time. Our scheme is based on the framework of axisymmetric modeling but has been extended to treat asymmetric structures, arbitrary moment-tensor point sources, anelastic attenuation, and the Earth center which is a singularity of wave equations in spherical coordinates. All these types of schemes which solve 3-D wavefields on a 2-D model cross-section are classified as 2.5-D modeling, so we have named our scheme the spherical 2.5-D FDM. In this study, we compare synthetic seismograms calculated using our FDM scheme with three-component observed long-period seismograms including data from stations newly installed in Antarctica in conjunction with the International Polar Year (IPY) 2007–2008. Seismic data from inland Antarctica are expected to reveal images of the Earth's deep interior with enhanced resolution because of the high signal-to-noise ratio and wide extent of this region, in addition to the rarity of sampling paths along the rotation axis of the Earth. We calculate synthetic seismograms through the preliminary reference earth model (PREM) including attenuation using a moment-tensor point source for the November 9, 2009 Fiji earthquake. Our results show quite good agreement between synthetic and observed seismograms, which indicates the accuracy of observations in the Antarctica, as well as the feasibility of the spherical 2.5-D modeling scheme.     

Effect of a simple mountain range on underground seismic motion

Summary. The influence of a simple mountain range on seismic ground motion is studied. A two-dimensional model of the medium and vertically incident plane SH -waves are considered. Attention is devoted not only to the wavefield along the Earth's surface, but also within the medium. The wavefield is computed in two steps: (1) the computation of the impulse response by the finite-difference method, (2) the computation of the response to a time variation of the incident wave. approximately corresponding to a 'hundred-year' local earthquake at the site. Numerical results (the impulse response, the transfer function, the accelerograms, and their spectra) indicate strong spatial variability of the wavefield due to the topographical anomaly. The differences as large as 100 per cent in the peak amplitudes of the accelerograms, and of the order of 5 Hz in their predominant frequencies, at closely (∼ 200 m) spaced internal points of the medium have been found. Attention is focused also on the effect of causal absorption. Even unrealistically strong absorption seems to be unable to reduce significantly the spatial variability of the ground motion, caused by the topography. A variability like this. implying the occurrence of underground differential strains, might be of engineering importance in the antiseismic design of underground structures (tunnels, for example) in mountainous seismic regions. The ground-motion variability along the surface of the mountain is considerably smaller than within the medium.     

Surface wave tomography: global membrane waves and adjoint methods

We implement the wave equation on a spherical membrane, with a finite-difference algorithm that accounts for finite-frequency effects in the smooth-Earth approximation, and use the resulting 'membrane waves' as an analogue for surface wave propagation in the Earth. In this formulation, we derive fully numerical 2-D sensitivity kernels for phase anomaly measurements, and employ them in a preliminary tomographic application. To speed up the computation of kernels, so that it is practical to formulate the inverse problem also with respect to a laterally heterogeneous starting model, we calculate them via the adjoint method, based on backpropagation, and parallelize our software on a Linux cluster. Our method is a step forward from ray theory, as it surpasses the inherent infinite-frequency approximation. It differs from analytical Born theory in that it does not involve a far-field approximation, and accounts, in principle, for non-linear effects like multiple scattering and wave front healing. It is much cheaper than the more accurate, fully 3-D numerical solution of the Earth's equations of motion, which has not yet been applied to large-scale tomography. Our tomographic results and trade-off analysis are compatible with those found in the ray- and analytical-Born-theory approaches.     

地表能量平衡方程在新疆地表温度模拟中的应用

地表温度是影响地表能量收支平衡和气候变化的重要因子,也是评价陆面模式性能的主要指标。新疆下垫面包括沙漠、绿洲、冰川、积雪等多种类型,其地表温度变化对生态环境和天气气候变化有重要影响。为了找到影响第三代陆面模式NCAR CLM4.5模拟新疆地表温度效果的主要因子,促进该模式发展,在对ERA-5再分析资料评估基础上,以其作为实况资料,首先对CLM4.5长时段模拟新疆地表温度的效果开展评估,然后利用地表能量平衡方程得到的地表温度变化影响因子分离方法,揭示出引起新疆地表温度模拟偏差的主要因子。结果表明：ERA-5地表温度的空间分布特征、量级和时间变化与观测结果一致,可用来研究该区域地表温度变化。CLM4.5对新疆地表温度模拟冷偏差值较大,但相关性显著,且误差与气候背景和下垫面条件有关。CLM4.5模拟的潜热值偏大是地表温度模拟值偏低的主要原因,CLM4.5模拟的潜热值在新疆大部分区域比ERA-5的多90%以上。CLM4.5对新疆0—10 cm土壤湿度模拟值偏大,蒸散过强,导致潜热通量偏大,因此改进土壤湿度计算方案是提高该陆面模式对新疆地表温度模拟水平的有效途径。     

Acceleration of buoyancy-driven fractures and magmatic dikes beneath the free surface

By performing buoyancy-driven fracture experiments in brittle gelatin we observe that the ascent velocity of a fracture containing a finite volume of fluid increases when approaching the free surface. We theoretically describe this free-surface effect and quantify it by introducing an effective depth-dependent fracture toughness and developing an ascent model on the basis of linear fracture mechanics. We develop a successful inversion approach and resolve the actual and critical fracture length and the ascent velocity far away from the free surface from the observation of the fracture tip migration alone. Other parameters, as the fluid volume included in the fracture and the in situ fracture toughness, can be derived. Applying the model and inversion to the 1998 eruption at Piton de la Fournaise, Reunion Island, reveals estimates of the length and critical length of the feeding magma batch, the magma batch volume and the in situ fracture toughness. It further indicates that the ascent velocity of the magma batch was probably much smaller at greater depths and that the batch might have been initiated several months or years before the eruption.     

基于地理特征的CBERS-02 CCD非正射影像镶嵌研究

非正射影像镶嵌效果与控制点的布局及镶嵌线的选取密切相关。基于地理特征分析,提出一种针对非正射影像的镶嵌方案:将几何精纠正的控制点与影像镶嵌线布设在主体高程面内。遵循此原则,能降低利用推扫式中心投影成像原理获取地形起伏区遥感影像投影差对影像镶嵌的影响。对太原市2006年3月一组CBERS-02 CCD影像镶嵌的实验表明,该方法能够有效削弱甚至消除非正射影像镶嵌中的几何拼接缝。     

Winter estimation of surface roughness length over eastern Qinghai-Tibetan Plateau

Based on the Monin-Obukhov similarity theory, a scheme was developed to calculate surface roughness length. Surface roughness length over the eastern Qinghai-Tibetan Plateau during the winter season was then estimated using the scheme and eddy covariance measurement data. Comparisons of estimated and measured wind speeds show that the scheme is feasible to calculate surface roughness length. The estimated roughness lengths at the measurement site during unfrozen,frozen and melted periods are 3.23×10-3, 2.27×10-3 and 1.92×10-3 m, respectively. Surface roughness length demonstrates a deceasing trend with time during the winter season. Thereby, setting the roughness length to be a constant value in numerical models could lead to certain degree of simulation errors. The variation of surface roughness length may be caused by the change in land surface characteristic.