Complete Synthetic Seismograms Research Articles

The algorithms for calculating synthetic seismograms from a dipole source using the derivatives of Green’s function

The problem of calculating complete synthetic seismograms from a point dipole with an arbitrary seismic moment tensor in a plane parallel medium composed of homogeneous elastic isotropic layers is considered. It is established that the solutions of the system of ordinary differential equations for the motion–stress vector have a reciprocity property, which allows obtaining a compact formula for the derivative of the motion vector with respect to the source depth. The reciprocity theorem for Green’s functions with respect to the interchange of the source and receiver is obtained for a medium with cylindrical boundary. The differentiation of Green’s functions with respect to the coordinates of the source leads to the same calculation formulas as the algorithm developed in the previous work (Pavlov, 2013). A new algorithm appears when the derivatives with respect to the horizontal coordinates of the source is replaced by the derivatives with respect to the horizontal coordinates of the receiver (with the minus sign). This algorithm is more transparent, compact, and economic than the previous one. It requires calculating the wavenumbers associated with Bessel function’s roots of order 0 and order 1, whereas the previous algorithm additionally requires the second order roots.

Complete synthetic seismograms based on a spherical self-gravitating Earth model with an atmosphere–ocean–mantle–core structure

Complete synthetic seismograms based on a spherical self-gravitating Earth model with an atmosphere–ocean–mantle–core structure

2(3)D convolution modelling of complex geological targets beyond – 1D convolution

Seismic modelling is necessary to understand elasticwave propagation in the subsurface. Modelling is costeffective and insightful, as long as adequate methods are used. An ideal seismic-modelling strategy is to generate complete synthetic seismograms for realistic earth models, then process them as performed with real seismic data. These complete seismograms are best obtained by full-wavefield (FW) approaches. FW methods are therefore used in extensive benchmarking studies, which may require the joint effort of several institutions owing to high resource costs. Though computer power continues to increase, we are still a long way from applying ideal modelling to all cases where synthetic data is necessary to constrain results. This is especially true in interpretation and sensitivity analysis, where the influence of multiple parameters must be assessed. In many situations, ray-based (RB) methods are suitable alternatives, especially when rays, traveltimes, etc., are useful information for the problem at hand (Gjoystdal et al., 2007). However, standard RB methods do not allow modelling of detailed target structures owing to smoothness requirements. Geoscientists needing seismic modelling of such targets will then often resort to 1D convolution (Lecomte et al., 2015). 1D convolution has been successfully used for decades, and is still the method behind most standard well calibration, seismic inversion, etc. However, its conceptual validity is in reality very limited.

Teleseismic inversion of the 2004 Sumatra-Andaman earthquake rupture process using complete Green’s functions

Spatial and temporal variations in coseismic slip distribution are often obtained by rupture process analyses using teleseismic body waves. Many analyses using teleseismic body waves were based on the ray theory because of the very efficiently computable direct P-, S-, and major reflected waves near the source. The 2004 Sumatra-Andaman earthquake was one of the largest earthquakes recorded in history, and the data that are required for the entire rupture process analysis include later phases such as PP waves and a very long period phase called a W phase. However, calculating these later phases using the conventional ray theoretical method is difficult. Here we investigate the rupture process of the 2004 Sumatra-Andaman earthquake using complete Green’s functions, including all later phases such as PP waves and W phase. We use the direct solution method, which computes complete synthetic seismograms up to 2 Hz for transversely isotropic spherically symmetric media, to calculate the Green’s functions. The obtained synthetic seismograms generated fit the observed seismograms quite well from a short period to a long period. The estimated slip distribution consists of four large slip areas: the largest slip occurred in the shallow part off the Sumatra west coast with a maximum slip of approximately 29 m, the second and third largest slips occurred in the shallow and deep parts of the Nicobar region with maximum slips of approximately 8 m and 7 m, respectively, and the fourth slip occurred in the middle Andaman region with a maximum slip of approximately 6 m. The estimated average rupture velocity is 2.8 km/s, but the rupture may have slowed between the Sumatra and the shallow Nicobar slip areas, and between the Nicobar and the middle Andaman slip areas. The delayed initiation of the shallow slips in the Nicobar region may possibly have been triggered by the deeper slip in the Nicobar region. There were no distinct depth-varying properties for the shallow and deep slips in the Nicobar region, as were reported for the 2011 Tohoku-Oki and 2010 Chile earthquake.

Focal mechanism and depth of the 1956 Amorgos twin earthquakes from waveform matching of analogue seismograms

Abstract. Historic analogue seismograms of the large 1956 Amorgos twin earthquakes which occurred in the volcanic arc of the Hellenic subduction zone (HSZ) were collected, digitized and reanalyzed to obtain refined estimates of their depth and focal mechanism. In total, 80 records of the events from 29 European stations were collected and, if possible, digitized. In addition, bulletins were searched for instrument parameters required to calculate transfer functions for instrument correction. A grid search based on matching the digitized historic waveforms to complete synthetic seismograms was then carried out to infer optimal estimates for depth and focal mechanism. Owing to incomplete or unreliable information on instrument parameters and frequently occurring technical problems during recording, such as writing needles jumping off mechanical recording systems, much less seismograms than collected proved suitable for waveform matching. For the first earthquake, only seven seismograms from three different stations at Stuttgart (STU), Göttingen (GTT) and Copenhagen (COP) could be used. Nevertheless, the waveform matching grid search yields two stable misfit minima for source depths of 25 and 50 km. Compatible fault plane solutions are either of normal faulting or thrusting type. A separate analysis of 42 impulsive first-motion polarities taken from the International Seismological Summary (ISS bulletin) excludes the thrusting mechanism and clearly favors a normal faulting solution with at least one of the potential fault planes striking in SW–NE direction. This finding is consistent with the local structure and microseismic activity of the Santorini–Amorgos graben. Since crustal thickness in the Amorgos area is generally less than 30 km, a source depth of 25 km appears to be more realistic. The second earthquake exhibits a conspicuously high ratio of body wave to surface wave amplitudes suggesting an intermediate-depth event located in the Hellenic Wadati–Benioff zone. This hypothesis is supported by a focal mechanism analysis based on first-motion polarities, which indicates a mechanism very different from that of the first event. A waveform matching grid search done to support the intermediate-depth hypothesis proved not to be fruitful because the body wave phases are overlain by strong surface wave coda of the first event inhibiting a waveform match. However, body to surface wave amplitude ratios of a modern intermediate-depth event with an epicenter close to the island of Milos observed at stations of the German Regional Seismic Network (GRSN) exhibit a pattern similar to the one observed for the second event with high values in a frequency band between 0.05 Hz and 0.3 Hz. In contrast, a shallow event with an epicenter in western Crete and nearly identical source mechanism and magnitude, shows very low ratios of body and surface wave amplitude up to 0.17 Hz and higher ratios only beyond that frequency. Based on this comparison with a modern event, we estimate the source depth of the second event to be greater than 100 km. The proximity in time and space of the two events suggests a triggering of the second, potentially deep event by the shallow first one.

Algorithm for calculating synthetic seismograms in a layered half-space with application of matrix impedance

A new algorithm is proposed for calculating the complete synthetic seismograms from a point source in the form of the sum of a single force and a dipole with an arbitrary seismic moment tensor in a plane layered medium composed of homogenous elastic isotropic layers. Following the idea of (Alekseev and Mikhailenko, 1978), an artificial cylindrical boundary is introduced, on which the boundary conditions are specified. For this modified problem, the exact solution (in terms of the displacements and stresses on the horizontal plane areal element) in the frequency domain is derived and substantiated. The unknown depth-dependent coefficients form the motion-stress vector, whose components satisfy the known system of ordinary differential equations. This system is solved by the method that involves the matrix impedance and propagator for the vector of motion, as previously suggested by the author in (Pavlov, 2009). In relation to the initial problem, the reflections from the artificial boundary are noise, which, to a certain degree, can be suppressed by selecting a long enough distance to this boundary and owing to the presence of a purely imaginary addition to the frequency. The algorithm is not constrained by the thickness of the layers, is applicable for any frequency range, and is suitable for computing the static offset.

Numerically stable computation of complete synthetic seismograms including the static displacement in plane layered media

SUMMARY We present a new method for calculating dynamic and static displacements arising from a buried point source in a plane layered, isotropic, elastic half-space. We avoid the problem of the loss of precision in solutions for P-SV wave motion by integrating a system of minor vector equations with solutions constructed from pairs of solutions to the equations of motion in the frequency–wavenumber domain. The resulting algorithm is efficient, and numerically stable at zero frequency and at high frequencies, and thus can be used to compute complete synthetic seismograms that include the static offset. For the special case of the static deformation of a homogeneous half-space, we show that our results are equivalent to formulae that are often currently used in the geodetic community. An advantage of our method is that a single algorithm can be used to model both seismological and geodetic measurements of ground motion. We illustrate one seismological application of our method using continuous Global Positioning System (GPS) data obtained for a recent earthquake; seismograms calculated in a realistic layered crustal model emulate the observed dynamic and static displacements well. If the calculation is only carried out at 0 Hz, we obtain an image of the static surface deformation, and thus our method also has potential applications in the modelling of coseismic surface displacements measured with, for example, Synthetic Aperture Radar Interferometry.

First Comparison of Array-Derived Rotational Ground Motions with Direct Ring Laser Measurements

Recently, ring laser technology has provided the first consistent obser- vations of rotational ground motions around a vertical axis induced by earthquakes. Consistent, in this context, implies that the observed waveforms and amplitudes are compatible with collocated recordings of translational ground motions. In partic- ular, transverse accelerations should be in phase with rotation rate and their ratio proportional to local horizontal phase velocity assuming plane-wave propagation. The ring laser installed at the Fundamentalstation Wettzell in the Bavarian Forest, Southeast Germany, is recording the rotation rate around a vertical axis, theoretically a linear combination of the space derivatives of the horizontal components of motion. This suggests that, in principle, rotation can be derived from seismic-array experi- ments by finite differencing. This has been attempted previously in several studies; however, the accuracy of these observations could never be tested in the absence of direct measurements. We installed a double cross-shaped array of nine stations from December 2003 to March 2004 around the ring laser instrument and observed several large earthquakes on both the ring laser and the seismic array. Here we present for the first time a comparison of array-derived rotations with direct measurements of rotations for ground motions induced by the M 6.3 Al Hoceima, Morocco, earthquake of 24 February 2004. With complete 3D synthetic seismograms calculated for this event we show that even low levels of noise may considerably influence the accuracy of the array-derived rotations when the minimum number of required stations (three) is used. Nevertheless, when using all nine stations, the overall fit between direct and array-derived measurements is surprisingly good (maximum correlation coefficient of 0.94).

Complete synthetic seismograms up to 2 Hz for transversely isotropic spherically symmetric media

SUMMARY We use the direct solution method (DSM) with optimally accurate numerical operators to calculate complete (including both body and surface waves) three-component synthetic seismograms for transversely isotropic (TI), spherically symmetric media, up to 2 Hz. We present examples of calculations for both deep (600 km) and shallow (5 km) sources. Such synthetics should be useful in forward and inverse studies of earth structure. In order to make these calculations accurately and efficiently the vertical grid spacing, maximum angular order, and cut-off depth must be carefully and systematically chosen.

Beno Gutenberg contribution to seismic hazard assessment and recent progress in the European–Mediterranean region

The fundamental work of Beno Gutenberg has inspired and guided an appreciable part of research in modern seismology, both from the experimental and the theoretical point of view. Among the several topics of seismology that have benefited from the fundamental contribution of Gutenberg, we consider particularly relevant the description of the asthenospheric low-velocity channel, the definition of the surface waves magnitude and the Gutenberg–Richter law, since they are pivotal tools for seismic hazard assessment. The quite revolutionary model for the lithosphere–asthenosphere system in the European area predicts the existence of almost aseismic lithospheric roots. These roots are located in correspondence of most of the orogenic belts and interrupt the asthenosphere low velocity channel that has been identified by Beno Gutenberg in 1948. The model of the European upper mantle, proposed for the first time in 1979 and subsequently refined, has stimulated a considerable amount of research, which has nicely confirmed the major innovative features of the early model. At present, the subduction of the lithosphere at continent–continent collisions, supported not only by seismological data, is a widely accepted concept within the community of Earth scientists, even if it contradicts one of the basic dogmas of the original formulation of plate tectonics. The proposed model for the Alpine–Apennines area supplies a new and unifying framework for the interpretation of the Quaternary magmatism, at present generally accepted by petrologists and geochemists. The theoretical basis for the Gutenberg's surface-wave magnitude calibration function has been supplied by the use of complete synthetic seismograms, and thus it has been possible to formulate the theoretical M S depth correction. The introduction of the depth correction for M S enables the computation of surface wave magnitude for all earthquakes, regardless of their focal depth. This is especially important for the quantification of subcrustal historical earthquakes, for which the seismic moment may be difficult to estimate from recordings of early mechanical seismographs. The new M S calibrating function yield both distance- and depth-independent magnitude estimates. The analysis of the global seismicity, using the seismotectonic regionalization in subduction zones, mid oceanic ridge zones, island arcs, shows that a single Gutenberg–Richter (GR) relation is not universally applicable and that a multiscale seismicity model can reconcile two apparently conflicting paradigms: the Self-Organized Criticality mechanism and the Characteristic Earthquake concept. The multiscale representation has been applied to Italy, where the zones at the space scale of 400–500 km quite well reproduce the shapes of the regions used to apply the, globally tested, CN intermediate term earthquake prediction algorithm.

Finite difference modelling of P-wave scattering in the upper mantle

Summary Scattering within the upper mantle is prominently documented in the so-called high-frequency teleseismic Pn phase, generated by velocity fluctuations in the upper mantle. This phase is seen in many data sets from active and passive seismology. Whenever the source-generated frequencies are high enough (5–10 Hz) to allow their observation at large distances (>2000 km), this phase appears in the recordings of the Russian PNE programme, in other long-range refraction experiments, and in numerous recordings of earthquakes. Using a 2-D finite difference scheme we calculate complete synthetic seismograms in order to relate the scales of the fluctuation of the elastic parameters to wavefield properties. We discuss the influence of critical parameters on the scattering properties of the upper mantle and on the propagation mechanism of the teleseismic Pn: variations in vertical and horizontal correlation lengths, rms velocity fluctuations, thickness of the heterogeneous layer, and cross-correlation of P- and S-wave velocities. A teleseismic Pn phase arises only if the velocity fluctuations are strong enough to cause multiple scattering and if the horizontal correlation length exceeds the vertical by an order of magnitude. Our favourite model for the Russian PNE profiles includes a 100–150-km-thick zone of scatterers below the Moho, containing heterogeneities with a horizontal correlation length of 20 km, a vertical correlation length of 0.5 km and an rms velocity perturbation of 2 per cent, while keeping the P- and S-wave velocities correlated. As the teleseismic Pn is often observed, this scale of velocity fluctuations seems to be characteristic of the uppermost mantle. Anisotropy is likely to contribute to the inferred velocity fluctuations.

Complete synthetic seismograms for 3-D heterogeneous Earth models computed using modified DSM operators and their applicability to inversion for Earth structure

We compute complete (including both body and surface waves) synthetic seismograms for laterally and vertically heterogeneous Earth models using the Direct Solution Method (DSM). We use the optimally accurate modified operators derived by Geller and Takeuchi [Geller, R.J., Takeuchi, N., 1995. A new method for computing highly accurate DSM synthetic seismograms. Geophys. J. Int. 123, 449–470] and extended to spherical coordinates by Takeuchi et al. [Takeuchi, N., Geller, R.J., Cummins, P.R., 1996. Highly accurate P-SV complete synthetic seismograms using modified DSM operators. Geophys. Res. Lett. 23, 1175–1178] and Cummins et al. [Cummins, P.R., Takeuchi, N., Geller, R.J., 1997. Computation of complete synthetic seismograms for laterally heterogenous models using the Direct Solution Method. Geophys. J. Int. 130, 1–16] for 1- and 3-D models, respectively. In this study we greatly reduce the CPU time by treating the laterally heterogeneous structure as a perturbation to a spherically symmetric model (i.e., using the Born approximation). Note, however, that (1) our methods do not require the use of the Born approximation and (2) the reference model for the Born approximation is not required to be spherically symmetric. The synthetic seismograms in this paper are computed using the first-order Born approximation. However, accuracy can be greatly improved by using higher order terms of the Born series; theoretical results are presented in this paper, and some preliminary numerical examples are presented in this volume by Igel et al. [Igel, H., Takeuchi, N., Geller, R.J., Megnin, C., Bunge, H.P., Clévédé, E., Dalkolmo, J., Romanowicz, B., 1998. The COSY project: verification of global seismic modeling algorithms, Phys. Earth Planet. Inter., this issue].

Seismic Zoning of Slovenia Based on Deterministic Hazard Computations

Seismic hazard of the territory of Slovenia is estimated using a deterministic approach based on the computation of complete synthetic seismograms. The input data are the catalogues of earthquakes and fault plane solutions for Slovenia and surrounding regions. Structural models are defined based on available seismological and geophysical information, but are mainly constrained by surface-wave dispersion and 3-D tomographic modelling of the upper crust. Seismogenic zones are delineated considering geotectonic characteristics, fault plane solutions and distribution of earthquake hypocentres. Outside Slovenia seismogenic zones are extended up to distances from which they can considerably influence seismic hazard estimates.¶Synthetic seismograms are computed using the “receiver” structure along the entire path by normal mode summation (up to 1 Hz) for receiver sites on a 0.2 × 0.2 degrees grid and scaled to the magnitude of the earthquake allowing for spectral falloff. At each site the maximum value of horizontal velocity, horizontal displacement and design ground acceleration are considered as hazard parameter. The highest values are obtained for western Slovenia where the hazard is controlled by the strongest earthquake in the catalogue, the “Idrija” event of March 26, 1511.

Seismogenic potential and earthquake hazard assessment in the Tell Atlas of Algeria

Seismotectonic zonation studies in the Tell Atlas of Algeria, a branch of the Africa-Eurasia plate boundary, provide a valuable input for deterministic seismic hazard calculations. We delineate a number of seismogenic zones from causal relationships established between geological structures and earthquakes and compile a working seismic catalogue mainly from readily available sources. To this catalogue, for a most rational and best-justified hazard analysis, we add estimates of earthquake size translated from active faulting characteristics. We assess the regional seismic hazard using a deterministic procedure based on the computation of complete synthetic seismograms (up to 1 Hz) by the modal summation technique. As a result, we generate seismic hazard maps of maximum velocity, maximum displacement, and design ground acceleration that blend information from geology, historical seismicity and observational seismology, leading to better estimates of the earthquake hazard throughout northern Algeria. Our analysis and the resulting maps illustrate how different the estimate of seismic hazard is based primarily on combined geologic and seismological data with respect to the one for which only information from earthquake catalogues has been used.

Seismic Hazard of Romania: Deterministic Approach

—The seismic hazard of Romania is estimated in terms of peak-ground motion values—displacement, velocity, design ground acceleration (DGA)—computing complete synthetic seismograms, which are considered to be representative of the different seismogenic and structural zones of the country. The deterministic method addresses issues largely neglected in probabilistic hazard analysis, e.g., how crustal properties affect attenuation, since the ground motion parameters are not derived from overly simplified attenuation “functions,” but rather from synthetic time histories. The synthesis of the hazard is divided into two parts, one that of shallow-focus earthquakes, and the other, that of intermediate-focus events of the Vrancea region.¶The previous hazard maps of Romania completely ignore the seismic activity in the southeastern part of the country (due to the seismic source of Shabla zone). For the Vrancea intermediate-depth earthquakes, which control the seismic hazard level over most of the territory, the comparison of the numerical results with the historically-based intensity map show significant differences. They could be due to possible structural or source properties not captured by our modeling, or to differences in the distribution of damageable buildings over the territory (meaning that future earthquakes can be more spectacularly damaging in regions other than those regions experiencing damage in the past). Since the deterministic modeling is highly sensitive to the source and path effects, it can be used to improve the seismological parameters of the historical events.

Seismic characterization of the shore sand at Catania

In the framework of the Catania Project supported by GNDT (Gruppo Nazionale per la Difesa dai Terremoti) seismic measurements have been carried out in the shore sands of Catania to define their detailed shear wave velocity (V) profiles and attenuation properties (Q). The method used, FTAN, is based on the analysis of the dispersion curve extracted from signals generated by controlled sources. The computation of complete synthetic seismograms enables the estimation of the attenuation properties of the surface layers. Hereinafter we analyse the shear seismic characteristics of the shore sands as indicators of potential liquefaction phenomena.

Computation of complete synthetic seismograms for laterally heterogeneous models using the Direct Solution Method

SUMMARY We use the Direct Solution Method (DSM) together with the modified operators derived by Geller & Takeuchi (1995) and Takeuchi, Geller & Cummins (1996) to compute complete synthetic seismograms and their partial derivatives for laterally heterogeneous models in spherical coordinates. The methods presented in this paper are well suited to conducting waveform inversion for 3-D Earth structure. No assumptions of weak perturbation are necessary, although such approximations greatly improve computational efficiency when their use is appropriate. An example calculation is presented in which the toroidal wavefield is calculated for an axisymmetric model for which velocity is dependent on depth and latitude but not longitude. The wavefield calculated using the DSM agrees well with wavefronts calculated by tracing rays. To demonstrate that our algorithm is not limited to weak, aspherical perturbations to a spherically symmetric structure, we consider a model for which the latitude-dependent part of the velocity structure is very strong.

Highly accurate P‐SV complete synthetic seismograms using modified DSM operators

In previous papers [Cummins et al., 1994ab] (hereafter referred to as DSMI and DSMII respectively), we presented accurate methods for computing complete synthetic seismograms for SH and P‐SV respectively in a spherical earth model. The SH calculations used computationally efficient modified matrix operators, but the P‐SV synthetics were computationally intensive. Geller and Takeuchi [1995] (hereafter referred to as GT95) presented a general theory for deriving modified operators and gave the explicit form of the modified operators for the P‐SV case in cylindrical or cartesian coordinates. In this paper we extend GT95's results to derive modified operators for the P‐SV case in spherical coordinates. The use of the modified operators reduces the CPU time by a factor of about 5 without a loss of accuracy. 10 CPU min on a SPARC‐20 workstation with one CPU are required to compute a profile of synthetic seismograms from DC to 20 sec period.

Complete synthetic seismograms for a spherically symmetric earth by a numerical computation of the Green's function in the frequency domain

We present a new method to calculate complete synthetic seismograms for a spherically symmetric earth model which uses neither eigenfrequencies and eigenfunctions nor an earth-flattening transformation. The response of the earth to a moment tensor point source is evaluated in the frequency domain for both spheroidal and toroidal motion by numerical integration of the appropriate system of ordinary differential equations with source term and summation over vector spherical harmonics. Attenuation is included by using complex elastic moduli. Owing to the discrete sampling of the response in the frequency domain, the numerical effort is proportional to the length of the desired time series for a fixed maximum frequency. This makes the method much more efficient than normal-mode calculations for higher frequency applications, where often seismogram lengths of 20 to 40 min are sufficient. Since the angular degree of the spherical harmonics provide a natural discretization in the wavenumber domain, spatial aliasing is unimportant. Time aliasing is suppressed by evaluating the response at complex frequencies with constant imaginary part. We have compared synthetic seismograms obtained by the new method with normal-mode seismograms up to a frequency of 20 mHz and achieve excellent agreement for all three components. The accuracy of the method is further corroborated by comparisons with real data up to a frequency of 200 mHz. We tested the numerical scheme up to frequencies of 1 Hz and harmonic degrees of 12 000 and did not find any numerical instabilities. Incidentally, the approach sheds some light on how normal modes make up body waves.

Efficient calculation of complete differential seismograms for laterally homogeneous earth models

SUMMARY Differential seismograms, the sensitivity of the seismogram to perturbations of parameters of the earth model, are essential for waveform modelling with iterative non-linear optimization techniques. The growing number of broad-band highdynamic-range fixed and portable stations provide a rich data set of regional distance events with accurate locations and source mechanisms that are ideal for modelling lithospheric structure. An efficient technique for calculating differential seismograms for complete synthetic seismograms in laterally homogeneous earth models will significantly decrease the computational burden of seismic-waveform inverse modelling. The technique we present here is equivalent to computing three synthetic seismograms in multilayered models, contrasted with N + 1 seismograms for a typical brute-force approach, which reduces the work to 3/(N + 1) or about a 90 per cent saving for 30 layers. The lateral-homogeneity requirement will yield path-averaged structures which may form the basis for later tomographic studies to resolve lateral variations.