Computation Of Seismograms Research Articles

Computing Theoretical Seismograms from a Point Source in a Spherical Multilayered Medium

Abstract Computing theoretical seismograms from a point source in a given Earth model is essential for modeling and inversion of observed seismic waveforms for Earth’s structure and earthquake source parameters. Here, we derived the propagator matrices and source terms for a spherical multilayered Earth model using the exact earth flattening transformation. We found that their differences from their counterparts in horizontal layered media are inversely proportional to the nondimensional horizontal wavenumber and its higher order. In addition, all the source terms in a spherical layered model have a source-depth dependent scaling factor that differs from in a horizontal layered model by up to 6% for deep earthquakes. The surface displacement produced by a point source can be obtained in a similar form as in horizontal layered media. Computation of theoretical seismograms was implemented using the generalized reflection and transmission coefficients method. Numerical tests show that our formulae and implementation are correct and efficient for computing full-wave seismograms, including the permanent displacements, at teleseismic distances up to 100°. Individual seismic phases can be isolated and analyzed semianalytically because the generalized reflection and transmission method is used. Furthermore, our analytic expression of displacement in terms of the propagator matrices and source terms can be used to derive analytic derivatives of seismograms for full-wave waveform inversion.

Interacting Seismic Waveguides: Multimode Surface Waves and Leaking Modes

Recent developments in seismic recording provide dense sampling of the seismic wavefield that allows the extraction of higher surface wave modes as well as the fundamental, and in some circumstances also P-dominated modes. The character of the modal branches and their dispersion with frequency and phase speed depends on interactions between different aspects of the seismic structure and between wavetypes. The influence of P waves becomes quite strong in the very near surface when there are strong wavespeed gradients. An effective tool for understanding the nature of the modal interactions due to structure is provided by the seismic response in frequency-phase speed space, the kernel for seismogram calculation. Such displays for all three-components of the surface response extract the full modal response, including leaking modes effects, for Rayleigh and Love waves and an indication of the way that the modes are excited.

Sliding speed calculation of the Wenchuan earthquake fault using the Doppler effect

At present, the fault sliding speed functions, or source time functions that used in theoretical seismogram calculation are all theoretical results. That is, the fault speed time functions are the results of the theoretical study, and they are not proved by the measured data. The method that used to calculate fault sliding speed by using Doppler effect has been greatly improved in this paper from the following four perspectives: Firstly, the paper proves theoretically how to confirm the seismic waves that received in some frequency bands by two receivers are emitted by the same source. Secondly, the paper puts forward the method to make sure whether seismic waves received by two receivers in some frequency bands are similar by using relative change of Fourier amplitude standard deviation in some frequency bands of two receivers, and similar seismic waves are emitted by the same frequency bands. Thirdly, to eliminate the interference of reflection and refraction waves, S wave records in fault sliding time are selected for data processing. Finally, long time Fourier transform is replaced by short time Fourier transform (STFT) to enhance fault sliding time positioning accuracy. On the basis of the work above, a general method to calculate the fault sliding speed by using Doppler effect is summarized systematically. The fault sliding speed of Wenchuan earthquake is calculated by the method mentioned above. The calculations show that Wenchuan earthquake fault sliding speeds are basically consistent with the seismic moment changes. It has proved that the sliding speed of Wenchuan earthquake fault has the characteristics of abrupt change, that is, the sliding speed increases suddenly and decreases rapidly. For most of the time, the sliding speed is not large, and sometimes the fault does not slide. There are obvious differences from sliding speed functions that are applied currently, such as Haskell function, Bell shaped functions, Exponential function, Triangle function, etc. To compute fault sliding speed by using Doppler effect, instead of grasping so-far unknown crust parameters, we only need to know the earthquake records, the locations of the epicenter and the receivers. In a word, the calculation method has clear physical meaning and the parameters required are easier to be obtained.

Realistic Shakemap M6.5 Pidie Jaya Earthquake 7 December 2016 Based on Modal Summation Technique

The shakemap were usually generated based on strong ground motion calculations on attenuation equations of the Ground Motion Prediction Equation (GMPE). However, the GMPE is a general equation which does not take into account the structural model of each layer of the subsurface. It considers only the structure of the surface layer known as Vs30 calculated as amplification in shakemap calculation. In this paper, the realistic calculation of shakemap for Pidie Jaya earthquake of 7 December 2016 with magnitude M6.5 is performed. Shakemap is generated from a number of synthetics seismogram calculations from source to site-based capital calculation summation Technique. Neo Deterministic Seismic Hazard Analysis (NDSHA) has included in the calculation of the synthetics seismograph. Before engaging the Vs30 effect, the realistic shakemap capable of producing effects of the angular distribution of a double-couple model in a realistic calculation. This could give an advantage of realistic shakemap while the GMPE based shakemap had not taken into account the double couple model. The realistic calculation of modal summation technique for shakemap can be further developed by involving local structural models in a more comprehensive manner.

Constraints on D″ beneath the North Atlantic region fromPandStraveltimes and amplitudes

In this study, we discuss possible origins of the D″ reflector beneath the North Atlantic region based on a combined analysis of P and S wave data. We use over 700 USArray station recordings of the M_w 6.3 earthquake that occurred in April 2010 in Spain. In order to investigate the D″ layer we look for waves reflected off the top of it, namely PdP and SdS waves, and compare them to the core–mantle boundary (CMB) reflections used as reference phases. The differences in traveltimes and amplitudes are sensitive to D″ properties. Because the USArray installation generates a dense array, we are able to provide an almost continuous map of the detection or absence of PdP and SdS waves in the North Atlantic region. We use a Bayesian inversion for traveltimes, together with synthetic seismogram calculations, to find the best-fitting D″ properties, (Vp, Vs) jumps across the D″ interface and D″ thickness. We find that the best-fitting models are for a D″ layer of about 300 km thick, with or without a velocity gradient of about 30 km at the top of it. Regardless of the model type, positive and similar velocity increases in both P and S velocities at the D″ interface, ranging from 2.7 to 3.8 per cent, are required to fit the data well. Our data rule out velocity decreases in P and S waves at the D″ interface as well as no velocity reduction above the CMB. There are also regions where we do not observe PdP and SdS waves. Collectively, these observations suggest lateral variations in both chemistry and temperature, combined with phase transitions. For instance, ancient oceanic basalt debris from the Farallon slab could be modulating the detection of the D″ reflector in this region.

GMPE based Shakemap Generation of Peak Ground Motion and Intensity Maps for Pidie Jaya Earthquake

An algorithm based on GMPE (Ground Motion Prediction Equation) is used for Generating standard and officially Shakemap as early information for disaster mitigation quick respond. We investigate the method developed by USGS (US Geological Survey) namely ShakeMap for generating Peak Ground Motion and Intensity Maps. The 2016 Mw 6.5Pidie Jaya Earthquake which strike-slip mechanism was selected for the scenario earthquake. The appropriate GMPE (i.e.Zhao06_crustal) was selected for this earthquake, which is located at the crust with focal depth of13km and strike-slip mechanism. Strong ground motion on bedrock, without sediment amplification, generate PGA (Peak Ground Acceleration) and PGV (Peak Ground Velocity) spherically uniform. The algorithm able to consider the site effect only for thin sedimentary layer Vs30 without considering deeper structural layer. The amplification effect is calculated based on Vs30 without considering the nonlinear effect and independently each site-class. Because of unavailable spatial distribution Vs30 measurement, the Vs30 was estimated from the topographic gradient and the method is commonly used. Strong ground motion on sediment layer will be amplified following the Vs30 distribution. The intensity map, which is easier to read, shown in MMI (Modified Mercalli Intensity) and SIG-BMKG (Indonesia Earthquake Intensity Scale). The GMPE based shakemap has bias error because of the general equation has not considering complete structural model and radiation pattern of a focal mechanism. An advance method for thefuture development will be much more realistic with the shakemap based on synthetics seismogram calculation.

Verification Modal Summation Technique for Synthetic and Observation Seismogram for Pidie Jaya Earthquake M6.5

Neo-Deterministic Seismic Hazard Assessment (NDSHA) method is a seismic hazard assessment method that has an advantage on realistic physical simulation of the source, propagation, and geological-geophysical structure. This simulation is capable on generating the synthetics seismograms at the sites that being observed. At the regional NDSHA scale, calculation of the strong ground motion is based on 1D modal summation technique because it is more efficient in computation. In this article, we verify the result of synthetic seismogram calculations with the result of field observations when Pidie Jaya earthquake on 7 December 2016 occurred with the moment magnitude of M6.5. Those data were recorded by broadband seismometers installed by BMKG (Indonesian Agency for Meteorology, Climatology and Geophysics). The result of the synthetic seismogram calculations verifies that some stations well show the suitability with observation while some other stations show the discrepancies with observation results. Based on the results of the observation of some stations, evidently 1D modal summation technique method has been well verified for thin sediment region (near the pre-tertiary basement), but less suitable for thick sediment region. The reason is that the 1D modal summation technique excludes the amplification effect of seismic wave occurring within thick sediment region. So, another approach is needed, e.g., 2D finite difference hybrid method, which is a part of local scale NDSHA method.

Forward Modeling Time-Lapse Seismic based on Reservoir Simulation Result on The CCS Project at Gundih Field, Indonesia

This paper contains reservoir simulation study of carbon storage at Gundih field in Central Java Island, Indonesia. Two different cases of injection simulation were performed and analyzed in this paper. The cases represent the conditions when the smallest and largest volumes of CO2areinjected into the subsurface to see the changes of reservoir that happen after the injection processes. The simulation result shows that when a larger amount of CO2 is injected into the targeted reservoir, it will migrate to the peak of anticline structure located in the southeast of CO2 injection well. The displacement of CO2 in the simulation progress shows that it will not reach the fault location. The geological model for synthetic seismogram calculation is then built based on the simulation reservoir result. The furthest displacement of CO2 is calculated on each case and described as the saturated CO2 layers. Forward modeling is performed to create synthetic seismic gather which will be processed to construct seismic section. The difference between the initial seismic section before the injection process and seismic section including saturated CO2 layer after the injection process will be evaluated by the potential of injected CO2 monitoring using time-lapse seismic survey in the Gundih field.

Calculation of differential seismograms using analytic partial derivatives—II Regional waveforms

Calculation of differential seismograms using analytic partial derivatives—II Regional waveforms

Calculation of differential seismograms using analytic partial derivatives – I: teleseismic receiver functions

Calculation of differential seismograms using analytic partial derivatives – I: teleseismic receiver functions

Amplitude variation with offset inversion using the reflectivity method

Most existing amplitude variation with offset (AVO) inversion methods are based on the Zoeppritz’s equation or its approximations. These methods assume that the amplitude of seismic data depends only on the reflection coefficients, which means that the wave-propagation effects, such as geometric spreading, attenuation, transmission loss, and multiples, have been fully corrected or attenuated before inversion. However, these requirements are very strict and can hardly be satisfied. Under a 1D assumption, reflectivity-method-based inversions are able to handle transmission losses and internal multiples. Applications of these inversions, however, are still time-consuming and complex in computation of differential seismograms. We have evaluated an inversion methodology based on the vectorized reflectivity method, in which the differential seismograms can be calculated from analytical expressions. It is computationally efficient. A modification is implemented to transform the inversion from the intercept time and ray-parameter domain to the angle-gather domain. AVO inversion is always an ill-posed problem. Following a Bayesian approach, the inversion is stabilized by including the correlation of the P-wave velocity, S-wave velocity, and density. Comparing reflectivity-method-based inversion with Zoeppritz-based inversion on a synthetic data and a real data set, we have concluded that reflectivity-method-based inversion is more accurate when the propagation effects of transmission losses and internal multiples are not corrected. Model testing has revealed that the method is robust at high noise levels.

Inverse AVO problem for a stack of layers

The problem of estimating thin layered model parameters by amplitude variation with offset (AVO) inversion has been studied. The motivation was resolving of the thin layers in inverted prestack seismic data as it contains more information on elastic properties of the subsurface than poststack seismic data.In this paper, an algorithm for solving the prestack inverse AVO problem in the case of multilayered media is derived. This algorithm is based on iterative corrections to the parameters of the initial model which tend to minimise the misfits between observed and synthetic seismograms. The synthetic seismograms are calculated using the reflection–transmission (RT)-matrices method, assuming a plane-wave with respect to the source position.A regularised Gauss-type algorithm for the inversion of prestack seismic data has been used. A differential seismogram computation algorithm to characterise the sensitivity of the seismic signal to the variations of a model parameter was used. The derived solution of the inverse problem is constructed in the time domain. This gives a slight advantage because it allows for visual control of the solution process. One can monitor the amplitude reduction of the data residual (difference between observed and synthetic seismograms) during the iteration process. Numerical examples show the accuracy and efficiency of the method.

Instaseis: instant global seismograms based on a broadband waveform database

Abstract. We present a new method and implementation (Instaseis) to store global Green's functions in a database which allows for near-instantaneous (on the order of milliseconds) extraction of arbitrary seismograms. Using the axisymmetric spectral element method (AxiSEM), the generation of these databases, based on reciprocity of the Green's functions, is very efficient and is approximately half as expensive as a single AxiSEM forward run. Thus, this enables the computation of full databases at half the cost of the computation of seismograms for a single source in the previous scheme and allows to compute databases at the highest frequencies globally observed. By storing the basis coefficients of the numerical scheme (Lagrange polynomials), the Green's functions are 4th order accurate in space and the spatial discretization respects discontinuities in the velocity model exactly. High-order temporal interpolation using Lanczos resampling allows to retrieve seismograms at any sampling rate. AxiSEM is easily adaptable to arbitrary spherically symmetric models of Earth as well as other planets. In this paper, we present the basic rationale and details of the method as well as benchmarks and illustrate a variety of applications. The code is open source and available with extensive documentation at www.instaseis.net .

Gridification of regional Fault Plain Solution (FPS)

The main goal of the Fault Plain Solution (FPS) application is the determination of the faulting mechanism that characterizes earthquake parameters (dip, strike, rake). For this purpose, regional Moment Tensor Inversion (MTI) method that calculates faulting parameters using three-dimensional broad-band waveform seismic data was used. FPS application is composed of two main software components: the first one is Time Domain Moment Tensor Inverse Code (TDMT_INVC with synthetic seismogram calculation) and the second one is Seismic Data Server Application Service (SDSAS) that is used for retrieving waveform data from storage elements. The most time consuming part of the FPS application involves the calculation of Green’s functions. Therefore, we parallelized calculation of the Greens functions. A Job Description Language (JDL) workflow is automatically created based on the initial parameters of individual crustal models that are necessary for producing Green functions. Timing results obtained show that our FPS application did achieve speed-up as a result of parallelization.

Matrix impedance in the problem of the calculation of synthetic seismograms for a layered-homogeneous isotropic elastic medium

A new method is proposed for calculating synthetic seismograms caused by a force in a plane-parallel medium consisting of homogeneous elastic isotropic layers. The matrix impedance, i.e., the matrix function of depth, by which motion vector must be multiplied in order to obtain the stress vector, is introduced for solving a system of ordinary differential equations with respect to the motion-stress vector, which appears during the separation of variables. An independent nonlinear equation is obtained for the impedance. The propagator for the motion vector is constructed with the aid of the impedance. The closed analytical formulas, which do not contain any exponents with positive indices, are obtained both for the impedance and for the motionvector propagator. The algorithm for the calculation of seismograms, free of limitations on the number and thickness of layers, as well as on the frequency range of interest, is constructed on the basis of these formulas. The algorithm is tested with the aid of an analytical solution.

The Earth Flattening Transformation in Body Wave Theory

Summary A simple class of power law transformations between spherical and plane models is investigated to determine the effect of the Earth flattening approximation in body wave studies. The transformations lead to the correct kinematic properties and should be chosen to reduce the inaccuracies in the dynamic properties. It is concluded that except in a fluid, the Love wave transformation of Biswas & Knopoff is the most satisfactory for SH, SV and P ray studies, Biswas's transformation for Rayleigh waves is not suitable for body wave studies, and the transformation used by Helmberger is optimum in a fluid. htrodnction The computation of theoretical seismograms in an inhomogeneous, spherically symmetric, Earth model has received considerable attention recently. Gilbert & Helmberger (1972) have applied the Cagniard-de Hoop generalized ray method to a layered sphere and the method has been used with considerable success by Helmberger & Wiggins (1971) in upper mantle studies. Gilbert & Helmberger (1972) describe their method as a direct approximation to the spherical wave functions. These approximations are invalid in layers where the ray is travelling nearly tangentially, e.g. near the turning point of the ray. If an Earth flattening approximation is applied first, the Cagniard-de Hoop method can be applied exactly to the approximate model of homogeneous plane layers. This method has been used by Miiller (1970) and is in fact the method used by Helmberger & Wiggins (1971) (Wiggins, private communication). Helmberger (1973) has used an identical transformation for velocity and density, and compared the exact solution in a homogeneous model with the Cagniardde Hoop solution. In this paper we shall discuss the nature of the approximation made in the Earth flattening approximation. In a recent paper, Hill (1972) has considered the approximations to the spherical wave functions in more detail. He obtains phase corrections to the plane wave approximations and frequency dependent corrections to the plane wave reflection coefficients. The combined effect of these approximations on the dynamic properties of the body wave and the choice of density transformation is not investigated. Chapman & Phinney (1972) have described an alternative spectral method for computing body wave seismograms. Although their method includes a solution of the fourth order spherical wave equations, it is necessary to match this to an approximate WKBJ solution. The results in this paper can also be used to indicate the approximation made here.

ACE—A FORTRAN subroutine for analytical computation of effective grid parameters for finite-difference seismic waveform modeling with standard Earth models

Despite the broad use of the heterogeneous finite-difference (FD) method for seismic waveform modeling, accurate treatment of material discontinuities inside the grid cells has been a serious problem for many years. One possible way to solve this problem is to introduce effective grid elastic moduli and densities (effective parameters) calculated by the volume harmonic averaging of elastic moduli and volume arithmetic averaging of densities in grid cells. This scheme enables us to place a material discontinuity in an arbitrary position in the spatial grids. Standard Earth models have made a significant contribution to synthetic seismogram calculations with a variety of numerical procedures such as the FD method. For the FD computation of seismic waveform with these models, we must first ensure accurate treatment of material discontinuities in the radius (or depth). The present paper introduces a FORTRAN subroutine ACE which calculates effective parameters analytically for an arbitrary spatial region in either the radius or depth direction for four major standard Earth models, namely, the PREM, IASP91, SP6, and AK135. This program is intended for all FD users who are concerned with seismic wave simulation for these models.

Method for rapid high-frequency seismogram calculation

We present a method for rapid, high-frequency seismogram calculation that makes use of an algorithm to automatically generate an exhaustive set of seismic phases with an appreciable amplitude on the seismogram. The method uses a hierarchical order of ray and seismic-phase generation, taking into account some existing constraints for ray paths and some physical constraints. To compute synthetic seismograms, the COMRAD code (from the Italian: “COdice Multifase per il RAy-tracing Dinamico”) uses as core a dynamic ray-tracing code. To validate the code, we have computed in a layered medium synthetic seismograms using both COMRAD and a code that computes the complete wave field by the discrete wave number method. The seismograms are compared according to a time–frequency misfit criteria based on the continuous wavelet transform of the signals. Although the number of phases is considerably reduced by the selection criteria, the results show that the loss in amplitude on the whole seismogram is negligible. Moreover, the time for the computing of the synthetics using the COMRAD code (truncating the ray series at the 10th generation) is 3–4-fold less than that needed for the AXITRA code (up to a frequency of 25 Hz).

Deep structure at northern margin of Tarim Basin

In this paper, a 2D velocity structure of the crust and the upper mantle of the northern margin of the Tarim Basin (TB) has been obtained by ray tracing and theoretical seismogram calculation under the condition of 2D lateral inhomogeneous medium using the data of seismic wide angle reflection/refraction profile from Baicheng to Da Qaidam crossing the Kuqa Depression (KD) and Tabei Uplift (TU). And along the Baicheng to Da Qaidam profile, 4 of the 10 shot points are located in the northern margin of the TB. The results show that the character of the crust is uniform on the whole between the KD and TU, but the depth of the layers, thickness of the crust and the velocity obviously vary along the profile. Thereinto, the variation of the crust thickness mainly occurs in the middle and lower crust. The Moho has an uplifting trend near the Baicheng shot point in KD and Luntai shot point in TU, and the thickness of the crust reduces to 42 km and 47 km in these two areas, respectively. The transition zone between the KD and TU has a thickest crust, up to 52 km. In this transition zone, there are high velocity anomalies in the upper crust, and low velocity anomalies in the lower crust, these velocity anomalies zone is near vertical, and the sediment above them is thicker than the other areas. According to the velocity distributions, the profile can be divided into three sections: KD, TU and transition zone between them. Each section has a special velocity structural feature, the form of the crystalline basement and the relationship between the deep structure and the shallow one. The differences of velocity and tectonic between eastern and western profile in the northern margin of the Tarim Basin (NMTB) may suggest different speed and intensity of the subduction from the Tarim basin to the Tianshan orogenic belt (TOB).

Inversion of far‐regional broadband P waves for the estimation of source parameters from shallow depth earthquakes

Attempts in separating the earthquake source effect from complex upper mantle wave propagation effects have been unsuccessful due to band‐limited World‐Wide Standardized Seismograph Network data and limitations in the computation of accurate seismograms. In this study, we examine the potential for estimating source parameters using broadband P waves observed between distances of 10° and 30°. We used the reflectivity method, which retains all of the rays, to compute synthetic Green's functions from the Earth‐flattened global AK135 and regional T7 isotropic velocity models. We applied the moment tensor inversion and iterative dislocation grid search methods to estimate seismic moment, moment tensor, focal mechanism, and focal depth for three example seismic events, the 1995 (Mw 5.7) west Texas, 1999 (Mw 5.6) Scotty's Junction, Nevada, and 1999 southern Iran (Mw 4.8) earthquakes. The west Texas and Scotty's Junction earthquakes were well recorded, and they had independently established source parameters. Independent estimates from different data sets and agencies were used to validate the inversion results. We selected the southern Iran earthquake to show the utility of the method at lower magnitudes. Selected vertical and radial component seismograms were windowed about the P wave arrival, corrected for instrument response, and converted to velocity, preserving the absolute amplitudes and periods within approximately 0.5 to 100 s. Both inverse methods were successful in separating the source effect with the use of the AK135 velocity model. The inversion results using the T7 model appear slightly biased and had worse waveform fits compared to AK135. More importantly, the use of two different velocity models resulted in biases in focal depths that were less than ±2 km, thereby providing a robust inverse method for focal depth. Utilizing far‐regional P waves would be most useful for remote earthquakes not adequately recorded locally and too small in magnitude to be recorded beyond epicentral distances of 30° for teleseismic analysis.