ALBERT

Description: This paper presents a generalized wave equation which unifies viscoelastic and pure elastic cases into a single wave equation. In the generalized wave equation, the degree of viscoelasticity varies between zero and unity, and is defined by a controlling parameter. When this viscoelastic controlling parameter equals to 0, the viscous property vanishes and the generalized wave equation becomes a pure elastic wave equation. When this viscoelastic controlling parameter equals to 1, it is the Stokes equation made up of a stack of pure elastic and Newtonian viscous models. Given this generalized wave equation, an analytical solution is derived explicitly in terms of the attenuation and the velocity dispersion. It is proved that, for any given value of the viscoelastic controlling parameter, the attenuation component of this generalized wave equation perfectly satisfies the power laws of frequency. Since the power laws are the fundamental characteristics in physical observations, this generalized wave equation can well represent seismic wave propagation through subsurface media.

Description: We investigate the relationship between subduction processes and related seismicity for the Lesser Antilles Arc using the Gutenberg–Richter law. This power law describes the earthquake-magnitude distribution, with the gradient of the cumulative magnitude distribution being commonly known as the b -value. The Lesser Antilles Arc was chosen because of its along-strike variability in sediment subduction and the transition from subduction to strike-slip movement towards its northern and southern ends. The data are derived from the seismicity catalogues from the Seismic Research Centre of The University of the West Indies and the Observatoires Volcanologiques et Sismologiques of the Institut de Physique du Globe de Paris and consist of subcrustal events primarily from the slab interface. The b -value is found using a Kolmogorov–Smirnov test for a maximum-likelihood straight line-fitting routine. We investigate spatial variations in b -values using a grid-search with circular cells as well as an along-arc projection. Tests with different algorithms and the two independent earthquake cataloges provide confidence in the robustness of our results. We observe a strong spatial variability of the b -value that cannot be explained by the uncertainties. Rather than obtaining a simple north–south b -value distribution suggestive of the dominant control on earthquake triggering being water released from the sedimentary cover on the incoming American Plates, or a b -value distribution that correlates with on the obliquity of subduction, we obtain a series of discrete, high b -value ‘bull's-eyes’ along strike. These bull's-eyes, which indicate stress release through a higher fraction of small earthquakes, coincide with the locations of known incoming oceanic fracture zones on the American Plates. We interpret the results in terms of water being delivered to the Lesser Antilles subduction zone in the vicinity of fracture zones providing lubrication and thus changing the character of the related seismicity. Our results suggest serpentinization around mid-ocean ridge transform faults, which go on to become fracture zones on the incoming plate, plays a significant role in the delivery of water into the mantle at subduction zones.

Description: Given a moment tensor m inferred from seismic data for an earthquake, we define ${\scr P}(V)$ to be the probability that the true moment tensor for the earthquake lies in the neighbourhood of m that has fractional volume V . The average value of ${\scr P}(V)$ is then a measure of our confidence in  m . The calculation of ${\scr P}(V)$ requires knowing both the probability $\skew4\hat{P}(\omega )$ and the fractional volume $\skew4\hat{V}(\omega )$ of the set of moment tensors within a given angular radius of  m . We explain how to construct $\skew4\hat{P}(\omega )$ from a misfit function derived from seismic data, and we show how to calculate $\skew4\hat{V}(\omega )$ , which depends on the set $\mathbb {M}$ of moment tensors under consideration. The two most important instances of $\mathbb {M}$ are where $\mathbb {M}$ is the set of all moment tensors of fixed norm, and where $\mathbb {M}$ is the set of all double couples of fixed norm.

Description: The M s ~ 7.7 Sarez-Pamir earthquake of 1911 February 18 is the largest instrumentally recorded earthquake in the Pamir region. It triggered one of the largest landslides of the past century, building a giant natural dam and forming Lake Sarez. As for many strong earthquakes from that time, information about source parameters of the Sarez-Pamir earthquake is limited due to the sparse observations. Here, we present the analysis of analogue seismic records of the Sarez-Pamir earthquake. We have collected, scanned and digitized 26 seismic records from 13 stations worldwide to relocate the epicentre and determine the event's depth (~26 km) and magnitude ( m B 7.3 and M s 7.7). The unusually good quality of the digitized waveforms allowed their modelling, revealing an NE-striking sinistral strike-slip focal mechanism in accordance with regional tectonics. The shallow depth and magnitude ( M w 7.3) of the earthquake were confirmed. Additionally, we investigated the possible contribution of the landslide to the waveforms and present an alternative source model assuming the landslide and earthquake occurred in close sequence.

Description: Array measurements of surface wave phase velocity can be biased by multipath arrivals. A two-plane-wave (TPW) inversion method, in which the incoming wavefield is represented by the interference of two plane waves, is able to account for the multipath effect and solve for laterally varying phase velocity. Despite broad applications of the TPW method, its usage has been limited to Rayleigh waves. In this study, we have modified the TPW approach and applied it to Love waves. Main modifications include decomposing Love wave amplitude on the transverse component to x and y components in a local Cartesian system for each earthquake and using both components in the inversion. Such decomposition is also applied to the two plane waves to predict the incoming wavefield of an earthquake. We utilize fundamental mode Love wave data recorded at 85 broad-band stations from 69 distant earthquakes and solved for phase velocity in nine frequency bands with centre periods ranging from 34 to 100 s. The average phase velocity in southern Africa increases from 4.30 km s –1 at 34 s to 4.87 km s –1 at 100 s. Compared with predicted Love wave phase velocities from the published 1-D SV velocity model and radial anisotropy model in the region, these values are compatible from 34 to 50 s and slightly higher beyond 50 s, indicating radial anisotropy of V SH 〉 V SV in the shallow upper mantle. A high Love wave velocity anomaly is imaged in the central and southern Kaapvaal craton at all periods, reflecting a cold and depleted cratonic lithosphere. A low velocity anomaly appears in the Bushveld Complex from 34 to 50 s, which can be interpreted as being caused by high iron content from an intracratonic magma intrusion. The modified TPW method provides a new way to measure Love wave phase velocities in a regional array, which are essential in developing radial anisotropic models and understanding the Earth structure in the crust and upper mantle.

Description: Ocean waves activity is a major source of microvibrations that travel through the solid Earth, known as microseismic noise and recorded worldwide by broadband seismometers. Analysis of microseismic noise in continuous seismic records can be used to investigate noise sources in the oceans such as storms, and their variations in space and time, making possible the regional and global-scale monitoring of the wave climate. In order to complete the knowledge of the Atlantic and Pacific oceans microseismic noise sources, we analyse 1 yr of continuous data recorded by permanent seismic stations located in the Indian Ocean basin. We primarily focus on secondary microseisms (SM) that are dominated by Rayleigh waves between 6 and 11 s of period. Continuous polarization analyses in this frequency band at 15 individual seismic stations allow us to quantify the number of polarized signal corresponding to Rayleigh waves, and to retrieve their backazimuths ( BAZ ) in the time–frequency domain. We observe clear seasonal variations in the number of polarized signals and in their frequencies, but not in their BAZ that consistently point towards the Southern part of the basin throughout the year. This property is very peculiar to the Indian Ocean that is closed on its Northern side, and therefore not affected by large ocean storms during Northern Hemisphere winters. We show that the noise amplitude seasonal variations and the backazimuth directions are consistent with the source areas computed from ocean wave models.

Description: A method is presented to simulate the propagation of elastic waves in a 2-D piecewise homogeneous domain composed of several regions of arbitrary shapes. The regions can be isolated within an elastic background or also embedded among themselves. The method is based on the indirect boundary element method (IBEM) in which the diffracted field is obtained as the contributions of distributed sources along the regions boundaries. This is somehow reminiscent of Huygens's Principle for wave fronts. The fields emitted by the sources are computed in frequency domain using analytical expressions of Green's functions for displacements and tractions. This approach is validated, in the antiplane case, by comparing against a semi-analytical solution for a basin composed of concentric semi-circular layers and embedded in a half-space. Finally, an example of in-plane wave propagation is given for a complex configuration involving five interlocked regions.

Description: A set of uniaxial compression tests of granite specimens taken from five localities across Japan was conducted to identify the factors controlling the quantity of radon (Rn) emission (sum of 222 Rn and 220 Rn) during compression and failure. An α-scintillation detector and a gas flow unit were installed with a testing machine to enable continuous measurement of Rn emissions. Common to all specimens, Rn emissions remained at or slightly declined from the background level after the start of loading; this is similar to the natural phenomenon of decline in groundwater-dissolved Rn before an earthquake. Closure of original microcracks is the most likely cause of the initial Rn decline. Then, Rn emissions begin to increase at 46–57 per cent stress level to the uniaxial compressive strength, and continue to increase even after the failure of specimen. This commencement stress level is close to the general stress level at outbreak of acoustic emissions caused by the development and connection of microcracks. The Rn increase after failure is similar to a phenomenon observed in aftershocks, which may originate from the enhancement of Rn emanations from grains due to the large increase in total surface area and stress release. In addition to the initial radioelement content in rock, the failure pattern (conjugate shear versus longitudinal tensile type), compressive strength, and grain size are possible control factors of the maximum quantity of Rn emissions induced by failure. This maximum may also be affected by the development velocity of the emanation area, which is related to the Rn emanation fraction, associated with the fragmentation. In addition to the magnitude of an earthquake and its hypocentre distance to Rn detectors, the magnitude of increase in Rn concentration in soil gas and groundwater before, during, and after an earthquake in crystalline rocks depends on the intrinsic radioelement content, the mineral texture, and the mechanical properties of rocks. Rock fracturing and failure do not necessarily induce increase in Rn emission due to these rock properties, which can be used to understand the sensitivity of Rn concentration in soil gas or groundwater in connection with an earthquake.

Description: The Gaussian wavelet is widely used as a shaping wavelet for scattered wave imaging with P receiver functions due to widespread use of the iterative deconvolution method. We show the Gaussian wavelet degrades the resolution of plane wave migration by comparing results from the latest USArray data shaped with Gaussian and Ricker wavelets. We use simulations of primary conversions from the 410 and 660 km discontinuity to show this is a property of the algorithm and not the data. Simulations also show the more conventional common conversion point (CCP) method is not subject to this behaviour for flat horizons, but the CCP method penalizes dipping horizons focusing only nearly horizontal features for any choice of shaping wavelet. We explain these results using the concept of migration impulse response for an individual data sample. Applications to data from USArray show dramatic improvements in the resolution of plane wave migration images produced using Ricker wavelet in comparison to a comparable resolution a Gaussian shaping wavelet. The 410 and 660 discontinuities are resolved to higher precision, and we find the upper mantle and transition zone are full of previously unresolved dipping horizons that remain to be interpreted.

Description: The space–time Epidemic-Type Aftershock Sequence (ETAS) model is extended by incorporating the depth component of earthquake hypocentres. The depths of the direct offspring produced by an earthquake are assumed to be independent of the epicentre locations and to follow a beta distribution, whose shape parameter is determined by the depth of the parent event. This new model is verified by applying it to the Southern California earthquake catalogue. The results show that the new model fits data better than the original epicentre ETAS model and that it provides the potential for modelling and forecasting seismicity with higher resolutions.

Description: We perform a systematic investigation of along-strike rupture directivity of 70 earthquakes (3.0 ≤ M w ≤ 6.1) of the 2009 L'Aquila seismic sequence by analysing azimuthal variations of broad-band seismograms recorded in the distance range 60 km 〈 R 〈 230 km. We use reference spectra of events with little directivity (similar to the empirical Green's function method) to deconvolve propagation-site effects and focus on source properties. A directivity index (0 ≤ I DIR ≤ 1) calculated for each earthquake quantifies the spectral separation above the corner frequency of the target event at opposite along-strike directions. A large number (73 per cent) of events including the M w 6.1 main shock show high (〉0.7) I DIR values indicating predominantly unilateral rupture propagation. The preferred rupture propagation direction is generally to the southeast with no dependence on the earthquake magnitude or occurrence time. Events on two main faults (L'Aquila and Campotosto) show somewhat different behaviour. Almost all earthquakes on the L'Aquila fault have strong unilateral directivity to the southeast, whereas earthquakes on the Campotosto fault show more diverse behaviour. However, there is a predominance of unilateral ruptures (14 out of 22) also on the Campotosto fault, and the few (five) earthquakes with ruptures to the northwest are limited to the most northwestern segment of the fault. The spectral results are consistent with time-domain analysis when the latter samples adequately the frequency band above corner frequency. The preferred rupture direction may be produced at least in part by a velocity contrast across the fault. The results provide important input for estimates of seismic motion and physics of earthquake ruptures.

Description: We show a case study of full-waveform inversion (FWI) applied to a regional 240-km long seismic profile (called POLCRUST-01) located in southeast Poland. The experiment is unique because it portrays different tectonic deformation styles over a relatively short distance. It was designed for deep reflection imaging, but we decided to use FWI to provide a high-resolution P -wave velocity model of the shallow (down to 3 km) subsurface. The acquisition parameters, namely the use of 10-Hz geophones and Vibroseis sweeps starting at 6 Hz, made it necessary to design a non-standard data preconditioning workflow for enhancing low frequencies including match filtering and curvelet denoising. Final model was validated using multiple procedures, such as comparison with borehole data and independently processed Kirchhoff pre-stack depth migration. Although in some deeper areas of the model we identified artefacts that are most likely a result of cycle skip, the velocity profiles from our model match the borehole check-shot data and the velocity anomalies coincide very well with the reflectivity in the migrated data. The resolution of the model allows us to delineate the internal structure of the Carpathian thrust system and presumably small natural gas reservoirs in the Miocene sediments of the Carpathian Foredeep.

Description: Surface waves contain fundamental mode and higher modes, which could interfere with each other. If different modes are not properly separated, the inverted Earth structures using surface waves could be biased. In this study, we apply linear radon transform (LRT) to synthetic seismograms and real seismograms from the USArray to demonstrate the effectiveness of LRT in separating fundamental-mode Love waves from higher modes. Analysis on synthetic seismograms shows that two-station measurements on reconstructed data obtained after mode separation can completely retrieve the fundamental-mode Love-wave phase velocities. Results on USArray data show that higher mode contamination effects reach up to ~10 per cent for two-station measurements of Love waves, while two-station measurements on mode-separated data obtained by LRT are very close to the predicted values from a global dispersion model of GDM52, demonstrating that the contamination of overtones on fundamental-mode Love-wave phase velocity measurements is effectively mitigated by the LRT method and accurate fundamental-mode Love-wave phase velocities can be measured.

Description: Various factors need to be considered when inverting for surface wave azimuthal anisotropic structure. This paper focuses on the 2 terms for Rayleigh wave azimuthal anisotropy and shows that the uncertainties of earthquake locations also have significant impacts on the resulting anisotropic structure. We use the global Rayleigh wave phase velocity data set collected in a previous study to demonstrate this effect. The differences between azimuthal anisotropic patterns with and without source relocations are greatest near plate boundaries. Large differences around the South American plate are also identified. Although most of the earthquakes are shifted by less than 15 km from the CMT locations, earthquakes near the Andes can be systematically shifted by more than 30 km. Our final epicentres for earthquakes on ridge-transform fault systems better match the plate boundaries.

Description: The finite-difference method is among the most popular methods for modelling seismic wave propagation. Although the method has enjoyed huge success for its ability to produce full wavefield seismograms in complex models, it has one major limitation which is of critical importance for many modelling applications; to naturally output up- and downgoing and P - and S -wave constituents of synthesized seismograms. In this paper, we show how such wavefield constituents can be isolated in finite-difference-computed synthetics in complex models with high numerical precision by means of a simple algorithm. The description focuses on up- and downgoing and P - and S -wave separation of data generated using an isotropic elastic finite-difference modelling method. However, the same principles can also be applied to acoustic, electromagnetic and other wave equations.

Description: A second-order Born approximation is used to formulate a self-consistent theory for the effective elastic parameters of stochastic media with ellipsoidal distributions of small-scale heterogeneity. The covariance of the stiffness tensor is represented as the product of a one-point tensor variance and a two-point scalar correlation function with ellipsoidal symmetry, which separates the statistical properties of the local anisotropy from those of the geometric anisotropy. The spatial variations can then be rescaled to an isotropic distribution by a simple metric transformation; the spherical average of the strain Green's function in the transformed space reduces to a constant Kneer tensor, and the second-order corrections to the effective elastic parameters are given by the contraction of the rescaled Kneer tensor against the single-point variance of the stiffness tensor. Explicit results are derived for stochastic models in which the heterogeneity is transversely isotropic and its second moments are characterized by a horizontal-to-vertical aspect ratio . If medium is locally isotropic, the expressions for the anisotropic effective moduli reduce in the limit -〉 to Backus's second-order expressions for a 1-D stochastic laminate. Comparisons with the exact Backus theory show that the second-order approximation predicts the effective anisotropy for non-Gaussian media fairly well for relative rms fluctuations in the moduli smaller than about 30 per cent. A locally anisotropic model is formulated in which the local elastic properties have hexagonal symmetry, guided by a Gaussian random vector field that is transversely isotropic and specified by a horizontal-to-vertical orientation ratio . The self-consistent theory provides closed-form expressions for the dependence of the effective moduli on 0 〈 〈 and 0 〈 〈 . The effective-medium parametrizations described here appear to be suitable for incorporation into tomographic modelling.

Description: Deviation of seismic surface waves from the great-circle between source and receiver is illustrated by the anomalies in the arrival angle, that is the difference between the observed backazimuth of the incident waves and the great-circle. Such arrival angle anomalies have been known for decades, but observations remain scattered. We present a systematic study of arrival angle anomalies of fundamental mode Rayleigh waves (20–100 s period interval) from 289 earthquakes and recorded by a broadband network LAPNET, located in northern Finland. These observations are compared with those of full waveform synthetic seismograms for the same events, calculated in a 3-D Earth and also compared with those of seismograms obtained by ambient noise correlation. The arrival angle anomalies for individual events are complex, and have significant variations with period. On average, the mean absolute deviation decreases from ~9° at 20 s period to ~3° at 100 s period. The synthetic seismograms show the same evolution, albeit with somewhat smaller deviations. While the arrival angle anomalies are fairly well simulated at long periods, the deviations at short periods are very poorly modelled, demonstrating the importance of the continuous improvement of global crustal models. At 20–30 s period, both event data and numerical simulations have strong multipathing, and relative amplitude changes between different waves will induced differences in deviations between very closely located events. The source mechanism has only limited influence on the deviations, demonstrating that they are directly linked to propagation effects, including near-field effects in the source area. This observation is confirmed by the comparison with seismic noise correlation records, that is where the surface waves correspond to those emitted by a point source at the surface, as the two types of observations are remarkably similar in the cases where earthquakes are located close to seismic stations. This agreement additionally confirms that the noise correlations capture the complex surface wave propagation.

Description: We present a method to construct non-stationary and anisotropic second-order random model realizations that can be used for numerical wave propagation simulations in various geometries. Models are generated directly from a given covariance matrix using its eigenvector decomposition (principal component or Karhunen-Loève method). Because this decomposition is very expensive computationally in 3-D, we use model symmetries to reduce the size of the covariance matrix to its non-stationary components. Stationary components can then be described through their power spectrum, such that models with axisymmetric or spherically symmetric statistics can be generated from a 1-D covariance matrix. We focus in particular on models with spherically symmetric statistics that are important to study wave propagation in the Earth. We use this method to show the influence of hypothetical small-scale structure in the Earth's mantle on the elastic wavefield. To this end, we extend tomographic models beyond their spatial resolution limit with different distributions of small-scale scatterers that generate a coda and attenuate direct waves (scattering attenuation). We observe that scattering attenuation of fundamental mode Rayleigh waves is small (0.5–2 per cent of the total attenuation), if the elastic mantle structure does not become significantly stronger at smaller scales. At the examined heterogeneity strengths, scattering attenuation scales linearly with the model variance. The long-period fundamental mode Rayleigh wave coda is difficult to measure because it is weak and overlaps with other signals. However, it can be shown that its intensity also scales linearly with model power, and that it depends strongly on the spherical geometry of the Earth. It can therefore be used to distinguish between models with different small-scale power. We show qualitatively that the coda generated by the type of random models we consider can explain observed scattered energy at long periods (100 s).

Description: We analysed the ground deformation produced by the M w = 6.1 2014 January 26 and M w = 6.0 2014 February 3 Cephalonia earthquakes, western Greece. Campaign GPS measurements and RADARSAT-2 synthetic aperture radar (SAR) interferometry provide constraints on the overall deformation produced by the sequence. TerraSAR-X and COSMO-SkyMed SAR interferometry provide constraints on the second earthquake separately. Two permanent GPS stations captured the two coseismic offsets and show no pre- or post-seismic transients. Most of the deformation is concentrated in the Paliki peninsula which is consistent with the location of the seismicity and the damages. Both GPS and SAR interferometry indicate areas with large deformation gradients probably due to shallow effects. Given the limitations on the data and on the knowledge of the structure and rheology of the crust, we used a simple elastic model to fit the ground displacements. Although such model cannot fit all the detail of the deformation, it is expected to provide a robust estimate of the overall geometry and slip of the fault. The good data coverage in azimuth and distance contributes to the robustness of the model. The entire sequence is modelled with a strike slip fault dipping 70° east and cutting most of the brittle crust beneath Paliki, with an upper edge located at 2.5 km depth and a deeper edge at 8.5 km. This fault is oriented N14° which corresponds to the azimuth of the Cephalonia Transform Fault (CTF). The fit to the data is significantly improved by adding a secondary shallow strike-slip fault with low dip angle (30°) with a component of reverse faulting on that shallow fault. The modelling of the February 3 event indicates that the faulting is shallow in the north of Paliki, with a centroid depth of ~3.2 km. The fit is improved when a single planar fault is replaced by a bent fault dipping ~30° in the uppermost 2 km and ~70° below. The fault of the January 26 earthquake, inferred from the difference between the two above models, is located south and beneath the February 3 fault, with a centroid depth of ~6.4 km. We interpret the 2014 fault zone as an east segment of the CTF located ~7 km east of the main axis of the CTF, which location is constrained by the elastic modelling of the interseismic GPS velocities. The aftershock sequence is mostly located between the January 26 fault and the axis of the CTF. According to our analysis, the Paliki peninsula is partly dragged north with the Apulian platform with ~7 mm yr –1 of shear accommodated offshore to the west. During the last 30 yr three main sequences occurred along the CTF, in 1983, 2003 and 2014 breaking a large part of the fault, with a gap of 20–40 km left between Cephalonia and Lefkada.

Description: The mathematical framework used in the Bayesian Infrasonic Source Localization (BISL) methodology is examined and simplified providing a generalized method of estimating the source location and time for an infrasonic event. The likelihood function describing an infrasonic detection used in BISL has been redefined to include the von Mises distribution developed in directional statistics and propagation-based, physically derived celerity-range and azimuth deviation models. Frameworks for constructing propagation-based celerity-range and azimuth deviation statistics are presented to demonstrate how stochastic propagation modelling methods can be used to improve the precision and accuracy of the posterior probability density function describing the source localization. Infrasonic signals recorded at a number of arrays in the western United States produced by rocket motor detonations at the Utah Test and Training Range are used to demonstrate the application of the new mathematical framework and to quantify the improvement obtained by using the stochastic propagation modelling methods. Using propagation-based priors, the spatial and temporal confidence bounds of the source decreased by more than 40 per cent in all cases and by as much as 80 per cent in one case. Further, the accuracy of the estimates remained high, keeping the ground truth within the 99 per cent confidence bounds for all cases.

Description: The Borborema province in NE Brazil is characterized by seismic sequences with small earthquakes that can last 10 yr or more. The seismicity in this region is concentrated in three main seismic zones. In this work, we investigate the stress field in one of these zones, the Acaraú Seismic Zone, which is located in the NW part of the Borborema province. This seismic zone exhibits earthquake sequences that contain repeated earthquakes with similar waveforms and a shallow depth. Using a local network, we investigated a seismic sequence close to the town of Santana do Acaraú from December 2009 to December 2010, and we present detailed results (velocity model, hypocentres and focal mechanism) from this network. In addition, we inverted seven focal mechanisms, including six that were used in previous studies, and determined the directions of the three main axes of the regional stress field. Selecting a very precise set of 12 earthquakes, we found an active seismic zone with a depth between 3.5 and 4.8 km and with a horizontal dimension of approximately 2.5 km in the NW–SE direction (azimuth of 118°) and a strike-slip focal mechanism. The new seismic fault and some of the previous seismic faults determined in previous studies occur near the continental-scale Transbrasiliano lineament, but they exhibit no direct relationship with that ancient structure. The stress field is characterized by NW–SE trending compression and NE–SW trending extension. This result suggests that the rheological contrast between the continental–oceanic crusts created flexural stresses with maximum horizontal compression parallel to the continental margin. This stress pattern occurs along the Potiguar basin and continues west as far as the Amazon fan along the Equatorial margin of Brazil. This stress field and related seismicity may be a characteristic of this type of passive margin that is generated during the transform shearing between the South America and Africa plates and that exhibits an abrupt oceanic–continent transition, steep continental slopes and high bathymetric gradients.

Description: A new method of fully nonlinear waveform fitting to measure interstation phase speeds and amplitude ratios is developed and applied to USArray. The Neighbourhood Algorithm is used as a global optimizer, which efficiently searches for model parameters that fit two observed waveforms on a common great-circle path by modulating the phase and amplitude terms of the fundamental-mode surface waves. We introduce the reliability parameter that represents how well the waveforms at two stations can be fitted in a time–frequency domain, which is used as a data selection criterion. The method is applied to observed waveforms of USArray for seismic events in the period from 2007 to 2010 with moment magnitude greater than 6.0. We collect a large number of phase speed data (about 75 000 for Rayleigh and 20 000 for Love) and amplitude ratio data (about 15 000 for Rayleigh waves) in a period range from 30 to 130 s. The majority of the interstation distances of measured dispersion data is less than 1000 km, which is much shorter than the typical average path-length of the conventional single-station measurements for source-receiver pairs. The phase speed models for Rayleigh and Love waves show good correlations on large scales with the recent tomographic maps derived from different approaches for phase speed mapping; for example, significant slow anomalies in volcanic regions in the western Unites States and fast anomalies in the cratonic region. Local-scale phase speed anomalies corresponding to the major tectonic features in the western United States, such as Snake River Plains, Basin and Range, Colorado Plateau and Rio Grande Rift have also been identified clearly in the phase speed models. The short-path information derived from our interstation measurements helps to increase the achievable horizontal resolution. We have also performed joint inversions for phase speed maps using the measured phase and amplitude ratio data of vertical component Rayleigh waves. These maps exhibit better recovery of phase speed perturbations, particularly where the strong lateral velocity gradient exists in which the effects of elastic focussing can be significant; that is, the Yellowstone hotspot, Snake River Plains, and Rio Grande Rift. The enhanced resolution of the phase speed models derived from the interstation phase and amplitude measurements will be of use for the better seismological constraint on the lithospheric structure, in combination with dense broad-band seismic arrays.

Description: Motivated by the need to derive and characterize increasingly sophisticated seismic data analysis and inversion methods incorporating wave dissipation, we consider the problem of scattering of homogeneous and inhomogeneous waves from perturbations in five viscoelastic parameters (density, P - and S -wave velocities, and P - and S -wave quality factors), as formulated in the context of the Born approximation. Within this approximation the total wave field is the superposition of an incident plane wave and a scattered wave, the latter being a spherical wave weighted by a function of solid angle called the scattering potential. In elastic media the scattering potential is real, but if dissipation is included through a viscoelastic model, the potential becomes complex and thus impacts the amplitude and phase of the outgoing wave. The isotropic-elastic scattering framework of Stolt and Weglein, extended to admit viscoelastic media, exposes these amplitude and phase phenomena to study, and in particular allows certain well-known layered-medium viscoelastic results due to Borcherdt to be re-considered in an arbitrary heterogeneous Earth. The main theoretical challenge in doing this involves the choice of coordinate system over which to evaluate and analyse the waves, which in the viscoelastic case must be based on complex vector analysis. We present a candidate system within which several of Borcherdt's key results carry over; for instance, we show that elliptically polarized P and SI waves cannot be scattered into linearly polarized SII waves. Furthermore, the elastic formulation is straightforwardly recovered in the limit as P - and S -wave quality factors tend to infinity.

Description: We present an extension of the nodal discontinuous Galerkin method for elastic wave propagation to high interpolation orders and arbitrary heterogeneous media. The high-order lagrangian interpolation is based on a set of nodes with excellent interpolation properties in the standard triangular element. In order to take into account highly variable geological media, another set of suitable quadrature points is used where the physical and mechanical properties of the medium are defined. We implement the methodology in a 2-D discontinuous Galerkin solver. First, a convergence study confirms the hp -convergence of the method in a smoothly varying elastic medium. Then, we show the advantages of the present methodology, compared to the classical one with constant properties within the elements, in terms of the complexity of the mesh generation process by analysing the seismic amplification of a soft layer over an elastic half-space. Finally, to verify the proposed methodology in a more complex and realistic configuration, we compare the simulation results with the ones obtained by the spectral element method for a sedimentary basin with a realistic gradient velocity profile. Satisfactory results are obtained even for the case where the computational mesh does not honour the strong impedance contrast between the basin bottom and the bedrock.

Description: To obtain the synthetic seismogram using the Cagniard-de Hoop method, one needs to calculate the integral over slowness. When the source is shallow and the slowness is near the zero of the Rayleigh function, the integrand behaves like a sharp pulse. In this study, we attempt to study this pulse with an asymptotic approach, and conclude that the Rayleigh wave in the time domain originates from this pulse in the slowness domain. We therefore offer an explanation of the excitation of the Rayleigh wave in a mathematical point of view. In addition, we propose a method to improve the efficiency of the numerical quadrature in the calculation of the synthetic seismogram.

Description: Many studies have sought to seismically image plumes rising from the deep mantle in order to settle the debate about their presence and role in mantle dynamics, yet the predicted seismic signature of realistic plumes remains poorly understood. By combining numerical simulations of flow, mineral-physics constraints on the relationships between thermal anomalies and wave speeds, and spectral-element method based computations of seismograms, we estimate the delay times of teleseismic S and P waves caused by thermal plumes. Wave front healing is incomplete for seismic periods ranging from 10 s (relevant in traveltime tomography) to 40 s (relevant in waveform tomography). We estimate P -wave delays to be immeasurably small (〈0.3 s). S -wave delays are larger than 0.4 s even for S waves crossing the conduits of the thinnest thermal plumes in our geodynamic models. At longer periods (〉20 s), measurements of instantaneous phase misfit may be more useful in resolving narrow plume conduits. To detect S -wave delays of 0.4–0.8 s and the diagnostic frequency dependence imparted by plumes, it is key to minimize the influence of the heterogeneous crust and upper mantle. We argue that seismic imaging of plumes will advance significantly if data from wide-aperture ocean-bottom networks were available since, compared to continents, the oceanic crust and upper mantle are relatively simple.

Description: Field experiments are used to unequivocally demonstrate seismic superresolution imaging of subwavelength objects in the near-field region of the source. The field test is for a conventional hammer source striking a metal plate near subwavelength scatterers and the seismic data are recorded by vertical-component geophones in the far-field locations of the sources. Time-reversal mirrors (TRMs) are then used to refocus the scattered energy with subwavelength resolution to the position of the original source. A spatial resolution of /10, where is the dominant wavelength associated with the data, is seen in the field tests that exceeds the Abbe resolution limit of /2.

Description: We propose an adaptive root-determining strategy that is very useful when dealing with trapped modes or Stoneley modes whose energies become very insignificant on the free surface in the presence of low-velocity layers or fluid layers in the model. Loss of modes in these cases or inaccuracy in the calculation of these modes may then be easily avoided. Built upon the generalized reflection/transmission coefficients, the concept of ‘family of secular functions’ that we herein call ‘adaptive mode observers’ is thus naturally introduced to implement this strategy, the underlying idea of which has been distinctly noted for the first time and may be generalized to other applications such as free oscillations or applied to other methods in use when these cases are encountered. Additionally, we have made further improvements upon the generalized reflection/transmission coefficient method; mode observers associated with only the free surface and low-velocity layers (and the fluid/solid interface if the model contains fluid layers) are adequate to guarantee no loss and high precision at the same time of any physically existent modes without excessive calculations. Finally, the conventional definition of the fundamental mode is reconsidered, which is entailed in the cases under study. Some computational aspects are remarked on. With the additional help afforded by our superior root-searching scheme and the possibility of speeding calculation using a less number of layers aided by the concept of ‘turning point’, our algorithm is remarkably efficient as well as stable and accurate and can be used as a powerful tool for widely related applications.

Description: Numerical solvers of wave equations have been widely used to simulate global seismic waves including PP waves for modelling 410/660 km discontinuity and Rayleigh waves for imaging crustal structure. In order to avoid extra computation cost due to ocean water effects, these numerical solvers usually adopt water column approximation, whose accuracy depends on frequency and needs to be investigated quantitatively. In this paper, we describe a unified representation of accurate and approximate forms of the equivalent water column boundary condition as well as the free boundary condition. Then we derive an analytical form of the PP -wave reflection coefficient with the unified boundary condition, and quantify the effects of water column approximation on amplitude and phase shift of the PP waves. We also study the effects of water column approximation on phase velocity dispersion of the fundamental mode Rayleigh wave with a propagation matrix method. We find that with the water column approximation: (1) The error of PP amplitude and phase shift is less than 5 per cent and 9° at periods greater than 25 s for most oceanic regions. But at periods of 15 s or less, PP is inaccurate up to 10 per cent in amplitude and a few seconds in time shift for deep oceans. (2) The error in Rayleigh wave phase velocity is less than 1 per cent at periods greater than 30 s in most oceanic regions, but the error is up to 2 per cent for deep oceans at periods of 20 s or less. This study confirms that the water column approximation is only accurate at long periods and it needs to be improved at shorter periods.

Description: Seismic wave resonance in sedimentary basins is a well-recognized seismic hazard; however, concentrated areas of earthquake damage have been observed near basin edges, where wave propagation is particularly complex and difficult to understand with sparse observations. The Tokyo metropolitan area is densely populated, subject to strong shaking from a diversity of earthquake sources, and sits atop the deep Kanto sedimentary basin. It is also instrumented with two seismic arrays: the dense MEtropolitan Seismic Observation network (MeSO-net) within the basin, and the High sensitivity seismograph network (Hi-net) surrounding it. In this study, we explore the 3-D seismic wavefield within and throughout the Kanto basin, including near and across basin boundaries, using cross-correlations of all components of ambient seismic field between the stations of these two arrays. Dense observations allow us to observe clearly the propagation of three modes of both Rayleigh and Love waves. They also show how the wavefield behaves in the vicinity of sharp basin edges with reflected/converted waves and excitation of higher modes.

Description: We present an analytical approach to compute the curvature effect by the new analytical solutions of coseismic deformation derived for the homogeneous sphere model. We consider two spheres with different radii: one is the same as earth and the other with a larger radius can approximate a half-space model. Then, we calculate the coseismic displacements for the two spheres and define the relative percentage of the displacements as the curvature effect. The near-field curvature effect is defined relative to the maximum coseismic displacement. The results show that the maximum curvature effect is about 4 per cent for source depths of less than 100 km, and about 30 per cent for source depths of less than 600 km. For the far-field curvature effect, we define it relative to the observing point. The curvature effect is extremely large and sometimes exceeds 100 per cent. Moreover, this new approach can be used to estimate any planet's curvature effect quantitatively. For a smaller sphere, such as the Moon, the curvature effect is much larger than that of the Earth, with an inverse ratio to the earth's radius.

Description: In Antarctica, locally grounded ice, such as ice rises bordering floating ice shelves, plays a major role in the ice mass balance as it stabilizes the ice sheet flow from the hinterland. When in direct contact with the ocean, the ice rise buttressing effect may be altered in response of changing ocean forcing. To investigate this vulnerable zone, four sites near the boundary of an ice shelf with an ice rise promontory in Dronning Maud Land, East-Antarctica were monitored for a month in early 2014 with new instruments that include both seismic and GPS sensors. Our study indicated that this transition zone experiences periodic seismic activity resulting from surface crevassing during oceanic tide-induced flexure of the ice shelf. The most significant finding is the observation of apparent fortnightly tide-modulated low-frequency, long-duration seismic events at the seaward front of the ice rise promontory. A basal origin of these events is postulated with the ocean water surge at each new spring tide triggering basal crevassing or basal slip on a local bedrock asperity. Detection and monitoring of such seismicity may help identifying ice rise zones vulnerable to intensified ocean forcing.

Description: SKS arrivals from ocean bottom seismometer (OBS) data from an offshore southern California deployment are analysed for shear wave splitting. The project involved 34 OBSs deployed for 12 months in a region extending up to 500 km west of the coastline into the oceanic Pacific plate. The measurement process consisted of removing the effects of anisotropy using a range of values for splitting fast directions and delay times to minimize energy along the transverse seismometer axis. Computed splitting parameters are unexpectedly similar to onland parameters, exhibiting WSW–ENE fast polarization directions and delays between 0.8 and 1.8 s, even for oceanic plate sites. This is the first SKS splitting study to extend across the entire boundary between the North America and Pacific plates, into the oceanic part of the Pacific plate. The splitting results show that the fast direction of anisotropy on the Pacific plate does not align with absolute plate motion (APM), and they extend the trend of anisotropy in southern California an additional 500 km west, well onto the oceanic Pacific plate. We model the finite strain and anisotropy within the asthenosphere associated with density–buoyancy driven mantle flow and the effects of APM. In the absence of plate motion effects, such buoyancy driven mantle flow would be NE-directed beneath the Pacific plate observations. The best-fit patterns of mantle flow are inferred from the tomography-based models that show primary influences from foundering higher-density zones associated with the history of subduction beneath North America. The new offshore SKS measurements, when combined with measurements onshore within the plate boundary zone, indicate that dramatic lateral variations in density-driven upper-mantle flow are required from offshore California into the plate boundary zone in California and western Basin and Range.

Description: The islands composing the Maltese archipelago (Central Mediterranean) are characterized by a four-layer sequence of limestones and clays. A common feature found in the western half of the archipelago is Upper Coralline Limestone (UCL) plateaus and hillcaps covering a soft Blue Clay (BC) layer which can be up to 75 m thick. The BC layer introduces a velocity inversion in the stratigraphy, implying that the V S 30 (traveltime average sear wave velocity ( V S ) in the upper 30 m) parameter is not always suitable for seismic microzonation purposes. Such a layer may produce amplification effects, however might not be included in the V S 30 calculations. In this investigation, V S profiles at seven sites characterized by such a lithological sequence are obtained by a joint inversion of the single-station Horizontal-to-Vertical Spectral Ratios (H/V or HVSR) and effective dispersion curves from array measurements analysed using the Extended Spatial Auto-Correlation technique. The lithological sequence gives rise to a ubiquitous H/V peak between 1 and 2 Hz. All the effective dispersion curves obtained exhibit a ‘normal’ dispersive trend at low frequencies, followed by an inverse dispersive trend at higher frequencies. This shape is tentatively explained in terms of the presence of higher mode Rayleigh waves, which are commonly present in such scenarios. Comparisons made with the results obtained at the only site in Malta where the BC is missing below the UCL suggest that the characteristics observed at the other seven sites are due to the presence of the soft layer. The final profiles reveal a variation in the V S of the clay layer with respect to the depth of burial and some regional variations in the UCL layer. This study presents a step towards a holistic seismic risk assessment that includes the implications on the site effects induced by the buried clay layer. Such assessments have not yet been done for Malta.

Description: Noise is a persistent feature in seismic data and so poses challenges in extracting increased accuracy in seismic images and physical interpretation of the subsurface. In this paper, we analyse passive seismic data from the Aquistore carbon capture and storage pilot project permanent seismic array to characterise, classify and model seismic noise. We perform noise analysis for a three-month subset of passive seismic data from the array and provide conclusive evidence that the noise field is not white, stationary, or Gaussian; characteristics commonly yet erroneously assumed in most conventional noise models. We introduce a novel noise modelling method that provides a significantly more accurate characterisation of real seismic noise compared to conventional methods, which is quantified using the Mann–Whitney–White statistical test. This method is based on a statistical covariance modelling approach created through the modelling of individual noise signals. The identification of individual noise signals, broadly classified as stationary, pseudo-stationary and non-stationary, provides a basis on which to build an appropriate spatial and temporal noise field model. Furthermore, we have developed a workflow to incorporate realistic noise models within synthetic seismic data sets providing an opportunity to test and analyse detection and imaging algorithms under realistic noise conditions.

Description: Plate-scale deformation is expected to impart seismic anisotropic fabrics on the lithosphere. Determination of the fast shear wave orientation ( ) and the delay time between the fast and slow split shear waves ( t ) via SKS splitting can help place spatial and temporal constraints on lithospheric deformation. The Canadian Appalachians experienced multiple episodes of deformation during the Phanerozoic: accretionary collisions during the Palaeozoic prior to the collision between Laurentia and Gondwana, and rifting related to the Mesozoic opening of the North Atlantic. However, the extent to which extensional events have overprinted older orogenic trends is uncertain. We address this issue through measurements of seismic anisotropy beneath the Canadian Appalachians, computing shear wave splitting parameters ( , t ) for new and existing seismic stations in Nova Scotia and New Brunswick. Average t values of 1.2 s, relatively short length scale (≥100 km) splitting parameter variations, and a lack of correlation with absolute plate motion direction and mantle flow models, demonstrate that fossil lithospheric anisotropic fabrics dominate our results. Most fast directions parallel Appalachian orogenic trends observed at the surface, while t values point towards coherent deformation of the crust and mantle lithosphere. Mesozoic rifting had minimal impact on our study area, except locally within the Bay of Fundy and in southern Nova Scotia, where fast directions are subparallel to the opening direction of Mesozoic rifting; associated t values of 〉1 s require an anisotropic layer that spans both the crust and mantle, meaning the formation of the Bay of Fundy was not merely a thin-skinned tectonic event.

Description: A new method for discrete sampling of signals is presented with specific applications to the reconstruction of recorded interfering wavefields from two or more sources excited simultaneously at discrete positions along lines. By utilizing a periodic sequence of source signatures along one of the source lines, the corresponding wavefield becomes separately visible in a part of the spectral domain where it can be isolated, processed and subtracted from the interfering wavefields. As a result, interfering wavefields from multiple sources recorded at a single location can be fully separated from each other. The concept is referred to as signal apparition which we suggest refers to ‘the act of becoming visible’. It may find applications in a wide range of disciplines relying on wave experimentation, such as acoustic, seismic and electromagnetic imaging of the Earth's interior for instance to significantly enhance resolution of subsurface images.

Description: In the context of seismic monitoring, recent studies made successful use of seismic coda waves to locate medium changes on the horizontal plane. Locating the depth of the changes, however, remains a challenge. In this paper, we use 3-D wavefield simulations to address two problems: first, we evaluate the contribution of surface- and body-wave sensitivity to a change at depth. We introduce a thin layer with a perturbed velocity at different depths and measure the apparent relative velocity changes due to this layer at different times in the coda and for different degrees of heterogeneity of the model. We show that the depth sensitivity can be modelled as a linear combination of body- and surface-wave sensitivity. The lapse-time-dependent sensitivity ratio of body waves and surface waves can be used to build 3-D sensitivity kernels for imaging purposes. Second, we compare the lapse-time behaviour in the presence of a perturbation in horizontal and vertical slabs to address, for instance, the origin of the velocity changes detected after large earthquakes.

Description: Traveltime tomography is the main method by which the Earth's seismic velocity is determined on all scales, from the near-surface (〈100 m) to the core. Usually traveltime tomography uses ray theory, an infinite-frequency approximation of wave propagation. A theory developed in global seismology to account for the finite-frequency nature of seismic data, known as finite-frequency traveltime tomography (FFTT), can theoretically provide a more accurate estimation of velocity. But the FFTT theory is generally not applicable to near-surface data because there is no reference velocity model known in advance that is capable of yielding synthetic waveforms that are close enough to the recorded seismograms to yield a reliable delay time. Also, there is usually no reference model for which the unknown velocity model represents a small (linear) perturbation from the reference model. This paper presents a frequency dependent form of non-linear traveltime tomography specifically designed for near-surface seismic data in which a starting model, iterative approach with recalculated travel paths at each iteration, and the calculation of a frequency-dependent total traveltime, as opposed to a delay time, are used. Frequency-dependent traveltime tomography (FDTT) involves two modifications to conventional traveltime tomography: (1) the calculation of frequency-dependent traveltimes using wavelength-dependent velocity smoothing (WDVS) and (2) the corresponding sensitivity kernels that arise from using WDVS. Results show that the former modification is essential to achieve significant benefits from FDTT, whereas the latter is optional in that similar results can be achieved using infinite-frequency kernels. The long seismic wavelengths relative to the total path lengths and the size of subsurface heterogeneities of typical near-surface data means the improvements over ray theory tomography are significant. The benefits of FDTT are demonstrated using conventional minimum-structure regularization techniques to address the issue of model non-uniqueness. For synthetic data, the estimated FDTT models are shown to be more accurate than the corresponding infinite-frequency-derived models. Both 2-D and 3-D applications of FDTT to real data from a near-surface study yield estimated models that contain more structure than the corresponding infinite-frequency-derived models. Applications of FDTT without regularization demonstrate the potential of the WDVS-derived sensitivity kernels to provide a natural smoothing of the velocity model and thereby allow the data alone to determine the final model structure.

Description: We analysed a catalogue of Italian earthquakes, covering 55 yr of data from 1960 to 2014 with magnitudes homogeneously converted to M w , to compute the time-dependent relative frequencies with which strong seismic shocks (4.0 ≤ M w 〈 5.0), widely felt by the population, have been followed by main shocks ( M w ≥ 5.0) that threatened the health and the properties of the persons living in the epicentral area. Assuming the stationarity of the seismic release properties, such frequencies are estimates of the probabilities of potentially destructive shocks after the occurrence of future strong shocks. We compared them with the time-independent relative frequencies of random occurrence in terms of the frequency gain that is the ratio between the time-dependent and time-independent relative frequencies. The time-dependent relative frequencies vary from less than 1 per cent to about 20 per cent, depending on the magnitudes of the shocks and the time windows considered (ranging from minutes to years). They remain almost constant for a few hours after the strong shock and then decrease with time logarithmically. Strong earthquakes (with M w ≥ 6.0) mainly occurred within two or three months of the strong shock. The frequency gains vary from about 10 000 for very short time intervals to less than 10 for a time interval of 2 yr. Only about 1/3 of main shocks were preceded by at least a strong shock in the previous day and about 1/2 in the previous month.

Description: We present the first high-resolution Rayleigh-wave phase-velocity azimuthal anisotropy tomography of the Japan subduction zone at periods of 20–150 s, which is determined using a large number of high-quality amplitude and phase data of teleseismic fundamental-mode Rayleigh waves. The obtained 2-D anisotropic phase-velocity models are then inverted for a 3-D shear-wave velocity azimuthal anisotropy tomography down to a depth of ~300 km beneath Japan. The subducting Pacific slab is imaged as a dipping high-velocity zone with trench-parallel fast-velocity directions (FVDs) which may indicate the anisotropy arising from the normal faults produced at the outer-rise area near the Japan trench axis, overprinting the slab fossil fabric, whereas the mantle wedge generally exhibits lower velocities with trench-normal FVDs which reflect subduction-driven corner flow and anisotropy. Depth variations of azimuthal anisotropy are revealed in the big mantle wedge beneath the Japan Sea, which may reflect past deformations in the Eurasian lithosphere related to backarc spreading during 21 to 15 Ma and complex current convection in the asthenosphere induced by active subductions of both the Pacific and Philippine Sea plates.

Description: Markov chain Monte Carlo sampling methods are widely used for non-linear Bayesian inversion where no analytical expression for the forward relation between data and model parameters is available. Contrary to the linear(ized) approaches, they naturally allow to evaluate the uncertainties on the model found. Nevertheless their use is problematic in high-dimensional model spaces especially when the computational cost of the forward problem is significant and/or the a posteriori distribution is multimodal. In this case, the chain can stay stuck in one of the modes and hence not provide an exhaustive sampling of the distribution of interest. We present here a still relatively unknown algorithm that allows interaction between several Markov chains at different temperatures. These interactions (based on importance resampling) ensure a robust sampling of any posterior distribution and thus provide a way to efficiently tackle complex fully non-linear inverse problems. The algorithm is easy to implement and is well adapted to run on parallel supercomputers. In this paper, the algorithm is first introduced and applied to a synthetic multimodal distribution in order to demonstrate its robustness and efficiency compared to a simulated annealing method. It is then applied in the framework of first arrival traveltime seismic tomography on real data recorded in the context of hydraulic fracturing. To carry out this study a wavelet-based adaptive model parametrization has been used. This allows to integrate the a priori information provided by sonic logs and to reduce optimally the dimension of the problem.

Description: Improving both accuracy and computational performance of numerical tools is a major challenge for seismic imaging and generally requires specialized implementations to make full use of modern parallel architectures. We present a computational strategy for reverse-time migration (RTM) with accelerator-aided clusters. A new imaging condition computed from the pressure and velocity fields is introduced. The model solver is based on a high-order discontinuous Galerkin time-domain (DGTD) method for the pressure–velocity system with unstructured meshes and multirate local time stepping. We adopted the MPI+X approach for distributed programming where X is a threaded programming model. In this work we chose OCCA, a unified framework that makes use of major multithreading languages (e.g. CUDA and OpenCL) and offers the flexibility to run on several hardware architectures. DGTD schemes are suitable for efficient computations with accelerators thanks to localized element-to-element coupling and the dense algebraic operations required for each element. Moreover, compared to high-order finite-difference schemes, the thin halo inherent to DGTD method reduces the amount of data to be exchanged between MPI processes and storage requirements for RTM procedures. The amount of data to be recorded during simulation is reduced by storing only boundary values in memory rather than on disk and recreating the forward wavefields. Computational results are presented that indicate that these methods are strong scalable up to at least 32 GPUs for a three-dimensional RTM case.

Description: The study of acoustic emissions (AEs) is of paramount importance to understand rock deformation processes. AE recorded during laboratory experiments mimics, in a controlled geometry and environment, natural and induced seismicity. However, these experiments are destructive, time consuming and require a significant amount of resources. Lately, significant progresses have been made in numerical simulations of rock failure processes, providing detailed insights into AE. We utilized the 2-D combined finite-discrete element method to simulate the deformation of Stanstead Granite under varying confining pressure ( P c ) and demonstrated that the increase of confining pressure, P c , (i) shifts failures from tensile towards shear dominated and (ii) enhance the macroscopic ductility. We quantitatively describe the AE activity associated with the fracturing process by assessing the spatial fractal dimension (D-value), the temporal distribution (AE rate) and the slope of the frequency–magnitude distribution (b-value). Based on the evaluation of D-value and AE rate, we defined two distinct deformation phases: Phase I and Phase II. The influence of P c on the spatial distribution of AE varies according to the deformation phase: for increasing P c , D-value decreases and increases during Phases I and II, respectively. In addition, b-value decreases with increasing P c during the entire experiment. Our numerical results show for the first time that variations of D- and b-values as a function of in situ stress can be simulated using the combined finite-discrete element approach. We demonstrate that the examination of seismicity should be carried out carefully, taking into consideration the deformation phase and in situ stress conditions.

Description: We present an approach for subspace detection of small seismic events that includes methods for estimating magnitudes and associating detections from multiple stations into unique events. The process is used to identify mining related seismicity from a surface coal mine and an underground coal mining district, both located in the Western U.S. Using a blasting log and a locally derived seismic catalogue as ground truth, we assess detector performance in terms of verified detections, false positives and failed detections. We are able to correctly identify over 95 per cent of the surface coal mine blasts and about 33 per cent of the events from the underground mining district, while keeping the number of potential false positives relatively low by requiring all detections to occur on two stations. We find that most of the potential false detections for the underground coal district are genuine events missed by the local seismic network, demonstrating the usefulness of regional subspace detectors in augmenting local catalogues. We note a trade-off in detection performance between stations at smaller source–receiver distances, which have increased signal-to-noise ratio, and stations at larger distances, which have greater waveform similarity. We also explore the increased detection capabilities of a single higher dimension subspace detector, compared to multiple lower dimension detectors, in identifying events that can be described as linear combinations of training events. We find, in our data set, that such an advantage can be significant, justifying the use of a subspace detection scheme over conventional correlation methods.

Description: Near-surface seismic field experiments using both P - and S -wave sources were carried out to image fractured limestones at two sites in southwest England. We measured P - and S -wave seismic velocities at multiple azimuths to aligned fracture sets, allowing us to determine the seismic anisotropy generated by these fractures. The effect of aligned fractures on seismic anisotropy is commonly modelled in terms of the additional compliance introduced by the fractures. Therefore, an understanding of fracture compliance is crucial both in terms of interpreting observations of anisotropy in the field and in forward modelling the effects of fractures on seismic wave propagation. Of particular concern is (1) the scaling of fracture compliance with fracture length scale, and (2) the controls on the ratio of normal to tangential compliance of the fractures (  =  Z N /Z T ). Our experimental design allows us to image both, and we find that  = 0.37 ± 0.06 and  = 0.75 ± 0.10 for our two study sites, while the absolute values of the tangential compliance range from 0.66  x  10 –11 to 5.0  x  10 –11 Pa –1  m.

Description: Global-scale tomographic models should aim at satisfying the full seismic spectrum. For this purpose, and to better constrain isotropic 3-D variations of shear velocities in the mantle, we tackle a joint inversion of spheroidal normal-mode structure coefficients and multiple-frequency S -wave delay times. In all previous studies for which normal modes were jointly inverted for, with body and/or surface waves, the mantle was laterally parametrized with uniform basis functions, such as spherical harmonics, equal-area blocks and evenly spaced spherical splines. In particular, spherical harmonics naturally appear when considering the Earth's free oscillations. However, progress towards higher resolution joint tomography requires a movement away from such uniform parametrization to overcome its computational inefficiency to adapt to local variations in resolution. The main goal of this study is to include normal modes into a joint inversion based upon a non-uniform parametrization that is adapted to the spatially varying smallest resolving length of the data. Thus, we perform the first joint inversion of normal-mode and body-wave data using an irregular tomographic grid, optimized according to ray density. We show how to compute the projection of 3-D sensitivity kernels for both data sets onto our parametrization made up of spherical layers spanned with irregular Delaunay triangulations. This approach, computationally efficient, allows us to map into the joint model multiscale structural informations from data including periods in the 10–51 s range for body waves and 332–2134 s for normal modes. Tomographic results are focused on the 400–2110 km depth range, where our data coverage is the most relevant. We discuss the potential of a better resolution where the grid is fine, compared to spherical harmonics up to degree 40, as the number of model parameters is similar. Our joint model seems to contain coherent structural components beyond degree 40, such as those related to the Farallon subduction. Assessing their robustness is postponed to a future work. A wider application of this tomographic workflow, holding promise to better understand mantle dynamics at various spatial scales, should primarily consist in adding surface-wave data and extending our sets of normal-mode and body-wave data.

Description: The Shanxi rift system is one of the most active intraplate tectonic zones in the North China Block, resulting in devastating seismicity. Since 1303, the rift has experienced fifteen M s ≥ 6.5 earthquakes. Aiming at a better understanding of Coulomb stress evolution and its relationship with the seismicity in the rift system, we investigated the Coulomb stress changes due to coseismic slip and post-seismic relaxation processes following strong earthquakes as well as the interseismic tectonic loading since the 1303 Hongdong M s = 8.0 earthquake. Our calculation applies a specified regional medium model, takes the gravity effect into account and uses the fault geometry of the next event as the receiver fault in a given calculation. Our results show that nine out of 12 M s ≥ 6.5 earthquakes since the 1303 Hongdong earthquake and more than 82 per cent of small-medium instrumental events after the 1989 Datong-Yanggao M s = 6.1 earthquake fall into the total stress increased areas. Our results also reveal the different roles of the coseismic, post-seismic and interseismic Coulomb stress changes in the earthquake triggering process in the Shanxi rift system. In a short period after a strong event, the stress field changes are dominated by coseismic Coulomb stress due to sudden slip of the ruptured fault, while in the long term, the stress field is mainly dominated by the accumulation of interseismic tectonic loading. Post-seismic stress changes play an important role by further modifying the distribution of stress and therefore cannot be ignored. Based on the current stress status in the Shanxi rift system, the Linfen basin, southern and northern Taiyuan basin, Xinding basin and the north part of the rift system are identified as the most likely locations of large events in the future. The results of this study can provide important clues for the further understanding of seismic hazard in the Shanxi rift system and thus help guiding earthquake risk mitigation efforts in the future.

Description: We present spectra concentrating on the lowest-frequency normal modes of the Earth obtained from records of the invar-wire strainmeters and STS-1 broad-band seismometers located in the Black Forest Observatory, Germany after the disastrous earthquakes off the NW coast of Sumatra in 2004 and off the coast near Tohoku, Japan in 2011. We compare the spectra to ones obtained from synthetic seismograms computed using a mode summation technique for an anelastic, elliptical, rotating, spherically symmetric Earth model. The synthetics include strain–strain-coupling effects by using coupling coefficients obtained from comparisons between Earth tide signals recorded by the strainmeters and synthetic tidal records. We show that for the low-frequency toroidal and spheroidal modes up to 1 mHz, the strainmeters produce better signal-to-noise ratios than the broad-band horizontal seismometers. Overall, the comparison with the synthetics is satisfactory but not as good as for vertical accelerations. In particular, we demonstrate the high quality of the strainmeter data by showing the Coriolis splitting of toroidal modes for the first time in individual records, the first clear observation of the singlet $_2S_1^0$ and the detection of the fundamental radial mode 0 S 0 with good signal-to-noise ratio and with a strain amplitude of 10 –11 . We also identify the latter mode in a record of the Isabella strainmeter after the great Chilean quake in 1960, the detection of which was missed by the original studies.

Description: Cross-correlation of ambient seismic noise recorded by two seismic stations may result in an estimate of the Green's function between those two receivers. Several authors have recently attempted to measure attenuation based on these interferometric, receiver–receiver surface waves. By now, however, it is well established that the loss of coherence of the cross-correlation as a function of space depends strongly on the excitation of the medium. In fact, in a homogeneous dissipative medium, uniform excitation is required to correctly recover attenuation. Applied to fundamental-mode ambient seismic surface waves, this implies that the cross-correlation will decay at the local attenuation rate only if noise sources are distributed uniformly on the Earth's surface. In this study we show that this constraint can be relaxed in case the observed loss of coherence is due to multiple scattering instead of dissipation of energy. We describe the scattering medium as an effective medium whose phase velocity and rate of attenuation are a function of the scatterer density and the average strength of the scatterers. We find that the decay of the cross-correlation in the effective medium coincides with the local attenuation of the effective medium in case the scattering medium is illuminated uniformly from all angles. Consequently, uniform excitation is not a necessary condition for the correct retrieval of scattering attenuation. We exemplify the implications of this finding for studies using the spectrally whitened cross-correlation to infer subsurface attenuation.

Description: This paper presents direct-seismogram inversion (DSI) for receiver-side structure which treats the source signal incident from below (the effective source–time function—STF) as a vector of unknown parameters in a Bayesian framework. As a result, the DSI method developed here does not require deconvolution by observed seismogram components as typically applied in receiver-function inversion and avoids the problematic issue of choosing subjective tuning parameters in this deconvolution. This results in more meaningful inversion results and uncertainty estimation compared to classic receiver-function inversion. A rigorous derivation is presented of the likelihood function required for unbiased inversion results. The STF is efficiently inferred by a maximum-likelihood closed-form expression that does not require deconvolution by noisy waveforms. Rather, deconvolution is only by predicted impulse responses for the unknown environment (considered to be a 1-D, horizontally stratified medium). For a given realization of the parameter vector which describes the medium below the station, data predictions are computed as the convolution of the impulse response and the maximum-likelihood source estimate for that medium. Therefore, the assumption of a Gaussian pulse with specified parameters, typical for the prediction of receiver functions, is not required. Directly inverting seismogram components has important consequences for the noise on the data. Since the signal processing does not require filtering and deconvolution, data errors are less correlated and more straightforward to model than those for receiver functions. This results in better inversion results (parameter values and uncertainties), since assumptions made in the derivation of the likelihood function are more likely to be met by the inversion process. The DSI method is demonstrated for simulated waveforms and then applied to data for station Hyderabad on the Indian craton. The measured data are inverted with both the new DSI and traditional receiver-function inversion. All inversions are carried out for a trans-dimensional model that treats the number of layers in the model as unknown. Results for DSI are consistent with previous studies for the same location. The DSI has clear advantages in trans-dimensional inversion. Uncertainty estimates appear more realistic (larger) in both model complexity (number of layers) and in terms of seismic velocity profiles. Receiver-function inversion results in more complex profiles (highly-layered structure) and suggests unreasonably small uncertainties. This effect is likely also significant when the parametrization is considered to be fixed but exacerbated for the trans-dimensional model: If hierarchical errors are poorly estimated, trans-dimensional models overestimate the structure which produces unfavourable results for the receiver-function inversion.

Description: Ultra-low velocity zones (ULVZs) are small-scale structures in the Earth's lowermost mantle inferred from the analysis of seismological observations. These structures exhibit a strong decrease in compressional ( P )-wave velocity, shear ( S )-wave velocity, and an increase in density. Quantifying the elastic properties of ULVZs is crucial for understanding their physical origin, which has been hypothesized either as partial melting, iron enrichment, or a combination of the two. Possible disambiguation of these hypotheses can lead to a better understanding of the dynamic processes of the lowermost mantle, such as, percolation, stirring and thermochemical convection. To date, ULVZs have been predominantly studied by forward waveform modelling of seismic waves that sample the core–mantle boundary region. However, ULVZ parameters (i.e. velocity, density, and vertical and lateral extent) obtained through forward modelling are poorly constrained because inferring Earth structure from seismic observations is a non-linear inverse problem with inherent non-uniqueness. To address these issues, we developed a trans-dimensional hierarchical Bayesian inversion that enables rigorous estimation of ULVZ parameter values and their uncertainties, including the effects of model selection. The model selection includes treating the number of layers and the vertical extent of the ULVZ as unknowns. The posterior probability density (solution to the inverse problem) of the ULVZ parameters is estimated by reversible jump Markov chain Monte Carlo sampling that employs parallel tempering to improve efficiency/convergence. First, we apply our method to study the resolution of complex ULVZ structure (including gradually varying structure) by probabilistically inverting simulated noisy waveforms. Then, two data sets sampling the CMB beneath the Philippine and Tasman Seas are considered in the inversion. Our results indicate that both ULVZs are more complex than previously suggested. For the Philippine Sea data, we find a strong decrease in S -wave velocity, which indicates the presence of iron-rich material, albeit this result is accompanied with larger parameter uncertainties than in a previous study. For the Tasman Sea data, our analysis yields a well-constrained S -wave velocity that gradually decreases with depth. We conclude that this ULVZ represents a partial melt of iron-enriched material with higher melt content near its bottom.

Description: In this study, frequency-dependent seismic scattering and intrinsic attenuation parameters for the crustal structure beneath the W-Bohemia/Vogtland swarm earthquake region close to the border of Czech Republic and Germany are estimated. Synthetic seismogram envelopes are modelled using elastic and acoustic radiative transfer theory. Scattering and absorption parameters are determined by fitting these synthetic envelopes to observed seismogram envelopes from 14 shallow local events from the October 2008 W-Bohemia/Vogtland earthquake swarm. The two different simulation approaches yield similar results for the estimated crustal parameters and show a comparable frequency dependence of both transport mean free path and intrinsic absorption path length. Both methods suggest that intrinsic attenuation is dominant over scattering attenuation in the W-Bohemia/Vogtland region for the investigated epicentral distance range and frequency bands from 3 to 24 Hz. Elastic simulations of seismogram envelopes suggest that forward scattering is required to explain the data, however, the degree of forward scattering is not resolvable. Errors in the parameter estimation are smaller in the elastic case compared to results from the acoustic simulations. The frequency decay of the transport mean free path suggests a random medium described by a nearly exponential autocorrelation function. The fluctuation strength and correlation length of the random medium cannot be estimated independently, but only a combination of the parameters related to the transport mean free path of the medium can be computed. Furthermore, our elastic simulations show, that using our numerical method, it is not possible to resolve the value of the mean free path of the random medium.

Description: In this manuscript, we introduce a framework for developing earthquake forecasts using Virtual Quake (VQ), the generalized successor to the perhaps better known Virtual California (VC) earthquake simulator. We discuss the basic merits and mechanics of the simulator, and we present several statistics of interest for earthquake forecasting. We also show that, though the system as a whole (in aggregate) behaves quite randomly, (simulated) earthquake sequences limited to specific fault sections exhibit measurable predictability in the form of increasing seismicity precursory to large m 〉 7 earthquakes. In order to quantify this, we develop an alert-based forecasting metric, and show that it exhibits significant information gain compared to random forecasts. We also discuss the long-standing question of activation versus quiescent type earthquake triggering. We show that VQ exhibits both behaviours separately for independent fault sections; some fault sections exhibit activation type triggering, while others are better characterized by quiescent type triggering. We discuss these aspects of VQ specifically with respect to faults in the Salton Basin and near the El Mayor-Cucapah region in southern California, USA and northern Baja California Norte, Mexico.

Description: Mud volcanoes are geological systems often characterized by elevated fluid pressures at depth deviating from hydrostatic conditions. This near-critical state makes mud volcanoes particularly sensitive to external forcing induced by natural or man-made perturbations. We used the Nirano mud volcanic field as a natural laboratory to test pre- and post-seismic effects generated by distant earthquakes. We first characterized the subsurface structure of the Nirano mud volcanic field with a geoelectrical study. Next, we deployed a broad-band seismic station in the area to understand the typical seismic signal generated by the mud volcano. Seismic records show a background noise below 2 s, sometimes interrupted by pulses of drumbeat-like high-frequency signals lasting from several minutes to hours. To date this is the first observation of drumbeat signal observed in mud volcanoes. In 2013 June we recorded a M4.7 earthquake, that occurred approximately 60 km far from our seismic station. According to empirical estimations the Nirano mud volcanic field should not have been affected by the M4.7 earthquake. Yet, before the seismic event we recorded an increasing amplitude of the signal in the 10–20 Hz frequency band. The signal emerged approximately two hours before the earthquake and lasted for about three hours. Our statistical analysis suggests the presence of a possible precursory signal about 10 min before the earthquake.

Description: Assessing the similarities of seismic attributes can help identify tremor, low signal-to-noise (S/N) signals and converted or reflected phases, in addition to diagnosing site noise and sensor misalignment in arrays. Polarization analysis is a widely accepted method for studying the orientation and directional characteristics of seismic phases via computed attributes, but similarity is ordinarily discussed using qualitative comparisons with reference values or known seismic sources. Here we introduce a technique for quantitative polarization similarity that uses weighted histograms computed in short, overlapping time windows, drawing on methods adapted from the image processing and computer vision literature. Our method accounts for ambiguity in azimuth and incidence angle and variations in S/N ratio. Measuring polarization similarity allows easy identification of site noise and sensor misalignment and can help identify coherent noise and emergent or low S/N phase arrivals. Dissimilar azimuths during phase arrivals indicate misaligned horizontal components, dissimilar incidence angles during phase arrivals indicate misaligned vertical components and dissimilar linear polarization may indicate a secondary noise source. Using records of the M w = 8.3 Sea of Okhotsk earthquake, from Canadian National Seismic Network broad-band sensors in British Columbia and Yukon Territory, Canada, and a vertical borehole array at Hoadley gas field, central Alberta, Canada, we demonstrate that our method is robust to station spacing. Discrete wavelet analysis extends polarization similarity to the time–frequency domain in a straightforward way. Time–frequency polarization similarities of borehole data suggest that a coherent noise source may have persisted above 8 Hz several months after peak resource extraction from a ‘flowback’ type hydraulic fracture.

Description: Differential traveltimes of P4KP and PKP waves are used in an attempt to constrain the topography of the core–mantle boundary (CMB) in localized areas. For epicentral distances of approximately 150°, PKP and P4KP follow the same path through the crust and mantle and differ only in the core, since P4KP reflects three times at the underside of the CMB. Hence traveltime differences should only be due to the path differences in the core. We use array techniques in order to enhance the signal to noise ratio and confirm the slowness and backazimuth of the P4KP waves. Data from the Gräfenberg Array in southern Germany yield a total of 27 events with observable P4KP waves. We find most traveltime residuals between theoretical and observed P4KP minus PKP arrival time generally within a range of ±7 s. A ray tracing method is used to model the deviations of traveltimes due to changes of CMB depth at the reflection points of the P4KP waves. It is not possible to explain the observed traveltime residuals with changes in depth of the CMB only; in several cases, the residuals are too large to be caused by CMB topography alone. Due to the expected coinciding ray paths of PKP and P4KP, the remaining variations are unlikely to be produced by heterogeneities in the lowermost mantle, and so remain for the moment unaccounted for. We find that the reflection coefficient of P4KP undergoes a polarity reversal for some scenarios of lower mantle velocities and densities and the observed P4KP polarity at this distance range can therefore help to investigate the lowermost mantle in addition to CMB topography.

Description: We present the new model SP12RTS of isotropic shear-wave ( V S ) and compressional-wave ( V P ) velocity variations in the Earth's mantle. SP12RTS is derived using the same methods as employed in the construction of the shear-wave velocity models S20RTS and S40RTS, and the same data types. SP12RTS includes additional traveltime measurements of P -waves and new splitting measurements: 33 normal modes with sensitivity to the compressional-wave velocity and 9 Stoneley modes with sensitivity primarily to the lowermost mantle. Contrary to S20RTS and S40RTS, variations in V S and V P are determined without invoking scaling relationships. Lateral velocity variations in SP12RTS are parametrised using spherical harmonics up to degree 12, to focus on long-wavelength features of V S and V P and their ratio R . Large-low-velocity provinces (LLVPs) are observed for both V S and V P . SP12RTS also features an increase of R up to 2500 km depth, followed by a decrease towards the core–mantle boundary. A negative correlation between the shear-wave and bulk-sound velocity variations is observed for both the LLVPs and the surrounding mantle. These characteristics can be explained by the presence of post-perovskite or large-scale chemical heterogeneity in the lower mantle.

Description: We present a new conceptual approach to scattering-integral-based seismic full waveform inversion (FWI) that allows a flexible, extendable, modular and both computationally and storage-efficient numerical implementation. To achieve maximum modularity and extendability, interactions between the three fundamental steps carried out sequentially in each iteration of the inversion procedure, namely, solving the forward problem, computing waveform sensitivity kernels and deriving a model update, are kept at an absolute minimum and are implemented by dedicated interfaces. To realize storage efficiency and maximum flexibility, the spatial discretization of the inverted earth model is allowed to be completely independent of the spatial discretization employed by the forward solver. For computational efficiency reasons, the inversion is done in the frequency domain. The benefits of our approach are as follows: (1) Each of the three stages of an iteration is realized by a stand-alone software program. In this way, we avoid the monolithic, unflexible and hard-to-modify codes that have often been written for solving inverse problems. (2) The solution of the forward problem, required for kernel computation, can be obtained by any wave propagation modelling code giving users maximum flexibility in choosing the forward modelling method. Both time-domain and frequency-domain approaches can be used. (3) Forward solvers typically demand spatial discretizations that are significantly denser than actually desired for the inverted model. Exploiting this fact by pre-integrating the kernels allows a dramatic reduction of disk space and makes kernel storage feasible. No assumptions are made on the spatial discretization scheme employed by the forward solver. (4) In addition, working in the frequency domain effectively reduces the amount of data, the number of kernels to be computed and the number of equations to be solved. (5) Updating the model by solving a large equation system can be done using different mathematical approaches. Since kernels are stored on disk, it can be repeated many times for different regularization parameters without need to solve the forward problem, making the approach accessible to Occam's method. Changes of choice of misfit functional, weighting of data and selection of data subsets are still possible at this stage. We have coded our approach to FWI into a program package called ASKI (Analysis of Sensitivity and Kernel Inversion) which can be applied to inverse problems at various spatial scales in both Cartesian and spherical geometries. It is written in modern FORTRAN language using object-oriented concepts that reflect the modular structure of the inversion procedure. We validate our FWI method by a small-scale synthetic study and present first results of its application to high-quality seismological data acquired in the southern Aegean.

Description: This paper describes a generalization of the iterative deconvolution method commonly used as a component of passive array wavefield imaging. We show that the iterative method should be thought of as a sparse output deconvolution method with the number of terms retained dependent on the convergence criteria. The generalized method we introduce uses an inverse operator to shape the assumed wavelet to a peaked function at zero lag. We show that the conventional method is equivalent to using a damped least-squares spiking filter with extremely large damping and proper scaling. In that case, the inverse operator used in the generalized method reduces to the cross-correlation operator. The theoretical insight of realizing the output is a sparse series provides a basis for the second important addition of the generalized method—an output shaping wavelet. A constant output shaping wavelet is a critical component in scattered wave imaging to avoid mixing data of variable bandwidth. We demonstrate the new approach can improve resolution by using an inverse operator tuned to maximize resolution. We also show that the signal-to-noise ratio of the result can be improved by applying a different convergence criterion than the standard method, which measures the energy left after each iteration. The efficacy of the approach was evaluated with synthetic experiment in various signal and noise conditions. We further validated the approach with real data from the USArray. We compared our results with data from the EarthScope Automated Receiver Survey and found that our results show modest improvements in consistency measured by correlation coefficients with station stacks and a reduced number of outliers.

Description: Time-reversal modelling provides a simple and robust solution to source-imaging problems. However, for recovering a well-resolved image of the source, time-reversal requires a balanced illumination of the target from all angles. When acquisition is incomplete and a balanced illumination is not possible, the time-reversal solution may not be adequate. We present an inversion algorithm for computing the signals to be back propagated by an array of receivers with a given configuration in order to optimally image an unknown source in an elastic medium. This approach is based on minimizing the difference between the back-propagated wavefield and the time-reversed displacement field of the source in the near source region. The proposed method requires knowledge of the propagation medium and an estimate of the source location. Other information related to the source (e.g. source mechanism) is encoded within the data and not explicitly required by the method. The method is applicable in both elastic and acoustic media. We use synthetic examples to test the performance of our method, to analyse its sensitivities, potentials and limitations and to demonstrate its potential advantages over the time-reversal source-imaging technique.

Description: Using data from more than 2000 seismic stations from multiple networks arrayed throughout China (CEArray, China Array, NECESS, PASSCAL, GSN) and surrounding regions (Korean Seismic Network, F-Net, KNET), we perform ambient noise Rayleigh wave tomography across the entire region and earthquake tomography across parts of South China and Northeast China. We produce isotropic Rayleigh wave group and phase speed maps with uncertainty estimates from 8 to 50 s period across the entire region of study, and extend them to 70 s period where earthquake tomography is performed. Maps of azimuthal anisotropy are estimated simultaneously to minimize anisotropic bias in the isotropic maps, but are not discussed here. The 3D model is produced using a Bayesian Monte Carlo formalism covering all of China, extending eastwards through the Korean Peninsula, into the marginal seas, to Japan. We define the final model as the mean and standard deviation of the posterior distribution at each location on a 0.5° x 0.5° grid from the surface to 150 km depth. Surface wave dispersion data do not strongly constrain internal interfaces, but shear wave speeds between the discontinuities in the crystalline crust and uppermost mantle are well determined. We design the resulting model as a reference model, which is intended to be useful to other researchers as a starting model, to predict seismic wave fields and observables and to predict other types of data (e.g. topography, gravity). The model and the data on which it is based are available for download. In addition, the model displays a great variety and considerable richness of geological and tectonic features in the crust and in the uppermost mantle deserving of further focus and continued interpretation.

Description: Earthquake source inversion is highly dependent on location determination and velocity models. Uncertainties in both the model parameters and the observations need to be rigorously incorporated into an inversion approach. Here, we show a probabilistic Bayesian method that allows formal inclusion of the uncertainties in the moment tensor inversion. This method allows the combination of different sets of far-field observations, such as P -wave and S -wave polarities and amplitude ratios, into one inversion. Additional observations can be included by deriving a suitable likelihood function from the uncertainties. This inversion produces samples from the source posterior probability distribution, including a best-fitting solution for the source mechanism and associated probability. The inversion can be constrained to the double-couple space or allowed to explore the gamut of moment tensor solutions, allowing volumetric and other non-double-couple components. The posterior probability of the double-couple and full moment tensor source models can be evaluated from the Bayesian evidence, using samples from the likelihood distributions for the two source models, producing an estimate of whether or not a source is double-couple. Such an approach is ideally suited to microseismic studies where there are many sources of uncertainty and it is often difficult to produce reliability estimates of the source mechanism, although this can be true of many other cases. Using full-waveform synthetic seismograms, we also show the effects of noise, location, network distribution and velocity model uncertainty on the source probability density function. The noise has the largest effect on the results, especially as it can affect other parts of the event processing. This uncertainty can lead to erroneous non-double-couple source probability distributions, even when no other uncertainties exist. Although including amplitude ratios can improve the constraint on the source probability distribution, the measurements are often systematically affected by noise, leading to deviation from their noise-free true values and consequently adversely affecting the source probability distribution, especially for the full moment tensor model. As an example of the application of this method, four events from the Krafla volcano in Iceland are inverted, which show clear differentiation between non-double-couple and double-couple sources, reflected in the posterior probability distributions for the source models.

Description: Extracting accurate empirical Green's functions from the ambient seismic noise field requires the noise to be fully diffuse and that different frequency components are not correlated. Calculating a matrix of correlation coefficients of power spectral samples can be used to estimate deviations from a fully diffuse random noise field in the analysed frequency range. A fully diffuse field has correlations only in a narrow region around the diagonal of the matrix, with frequency resolution inversely proportional to length of the used time window. Analysis of low-frequency data (0.005–0.6 Hz) recorded by three broad-band stations of the southern California seismic network reveals three common types of correlations, manifested in the correlation coefficient matrix as square, diagonal halo and correlated stripes. Synthetic calculations show that these types of signatures in the correlation coefficient matrix can result from certain combinations of cross-frequency correlated random components and diffuse field. The analysis of observed data indicates that the secondary microseismic peak around 0.15 Hz is correlated with its neighbouring frequencies, while the primary peak around 0.06 Hz is more diffuse. This suggests that the primary and secondary peaks may be associated with somewhat different physical origins. In addition, significant correlation of frequencies below that of the primary microseismic peak suggests that the very low frequencies noise is less scattered during propagation. The power spectra recorded by a station close to the edge of the Los Angeles basin is higher compared to data recorded by stations outside the basin perhaps because of enhanced basin reverberations and/or closer proximity to the ocean. This and other regional variations should be tested further using data from many more stations.

Description: An improved description of scattering and inverse scattering processes in reflection seismology may be obtained on the basis of a scattering series solution to the Helmoltz equation, which allows one to separately model primary and multiple reflections. However, the popular scattering series of Born is of limited seismic modelling value, since it is only guaranteed to converge if the global contrast is relatively small. For frequency-domain waveform modelling of realistic contrasts, some kind of renormalization may be required. The concept of renormalization is normally associated with quantum field theory, where it is absolutely essential for the treatment of infinities in connection with observable quantities. However, the renormalization program is also highly relevant for classical systems, especially when there are interaction effects that act across different length scales. In the scattering series of De Wolf, a renormalization of the Green's functions is achieved by a split of the scattering potential operator into fore- and backscattering parts; which leads to an effective reorganization and partially re-summation of the different terms in the Born series, so that their order better reflects the physics of reflection seismology. It has been demonstrated that the leading (single return) term in the De Wolf series (DWS) gives much more accurate results than the corresponding Born approximation, especially for models with high contrasts that lead to a large accumulation of phase changes in the forward direction. However, the higher order terms in the DWS that are associated with internal multiples have not been studied numerically before. In this paper, we report from a systematic numerical investigation of the convergence properties of the DWS which is based on two new operator representations of the DWS. The first operator representation is relatively similar to the original scattering potential formulation, but more global and explicit in nature. The second representation is based on the T -operator formalism from quantum scattering theory, that offers a different perspective on the interaction between up- and downgoing waves, as well as significant computational advantages (e.g. domain decomposition and fast recursive methods for one-way propagators). Our numerical results demonstrate the convergence properties of the DWS are indeed superior to those of the Born series.

Description: In this study, we develop an efficient boundary integral equation method for estimation of seismic motion in a graded medium with multiple cavities under antiplane strain conditions. This inhomogeneous and heterogeneous medium is subjected to either time-harmonic incident shear seismic waves or to body waves radiating from a point seismic source. Three different types of soil material gradient are considered: (i) density and shear modulus vary proportionally as quadratic functions of depth, but the wave velocity remains constant; (ii) the soil material is viscoelastic, with a shear modulus and density that vary with respect to the spatial coordinates in an arbitrary fashion, so that the wave velocity is both frequency and position-dependent and (iii) the soil material has position-dependent shear modulus and constant density, yielding a linear profile for the wave velocity. Three different, frequency-dependent boundary integral equation schemes are respectively developed for the aforementioned three types of graded soil materials based on: (i) Green's function for the quadratically graded elastic half-plane; (ii) a fundamental solution for the viscoelastic full-plane with position-dependent wave speed profiles and (iii) a fundamental solution for an elastic full-plane with a linearly varying wave speed profile. Next, a number of cases involving geological media with position-dependent material properties and any number of cavities of various shapes and geometry are solved in the frequency domain. The numerical results reveal the dependency of the wave fields and zones of stress concentration on the following key factors: (i) type and properties of the soil material gradient; (ii) type and characteristics of the applied seismic load; (iii) shape, position and number of cavities and (iv) interaction phenomena between the cavities and the free surface.

Description: Seismic moment tensor is one of the most important source parameters defining the earthquake dimension and style of the activated fault. Geoscientists ordinarily use moment tensor catalogues, however, few attempts have been done to assess possible impacts of moment magnitude uncertainties upon their analysis. The 2012 May 20 Emilia main shock is a representative event since it is defined in literature with a moment magnitude value ( M w ) spanning between 5.63 and 6.12. A variability of ~0.5 units in magnitude leads to a controversial knowledge of the real size of the event and reveals how the solutions could be poorly constrained. In this work, we investigate the stability of the moment tensor solution for this earthquake, studying the effect of five different 1-D velocity models, the number and the distribution of the stations used in the inversion procedure. We also introduce a 3-D velocity model to account for structural heterogeneity. We finally estimate the uncertainties associated to the computed focal planes and the obtained M w . We conclude that our reliable source solutions provide a moment magnitude that ranges from 5.87, 1-D model, to 5.96, 3-D model, reducing the variability of the literature to ~0.1. We endorse that the estimate of seismic moment from moment tensor solutions, as well as the estimate of the other kinematic source parameters, requires coming out with disclosed assumptions and explicit processing workflows. Finally and, probably more important, when moment tensor solution is used for secondary analyses it has to be combined with the same main boundary conditions (e.g. wave-velocity propagation model) to avoid conflicting results.

Description: Trade-offs between velocity and anisotropy heterogeneity complicate the interpretation of differential traveltime data and have the potential to bias isotropic tomographic models. By constructing a simple parametrisation to describe an elastic tensor with hexagonal symmetry, we find analytic solutions to the Christoffel equations in terms of fast and slow horizontal velocities that allow us to simultaneously invert differential traveltime data and splitting data from teleseismic S arrivals to recover 3-D velocity and anisotropy structure. This technique provides a constraint on the depth-extent of shallow anisotropy, otherwise absent from interpretations based on SKS splitting alone. This approach is well suited to the young Woodlark Rift, where previous studies have found strong velocity variation and substantial SKS splitting in a continental rift with relatively simple geometry. This study images a low-velocity rift axis with ≤4 per cent spreading-parallel anisotropy at 50–100 km depth that separates regions of pre-existing lithospheric fabric, indicating the synchronous development of extensional crystallographic preferred orientation and lithospheric thinning. A high-velocity slab fragment north of the rift axis is associated with strike-parallel anisotropic fast axes, similar to that seen in the shallow mantle of some subduction zones. In addition to the insights provided by the anisotropy structure, the improvement in fit to the differential traveltime data demonstrates the merit to a joint inversion that accounts for anisotropy.

Description: Short-period small magnitude seismograms mainly comprise scattered waves in the form of coda waves (the tail part of the seismogram, starting after S waves and ending when the noise prevails), spanning more than 70 per cent of the whole seismogram duration. Corresponding coda envelopes provide important information about the earth inhomogeneity, which can be stochastically modeled in terms of distribution of scatterers in a random medium. In suitable experimental conditions (i.e. high earth heterogeneity), either the two parameters describing heterogeneity (scattering coefficient), intrinsic energy dissipation (coefficient of intrinsic attenuation) or a combination of them (extinction length and seismic albedo) can be used to image Earth structures. Once a set of such parameter couples has been measured in a given area and for a number of sources and receivers, imaging their space distribution with standard methods is straightforward. However, as for finite-frequency and full-waveform tomography, the essential problem for a correct imaging is the determination of the weighting function describing the spatial sensitivity of observable data to scattering and absorption anomalies. Due to the nature of coda waves, the measured parameter couple can be seen as a weighted space average of the real parameters characterizing the rock volumes illuminated by the scattered waves. This paper uses the Monte Carlo numerical solution of the Energy Transport Equation to find approximate but realistic 2-D space-weighting functions for coda waves. Separate images for scattering and absorption based on these sensitivity functions are then compared with those obtained with commonly used sensitivity functions in an application to data from an active seismic experiment carried out at Deception Island (Antarctica). Results show that these novel functions are based on a reliable and physically grounded method to image magnitude and shape of scattering and absorption anomalies. Their extension to 3-D holds promise to improve our ability to model volcanic structures using coda waves.

Description: We investigate the relationship between seismic moment M 0 and source duration t w of microearthquakes by using high-quality seismic data recorded with a vertical borehole array installed in central Taiwan. We apply a waveform cross-correlation method to the three-component records and identify several event clusters with high waveform similarity, with event magnitudes ranging from 0.3 to 2.0. Three clusters—Clusters A, B and C—contain 11, 8 and 6 events with similar waveforms, respectively. To determine how M 0 scales with t w , we remove path effects by using a path-averaged Q . The results indicate a nearly constant t w for events within each cluster, regardless of M 0 , with mean values of t w being 0.058, 0.056 and 0.034 s for Clusters A, B and C, respectively. Constant t w , independent of M 0 , violates the commonly used scaling relation ${t_w} \propto M_0^{1/3}$ . This constant duration may arise either because all events in a cluster are hosted on the same isolated seismogenic patch, or because the events are driven by external factors of constant duration, such as fluid injections into the fault zone. It may also be related to the earthquake nucleation size.

Description: Earthquake absolute location errors which can be encountered in an underground reservoir are investigated. In such an exploitation context, earthquake hypocentre errors can have an impact on the field development and economic consequences. The approach using the state-of-the-art techniques covers both the location uncertainty and the location inaccuracy—or bias—problematics. It consists, first, in creating a 3-D synthetic seismic cloud of events in the reservoir and calculating the seismic traveltimes to a monitoring network assuming certain propagation conditions. In a second phase, the earthquakes are relocated with assumptions different from the initial conditions. Finally, the initial and relocated hypocentres are compared. As a result, location errors driven by the seismic onset time picking uncertainties and inaccuracies are quantified in 3-D. Effects induced by erroneous assumptions associated with the velocity model are also modelled. In particular, 1-D velocity model uncertainties, a local 3-D perturbation of the velocity and a 3-D geostructural model are considered. The present approach is applied to the site of Rittershoffen (Alsace, France), which is one of the deep geothermal fields existing in the Upper Rhine Graben. This example allows setting realistic scenarios based on the knowledge of the site. In that case, the zone of interest, monitored by an existing seismic network, ranges between 1 and 5 km depth in a radius of 2 km around a geothermal well. Well log data provided a reference 1-D velocity model used for the synthetic earthquake relocation. The 3-D analysis highlights the role played by the seismic network coverage and the velocity model in the amplitude and orientation of the location uncertainties and inaccuracies at subsurface levels. The location errors are neither isotropic nor aleatoric in the zone of interest. This suggests that although location inaccuracies may be smaller than location uncertainties, both quantities can have a cumulative effect. Besides, small velocity uncertainties applied to the whole 1-D profile can lead to large increase of the location uncertainties. However, local variations of the velocity field around the well may have negligible effects that would make such a feature undetectable with an absolute location method. Although the reference 1-D velocity model was built from well log data, the results show that it is not a good representative of a more realistic 3-D model including a fault and its associated block shift. The amplitude and distribution of the induced location inaccuracies are such that the positioning and the orientation of features delineated by seismicity are distorted and may be difficult to correctly interpret.

Description: Compared with a single GPS system, GPS/GLONASS observations can improve the satellite visibility, optimize the spatial geometry and improve the precise positioning performance. Although the advantage over GPS-only methods in terms of positioning is clear, the potential contributions of GPS/GLONASS to co-seismic displacement determination and the subsequent seismic source inversion still require extensive study and validation. In this paper, we first extended a temporal point positioning model from GPS-only to GPS/GLONASS observations. Using this new model, the performance of the GPS/GLONASS method for obtaining co-seismic displacements was then validated via eight outdoor experiments on a shaking table. Our result reveals that the GPS/GLONASS method provides more accurate and robust co-seismic displacements than the GPS-only observations in a non-optimal observation environment. Furthermore, as a case study, observation data recorded during the September 2015 M w 8.3 Illapel earthquake in Chile were re-processed. At some stations, obvious biases were found between the co-seismic displacements derived from GPS-only and GPS/GLONASS observations. The subsequent slip distribution inversion on a curved fault confirms that the differences in the co-seismic displacements causes differences in the inversion results and that the slip distributions of the Illapel earthquake inferred from the GPS/GLONASS observations tend to be shallower and larger.

Description: Incorporation of attenuation into the normal mode sum representations of seismic signals is commonly effected by applying perturbation theory. This is fine for weak attenuation, but problematic for stronger attenuation. In this work, modes of the anelastic medium are represented as complex superpositions of elastic eigenfunctions. For the P – SV system, a generalized eigenvalue equation for the complex eigenwavenumbers and complex coefficients used to construct the anelastic eigenfunctions is derived. The generalized eigenvalue problem for the P – SV problem is exactly linear in the eigenwavenumber at the expense of doubling the dimension. The SH problem is exactly linear in the square of the eigenwavenumber. This is in contrast to a similar standing wave problem for the earth free oscillations. Attenuation is commonly incorporated into synthetic seismogram calculations by introduction of complex frequency-dependent elastic moduli. The moduli depend nonlinearly on the frequency. The independent variable in the standing wave free oscillation problem is the frequency, which makes the eigenvalue problem nonlinear. The choice of the wavenumber as the independent variable for the travelling wave problem leads to a linear problem. The Earth model may be transversely isotropic. Compressional waves and both polarizations of shear waves ( SV, SH ) are treated.

Description: We present a new technique for deriving detailed information on seismic velocities of the subsurface material from continuous ambient noise recorded by spatially dense seismic arrays. This method uses iterative double beamforming between various subarrays to extract surface wave contributions from the ambient-noise data in complex environments with unfavourable noise-source distributions. The iterative double beamforming extraction makes it possible to retrieve large amounts of Rayleigh wave traveltime information in a wide frequency band. The method is applied to data recorded by a highly dense Nodal array with 1108 vertical geophones, centred on the damage zone of the Clark branch of the San Jacinto Fault Zone south of Anza, California. The array covers a region of ~650 x 700 m 2 , with instrument spacing of 10–30 m, and continuous recording at 500 samples s –1 over 30 d in 2014. Using this iterative double beamforming on subarrays of 25 sensors and cross-correlations between all of the station pairs, we separate surface waves from body waves that are abundant in the raw cross-correlation data. Focusing solely on surface waves, maps of traveltimes are obtained at different frequencies with unprecedented accuracy at each point of a 15-m-spacing grid. Group velocity inversions at 2–4 Hz reveal depth and lateral variations in the structural properties within and around the San Jacinto Fault Zone in the study area. This method can be used over wider frequency ranges and can be combined with other imaging techniques, such as eikonal tomography, to provide unprecedented detailed structural images of the subsurface material.

Description: Knowledge of the quality factor of near-surface materials is of fundamental interest in various applications. Attenuation can be very strong close to the surface and thus needs to be properly assessed. In recent years, several researchers have studied the retrieval of attenuation coefficients from the cross correlation of ambient seismic noise. Yet, the determination of exact amplitude information from noise-correlation functions is, in contrast to the extraction of traveltimes, not trivial. Most of the studies estimated attenuation coefficients on the regional scale and within the microseism band. In this paper, we investigate the possibility to derive attenuation coefficients from seismic noise at much shallower depths and higher frequencies (〉1 Hz). The Euroseistest area in northern Greece offers ideal conditions to study quality factor retrieval from ambient noise for different rock types. Correlations are computed between the stations of a small scale array experiment (station spacings 〈2 km) that was carried out in the Euroseistest area in 2011. We employ the correlation of the coda of the correlation (C 3 ) method instead of simple cross correlations to mitigate the effect of uneven noise source distributions on the correlation amplitude. Transient removal and temporal flattening are applied instead of 1-bit normalization in order to retain relative amplitudes. The C 3 method leads to improved correlation results (higher signal-to-noise ratio and improved time symmetry) compared to simple cross correlations. The C 3 functions are rotated from the ZNE to the ZRT system and we focus on Love wave arrivals on the transverse component and on Love wave quality factors Q L . The analysis is performed for selected stations being either situated on soft soil or on weathered rock. Phase slowness is extracted using a slant-stack method. Attenuation parameters are inferred by inspecting the relative amplitude decay of Love waves with increasing interstation distance. We observe that the attenuation coefficient and Q L can be reliably extracted for stations situated on soft soil whereas the derivation of attenuation parameters is more problematic for stations that are located on weathered rock. The results are in acceptable conformance with theoretical Love wave attenuation curves that were computed using 1-D shear wave velocity and quality factor profiles from the Euroseistest area.

Description: Large earthquakes from the intermediate-depth Vrancea seismic zone are known to produce in Bucharest ground motion characterized by predominant long periods. This phenomenon has been interpreted as the combined effect of both seismic source properties and site response of the large sedimentary basin. The thickness of the unconsolidated Quaternary deposits beneath the city is more than 200 m, the total depth of sediments is more than 1000 m. Complex basin geometry and the low seismic wave velocities of the sediments are primarily responsible for the large amplification and long duration experienced during earthquakes. For a better understanding of the geological structure under Bucharest, a number of investigations using non-invasive methods have been carried out. With the goal to analyse and extract the polarization and dispersion characteristics of the surface waves, ambient vibrations and low-magnitude earthquakes have been investigated using single station and array techniques. Love and Rayleigh dispersion curves (including higher modes), Rayleigh waves ellipticity and SH -wave fundamental frequency of resonance ( f 0 SH ) have been inverted simultaneously to estimate the shear wave velocity structure under Bucharest down to a depth of about 8 km. Information from existing borehole logs was used as prior to reduce the non-uniqueness of the inversion and to constrain the shallow part of the velocity model (〈300 m). In this study, we use data from a 35-km diameter array (the URS experiment) installed by the National Institute for Earth Physics and by the Karlsruhe Institute of Technology during 10 months in the period 2003–2004. The array consisted of 32 three-component seismological stations, deployed in the urban area of Bucharest and adjacent zones. The large size of the array and the broad-band nature of the available sensors gave us the possibility to characterize the surface wave dispersion at very low frequencies (0.05–1 Hz) using frequency–wavenumber techniques. This is essential to explore and resolve the deeper portions of the basin. The horizontal to vertical spectral ratio (H/V) curves provide important additional information about the structure and are here characterized by two major peaks. The first is attributed to the fundamental frequency of the basin, while the second can be interpreted as a mixture of the second higher mode of Rayleigh waves and other types of waves such as SH waves. This hypothesis has been verified by comparing the H/V curves with the SH -wave transfer function from the retrieved velocity structure. We could also approximate the SH transfer function with H/V ratios of earthquake recordings, providing additional verification of the robustness of the proposed velocity model. The Cretaceous bedrock depth was then inverted at each URS station from the fundamental frequency of resonance and using this model. A 3-D geophysical model for Bucharest has been constructed based on the integration of the inverted velocity profiles and the available geological information using a geographic information system.

Description: Some modern seismicity in the magnitude range of M 4 and M 6 in California and eastern North America preferentially occurs at the edges of past large ruptures. Once a large earthquake rupture has occurred, stress is concentrated at the edges of the rupture, and apparently this stress concentration can trigger earthquakes at or near the rupture edges many decades or even longer after a main shock. Furthermore, the modern M ≥ 4 earthquakes in the vicinity of a past main shock usually have the same focal mechanism as the earlier main shock. There are a number of examples of this in California and Nevada, where there is a statistically significant correlation of the locations of M ≥ 4 earthquakes and the edges of 19th and 20th century fault ruptures in M w ≥ 6.5 earthquakes. In contrast, the M ≥ 4 earthquakes near the epicentres of future ruptures in California are randomly scattered around the fault with no concentration near the ends of the future fault rupture. The concentration of earthquakes near the ends of earlier large ruptures in California becomes progressively less pronounced as the smallest magnitude in the data set is reduced from M 4.0 to M 3.0. These observations also appear to be true for intraplate regions where aftershock sequences can last millennia. The identification of modern rupture-edge M ≥ 4 aftershocks can be used to help discover where and when past strong earthquakes took place, even if there is no historical record of the main shock. This is of great importance for seismic hazard studies.

Description: Surface wave methods provide a cost effective means of developing shear wave velocity ( Vs ) profiles for applications such as dynamic site characterization and seismic site response analyses. However, the inverse problem involved in obtaining a realistic layered earth model from surface wave dispersion data is inherently ill-posed, non-linear and mix-determined, without a unique solution. When available, a priori information such as geotechnical boreholes or geologic well logs should be used to aid in constraining site-specific inversion parameters. Unfortunately, a priori information is often unavailable, particularly at significant depths, and a ‘blind analysis’ must be performed. In these situations, the analyst must decide on an appropriate number of layers and ranges for their corresponding inversion parameters (i.e. trial number of layers and ranges in their respective thicknesses, shear wave velocities, compression wave velocities and mass densities). Selection of these parameters has been shown to significantly impact the results of an inversion. This paper presents a method for conducting multiple inversions utilizing systematically varied inversion layering parametrizations in order to identify and encompass the most reasonable layered earth models for a site. Each parametrization is defined by a unique layering ratio, which represents a multiplier that systemically increases the potential thickness of each layer in the inversion parametrization based on the potential thickness of the layer directly above it. The layering ratio method is demonstrated at two sites associated with the InterPacific Project, wherein it is shown to significantly aid in selecting reasonable Vs profiles that are close representations of the subsurface. While the goal of the layering ratio inversion methodology is not necessarily to find the ‘optimal’ or ‘best’ Vs profile for a site, it may be successful at doing so for certain sites/datasets. However, the primary reason for using the layering ratio method is to find Vs profiles that realistically represent the uncertainty in Vs resulting from surface wave inversion, and to avoid selection of Vs profiles that are unrealistic and adversely influenced by the choice of inversion parametrization.

Description: The appraisal of tomographic models, of fundamental importance towards better understanding the Earth's interior, consists in analysing their resolution and covariance. The discrete theory of Backus–Gilbert, solving all at once the linear problems of model estimation and appraisal, aims at evaluating weighted averages of the true model parameters. Contrary to damped least-squares techniques, one key advantage of Backus–Gilbert inversion is that no subjective regularization is needed to remove the non-uniqueness of the model solution. Indeed, it is often possible to identify unique linear combinations of the parameters even when the parameters themselves are not uniquely defined. In other words, the non-uniqueness can be broken by averaging rather than regularizing. Over the past few decades, many authors have considered that, in addition to a high computational cost, it could be a clumsy affair in the presence of data errors to practically implement the Backus–Gilbert approach to large-scale tomographic applications. In this study, we introduce and adapt to seismic tomography the Subtractive Optimally Localized Averages (SOLA) method, an alternative Backus–Gilbert formulation which retains all its advantages, but is more computationally efficient and versatile in the explicit construction of averaging kernels. As a leitmotiv, we focus on global-scale S -wave tomography and show that the SOLA method can successfully be applied to large-scale, linear and discrete tomographic problems.

Description: The envelope broadening and the peak delay of the S -wavelet of a small earthquake with increasing travel distance are results of scattering by random velocity inhomogeneities in the earth medium. As a simple mathematical model, Sato proposed a new stochastic synthesis of the scalar wavelet envelope in 3-D von Kármán type random media when the centre wavenumber of the wavelet is in the power-law spectral range of the random velocity fluctuation. The essential idea is to split the random medium spectrum into two components using the centre wavenumber as a reference: the long-scale (low-wavenumber spectral) component produces the peak delay and the envelope broadening by multiple scattering around the forward direction; the short-scale (high-wavenumber spectral) component attenuates wave amplitude by wide angle scattering. The former is calculated by the Markov approximation based on the parabolic approximation and the latter is calculated by the Born approximation. Here, we extend the theory for the envelope synthesis of a wavelet in 2-D random media, which makes it easy to compare with finite difference (FD) simulation results. The synthetic wavelet envelope is analytically written by using the random medium parameters in the angular frequency domain. For the case that the power spectral density function of the random velocity fluctuation has a steep roll-off at large wavenumbers, the envelope broadening is small and frequency independent, and scattering attenuation is weak. For the case of a small roll-off, however, the envelope broadening is large and increases with frequency, and the scattering attenuation is strong and increases with frequency. As a preliminary study, we compare synthetic wavelet envelopes with the average of FD simulation wavelet envelopes in 50 synthesized random media, which are characterized by the RMS fractional velocity fluctuation = 0.05, correlation scale a  = 5 km and the background wave velocity V 0 = 4 km s –1 . We use the radiation of a 2 Hz Ricker wavelet from a point source. For all the cases of von Kármán order  = 0.1, 0.5 and 1, we find the synthetic wavelet envelopes are a good match to the characteristics of FD simulation wavelet envelopes in a time window starting from the onset through the maximum peak to the time when the amplitude decreases to half the peak amplitude.

Description: We present a thermodynamically based formulation for modelling dynamic rupture processes in the brittle crust using a continuum damage-breakage rheology. The model combines aspects of a continuum viscoelastic damage framework for brittle solids with a continuum breakage mechanics for granular flow within dynamically generated slip zones. The formulation accounts for the density of distributed cracking and other internal flaws in damaged rocks with a scalar damage parameter, and addresses the grain size distribution of a granular phase in the slip zone with a breakage parameter. A dynamic brittle instability is associated with a critical level of damage in the solid, leading to loss of convexity of the solid strain energy, localization and transition to a granular phase associated with lower energy level. The continuum damage-breakage rheology model treats the localization to a slip zone at the onset of dynamic rupture and post-failure recovery process as phase transitions between solid and granular states. The model generates sub- and supershear rupture velocities and pulse-type ruptures seen also in frictional models, and additional important features such as strong dynamic changes of volumetric strain near the rupture front and diversity of nucleation mechanisms. The propagation of rupture front and slip accumulation at a point are correlated with sharp dynamic dilation followed by a gradual decay to a level associated with the final volumetric change associated with the granular phase transition in the slipping zone. The local brittle failure process associated with the solid–granular transition is expected to produce isotropic radiation in addition to the deviatoric terms. The framework significantly extends the ability to model brittle processes in complex geometrical structures and allows analysing the roles of gouge thickness and other parameters on nucleation, rupture and radiation characteristics.

Description: Numerical modeling of seismic waves in heterogeneous porous reservoir rocks is an important tool for the interpretation of seismic surveys in reservoir engineering. We apply globally optimal implicit staggered-grid finite differences (FD) to model 2-D wave propagation in heterogeneous poroelastic media at a low-frequency range (〈10 kHz). We validate the numerical solution by comparing it to an analytical-transient solution obtaining clear seismic wavefields including fast P and slow P and S waves (for a porous media saturated with fluid). The numerical dispersion and stability conditions are derived using von Neumann analysis, showing that over a wide range of porous materials the Courant condition governs the stability and this optimal implicit scheme improves the stability of explicit schemes. High-order explicit FD can be replaced by some lower order optimal implicit FD so computational cost will not be as expensive while maintaining the accuracy. Here, we compute weights for the optimal implicit FD scheme to attain an accuracy of = 10 –8 . The implicit spatial differentiation involves solving tridiagonal linear systems of equations through Thomas’ algorithm.

Description: Tectonic earthquake swarms challenge our understanding of earthquake processes since it is difficult to link observations to the underlying physical mechanisms and to assess the hazard they pose. Transient forcing is thought to initiate and drive the spatio-temporal release of energy during swarms. The nature of the transient forcing may vary across sequences and range from aseismic creeping or transient slip to diffusion of pore pressure pulses to fluid redistribution and migration within the seismogenic crust. Distinguishing between such forcing mechanisms may be critical to reduce epistemic uncertainties in the assessment of hazard due to seismic swarms, because it can provide information on the frequency–magnitude distribution of the earthquakes (often deviating from the assumed Gutenberg–Richter relation) and on the expected source parameters influencing the ground motion (for example the stress drop). Here we study the ongoing Pollino range (Southern Italy) seismic swarm, a long-lasting seismic sequence with more than five thousand events recorded and located since October 2010. The two largest shocks (magnitude M w = 4.2 and M w = 5.1) are among the largest earthquakes ever recorded in an area which represents a seismic gap in the Italian historical earthquake catalogue. We investigate the geometrical, mechanical and statistical characteristics of the largest earthquakes and of the entire swarm. We calculate the focal mechanisms of the M l 〉 3 events in the sequence and the transfer of Coulomb stress on nearby known faults and analyse the statistics of the earthquake catalogue. We find that only 25 per cent of the earthquakes in the sequence can be explained as aftershocks, and the remaining 75 per cent may be attributed to a transient forcing. The b -values change in time throughout the sequence, with low b -values correlated with the period of highest rate of activity and with the occurrence of the largest shock. In the light of recent studies on the palaeoseismic and historical activity in the Pollino area, we identify two scenarios consistent with the observations and our analysis: This and past seismic swarms may have been ‘passive’ features, with small fault patches failing on largely locked faults, or may have been accompanied by an ‘active’, largely aseismic, release of a large portion of the accumulated tectonic strain. Those scenarios have very different implications for the seismic hazard of the area.

Description: The electrical current density generated by the propagation of a seismic wave at the interface characterized by a drop in electrical, hydraulic or mechanical properties produces an electrical field of electrokinetic nature. This field can be measured remotely with a signal-to-noise ratio depending on the background noise and signal attenuation. The seismoelectric beamforming approach is an emerging imaging technique based on scanning a porous material using appropriately delayed seismic sources. The idea is to focus the hydromechanical energy on a regular spatial grid and measure the converted electric field remotely at each focus time. This method can be used to image heterogeneities with a high definition and to provide structural information to classical geophysical methods. A numerical experiment is performed to investigate the resolution of the seismoelectric beamforming approach with respect to the main wavelength of the seismic waves. The 2-D model consists of a fictitious water-filled bucket in which a cylindrical sandstone core sample is set up vertically. The hydrophones/seismic sources are located on a 50-cm diameter circle in the bucket and the seismic energy is focused on the grid points in order to scan the medium and determine the geometry of the porous plug using the output electric potential image. We observe that the resolution of the method is given by a density of eight scanning points per wavelength. Additional numerical tests were also performed to see the impact of a wrong velocity model upon the seismoelectric map displaying the heterogeneities of the material.

Description: The structure of seismic discontinuities in the mantle transition zone at depths of about 400 and 670 km provides important constraints on mantle convection as the associated mineral phase transformations are sensitive to thermal perturbations. Teleseismic P -to- S receiver functions have been widely used to map the depths of the two discontinuities. In this study, we investigate the resolution of receiver functions in imaging topographic variations of the 400-km and 670-km discontinuities based on wave propagation simulations using a Spectral Element Method (SEM). We investigate wave diffraction effects on direct P waves as well as P -to- S converted waves by varying the length scale of topography of the two discontinuities. We observe strong wave diffractional effects in both P waves and teleseismic receiver functions at periods of ~10 to 20 s. Ray theory overpredicts traveltime anomalies by a factor of 2–5 when the topography length scale is about 400 km. In addition, ray-theoretical predictions are out of phase with measurements which indicates that locations of small-scale topographic variations can not be resolved using ray theory. The observed traveltime anomalies further reduce to 10–20 per cent of ray-theoretical predictions when the topography length scale reduces to about 200 km. We calculate 2-D boundary sensitivity kernels for direct P waves as well as receiver functions. In general, calculations based Born sensitivity kernels fit the ‘ground-truth’ SEM measurements very well. They account for wave diffraction effects as well as phase interactions such as P and pP waves arriving in P -wave coda. 3-D wavespeed structure in the upper mantle beneath seismic stations may introduce significant traveltime anomalies on P waves and transition zone receiver functions. We show that traveltime corrections at periods of about 10 to 20 s are frequency dependent when the size of the anomalies becomes less than 500 km.

Description: The aim of this study is to improve the temporal resolution of seismic wave velocity variations measured using ambient noise correlations. We first reproduce the result obtained by Chen et al. using a network of 21 broad-band stations ideally located around the fault system activated during the Wenchuan earthquake.We measure a velocity drop of 0.07 per cent that was associated with the main shock, with a temporal resolution of 30 days. To determine whether this velocity drop is co-seismic or post-seismic, we attempt to increase the temporal resolution of our observations. By taking advantage of the properties of the curvelet transform, we increase the signal-to-noise ratio of the daily correlations computed between each station pair. It is then possible to measure the velocity drop associated with the Wenchuan earthquake with a temporal resolution of 1 day. This shows that the velocity drop started on 2008 May 12, which was the day of the earthquake, and the velocity reached its lowest value 2 days after the main shock. Moreover, there was a second velocity drop on 2008 May 27, which might relate to strong aftershocks.

Description: Several destructive earthquakes have occurred in Tien-Shan region at the beginning of 20th century. However, the detailed seismological characteristics, especially source parameters of those earthquakes are still poorly investigated. The Chon-Kemin earthquake is the strongest instrumentally recorded earthquake in the Tien-Shan region. This earthquake has produced an approximately 200 km long system of surface ruptures along Kemin-Chilik fault zone and killed about ~400 people. Several studies presented the different information on the earthquake epicentre location and magnitude, and two different focal mechanisms were also published. The reason for the limited knowledge of the source parameters for the Chon-Kemin earthquake is the complexity of old analogue records processing, digitization and analysis. In this study the data from 23 seismic stations worldwide were collected and digitized. The earthquake epicentre was relocated to 42.996N° and 77.367E°, the hypocentre depth is estimated between 10 and 20 km. The magnitude was recalculated to m B 8.05, M s 7.94 and M w 8.02. The focal mechanism, determined from amplitude ratios comparison of the observed and synthetic seismograms, was: str = 264°, dip = 52°, rake = 98°. The apparent source time duration was between ~45 and ~70 s, the maximum slip occurred 25 s after the beginning of the rupture. Two subevents were clearly detected from the waveforms with the scalar moment ratio between them of about $\frac{1}{3}$ , the third subevent was also detected with less certainty. Taking into account surface rupture information, the fault geometry model with three patches was proposed. Based on scaling relations we conclude that the total rupture length was between ~260 and 300 km and a maximum rupture width could reach ~70 km.

Description: New estimates of crustal thickness, Poisson's ratio and crustal shear wave velocity have been obtained for 39 stations in Angola, Botswana, the Democratic Republic of Congo, Malawi, Mozambique, Namibia, Rwanda, Tanzania and Zambia by modelling P -wave receiver functions using the H – stacking method and jointly inverting the receiver functions with Rayleigh-wave phase and group velocities. These estimates, combined with similar results from previous studies, have been examined for secular trends in Precambrian crustal structure within the southern African subcontinent. In both Archean and Proterozoic terranes we find similar Moho depths [38–39 ± 3 km SD (standard deviation)], crustal Poisson's ratio (0.26 ± 0.01 SD ), mean crustal shear wave velocity (3.7 ± 0.1 km s –1 SD ), and amounts of heterogeneity in the thickness of the mafic lower crust, as defined by shear wave velocities ≥4.0 km s –1 . In addition, the amount of variability in these crustal parameters is similar within each individual age grouping as between age groupings. Thus, the results provide little evidence for secular variation in Precambrian crustal structure, including between Meso- and Neoarchean crust. This finding suggests that (1) continental crustal has been generated by similar processes since the Mesoarchean or (2) plate tectonic processes have reworked and modified the crust through time, erasing variations in structure resulting from crustal genesis.

Description: The relative seismic velocity variations possibly associated to large earthquakes can be readily monitored via cross-correlation of seismic noise. In a recently published study, more than 2 yr of continuous seismic records have been analysed from three stations surrounding the epicentre of the 2009 April 6, M w 6.1 L'Aquila earthquake, observing a clear decrease of seismic velocities likely corresponding to the co-seismic shaking. Here, we extend the analysis in space, including seismic stations within a radius of 60 km from the main shock epicentre, and in time, collecting 5 yr of data for the six stations within 40 km of it. Our aim is to investigate how far the crustal damage is visible through this technique, and to detect a potential post-seismic recovery of velocity variations. We find that the co-seismic drop in velocity variations extends up to 40 km from the epicentre, with spatial distribution (maximum around the fault and in the north–east direction from it) in agreement with the horizontal co-seismic displacement detected by global positioning system (GPS). In the first few months after L'Aquila earthquake, the crust's perturbation in terms of velocity variations displays a very unstable behaviour, followed by a slow linear recovery towards pre-earthquake conditions; by almost 4 yr after the event, the co-seismic drop of seismic velocity is not yet fully recovered. The strong oscillations of the velocity changes in the first months after the earthquake prevent to detect the fast exponential recovery seen by GPS data. A test of differently parametrized fitting curves demonstrate that the post-seismic recovery is best explained by a sum of a logarithmic and a linear term, suggesting that processes like viscoelastic relaxation, frictional afterlip and poroelastic rebound may be acting concurrently.

Description: To improve the performance of Global Positioning System (GPS) in the earthquake/tsunami early warning and rapid response applications, minimizing the blind zone and increasing the stability and accuracy of both the rapid source and rupture inversion, the density of existing GPS networks must be increased in the areas at risk. For economic reasons, low-cost single-frequency receivers would be preferable to make the sparse dual-frequency GPS networks denser. When using single-frequency GPS receivers, the main problem that must be solved is the ionospheric delay, which is a critical factor when determining accurate coseismic displacements. In this study, we introduce a modified Satellite-specific Epoch-differenced Ionospheric Delay (MSEID) model to compensate for the effect of ionospheric error on single-frequency GPS receivers. In the MSEID model, the time-differenced ionospheric delays observed from a regional dual-frequency GPS network to a common satellite are fitted to a plane rather than part of a sphere, and the parameters of this plane are determined by using the coordinates of the stations. When the parameters are known, time-differenced ionospheric delays for a single-frequency GPS receiver could be derived from the observations of those dual-frequency receivers. Using these ionospheric delay corrections, coseismic displacements of a single-frequency GPS receiver can be accurately calculated based on time-differenced carrier-phase measurements in real time. The performance of the proposed approach is validated using 5 Hz GPS data collected during the 2012 Nicoya Peninsula Earthquake ( M w 7.6, 2012 September 5) in Costa Rica. This shows that the proposed approach improves the accuracy of the displacement of a single-frequency GPS station, and coseismic displacements with an accuracy of a few centimetres are achieved over a 10-min interval.

Description: We imaged the azimuthal anisotropy of Rayleigh wave phase velocity (10–60 s) in northeast North China Craton using the teleseismic data recorded by a dense temporary array, and then inverted for the 3-D azimuthal anisotropy of the crust and uppermost mantle (20–110 km). The results reveal that the azimuthal anisotropy varies both horizontally and vertically. Obvious stratified azimuthal anisotropy is shown in the Central Orogenic Belt, where the fast direction is NE–SW to NNE–SSW in the depth range of 20–40 km and changes to NW–SE to NWW–SEE in the depth range of 60–110 km. In the depth range of 30–40 km, a prominent low velocity belt is shown on the southwest of Zhangjiakou-Penglai fault zone (ZPFZ) and the fast direction is subparallel to the strike of the low velocity belt. Distinct lateral variations of azimuthal anisotropy are clearly shown at 110 km. Our results provide new evidence for the existence of upwelling asthenosphere beneath the Datong volcano and support the assumption that ZPFZ may act as the channel of upwelling asthenosphere. Historical strong earthquakes ( M ≥ 6.0) mainly occurred in the transition zone between low and high velocity anomalies in the upper and middle crust. The upwelling asthenosphere may prompt the generation of large earthquake.

Description: We explore with numerical simulations basic physical conditions leading to key observed features of non-volcanic tremor (NVT) in relation to slow slip events (SSEs) and earthquakes along the Guerrero segment of the Mexican subduction zone. To study the interactions between different modes of slip, and examine possible variations over timescales larger than the 15 year observational interval, we use a model with a planar interface governed by space-varying static/kinetic friction and dislocation creep in a 3-D elastic solid. A fault section with zero weakening during frictional slip fails in a mode corresponding to a ‘critical depinning transition’ that produces generically many observed features of NVT. A patch with elevated creep coefficients represents a section with SSE. Simulations with small added stress oscillations are used to examine triggering of NVT by large remote earthquakes. The results reproduce well the basic observed properties of NVT and SSE in the Guerrero area, while pointing to complex interactions between large earthquake cycles, quasi-period SSE and scale-invariant NVT behaviour. The model simulations provide additional information on expected frequency-magnitude statistics, slip distributions and space–time properties of the different event types that may be tested with accumulation of future data. Some earthquake and NVT events near the opposite sides of the SSE patch have significant separation between their hypocentres and centroids. The rates of these events are correlated with the creep evolution in the SSE section. The results also suggest that aseismic deformation in the area may have transients on timescales larger than the observational period.

Description: High-melt areas of glaciers and ice sheets foster a rich spectrum of ambient seismicity. These signals not only shed light on source mechanisms (e.g. englacial fracturing, water flow, iceberg detachment, basal motion) but also carry information about seismic wave propagation within glacier ice. Here, we present two approaches to measure and potentially monitor phase velocities of high-frequency seismic waves (≥1 Hz) using naturally occurring glacier seismicity. These two approaches were developed for data recorded by on-ice seasonal seismic networks on the Greenland Ice Sheet and a Swiss Alpine glacier. The Greenland data set consists of continuous seismograms, dominated by long-term tremor-like signals of englacial water flow, whereas the Alpine data were collected in triggered mode producing 1–2 s long records that include fracture events within the ice (‘icequakes’). We use a matched-field processing technique to retrieve frequency-dependent phase velocity measurements for the Greenland data. In principle, this phase dispersion relationship can be inverted for ice sheet thickness and bed properties. For these Greenland data, inversion of the dispersion curve yields a bedrock depth of 541 m, which may be too small by as much as 35 per cent. We suggest that the discrepancy is due to lateral changes in ice sheet depth and bed properties beneath the network, which may cause unaccounted mixing of surface wave modes in the dispersion curve. The Swiss Alpine icequake records, on the other hand, allow for reconstruction of the impulse response between two seismometers. The direct and scattered wave fields from the vast numbers of icequake records (tens of thousands per month) can be used to measure small changes in englacial velocities and thus monitor structural changes within the ice.

Description: Automatic classification of local seismic events which are only recorded at single stations poses great challenges because of weak hypocentre constraints. This study investigates how single-station event clusters relate to geographic hypocentre regions and common source processes. Typical applications arise in local seismic networks where reliable ground truth by a dense temporal network precedes or follows a sparse (permanent) installation. The seismic signals for this study comprise a 3-month subset from a field campaign to map subduction below northern Chile (PISCO ’94). Due to favourable ground noise conditions in the Atacama desert, the data set contains an abundance of shallow and deeper earthquakes, and many quarry explosions. Often event signatures overlap, posing a challenge to any signal processing scheme. Pattern recognition must work on reduced seismograms to restrict parameter dimensionality. Continuous parameter extraction based on noise-adapted spectrograms was chosen instead of discrete representation by, for example, amplitudes, onset times or spectral ratios to ensure consideration of potentially hidden features. Visualization of the derived feature vectors for human inspection and template matching algorithms was hereby possible. Because event classes shall comprise earthquake regions regardless of magnitude, clustering based on amplitudes is prevented by proper normalization of feature vectors. Principal component analysis is applied to further reduce the number of features used to train a self-organizing map (SOM). The SOM will topologically arrange prototypes of each event class in a 2-D map. Overcoming the restrictions of this black-box approach, the arranged prototypes could be transformed back to spectrograms to allow for visualization and interpretation of event classes. The final step relates prototypes to ground-truth information, confirming the potential of automated, coarse-grain hypocentre clustering based on single-station seismograms. The approach was tested by a twofold cross-validation whereby multiple sets of feature vectors from half the events are compared by a one-nearest neighbour classifier in combination with an Euclidean distance measure resulting in an overall correct geographic separation rate of 80.5 per cent.

Description: Finite-fault earthquake source inversion is an ill-posed inverse problem leading to non-unique solutions. In addition, various fault parametrizations and input data may have been used by different researchers for the same earthquake. Such variability leads to large intra-event variability in the inferred rupture models. One way to understand this problem is to develop robust metrics to quantify model variability. We propose a Multi Dimensional Scaling (MDS) approach to compare rupture models quantitatively. We consider normalized squared and grey-scale metrics that reflect the variability in the location, intensity and geometry of the source parameters. We test the approach on two-dimensional random fields generated using a von Kármán autocorrelation function and varying its spectral parameters. The spread of points in the MDS solution indicates different levels of model variability. We observe that the normalized squared metric is insensitive to variability of spectral parameters, whereas the grey-scale metric is sensitive to small-scale changes in geometry. From this benchmark, we formulate a similarity scale to rank the rupture models. As case studies, we examine inverted models from the Source Inversion Validation (SIV) exercise and published models of the 2011 Mw 9.0 Tohoku earthquake, allowing us to test our approach for a case with a known reference model and one with an unknown true solution. The normalized squared and grey-scale metrics are respectively sensitive to the overall intensity and the extension of the three classes of slip (very large, large, and low). Additionally, we observe that a three-dimensional MDS configuration is preferable for models with large variability. We also find that the models for the Tohoku earthquake derived from tsunami data and their corresponding predictions cluster with a systematic deviation from other models. We demonstrate the stability of the MDS point-cloud using a number of realizations and jackknife tests, for both the random field and the case studies.

Description: Recently, various methods have been proposed and applied for earthquake source imaging, and theoretical relationships among the methods have been studied. In this study, we make a follow-up theoretical study to better understand the meanings of earthquake source imaging. For imaging problems, the point spread function (PSF) is used to describe the degree of blurring and degradation in an obtained image of a target object as a response of an imaging system. In this study, we formulate PSFs for earthquake source imaging. By calculating the PSFs, we find that waveform source inversion methods remove the effect of the PSF and are free from artefacts. However, the other source imaging methods are affected by the PSF and suffer from the effect of blurring and degradation due to the restricted distribution of receivers. Consequently, careful treatment of the effect is necessary when using the source imaging methods other than waveform inversions. Moreover, the PSF for source imaging is found to have a link with seismic interferometry with the help of the source–receiver reciprocity of Green's functions. In particular, the PSF can be related to Green's function for cases in which receivers are distributed so as to completely surround the sources. Furthermore, the PSF acts as a low-pass filter. Given these considerations, the PSF is quite useful for understanding the physical meaning of earthquake source imaging.

Description: We introduce the Complex-Step-Finite-Difference method (CSFDM) as a generalization of the well-known Finite-Difference method (FDM) for solving the acoustic and elastic wave equations. We have found a direct relationship between modelling the second-order wave equation by the FDM and the first-order wave equation by the CSFDM in 1-D, 2-D and 3-D acoustic media. We present the numerical methodology in order to apply the introduced CSFDM and show an example for wave propagation in simple homogeneous and heterogeneous models. The CSFDM may be implemented as an extension into pre-existing numerical techniques in order to obtain fourth- or sixth-order accurate results with compact three time-level stencils. We compare advantages of imposing various types of initial motion conditions of the CSFDM and demonstrate its higher-order accuracy under the same computational cost and dispersion–dissipation properties. The introduced method can be naturally extended to solve different partial differential equations arising in other fields of science and engineering.

Description: A discontinuous grid finite-difference (FD) method with non-uniform time step Runge–Kutta scheme on curvilinear collocated-grid is developed for seismic wave simulation. We introduce two transition zones: a spatial transition zone and a temporal transition zone, to exchange wavefield across the spatial and temporal discontinuous interfaces. A Gaussian filter is applied to suppress artificial numerical noise caused by down-sampling the wavefield from the finer grid to the coarser grid. We adapt the non-uniform time step Runge–Kutta scheme to a discontinuous grid FD method for further increasing the computational efficiency without losing the accuracy of time marching through the whole simulation region. When the topography is included in the modelling, we carry out the discontinuous grid method on a curvilinear collocated-grid to obtain a sufficiently accurate free-surface boundary condition implementation. Numerical tests show that the proposed method can sufficiently accurately simulate the seismic wave propagation on such grids and significantly reduce the computational resources consumption with respect to regular grids.

Description: SS precursor observations are a powerful tool to study the topography and character of transition zone discontinuities, especially in regions such as ocean basins where few seismic stations exist, precluding other high resolution approaches. Still, the available coverage is limited by the distribution of sources and stations, but also by the level of noise and by the fact that, in some distance ranges, interfering seismic phases mask the weak signal from the SS precursors. We introduce an array data processing tool, the local slant-stack filter, to address these challenges and clean up the otherwise noisy SS precursor record sections. We show that these filters are a powerful tool for extracting the weak yet coherent SS precursor signals while removing interfering seismic phases as well as random noise, yielding robust precursor traveltime measurements with spatial resolution higher than what can be achieved by the conventional common midpoint stacking method. The effectiveness of the filters are demonstrated by application to synthetic and real data. We systematically apply this filtering method to an SS precursor data set recorded by the U.S. Transportable Array that samples a vast region of the Pacific Ocean and its northwest margin, and present maps of 410 and 660 discontinuity topography. We discuss correlations observed between our discontinuity images and several fine-scale heterogeneities revealed by mantle shear wave tomography in the vicinity of Hawaii and the Pacific Superswell.