ALBERT

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: S -wave receiver functions (SRFs) are used to investigate crustal and upper-mantle structure beneath several ice-covered areas of Antarctica. Moho S-to-P (Sp) arrivals are observed at ~6–8 s in SRF stacks for stations in the Gamburtsev Mountains (GAM) and Vostok Highlands (VHIG), ~5–6 s for stations in the Transantarctic Mountains (TAM) and the Wilkes Basin (WILK), and ~3–4 s for stations in the West Antarctic Rift System (WARS) and the Marie Byrd Land Dome (MBLD). A grid search is used to model the Moho Sp conversion time with Rayleigh wave phase velocities from 18 to 30 s period to estimate crustal thickness and mean crustal shear wave velocity. The Moho depths obtained are between 43 and 58 km for GAM, 36 and 47 km for VHIG, 39 and 46 km for WILK, 39 and 45 km for TAM, 19 and 29 km for WARS and 20 and 35 km for MBLD. SRF stacks for GAM, VHIG, WILK and TAM show little evidence of Sp arrivals coming from upper-mantle depths. SRF stacks for WARS and MBLD show Sp energy arriving from upper-mantle depths but arrival amplitudes do not rise above bootstrapped uncertainty bounds. The age and thickness of the crust is used as a heat flow proxy through comparison with other similar terrains where heat flow has been measured. Crustal structure in GAM, VHIG and WILK is similar to Precambrian terrains in other continents where heat flow ranges from ~41 to 58 mW m –2 , suggesting that heat flow across those areas of East Antarctica is not elevated. For the WARS, we use the Cretaceous Newfoundland–Iberia rifted margins and the Mesozoic-Tertiary North Sea rift as tectonic analogues. The low-to-moderate heat flow reported for the Newfoundland–Iberia margins (40–65 mW m –2 ) and North Sea rift (60–85 mW m –2 ) suggest that heat flow across the WARS also may not be elevated. However, the possibility of high heat flow associated with localized Cenozoic extension or Cenozoic-recent magmatic activity in some parts of the WARS cannot be ruled out.

Description: The 2014 April 1, M w 8.1 Iquique (Chile) earthquake struck in the Northern Chile seismic gap. With a rupture length of less than 200 km, it left unbroken large segments of the former gap. Early studies were able to model the main rupture features but results are ambiguous with respect to the role of aseismic slip and left open questions on the remaining hazard at the Northern Chile gap. A striking observation of the 2014 earthquake has been its extensive preparation phase, with more than 1300 events with magnitude above M L 3, occurring during the 15 months preceding the main shock. Increasing seismicity rates and observed peak magnitudes accompanied the last three weeks before the main shock. Thanks to the large data sets of regional recordings, we assess the precursor activity, compare foreshocks and aftershocks and model rupture preparation and rupture effects. To tackle inversion challenges for moderate events with an asymmetric network geometry, we use full waveforms techniques to locate events, map the seismicity rate and derive source parameters, obtaining moment tensors for more than 300 events (magnitudes M w 4.0–8.1) in the period 2013 January 1–2014 April 30. This unique data set of fore- and aftershocks is investigated to distinguish rupture process models and models of strain and stress rotation during an earthquake. Results indicate that the spatial distributions of foreshocks delineated the shallower part of the rupture areas of the main shock and its largest aftershock, well matching the spatial extension of the aftershocks cloud. Most moment tensors correspond to almost pure double couple thrust mechanisms, consistent with the slab orientation. Whereas no significant differences are observed among thrust mechanisms in different areas, nor among thrust foreshocks and aftershocks, the early aftershock sequence is characterized by the presence of normal fault mechanisms, striking parallel to the trench but dipping westward. These events likely occurred in the shallow wedge structure close to the slab interface and are consequence of the increased extensional stress in this region after the largest events. The overall stress inversion result suggests a minor stress rotation after the main shock, but a significant release of the deviatoric stress. The temporal change in the distribution of focal mechanisms can also be explained in terms of the spatial heterogeneity of the stress field: under such interpretation, the potential of a large megathrust earthquake breaking a larger segment offshore Northern Chile remains high.

Description: We analyse the performance of a higher order accurate staggered viscoelastic time-domain finite-difference method, in which the staggered Adams–Bashforth (ABS) third-order and fourth-order accurate time integrators are used for temporal discretization. ABS is a multistep method that uses previously calculated wavefields to increase the order of accuracy in time. The analysis shows that the numerical dispersion is much lower than that of the widely used second-order leapfrog method. Numerical dissipation is introduced by the ABS method which is significantly smaller for fourth-order than third-order accuracy. In 1-D and 3-D simulation experiments, we verify the convincing improvements of simulation accuracy of the fourth-order ABS method. In a realistic elastic 3-D scenario, the computing time reduces by a factor of approximately 2.4, whereas the memory requirements increase by approximately a factor of 2.2. The ABS method thus provides an alternative strategy to increase the simulation accuracy in time by investing computer memory instead of computing time.

Description: Synthetic dispersion curves are here computed in the frame of an ambient-vibration full-wavefield model, which relies on the description of both ambient-vibration ground displacement and its sources as stochastic fields defined on the Earth's surface, stationary in time and homogeneous in space. In this model, previously developed for computing synthetic Horizontal-to-Vertical Spectral Ratio curves, the power spectral density function and the spatial autocorrelation of the displacement are naturally described as functions of the power spectral density function of the generating forces and of the subsoil properties ( via the relevant Green's function), by also accounting for spatial correlation of these forces. Dispersion curves are computed from the displacement power spectral density function and from the spatial autocorrelation according with the well-known f-k and SPAC techniques, respectively. Two examples illustrate the way this new ambient-vibration model works, showing its possible use in better understanding the role of the surface waves in forming the dispersion curves, as well as its capability to capture some remarkable experimental findings.

Description: The Abruzzi region is located in the Central Apennines Neogene fold-and-thrust belt and has one of the highest seismogenic potential in Italy, with high and diffuse crustal seismicity related to NE–SW oriented extension. In this study, we investigate the detailed spatial variation in shear wave splitting providing high-resolution anisotropic structure beneath the L'Aquila region. To accomplish this, we performed a systematic analysis of crustal anisotropic parameters: fast polarization direction ( ) and delay time ( t ). We benefit from the dense coverage of seismic stations operating in the area and from a catalogue of several accurate earthquake locations of the 2009 L'Aquila seismic sequence, related to the M w 6.1 2009 L'Aquila main shock, to describe in detail the geometry of the anisotropic volume around the active faults that ruptured. The spatial variations both in and t suggest a complex anisotropic structure beneath the region caused by a combination of both structural- and stress-induced mechanisms. The average is NNW–SSE oriented (N141°), showing clear similarity both with the local fault strike and the SH max . In the central part of the study area fast axes are oriented NW–SE, while moving towards the northeastern and northwestern sectors the fast directions clearly diverge from the general trend of NW–SE and rotate accordingly to the local fault strikes. The above-mentioned fault-parallel distribution suggests that the observed anisotropy is mostly controlled by the local fault-related structure. Toward the southeast fast directions become orthogonal both to strike of the local mapped faults and to the SH max . Here, are predominantly oriented NE–SW; we interpret this orientation as due to the presence of a highly fractured and overpressurized rock volume which should be responsible of the 90° flips in and the increase in t . Another possible mechanism for NE–SW orientation of in the southeastern sector could be ascribed to the presence of a buried, deep NE–SW oriented fault system. t , both unnormalized and normalized, does not show any clear evidence of increasing with increasing depth indicating that the anisotropy is confined primarily to the shallower crustal layers (~10 km depth). Interpolating t show that higher values are found at the edges of the main patches of the rupture related to the 2009 main shock, while lower values are limited in the central part of the fault plane, where the coseismic slip was higher. We infer that in the areas surrounding the ruptured region, lateral variations in material properties caused overpressurized fluid conditions, while within the main shock ruptured area, high energy released produced an open crack system such that overpressurization was not possible.

Description: This paper presents a set of deterministic 3-D ground motion simulations for the greater metropolitan area of Medellín in the Aburrá Valley, an earthquake-prone region of the Colombian Andes that exhibits moderate-to-strong topographic irregularities. We created the velocity model of the Aburrá Valley region (version 1) using the geological structures as a basis for determining the shear wave velocity. The irregular surficial topography is considered by means of a fictitious domain strategy. The simulations cover a 50 x 50 x 25 km 3 volume, and four M w = 5 rupture scenarios along a segment of the Romeral fault, a significant source of seismic activity in Colombia. In order to examine the sensitivity of ground motion to the irregular topography and the 3-D effects of the valley, each earthquake scenario was simulated with three different models: (i) realistic 3-D velocity structure plus realistic topography, (ii) realistic 3-D velocity structure without topography, and (iii) homogeneous half-space with realistic topography. Our results show how surface topography affects the ground response. In particular, our findings highlight the importance of the combined interaction between source-effects, source-directivity, focusing, soft-soil conditions, and 3-D topography. We provide quantitative evidence of this interaction and show that topographic amplification factors can be as high as 500 per cent at some locations. In other areas within the valley, the topographic effects result in relative reductions, but these lie in the 0–150 per cent range.

Description: Wave-induced local fluid flow is known as a key mechanism to explain the intrinsic wave dissipation in fluid-saturated rocks. Understanding the relationship between the acoustic properties of rocks and fluid patch distributions is important to interpret the observed seismic wave phenomena. A triple-layer patchy (TLP) model is proposed to describe the P -wave dissipation process in a double-porosity media saturated with two immiscible fluids. The double-porosity rock consists of a solid matrix with unique host porosity and inclusions which contain the second type of pores. Two immiscible fluids are considered in concentric spherical patches, where the inner pocket and the outer sphere are saturated with different fluids. The kinetic and dissipation energy functions of local fluid flow (LFF) in the inner pocket are formulated through oscillations in spherical coordinates. The wave propagation equations of the TLP model are based on Biot's theory and the corresponding Lagrangian equations. The P -wave dispersion and attenuation caused by the Biot friction mechanism and the local fluid flow (related to the pore structure and the fluid distribution) are obtained by a plane-wave analysis from the Christoffel equations. Numerical examples and laboratory measurements indicate that P -wave dispersion and attenuation are significantly influenced by the spatial distributions of both, the solid heterogeneity and the fluid saturation distribution. The TLP model is in reasonably good agreement with White's and Johnson's models. However, differences in phase velocity suggest that the heterogeneities associated with double-porosity and dual-fluid distribution should be taken into account when describing the P -wave dispersion and attenuation in partially saturated rocks.

Description: Seismic ambient noise tomography is applied to central and southern Mozambique, located in the tip of the East African Rift (EAR). The deployment of MOZART seismic network, with a total of 30 broad-band stations continuously recording for 26 months, allowed us to carry out the first tomographic study of the crust under this region, which until now remained largely unexplored at this scale. From cross-correlations extracted from coherent noise we obtained Rayleigh wave group velocity dispersion curves for the period range 5–40 s. These dispersion relations were inverted to produce group velocity maps, and 1-D shear wave velocity profiles at selected points. High group velocities are observed at all periods on the eastern edge of the Kaapvaal and Zimbabwe cratons, in agreement with the findings of previous studies. Further east, a pronounced slow anomaly is observed in central and southern Mozambique, where the rifting between southern Africa and Antarctica created a passive margin in the Mesozoic, and further rifting is currently happening as a result of the southward propagation of the EAR. In this study, we also addressed the question concerning the nature of the crust (continental versus oceanic) in the Mozambique Coastal Plains (MCP), still in debate. Our data do not support previous suggestions that the MCP are floored by oceanic crust since a shallow Moho could not be detected, and we discuss an alternative explanation for its ocean-like magnetic signature. Our velocity maps suggest that the crystalline basement of the Zimbabwe craton may extend further east well into Mozambique underneath the sediment cover, contrary to what is usually assumed, while further south the Kaapval craton passes into slow rifted crust at the Lebombo monocline as expected. The sharp passage from fast crust to slow crust on the northern part of the study area coincides with the seismically active NNE-SSW Urema rift, while further south the Mazenga graben adopts an N-S direction parallel to the eastern limit of the Kaapvaal craton. We conclude that these two extensional structures herald the southward continuation of the EAR, and infer a structural control of the transition between the two types of crust on the ongoing deformation.

Description: Seismic wavefields interact with low-velocity fault damage zones. Waveforms of ballistic fault zone head waves, trapped waves, reflected waves and signatures of trapped noise can provide important information on structural and mechanical fault zone properties. Here we extend the class of observable fault zone waves and reconstruct in-fault reverberations or multiples in a strike-slip faulting environment. Manifestations of the reverberations are significant, consistent wave fronts in the coda of cross-correlation functions that are obtained from scattered earthquake waveforms and seismic noise recorded by a linear fault zone array. The physical reconstruction of Green's functions is evident from the high similarity between the signals obtained from the two different scattered wavefields. Modal partitioning of the reverberation wavefield can be tuned using different data normalization techniques. The results imply that fault zones create their own ambiance, and that the here reconstructed reverberations are a key seismic signature of wear zones. Using synthetic waveform modelling we show that reverberations can be used for the imaging of structural units by estimating the location, extend and magnitude of lateral velocity contrasts. The robust reconstruction of the reverberations from noise records suggests the possibility to resolve the response of the damage zone material to various external and internal loading mechanisms.

Description: The Lower Rhine Embayment in Central Europe hosts a rift system that has very low deformation rates. The faults in this area have slip rates of less than 0.1 mm yr –1 , which does not allow to investigate ongoing tectonic deformation with geodetic techniques, unless they cover very long time spans. Instrumental seismicity does only cover a small fraction of the very long earthquake recurrence intervals of several thousands of years. Palaeoseismological studies are needed to constrain slip rates and the earthquake history of such faults. Destructive earthquakes are rare in the study area, but did occur in historic times. In 1755/1756, a series of strong earthquakes caused significant destruction in the city of Düren (Germany) and the surrounding areas. In this study we document palaeoseismological data from the nearby Rurrand Fault. In contrast to earlier studies on the same fault, we found evidence for a surface rupturing earthquake in the Holocene, and we identified at least one more surface rupturing event. Our study shows that the Rurrand Fault currently accommodates deformation in earthquakes rather than by creeping. The coseismic offsets were determined to be between less than 0.5 m per event. We assign maximum possible magnitudes of M w 5.9–6.8 for the Rurrand Fault and a slip rate of at least 0.02–0.03 mm yr –1 for the last ~130–50 kyr. The surface ruptures did not occur at the main fault trace that has a clear morphological expression due to older tectonic motions, but on a younger fault strand in the hanging wall of the main fault. Terrain analyses based on 1 m resolution airborne LiDAR data have been used to image the subtle morphological expression of this young fault zone. Georadar and electric resistivity tomography were applied to image the fault zone at depth and to test if these shallow geophysical methods can be used to identify and trace the fault zone. Georadar failed to produce reliable results, but geoelectrics were successfully applied and allowed us to retrieve slip rate estimates. Our results indicate that the Düren 1755/1756 earthquakes did not produce surface ruptures at the Rurrand Fault, either because they did not rupture the surface at all, or because they occurred at another, neighbouring fault.

Description: We have measured radial anisotropy in the crust beneath the northern part of North China by jointly inverting Rayleigh and Love wave phase velocities at periods less than 35 s from 14 months of ambient noise data recorded by 222 broad-band seismic stations. We also estimate the azimuthal anisotropy of phase velocity from Rayleigh wave data. The fast direction of azimuthal anisotropy varies with periods, NE–SW orientation at short and intermediate periods (10–16 s) and NW–SE orientation at periods larger than 20 s. The NE–SW oriented fast direction of azimuthal anisotropy may be related to the fossilized structural fabrics due to the compression during the Indosinian orogeny from late Palaeozoic to middle Mesozoic. The NW–SE trend of anisotropic fabric in the lower crust and uppermost mantle is probably associated with the later lithospheric extension. The observed radial anisotropy also shows a two-layer feature, negative radial anisotropy (Vsh 〈 Vsv) in the upper crust and positive (Vsh 〉 Vsv) in the middle-lower crust. The compressional tectonics from late Palaeozoic to middle Mesozoic may cause crustal materials align vertically throughout the crust. This vertical fabric could make Vsh slower than Vsv. However, the lithospheric extension in the late Mesozoic to Cenozoic time could overprint the older fabric in the middle and lower crust by magma intrusion and underplating. Horizontal alignment of the material or intruded melt sills due to the extension probably produce the observed strong positive radial anisotropy in the middle and lower crust.

Description: Using active and passive seismology data we derive a shear ( S ) wave velocity model and a Poisson's ratio () model across the Chilean convergent margin along a profile at 38°15'S, where the M w 9.5 Valdivia earthquake occurred in 1960. The derived S -wave velocity model was constructed using three independently obtained velocity models that were merged together. In the upper part of the profile (0–2 km depth), controlled source data from explosions were used to obtain an S -wave traveltime tomogram. For the middle part (2–20 km depth), data from a temporary seismology array were used to carry out a dispersion analysis. The resulting dispersion curves were used to obtain a 3-D S -wave velocity model. In the lower part (20–75 km depth, depending on the longitude), an already existent local earthquake tomographic image was merged with the other two sections. This final S -wave velocity model and already existent compressional ( P ) wave velocity models along the same transect allowed us to obtain a Poisson's ratio model. The results of this study show that the velocities and Poisson's ratios in the continental crust of this part of the Chilean convergent margin are in agreement with geological features inferred from other studies and can be explained in terms of normal rock types. There is no requirement to call on the existence of measurable amounts of present-day fluids, in terms of seismic velocities, above the plate interface in the continental crust of the Coastal Cordillera and the Central Valley in this part of the Chilean convergent margin. This is in agreement with a recent model of water being transported down and released from the subduction zone.

Description: In the estimation of spatiotemporal slip models, kinematic source inversions using Akaike's Bayesian Information Criterion (ABIC) and the multiple-time-window method have often been used. However, there are cases in which conventional ABIC-based source inversions do not work well in the determination of hyperparameters when a non-negative slip constraint is used. In order to overcome this problem, a new source inversion method was developed in this study. The new method introduces a fully Bayesian method into the kinematic multiple-time-window source inversion. The multiple-time-window method is one common way of parametrizing a source time function and is highly flexible in terms of the shape of the source time function. The probability distributions of model parameters and hyperparameters can be directly obtained by using the Markov chain Monte Carlo method. These probability distributions are useful for simply evaluating the uniqueness and reliability of the derived model, which is another advantage of a fully Bayesian method. This newly developed source inversion method was applied to the 2011 Ibaraki-oki, Japan, earthquake ( M w 7.9) to demonstrate its usefulness. It was demonstrated that the problem with using the conventional ABIC-based source inversion method for hyperparameter determination appeared in the spatiotemporal source inversion of this event and that the newly developed source inversion could overcome this problem.

Description: In a set of problems ranging from 4-D seismic to salt boundary estimation, updates to the velocity model often have a highly localized nature. Numerical techniques for these applications such as full-waveform inversion (FWI) require an estimate of the wavefield to compute the model updates. When dealing with localized problems, it is wasteful to compute these updates in the global domain, when we only need them in our region of interest. This paper introduces a local solver that generates forward and adjoint wavefields which are, to machine precision, identical to those generated by a full-domain solver evaluated within the region of interest. This means that the local solver computes all interactions between model updates within the region of interest and the inhomogeneities in the background model outside. Because no approximations are made in the calculation of the forward and adjoint wavefields, the local solver can compute the identical gradient in the region of interest as would be computed by the more expensive full-domain solver. In this paper, the local solver is used to efficiently generate the FWI gradient at the boundary of a salt body. This gradient is then used in a level set method to automatically update the salt boundary.

Description: Full waveform inversion using the conventional L 2 distance to measure the misfit between seismograms is known to suffer from cycle skipping. An alternative strategy is proposed in this study, based on a measure of the misfit computed with an optimal transport distance. This measure allows to account for the lateral coherency of events within the seismograms, instead of considering each seismic trace independently, as is done generally in full waveform inversion. The computation of this optimal transport distance relies on a particular mathematical formulation allowing for the non-conservation of the total energy between seismograms. The numerical solution of the optimal transport problem is performed using proximal splitting techniques. Three synthetic case studies are investigated using this strategy: the Marmousi 2 model, the BP 2004 salt model, and the Chevron 2014 benchmark data. The results emphasize interesting properties of the optimal transport distance. The associated misfit function is less prone to cycle skipping. A workflow is designed to reconstruct accurately the salt structures in the BP 2004 model, starting from an initial model containing no information about these structures. A high-resolution P -wave velocity estimation is built from the Chevron 2014 benchmark data, following a frequency continuation strategy. This estimation explains accurately the data. Using the same workflow, full waveform inversion based on the L 2 distance converges towards a local minimum. These results yield encouraging perspectives regarding the use of the optimal transport distance for full waveform inversion: the sensitivity to the accuracy of the initial model is reduced, the reconstruction of complex salt structure is made possible, the method is robust to noise, and the interpretation of seismic data dominated by reflections is enhanced.

Description: This paper concerns the numerical modelling of time-domain mechanical waves in viscoelastic media based on a generalized Zener model. To do so, classically in the literature relaxation mechanisms are introduced, resulting in a set of the so-called memory variables and thus in large computational arrays that need to be stored. A challenge is thus to accurately mimic a given attenuation law using a minimal set of relaxation mechanisms. For this purpose, we replace the classical linear approach of Emmerich & Korn with a nonlinear optimization approach with constraints of positivity. We show that this technique is more accurate than the linear approach. Moreover, it ensures that physically meaningful relaxation times that always honour the constraint of decay of total energy with time are obtained. As a result, these relaxation times can always be used in a stable way in a modelling algorithm, even in the case of very strong attenuation for which the classical linear approach may provide some negative and thus unusable coefficients.

Description: Fractures are a source of extra compliance in the rock mass. The mechanical properties of the rock matrix as well as the propagation of seismic waves inside the rock medium are dependent on the magnitude of roughness and offset between the imperfect fracture interfaces. Fracture compliance can estimate the degree of contact between fracture faces, type of fluid filling the fracture and the fracture roughness. To characterize these fracture properties, compliance ratio, known by the ratio of normal-to-shear compliance, can be a potential tool in the subsurface studies to improve the well layout design. The focus of this study is to illustrate how the compliance ratio of a rough fracture, with or without the offset between the fracture faces, can diverge from the compliance ratio of a fracture with smooth interface. Quasi-static and dynamic methods are two common ways to calculate the compliance. The former calculates the compliance by measuring the change in the displacement with the applied stress, while the latter estimates the compliance through monitoring the changes in propagation of seismic waves. To compare the compliance ratios of fractures with imperfect and smooth interfaces in an infinite medium, a numerical finite-element model is built in commercial finite-element software. The imperfect interface of the fracture is modeled with saw-tooth-like structures where they can be partially or fully in contact. The defined saw-tooth-like structures of contact asperities impose an in-plane asymmetry in the shear direction. This asymmetry causes two different values for the compliance in shear direction, known as the soft and stiff shear compliance. Our numerical simulations suggest the increase in the degree of contact between the fracture faces increases the compliance ratio in the stiff direction more than the soft direction. The compliance ratio of the fracture with the imperfect interface is larger than the compliance ratio of the smooth fracture. We suggest that the interlocking and riding up effects at the fracture interface may explain our findings in this study.

Description: Ambient noise correlation method has been widely used to extract surface waves and tomography. The extraction of body waves has been very limited, but recent reports have suggested promises for deep incident waves. Here we report our first observations of triplicated PKP phases (important phases for studying the Earth's core) and confirm observations of other body-wave core phases from noise correlations. We use dense seismic arrays in South America and China Regional Seismic Networks at distances from 145° to the antipode. We can clearly observe different PKP branches (df, bc and ab) in stacks of the station–station correlations. Both ambient noise and earthquake coda contribute to PKP phases. However, the contributions vary with frequency and with body-wave phases. At shorter periods (5–20 s), three branches of PKP (df, bc and ab) can be extracted from ambient noise and the ab phase from earthquake coda. At longer periods (15–50 s), earthquake coda are effective in generating the df branch, but not the ab branch. The generation of the PKIKP phase (df branch) from earthquake coda does not depend on earthquake focal mechanisms or focal depths. However, earthquakes far from the stations contribute more than events closer by. The best coda window is around 10 000–40 000 s and the best magnitude threshold is M w greater than 6.8 or 6.9. The observation of triplicated PKP branches from noise correlations provides a new type of data for studying the Earth's deep interior, in particularly the inner core anisotropy, which overcomes some of the limitations of traditional earthquake-based studies (such as limited source distributions and source location errors).

Description: The project Seismic Hazard Harmonization in Europe (SHARE), completed in 2013, presents significant improvements over previous regional seismic hazard modeling efforts. The Global Strain Rate Map v2.1, sponsored by the Global Earthquake Model Foundation and built on a large set of self-consistent geodetic GPS velocities, was released in 2014. To check the SHARE seismic source models that were based mainly on historical earthquakes and active fault data, we first evaluate the SHARE historical earthquake catalogues and demonstrate that the earthquake magnitudes are acceptable. Then, we construct an earthquake potential model using the Global Strain Rate Map data. SHARE models provided parameters from which magnitude–frequency distributions can be specified for each of 437 seismic source zones covering most of Europe. Because we are interested in proposed magnitude limits, and the original zones had insufficient data for accurate estimates, we combine zones into five groups according to SHARE's estimates of maximum magnitude. Using the strain rates, we calculate tectonic moment rates for each group. Next, we infer seismicity rates from the tectonic moment rates and compare them with historical and SHARE seismicity rates. For two of the groups, the tectonic moment rates are higher than the seismic moment rates of the SHARE models. Consequently, the rates of large earthquakes forecast by the SHARE models are lower than those inferred from tectonic moment rate. In fact, the SHARE models forecast higher seismicity rates than the historical rates, which indicate that the authors of SHARE were aware of the potentially higher seismic activities in the zones. For one group, the tectonic moment rate is lower than the seismic moment rates forecast by the SHARE models. As a result, the rates of large earthquakes in that group forecast by the SHARE model are higher than those inferred from tectonic moment rate, but lower than what the historical data show. For the other two groups, the seismicity rates from tectonic moment rate, historical data and SHARE models are consistent. For four groups, the maximum magnitudes used by SHARE are fairly consistent with the probable maximum magnitudes inferred from tectonic strain rates. This study demonstrates that: (1) the strain-rate data are useful for constraining seismicity rates and magnitude limits; and (2) SHARE seismic source models and historical earthquake catalogues are satisfactory.

Description: We present a method for estimating ratios of P and S waves corner frequencies ( R cf ) and earthquake strain drops by joint analysis of P and S source spectra of neighbouring groups of events. The method is applied systematically to data generated by ~9000 earthquakes around the Karadere segment of the North Anatolian Fault Zone. The results indicate several regions that produce consistently R cf values higher (e.g. 〉2) than expected from classical earthquake source models. These are associated generally with fault sections having strong geometrical heterogeneities, shallow depth sections and/or locations without large pre-existing surface trace. Earthquake ruptures in such regions are likely to generate significant rock damage and tensile components of faulting. To assess whether the observed high R cf values are produced by enriched high frequency P waves, reduced high frequency S waves or both, we compare the associated P and S spectra with mean/median results. The analysis suggests that the high R cf values of shallow events (depth 〈4 km) are generated primarily by reduced high frequency S radiation, and that the contribution from elevated high frequency P radiation increases with depth and proximity to geometrical complexities. The results highlight the importance of considering carefully the existence of some volumetric source components in earthquake rupture processes.

Description: Fractures are a source of extra compliance in the rock mass. The mechanical properties of the rock matrix as well as the propagation of seismic waves inside the rock medium are dependent on the magnitude of roughness and offset between the imperfect fracture interfaces. Fracture compliance can estimate the degree of contact between fracture faces, type of fluid filling the fracture and the fracture roughness. To characterize these fracture properties, compliance ratio, known by the ratio of normal-to-shear compliance, can be a potential tool in the subsurface studies to improve the well layout design. The focus of this study is to illustrate how the compliance ratio of a rough fracture, with or without the offset between the fracture faces, can diverge from the compliance ratio of a fracture with smooth interface. Quasi-static and dynamic methods are two common ways to calculate the compliance. The former calculates the compliance by measuring the change in the displacement with the applied stress, while the latter estimates the compliance through monitoring the changes in propagation of seismic waves. To compare the compliance ratios of fractures with imperfect and smooth interfaces in an infinite medium, a numerical finite-element model is built in commercial finite-element software. The imperfect interface of the fracture is modeled with saw-tooth-like structures where they can be partially or fully in contact. The defined saw-tooth-like structures of contact asperities impose an in-plane asymmetry in the shear direction. This asymmetry causes two different values for the compliance in shear direction, known as the soft and stiff shear compliance. Our numerical simulations suggest the increase in the degree of contact between the fracture faces increases the compliance ratio in the stiff direction more than the soft direction. The compliance ratio of the fracture with the imperfect interface is larger than the compliance ratio of the smooth fracture. We suggest that the interlocking and riding up effects at the fracture interface may explain our findings in this study.

Description: Full waveform inversion using the conventional L 2 distance to measure the misfit between seismograms is known to suffer from cycle skipping. An alternative strategy is proposed in this study, based on a measure of the misfit computed with an optimal transport distance. This measure allows to account for the lateral coherency of events within the seismograms, instead of considering each seismic trace independently, as is done generally in full waveform inversion. The computation of this optimal transport distance relies on a particular mathematical formulation allowing for the non-conservation of the total energy between seismograms. The numerical solution of the optimal transport problem is performed using proximal splitting techniques. Three synthetic case studies are investigated using this strategy: the Marmousi 2 model, the BP 2004 salt model, and the Chevron 2014 benchmark data. The results emphasize interesting properties of the optimal transport distance. The associated misfit function is less prone to cycle skipping. A workflow is designed to reconstruct accurately the salt structures in the BP 2004 model, starting from an initial model containing no information about these structures. A high-resolution P -wave velocity estimation is built from the Chevron 2014 benchmark data, following a frequency continuation strategy. This estimation explains accurately the data. Using the same workflow, full waveform inversion based on the L 2 distance converges towards a local minimum. These results yield encouraging perspectives regarding the use of the optimal transport distance for full waveform inversion: the sensitivity to the accuracy of the initial model is reduced, the reconstruction of complex salt structure is made possible, the method is robust to noise, and the interpretation of seismic data dominated by reflections is enhanced.

Description: We present a novel earthquake location method using acoustic wave-equation-based traveltime inversion. The linear relationship between the location perturbation ( t 0 , x s ) and the resulting traveltime residual t of a particular seismic phase, represented by the traveltime sensitivity kernel K ( t 0 , x s ) with respect to the earthquake location ( t 0 , x s ), is theoretically derived based on the adjoint method. Traveltime sensitivity kernel K ( t 0 , x s ) is formulated as a convolution between the forward and adjoint wavefields, which are calculated by numerically solving two acoustic wave equations. The advantage of this newly derived traveltime kernel is that it not only takes into account the earthquake–receiver geometry but also accurately honours the complexity of the velocity model. The earthquake location is obtained by solving a regularized least-squares problem. In 3-D realistic applications, it is computationally expensive to conduct full wave simulations. Therefore, we propose a 2.5-D approach which assumes the forward and adjoint wave simulations within a 2-D vertical plane passing through the earthquake and receiver. Various synthetic examples show the accuracy of this acoustic wave-equation-based earthquake location method. The accuracy and efficiency of the 2.5-D approach for 3-D earthquake location are further verified by its application to the 2004 Big Bear earthquake in Southern California.

Description: Seismic interferometry applied to 120 hr of railroad traffic recorded by an array of vertical component seismographs along a railway within the Rio Grande rift has recovered surface and body waves characteristic of the geology beneath the railway. Linear and hyperbolic arrivals are retrieved that agree with surface (Rayleigh), direct and reflected P waves observed by nearby conventional seismic surveys. Train-generated Rayleigh waves span a range of frequencies significantly higher than those recovered from typical ambient noise interferometry studies. Direct P -wave arrivals have apparent velocities appropriate for the shallow geology of the survey area. Significant reflected P -wave energy is also present at relatively large offsets. A common midpoint stack produces a reflection image consistent with nearby conventional reflection data. We suggest that for sources at the free surface (e.g. trains) increasing the aperture of the array to record wide angle reflections, in addition to longer recording intervals, might allow the recovery of deeper geological structure from railroad traffic. Frequency–wavenumber analyses of these recordings indicate that the train source is symmetrical (i.e. approaching and receding) and that deeper refracted energy is present although not evident in the time-offset domain. These results confirm that train-generated vibrations represent a practical source of high-resolution subsurface information, with particular relevance to geotechnical and environmental applications.

Description: An approach to depth migration, based on an integral representation of seismic data, that is, wavefields recorded on the boundary, is presented in terms of Poincaré wavelets. Each wavelet is taken as a boundary datum for a high-frequency asymptotic solution of the wave equation. This solution, which we call the quasiphoton or the Gaussian wave packet, decreases in a Gaussian manner away from a point running along a ray that is launched from the surface. The deformation of the propagating packet is taken into account in the migration algorithm. A numerical example of zero-offset migration with synthetic seismograms calculated for the 2-D SEG/EAGE salt model is presented. The result, which uses only 3.9 per cent of the total number of coefficients, is a satisfactory image, with a threshold of 0.75 per cent.

Description: The motion of a self-gravitating hyperelastic body is described through a time-dependent mapping from a reference body into physical space, and its material properties are determined by a referential density and strain-energy function defined relative to the reference body. Points within the reference body do not have a direct physical meaning, but instead act as particle labels that could be assigned in different ways. We use Hamilton's principle to determine how the referential density and strain-energy functions transform when the particle labels are changed, and describe an associated ‘particle relabelling symmetry’. We apply these results to linearized elastic wave propagation and discuss their implications for seismological inverse problems. In particular, we show that the effects of boundary topography on elastic wave propagation can be mapped exactly into volumetric heterogeneity while preserving the form of the equations of motion. Several numerical calculations are presented to illustrate our results.

Description: Scattering due to small-scale heterogeneities in the rock formation surrounding a wellbore can significantly change the acoustic waveform from a logging measurement which in turn can be used to characterize the formation heterogeneities. This study simulates the elastic heterogeneity scattering in monopole and dipole acoustic logging and analyse the resulting effects on the waveforms. The results show that significant coda waves are generated in both monopole and dipole waveforms and the dipole coda is dominated by S -to- S scattering, which can be effectively utilized to diagnose the heterogeneity in the rock formation. The coda wave modelling and analysis were used to characterize dipole acoustic data logged before and after fracturing a reservoir interval, with significant coda wave in the after-fracturing data indicating fracturing-induced heterogeneous property change in the rock volume surrounding the borehole.

Description: We present a novel earthquake location method using acoustic wave-equation-based traveltime inversion. The linear relationship between the location perturbation ( t 0 , x s ) and the resulting traveltime residual t of a particular seismic phase, represented by the traveltime sensitivity kernel K ( t 0 , x s ) with respect to the earthquake location ( t 0 , x s ), is theoretically derived based on the adjoint method. Traveltime sensitivity kernel K ( t 0 , x s ) is formulated as a convolution between the forward and adjoint wavefields, which are calculated by numerically solving two acoustic wave equations. The advantage of this newly derived traveltime kernel is that it not only takes into account the earthquake–receiver geometry but also accurately honours the complexity of the velocity model. The earthquake location is obtained by solving a regularized least-squares problem. In 3-D realistic applications, it is computationally expensive to conduct full wave simulations. Therefore, we propose a 2.5-D approach which assumes the forward and adjoint wave simulations within a 2-D vertical plane passing through the earthquake and receiver. Various synthetic examples show the accuracy of this acoustic wave-equation-based earthquake location method. The accuracy and efficiency of the 2.5-D approach for 3-D earthquake location are further verified by its application to the 2004 Big Bear earthquake in Southern California.

Description: Observations of shear wave splitting in SKS seismic phase play a key role in the current efforts to understand kinematics and dynamics of mantle flow, but azimuthal anisotropy as a depth-localized phenomenon still is poorly known. Here we analyse stratification of seismic azimuthal anisotropy beneath central and northern Anatolia (a microplate within the Alpine belt) by inverting P -wave receiver functions jointly with shear wave splitting in SKS seismic phase. The analysis is based on recordings of stations of the North Anatolian Fault (NAF) passive seismic experiment. In the resulting model in a depth interval from 120 to 200 km fast direction of anisotropy is nearly parallel to the plate motion direction (~E–W), whilst a normal direction (close to S–N) is found in the low velocity zone (LVZ) between 60 and 90 km. Our preferred interpretation of these data suggests that the flow in upper mantle is nearly parallel to the Anatolian plate motion direction in the depth range from the LAB to 200 km, but in part of the LVZ fast direction of anisotropy is normal to the direction of shear in the mantle. This relation between anisotropy and shear is known from laboratory experiments with peridotite-type rock containing melt. A similar relation between anisotropy and flow in the LVZ is found in Fennoscandia. These findings may have far-reaching implications for interpreting mantle anisotropy elsewhere.

Description: On 2013 April 13, an inland earthquake of M w 5.8 occurred in Awaji Island, which forms the western boundary of the Osaka sedimentary basin in western Japan. The strong ground motion data were collected from more than 100 stations within the basin and it was found that in the Osaka Plain, the pseudo velocity response spectra at a period of around 6.5 s were significantly larger than at other stations of similar epicentral distance outside the basin. The ground motion lasted longer than 3 min in the Osaka Plain where its bedrock depth spatially varies from approximately 1 to 2 km. We modelled long-period (higher than 2 s) ground motions excited by this earthquake, using the finite difference method assuming a point source, to validate the present velocity structure model and to obtain better constraint of the attenuation factor of the sedimentary part of the basin. The effect of attenuation in the simulation was included in the form of Q ( f ) =  Q 0 ( f / f 0 ), where Q 0 at a reference frequency f 0 was given by a function of the S -wave velocity, Q 0 = α V S . We searched for appropriate Q 0 values by changing α for a fixed value of f 0 = 0.2 Hz. It was found that values of α from 0.2 to 0.5 fitted the observations reasonably well, but that the value of α  = 0.3 performed best. Good agreement between the observed and simulated velocity waveforms was obtained for most stations within the Osaka Basin in terms of both amplitude and ground motion duration. However, underestimation of the pseudo velocity response spectra in the period range of 5–7 s was recognized in the central part of the Osaka Plain, which was caused by the inadequate modelling of later phases or wave packets in this period range observed approximately 2 min after the direct S -wave arrival. We analysed this observed later phase and concluded that it was a Love wave originating from the direction of the east coast of Awaji Island.

Description: Geodetic slip inversions for three major ( M w 〉 7) strike-slip earthquakes (1992 Landers, 1999 Hector Mine and 2010 El Mayor–Cucapah) show a 15–60 per cent reduction in slip near the surface (depth 〈 2 km) relative to the slip at deeper depths (4–6 km). This significant difference between surface coseismic slip and slip at depth has been termed the shallow slip deficit (SSD). The large magnitude of this deficit has been an enigma since it cannot be explained by shallow creep during the interseismic period or by triggered slip from nearby earthquakes. One potential explanation for the SSD is that the previous geodetic inversions lack data coverage close to surface rupture such that the shallow portions of the slip models are poorly resolved and generally underestimated. In this study, we improve the static coseismic slip inversion for these three earthquakes, especially at shallow depths, by: (1) including data capturing the near-fault deformation from optical imagery and SAR azimuth offsets; (2) refining the interferometric synthetic aperture radar processing with non-boxcar phase filtering, model-dependent range corrections, more complete phase unwrapping by SNAPHU (Statistical Non-linear Approach for Phase Unwrapping) assuming a maximum discontinuity and an on-fault correlation mask; (3) using more detailed, geologically constrained fault geometries and (4) incorporating additional campaign global positioning system (GPS) data. The refined slip models result in much smaller SSDs of 3–19 per cent. We suspect that the remaining minor SSD for these earthquakes likely reflects a combination of our elastic model's inability to fully account for near-surface deformation, which will render our estimates of shallow slip minima, and potentially small amounts of interseismic fault creep or triggered slip, which could ‘make up’ a small percentages of the coseismic SSD during the interseismic period. Our results indicate that it is imperative that slip inversions include accurate measurements of near-fault surface deformation to reliably constrain spatial patterns of slip during major strike-slip earthquakes.

Description: We present a mathematically consistent numerical algorithm for the simulation of earthquake rupture with rate-and-state friction. Its main features are adaptive time stepping, a novel algebraic solution algorithm involving nonlinear multigrid and a fixed point iteration for the rate-and-state decoupling. The algorithm is applied to a laboratory scale subduction zone which allows us to compare our simulations with experimental results. Using physical parameters from the experiment, we find a good fit of recurrence time of slip events as well as their rupture width and peak slip. Computations in 3-D confirm efficiency and robustness of our algorithm.

Description: Seismic moment tensors can provide information on the size and orientation of fractures producing acoustic emissions (AEs) and on the stress conditions in the sample. The moment tensor inversion of AEs is, however, a demanding procedure requiring carefully calibrated sensors and accurate knowledge of the velocity model. In field observations, the velocity model is usually isotropic and time independent. In laboratory experiments, the velocity is often anisotropic and time dependent and attenuation might be significant due to opening or closure of microcracks in the sample during loading. In this paper, we study the sensitivity of the moment tensor inversion to anisotropy of P -wave velocities and attenuation. We show that retrieved moment tensors critically depend on anisotropy and attenuation and their neglect can lead to misinterpretations of the source mechanisms. The accuracy of the inversion also depends on the fracturing mode of AEs: tensile events are more sensitive to P -wave anisotropy and attenuation than shear events. We show that geometry of faulting in anisotropic rocks should be studied using the source tensors, since the P - and T -axes of the moment tensors are affected by velocity anisotropy and deviate from the true orientation of faulting. The stronger the anisotropy is, the larger the deviations are. Finally, we prove that the moment tensor inversion applied to a large dataset of AEs can be utilized to provide information on the attenuation parameters of the rock sample. The method is capable of measuring anisotropic attenuation in the sample and allows for detection of dilatant cracking according to the stress regime.

Description: We investigate the lithospheric structure of Cameroon inverting Rayleigh waves obtained from the cross-correlation of ambient seismic noise. We correlate seismic records between 32 broad-band stations and we obtain good quality Rayleigh waves for 310 interstation paths. We measure group velocity dispersion curves from the reconstructed Rayleigh waves in the period range 10–35 s and we invert the group velocities for tomographic images. After the tomography the group velocities are then inverted, together with longer period group velocity measurements from existing literature, to compute a 3-D S -wave velocity model of the Cameroon lithosphere down to 100 km depth. Our results provide an unprecedented mapping of the physical properties of the different crustal units and their correlations with surface geology, as well as with mantle lithospheric variations. The Cameroon Volcanic Line (CVL) appears as a segmented feature exhibiting different physical properties along strike. The active Mt Cameroon volcano is underlain by very low velocities, unlike the other segments of the CVL. The along-strike variations in crustal structure suggest that lateral heterogeneities in lithospheric thickness and physical properties have influenced the location and distribution of magmatism. The crust beneath the Central African Shear Zone exhibits a sizeable low velocity anomaly. The lithosphere beneath Cameroon is characterised by a heterogeneous crust with a relatively constant thickness and a low velocity uppermost mantle at the edge of the Congo Craton. Our results favour processes combining small-scale upwelling at the edge of a thick lithosphere and reactivation of Precambrian basement structures to explain the distribution of Holocene-Recent magmatism and plateau uplift. Our results also indicate that Mt Cameroon and surroundings areas are the most at risk zones for magmatic activity during this stage of CVL development.

Description: Although optimal, computing the moment tensor solution is not always a viable option for the calculation of the size of an earthquake, especially for small events (say, below M w 2.0). Here we show an alternative approach to the calculation of the moment-rate spectra of small earthquakes, and thus of their scalar moments, that uses a network-based calibration of crustal wave propagation. The method works best when applied to a relatively small crustal volume containing both the seismic sources and the recording sites. In this study we present the calibration of the crustal volume monitored by the High-Resolution Seismic Network (HRSN), along the San Andreas Fault (SAF) at Parkfield. After the quantification of the attenuation parameters within the crustal volume under investigation, we proceed to the spectral correction of the observed Fourier amplitude spectra for the 100 largest events in our data set. Multiple estimates of seismic moment for the all events (1811 events total) are obtained by calculating the ratio of rms-averaged spectral quantities based on the peak values of the ground velocity in the time domain, as they are observed in narrowband-filtered time-series. The mathematical operations allowing the described spectral ratios are obtained from Random Vibration Theory (RVT). Due to the optimal conditions of the HRSN, in terms of signal-to-noise ratios, our network-based calibration allows the accurate calculation of seismic moments down to M w 〈 0. However, because the HRSN is equipped only with borehole instruments, we define a frequency-dependent Generalized Free-Surface Effect (GFSE), to be used instead of the usual free-surface constant F = 2. Our spectral corrections at Parkfield need a different GFSE for each side of the SAF, which can be quantified by means of the analysis of synthetic seismograms. The importance of the GFSE of borehole instruments increases for decreasing earthquake's size because for smaller earthquakes the bandwidth available for our calculations is consistently shifted towards higher frequencies.

Description: As a result of the accumulation of high-resolution observation data, 3-D high-fidelity crustal structure data for large domains are becoming available. However, it has been difficult to use such data to perform elastic/viscoelastic crustal deformation analyses in large domains with quality assurance of the numerical simulation that guarantees convergence of the numerical solution with respect to the discretization size because the costs of analysis are significantly high. This paper proposes a method of constructing a high-fidelity crustal structure finite element (FE) model using high-fidelity crustal structure data and fast FE analysis to reduce the costs of analysis (based on automatic FE model generation for parallel computation, OpenMP/MPI hybrid parallel computation on distributed memory computers, a geometric multigrid, variable preconditioning and multiple precision arithmetic). Using the proposed methods, we construct 10 billion degree-of-freedom high-fidelity crustal structure FE models for the entire Japan, and conduct elastic/viscoelastic crustal deformation analysis using this model with enough high accuracy of the numerical simulation.

Description: The EarthScope USArray provides an opportunity to obtain detailed images of the continental upper mantle at an unprecedented scale. The majority of mantle models derived from USArray data to date contain spatial variations in seismic-wave speed; however, in many cases these data sets do not by themselves allow a non-unique interpretation. Joint interpretation of seismic attenuation and velocity models can improve upon the interpretations based only on velocity and provide important constraints on the temperature, composition, melt content, and volatile content of the mantle. The surface wave amplitudes that constrain upper-mantle attenuation are sensitive to factors in addition to attenuation, including the earthquake source excitation, focusing and defocusing by elastic structure, and local site amplification. Because of the difficulty of isolating attenuation from these other factors, little is known about the attenuation structure of the North American upper mantle. In this study, Rayleigh wave traveltime and amplitude in the period range 25–100 s are measured using an interstation cross-correlation technique, which takes advantage of waveform similarity at nearby stations. Several estimates of Rayleigh wave attenuation and site amplification are generated at each period, using different approaches to separate the effects of attenuation and local site amplification on amplitude. It is assumed that focusing and defocusing effects can be described by the Laplacian of the traveltime field. All approaches identify the same large-scale patterns in attenuation, including areas where the attenuation values are likely contaminated by unmodelled focusing and defocusing effects. Regionally averaged attenuation maps are constructed after removal of the contaminated attenuation values, and the variations in intrinsic shear attenuation that are suggested by these Rayleigh wave attenuation maps are explored.

Description: Surface wave magnitude M s for a compilation of 72 nuclear tests detonated in hard rock media for which yields and burial depths have been reported in the literature is shown to scale with yield W as a + b x log[ W ], where a = 2.50 ± 0.08 and b = 0.80 ± 0.05. While the exponent b is consistent with an M s scaling model for fully coupled, normal containment-depth explosions, the intercept a is offset 0.45 magnitude units lower than the model. The cause of offset is important to understand in terms of the explosion source. Hard rock explosions conducted in extensional and compressional stress regimes show similar offsets, an indication that the tectonic setting in which an explosion occurs plays no role causing the offset. The scaling model accounts for the effects of source medium material properties on the generation of 20-s period Rayleigh wave amplitudes. Aided by thorough characterizations of the explosion and tectonic release sources, an extensive analysis of the 1963 October 26 Shoal nuclear test detonated in granite 27 miles southeast of Fallon NV shows that the offset is consistent with the predictions of a material damage source model related to non-linear stress wave interactions with the free surface. This source emits Rayleigh waves with polarity opposite to waves emitted by the explosion. The Shoal results were extended to analyse surface waves from the 1962 February 15 Hardhat nuclear test, the 1988 September 14 Soviet Joint Verification Experiment, and the anomalous 1979 August 18 northeast Balapan explosion which exhibits opposite polarity, azimuth-independent source component U 1 compared to an explosion. Modelling these tests shows that Rayleigh wave amplitudes generated by the damage source are nearly as large as or larger than amplitudes from the explosion. As such, destructive interference can be drastic, introducing metastable conditions due to the sensitivity of reduced amplitudes to Rayleigh wave initial phase angles of the explosion and damage sources. This meta-stability is a likely source of scatter in M s -yield scaling observations. The agreement of observed scaling exponent b with the model suggests that the damage source strength does not vary much with yield, in contrast to explosions conducted in weak media where M s scaling rates are greater than the model predicts, and the yield dependence of the damage source strength is significant. This difference in scaling behaviour is a consequence of source medium material properties.

Description: We present a systematic study of seismic velocity changes associated with a megathrust and five strong crustal earthquakes in Japan. We perform both cross-correlation and single-station cross-correlation analysis for station pairs and stations, respectively. The correlation of ambient seismic noise allows us to reconstruct the Green's functions of the wave propagation. By relating the coda parts of the daily Green's functions with the long-term reference Green's functions, shear wave velocity changes are determined. We analyse data from four areas in Japan where large earthquakes occurred: Iwate-Miyagi (2008 M W 6.9 Iwate-Miyagi Nairiku earthquake), Niigata (2004 M W 6.6 Chūetsu, 2007 M W 6.6 Chūetsu-oki and 2011 M W 6.2 Nagano/Niigata earthquakes), Noto Peninsula (2007 M W 6.7 Noto Hantō earthquake) and Fukuoka (2005 M W 6.6 Fukuoka earthquake). In all areas, we analyse time-series which start before the respective earthquakes and last until after the 2011 M W 9.0 Tōhoku-oki earthquake. The analysis in five different frequency ranges between 0.125 and 4.0 Hz yields time-series of the velocity changes for the different station pairs or stations. At the time of the respective earthquakes, we observe coseismic velocity drops in all areas which are followed by a partial post-seismic recovery process. For the Tōhoku-oki earthquake, coseismic velocity drops can also be observed in all regions. There is a general trend of increasing coseismic velocity drops with frequency in all four areas. The largest coseismic drops are observed close to the fault zones. Over the observed time range, the post-seismic recovery is only partial and around half of the coseismic velocity drops do not recover. The characteristic recovery times for the recovering part are similar in all areas and frequency ranges, with an average value of 0.55 yr. We model the volumetric strain changes for the different earthquakes and find that the observed pattern of the coseismic velocity drops cannot be explained by these models. The coseismic velocity drops at the different stations are better related with the peak ground velocities and the associated dynamic strain than with the peak ground accelerations, but the correlation is still poor. This suggests that non-linear effects caused by the strong ground motion during the earthquake can explain at least part of the coseismic velocity drops.

Description: Subducting slabs create strong short wavelength seismic anomalies in the upper mantle where much of Earth's seismicity is located. As such, they have the potential to bias longer wavelength seismic tomography models. To evaluate the effect of subducting slabs in global tomography, we performed a series of inversions using a global synthetic shear wave traveltime data set for a theoretical slab model based on predicted thermal anomalies within slabs. The spectral element method was applied to predict the traveltime anomalies produced by the 3-D slab model for paths corresponding to our current data used in actual tomography models. Inversion tests have been conducted first using the raw traveltime anomalies to check how well the slabs can be imaged in global tomography without the effect of earthquake mislocation. Our results indicate that most of the slabs can be identified in the inversion result but with smoothed and reduced amplitude. The recovery of the total mass anomaly in slab regions is about 88 per cent. We then performed another inversion test to investigate the effect of mislocation caused by subducting slabs. We found that source mislocation largely removes slab signal and significantly degrades the imaging of subducting slabs—potentially reducing the recovery of mass anomalies in slab regions to only 41 per cent. We tested two source relocation procedures—an iterative relocation inversion and joint relocation inversion. Both methods partially recover the true source locations and improve the inversion results, but the joint inversion method worked significantly better than the iterative method. In all of our inversion tests, the amplitudes of artefact structures in the lower mantle caused by the incorrect imaging of slabs (up to ~0.5 per cent S velocity anomalies) are comparable to some large-scale lower-mantle heterogeneities seen in global tomography studies. Based on our inversion tests, we suggest including a-priori subducting slabs in the starting models in global tomography studies and use joint relocation in the inversion.

Description: The crust at the junction of Qinling Orogen, Yangtze Craton and NE Tibetan Plateau bears imprints of the Triassic collision and later intracontinental orogeny between the Qinling Orogen and the Yangtze Craton, and the Cenozoic growth of Tibetan Plateau. Investigating detailed crustal structures in this region helps to better understand these tectonic processes. In this study, we construct a 3-D crustal Vs model using seismic ambient noise data recorded at 321 seismic stations. Ambient noise tomography is performed to generate Rayleigh wave phase velocity maps at 8–50 s periods, which are then inverted for a 3-D isotropic Vs model using a Bayesian Monte Carlo method. Our 3-D model reveals deep-rooted high velocities beneath the Hannan-Micang and Shennong-Huangling Domes, which are located on the west and east sides of the Dabashan Orocline. Similar high velocities are observed in the upper/mid crust of the western Qinling Orogen. We suggest the crustal-scale bodies with high velocity beneath the two domes and the western Qinling Orogen may represent mechanically strong rocks, which not only assisted the formation of the major Dabashan Orocline during late Mesozoic intracontinental orogeny, but also have impeded the northeastward expansion of the Tibetan Plateau during the Cenozoic era.

Description: The collective excitation of city structures by a seismic wavefield and the subsequent multiple Structure-Soil-Structure Interactions (SSSIs) between the buildings are usually disregarded in conventional seismology and earthquake engineering practice. The objective here is to qualify and quantify these complex multiple SSSIs through the design of an elementary study case, which serves as a benchmark for theoretical, numerical and experimental crossed-analysis. The experimental specimen consists of an idealized site-city setup with up to 37 anisotropic resonant structures arranged at the top surface of an elastic layer and in co-resonance with it. The experimental data from shaking table measurements is compared with the theoretical and numerical results provided respectively by an equivalent city-impedance model derived analytically from homogenization in the long-wavelength approximation and a model based on boundary elements. The signatures of the site-city interactions are identified in the frequency, time and space domain, and in particular consist of a frequency-dependent free/rigid switch in the surface condition at the city resonance, beatings in the records and the depolarization of the wavefield. A parametric study on the city density shows that multiple SSSIs among the city structures (five are sufficient) can have significant effects on both the seismic response of its implantation site and that of the buildings. Key parameters are provided to assess site-city interactions in the low seismic frequency range: They involve the mass and rigidity of the city compared to those of the soil and the damping of the building.

Description: We present a novel 3-D pre-stack Kirchhoff depth migration (PKDM) method for teleseismic receiver functions. The proposed algorithm considers the effects of diffraction, scattering and traveltime alteration caused by 3-D volumetric heterogeneities. It is therefore particularly useful for imaging complex 3-D structures such as dipping discontinuities, which is hard to accomplish with traditional methods. The scheme is based on the acoustic wave migration principle, where at each time step of the receiver function, the energy is migrated back to the ensemble of potential conversion points in the image, given a smooth 3-D reference model. Traveltimes for P and S waves are computed with an efficient eikonal solver, the fast marching method. We also consider elastic scattering patterns, where the amplitude of converted S waves depends on the angle between the incident P wave and the scattered S wave. Synthetic experiments demonstrate the validity of the method for a variety of dipping angle discontinuities. Comparison with the widely used common conversion point (CCP) stacking method reveals that our migration shows considerable improvement. For example, the effect of multiple reflections that usually produce apparent discontinuities is avoided. The proposed approach is practical, computationally efficient, and is therefore a potentially powerful alternative to standard CCP methods for imaging large-scale continental structure under dense networks.

Description: Wavelet estimation is an essential step in qualitatively and quantitatively analysing and interpreting seismic data. Applications span from seismic data quality assessment to well ties and seismic inversion. Wavelet estimation methods can be roughly separated into two approaches, data driven inversion methods and analytical definitions. We present a new analytical wavelet definition, which is based on Hermite basis functions. This wavelet model contains four parameters, which correspond to wavelet magnitude, phase, wavelet length and bandwidth. One of our main motivations for this development was to define a compact wavelet representation and an intrinsic parameter uncertainty assessment workflow, which allows us to quantify uncertainties in estimated wavelets, as well as the generation of wavelet realizations to be used, for example, in statistical seismic amplitude inversions. We present a statistical workflow to estimate the model parameters and to explore their posterior uncertainties given well log data and seismic amplitude data. This includes sampling the posterior distribution of the four wavelet parameters using Markov Chain Monte Carlo methods. We then discuss the applicability, limitations and challenges of the approach with the help of synthetic data and a North Sea data set with well logs and processed seismic amplitudes, where we also compare our method to Bayesian least-squares and a commercial wavelet estimation routine. Realizations of wavelets based on the optimized parameters and their uncertainties appear to sample the wavelet space well with reasonable variations in wavelet length, phase and amplitude while not introducing random fluctuations or wavelet lobes. The results indicate that the compact wavelet representation allows for an efficient and rather stable wavelet estimation workflow that achieves useful results in the presence of noisy data.

Description: We investigate the upper crustal velocity structure beneath the Tehri region of the Garhwal Himalaya. The investigated region is situated within the 700-km-long central seismic gap of the Himalaya that has experienced three gap-filling earthquakes since 1991 including the recent 2015 Nepal earthquake ( M w 7.8). The local tomographic inversion is based on a data set of 1365 events collected from 2008 January to 2012 December by a 12-station local network that covers an area of about 100  x  80 km around Tehri Dam. We perform a simultaneous inversion for P- and S -wave velocity anomalies. Tomograms are interpreted in the backdrop of the regional geological and tectonic framework of the region. The spatial distribution of relocated events from the 3-D velocity model has shed new light on the pattern of seismicity in the vicinity of the Main Central thrust (MCT), and has elucidated the structure of the underthrusting Indian plate. Our model exhibits a significant negative velocity anomaly up to ~5 per cent beneath the central part of the Garhwal Inner Lesser Himalaya, and a P -wave low velocity anomaly near the Chamoli region. The seismicity zone around the Chamoli region may be attributed to the presence of fluid-filled rocks. Furthermore, an area with ~3–4 per cent positive velocity anomaly is delineated to the northwest of the Uttarkashi thrust in the vicinity of the MCT. Significant findings of the study include: a flat-ramp-flat-type subsurface geometry of the underthrusting Indian plate below the Garhwal Himalaya, high-velocity images representing the trend and configuration of Delhi–Haridwar ridge below the Sub Himalaya and Lesser Himalaya and a seismically active zone representing geometrical asperity on the basement thrust in the vicinity of the MCT.

Description: Slow-slip events (SSEs) are common at subduction zone faults where large mega earthquakes occur. We report here that one of the best-recorded moderate size continental earthquake, the 2009 April 6 moment magnitude ( M w ) 6.3 L'Aquila (Italy) earthquake, was preceded by a 5.9 M w SSE that originated from the decollement beneath the reactivated normal faulting system. The SSE is identified from a rigorous analysis of continuous GPS stations and occurred on the 12 February and lasted for almost two weeks. It coincided with a burst in the foreshock activity with small repeating earthquakes migrating towards the main-shock hypocentre as well as with a change in the elastic properties of rocks in the fault region. The SSE has caused substantial stress loading at seismogenic depths where the magnitude 4.0 foreshock and M w 6.3 main shock nucleated. This stress loading is also spatially correlated with the lateral extent of the aftershock sequence.

Description: For more than two decades, the number of data and model parameters in seismic tomography problems has exceeded the available computational resources required for application of direct computational methods, leaving iterative solvers the only option. One disadvantage of the iterative techniques is that the inverse of the matrix that defines the system is not explicitly formed, and as a consequence, the model resolution and covariance matrices cannot be computed. Despite the significant effort in finding computationally affordable approximations of these matrices, challenges remain, and methods such as the checkerboard resolution tests continue to be used. Based upon recent developments in sparse algorithms and high-performance computing resources, we show that direct methods are becoming feasible for large seismic tomography problems. We demonstrate the application of QR factorization in solving the regional P -wave structure and computing the full resolution matrix with 267 520 model parameters.

Description: Inhomogeneous noise sources surrounding stations produce asymmetric amplitudes in cross-correlation functions that yield preferential source directions. Here we show that preprocessing biases the dominant source direction estimate towards the source producing long-duration signals by down-weighting high-amplitude signals. Tests with both synthetic data and observations show that conventional preprocessing, where only earthquakes and local transients (e.g. trawling, fish impacts) are removed, is more sensitive to coherent energy, while one-bit preprocessing and running-absolute-mean preprocessing are more influenced by signal duration. Comparisons between different preprocessing methods are made on data from the Cascadia Initiative ocean bottom seismometer array, where we find that the total energy arriving from pelagic and coastal areas is similar. Moreover, pelagic-generated signals tend to be weaker but have longer duration, in contrast to coastal-generated signals that tend to be stronger but have shorter duration.

Description: Understanding earthquake clustering in space and time is important but also challenging because of complexities in earthquake patterns and the large and diverse nature of earthquake catalogues. Swarms are of particular interest because they likely result from physical changes in the crust, such as slow slip or fluid flow. Both swarms and clusters resulting from aftershock sequences can span a wide range of spatial and temporal scales. Here we test and implement a new method to identify seismicity clusters of varying sizes and discriminate them from randomly occurring background seismicity. Our method searches for the closest neighbouring earthquakes in space and time and compares the number of neighbours to the background events in larger space/time windows. Applying our method to California's San Jacinto Fault Zone (SJFZ), we find a total of 89 swarm-like groups. These groups range in size from 0.14 to 7.23 km and last from 15 min to 22 d. The most striking spatial pattern is the larger fraction of swarms at the northern and southern ends of the SJFZ than its central segment, which may be related to more normal-faulting events at the two ends. In order to explore possible driving mechanisms, we study the spatial migration of events in swarms containing at least 20 events by fitting with both linear and diffusion migration models. Our results suggest that SJFZ swarms are better explained by fluid flow because their estimated linear migration velocities are far smaller than those of typical creep events while large values of best-fitting hydraulic diffusivity are found.

Description: High-frequency seismic ground motion (10+ Hz), as needed for earthquake engineering design purposes, is largely controlled by the metre-scale structure of the earth's crust. However, the state-of-the-art velocity models poorly resolve small-scale features of the subsurface velocity and density variation. We invert 35 sonic logs (up to 3000 m in depth) in and near Los Angeles basin, CA, to obtain a statistical description of the small-scale heterogeneities of the basin. Assuming a von Karman autocorrelation function, our analysis finds that Hurst numbers, , between 0.0 and 0.2, vertical correlation lengths, a z , of 15–150 m and standard deviations of about 5 per cent characterize the variability in the borehole data. We report average parameters for Los Angeles basin of = 0.064 (0.058, 0.069) ± 0.01 (0.006, 0.012) and a z = 54 (51.1, 57.6) ± 5.9 (1.79, 9.53) m with 95 per cent confidence intervals listed in the parentheses. Despite the large depth range of the logs, there is no significant variation of the statistical parameters with depth. Our analysis of 371 depth-averaged shear wave velocities in the upper 30 m, V s 30 , provides only an upper bound of basin scale-length estimates due to the coarse sampling distance, with a Hurst number of about 0.3 and lateral correlation lengths, a x , of 5–10 km.

Description: A subsurface extended image gather obtained during controlled-source depth imaging yields a blurred kernel of an interface reflection operator. This reflectivity kernel or reflection function is comprised of the interface plane-wave reflection coefficients and so, in principle, the gather contains amplitude versus offset or angle information. We present a modelling theory for extended image gathers that accounts for variable illumination and blurring, under the assumption of a good migration-velocity model. The method involves forward modelling as well as migration or back propagation so as to define a receiver-side blurring function, which contains the effects of the detector array for a given shot. Composition with the modelled incident wave and summation over shots then yields an overall blurring function that relates the reflectivity to the extended image gather obtained from field data. The spatial evolution or instability of blurring functions is a key concept and there is generally not just spatial blurring in the apparent reflectivity, but also slowness or angle blurring. Gridded blurring functions can be estimated with, for example, a reverse-time migration modelling engine. A calibration step is required to account for ad hoc band limitedness in the modelling and the method also exploits blurring-function reciprocity. To demonstrate the concepts, we show numerical examples of various quantities using the well-known SIGSBEE test model and a simple salt-body overburden model, both for 2-D. The moderately strong slowness/angle blurring in the latter model suggests that the effect on amplitude versus offset or angle analysis should be considered in more realistic structures. Although the description and examples are for 2-D, the extension to 3-D is conceptually straightforward. The computational cost of overall blurring functions implies their targeted use for the foreseeable future, for example, in reservoir characterization. The description is for scalar waves, but the extension to elasticity is foreseeable and we emphasize the separation of the overburden and survey-geometry blurring effects from the nature of the target scatterer.

Description: We analyse the influence of the solid Earth tides and ocean loading on the occurrence time of Southern California earthquakes. For each earthquake, we calculate tidal Coulomb failure stress and stress rate on a fault plane that is assumed to be controlled by the orientation of the adjacent fault. To reduce bias when selecting data for testing the tide-earthquake relationship, we create four earthquake catalogues containing events within 1, 1.5, 2.5 and 5 km of nearest faults. We investigate the difference in seismicity rates at times of positive and negative tidal stresses/stress rates given three different cases. We consider seismicity rates during times of positive versus negative stress and stress rate, as well as 2 and 3 hr surrounding the local tidal stress extremes. We find that tidal influence on earthquake occurrence is found to be statistically non-random only in close proximity to tidal extremes meaning that magnitude of tidal stress plays an important role in tidal triggering. A non-random tidal signal is observed for the reverse events. Along with a significant increase in earthquake rates around tidal Coulomb stress maxima, the strength of tidal correlation is found to be closely related to the amplitude of the peak tidal Coulomb stress ( p ). The most effective tidal triggering is found for p ≥ 1 kPa, which is much smaller than thresholds suggested for static and dynamic triggering of aftershocks.

Description: We use seismic waves for a magnitude 6.3 intermediate-depth (160 km) earthquake in the Bucaramanga Nest, Colombia, to infer a complex rupture process with two distinct stages, characterized by different rupture velocities possibly controlled by the evolution of strength on the fault. Our integrated data processing permitted to precisely characterize the multistage rupture and the presence of a strong weakening event. The resulting seismic radiation is interpreted as resulting from an extreme weakening due to a cascading thermal shear runaway, with an initial inefficient radiation process followed by a fast and dynamic efficient rupture. Our results imply dynamic complexity of the seismic rupture deep inside the Earth, and may help to give some new insights about the physical mechanism of intermediate-depth earthquakes.

Description: We investigated the shallow structure of the Solfatara, a volcano within the Campi Flegrei caldera, southern Italy, using surface waves as a diagnostic tool. We analysed data collected during the RICEN campaign, where a 3-D active seismic experiment was performed on a dense regular grid of 90 m  x  115 m using a Vibroseis as the seismic source. After removal of the source time function, we analysed the surface wave contribution to the Green's function. Here, a 1-D approximation can hold for subgrids of 40 m  x  40 m. Moreover, we stacked all of the signals in the subgrid according to source–receiver distance bins, despite the absolute location of the source and the receiver, to reduce the small-scale variability in the data. We then analysed the resulting seismic sections in narrow frequency bands between 7 and 25 Hz. We obtained phase and group velocities from a grid search, and a cost function based on the spatial coherence of both the waveforms and their envelopes. We finally jointly inverted the dispersion curves of the phase and group velocities to retrieve a 1-D S -wave model local to the subgrid. Together, the models provided a 3-D description of the S -wave model in the area. We found that the maximum penetration depth is 15 m. In the first 4 m, we can associate the changes in the S -wave field to the temperature gradient, while at greater depths, the seismic images correlate with the resistivity maps, which indicate the water layer close to the Fangaia area and an abrupt variation moving towards the northeast.

Description: Microseisms in the period of 2–10 s are generated in deep oceans and near coastal regions. It is common for microseisms from multiple sources to arrive at the same time at a given seismometer. It is therefore desirable to be able to measure multiple slowness vectors accurately. Popular ways to estimate the direction of arrival of ocean induced microseisms are the conventional (fk) or adaptive (Capon) beamformer. These techniques give robust estimates, but are limited in their resolution capabilities and hence do not always detect all arrivals. One of the limiting factors in determining direction of arrival with seismic arrays is the array response, which can strongly influence the estimation of weaker sources. In this work, we aim to improve the resolution for weaker sources and evaluate the performance of two deconvolution algorithms, Richardson–Lucy deconvolution and a new implementation of CLEAN-PSF. The algorithms are tested with three arrays of different aperture (ASAR, WRA and NORSAR) using 1 month of real data each and compared with the conventional approaches. We find an improvement over conventional methods from both algorithms and the best performance with CLEAN-PSF. We then extend the CLEAN-PSF framework to three components (3C) and evaluate 1 yr of data from the Pilbara Seismic Array in northwest Australia. The 3C CLEAN-PSF analysis is capable in resolving a previously undetected Sn phase.

Description: We investigated a detailed spatial distribution of principal stress axis orientations in the source area of the 2003 M6.4 Northern Miyagi Prefecture earthquake that occurred in the forearc of northeastern Japan. Aftershock hypocentres were precisely relocated by applying the double difference method to arrival time data obtained at temporary stations as well as at surrounding routine stations. We picked many P -wave polarity data from seismograms at these stations, which enabled us to obtain 312 well-determined focal mechanism solutions. Stress tensor inversions were performed by using these focal mechanism data. The results show that quite a lot of focal mechanisms are difficult to explain by the uniform stress field, especially near the large slip area of the main-shock rupture. Stress tensor inversions at the location of individual earthquakes show that 1 axes are orientated mainly to WSW–ENE in the northern part of the source area, while they are oriented to NW–SE in the southern part. This spatial pattern is roughly similar to those of the static stress change by the main shock, which suggests that the observed spatially heterogeneous stress field was formed by the static stress change. If this is the case, the deviatoric stress magnitude before the main shock was very small. Another possibility is the heterogeneous stress field observed after the main shock had existed even before the main shock, although we do not know why it was formed. Unfavourable orientation of the main shock fault with respect to this stress field suggests that the fault is not strong in this case too.

Description: Gassmann's fluid substitution theory is commonly used to predict seismic velocity change upon change in saturation, and is hence essential for 4-D seismic and AVO studies. This paper addresses the basics assumptions of the Gassmann theory, in order to see how well they are fulfilled in controlled laboratory experiments. Our focus is to investigate the sensitivity of shear modulus to fluid saturation, and the predictability of Gassmann's fluid substitution theory for P -wave modulus. Ultrasonic P - and S -wave velocities in dry and saturated (3.5 wt per cent NaCl) unconsolidated clean sands (Ottawa and Columbia) were measured in an oedometer test system (uniaxial strain conditions) over a range of 0.5–10 MPa external vertical stress. This study shows shear modulus hardening upon brine saturation, which is consistent with previous data found in the literature. Analysis of the data shows that most of the hardening of the ultrasonic shear modulus may be explained by Biot dispersion. Isotropic Gassmann's fluid substitution is found to underestimate the P -wave modulus upon fluid saturation. However, adding the Biot dispersion effect improves the prediction. More work is required to obtain good measurements of parameters influencing dispersion, such as tortuosity, which is very ambiguous and challenging to measure accurately.

Description: Modern tectonic studies often use regional moment tensors (RMTs) to interpret the seismotectonic framework of an earthquake or earthquake sequence; however, despite extensive use, little existing work addresses RMT parameter uncertainty. Here, we quantify how network geometry and faulting style affect RMT sensitivity. We examine how data-model fits change with fault plane geometry (strike and dip) for varying station configurations. We calculate the relative data fit for incrementally varying geometries about a best-fitting solution, applying our workflow to real and synthetic seismograms for both real and hypothetical station distributions and earthquakes. Initially, we conduct purely observational tests, computing RMTs from synthetic seismograms for hypothetical earthquakes and a series of well-behaved network geometries. We then incorporate real data and station distributions from the International Maule Aftershock Deployment (IMAD), which recorded aftershocks of the 2010 M W 8.8 Maule earthquake, and a set of regional stations capturing the ongoing earthquake sequence in Oklahoma and southern Kansas. We consider RMTs computed under three scenarios: (1) real seismic records selected for high data quality; (2) synthetic seismic records with noise computed for the observed source-station pairings and (3) synthetic seismic records with noise computed for all possible station-source pairings. To assess RMT sensitivity for each test, we observe the ‘fit falloff’, which portrays how relative fit changes when strike or dip varies incrementally; we then derive the ranges of acceptable strikes and dips by identifying the span of solutions with relative fits larger than 90 per cent of the best fit. For the azimuthally incomplete IMAD network, Scenario 3 best constrains fault geometry, with average ranges of 45° and 31° for strike and dip, respectively. In Oklahoma, Scenario 3 best constrains fault dip with an average range of 46°; however, strike is best constrained by Scenario 1, with a range of 26°. We draw two main conclusions from this study. (1) Station distribution impacts our ability to constrain RMTs using waveform time-series; however, in some tectonic settings, faulting style also plays a significant role and (2) increasing station density and data quantity (both the number of stations and the number of individual channels) does not necessarily improve RMT constraint. These results may be useful when organizing future seismic deployments (e.g. by concentrating stations in alignment with anticipated nodal planes), and in computing RMTs, either by guiding a more rigorous data selection process for input data or informing variable weighting among the selected data (e.g. by eliminating the transverse component when strike-slip mechanisms are expected).

Description: The polarity of the first motion of a seismic signal from an earthquake is an important constraint in earthquake source inversion. Microseismic events often have low signal-to-noise ratios, which may lead to difficulties estimating the correct first-motion polarities of the arrivals. This paper describes a probabilistic approach to polarity picking that can be both automated and combined with manual picking. This approach includes a quantitative estimate of the uncertainty of the polarity, improving calculation of the polarity probability density function for source inversion. It is sufficiently fast to be incorporated into an automatic processing workflow. When used in source inversion, the results are consistent with those from manual observations. In some cases, they produce a clearer constraint on the range of high-probability source mechanisms, and are better constrained than source mechanisms determined using a uniform probability of an incorrect polarity pick.

Description: Acoustic and gravity waves propagating in planetary atmospheres have been studied intensively as markers of specific phenomena such as tectonic events or explosions or as contributors to atmosphere dynamics. To get a better understanding of the physics behind these dynamic processes, both acoustic and gravity waves propagation should be modelled in a 3-D attenuating and windy atmosphere extending from the ground to the upper thermosphere. Thus, in order to provide an efficient numerical tool at the regional or global scale, we introduce a finite difference in the time domain (FDTD) approach that relies on the linearized compressible Navier–Stokes equations with a background flow (wind). One significant benefit of such a method is its versatility because it handles both acoustic and gravity waves in the same simulation, which enables one to observe interactions between them. Simulations can be performed for 2-D or 3-D realistic cases such as tsunamis in a full MSISE-00 atmosphere or gravity-wave generation by atmospheric explosions. We validate the computations by comparing them to analytical solutions based on dispersion relations in specific benchmark cases: an atmospheric explosion, and a ground displacement forcing.

Description: We measure ellipticity of teleseismic Rayleigh waves at 95 seismic stations in Northern Italy, for wave period between 10 and 110 s, using an automatic technique and a large volume of high-quality seismic recordings from over 500 global earthquakes that occurred in 2008–2014. Northern Italy includes a wide range of crustal structures, from the wide and deep Po Plain sedimentary basin to outcropping sedimentary and crystalline rocks in the Northern Apennines and Alps. It thus provides an excellent case for studying the influence of shallow earth structure on polarization of surface waves. The ellipticity measurements show excellent spatial correlation with geological features in the region, such as high ellipticity associated with regions of low seismic velocity in the Po Plain and low ellipticity values in faster, hard rock regions in the Alps and Apennine mountains. Moreover, the observed ellipticity values also relate to the thickness of the basement, as highlighted by observed differences beneath the Alps and the Apennines. Comparison between observations and predicted ellipticity from a reference crustal model of the region show substantial fit, particularly for T ~ 38 s data. Discrepancy for shorter wave period suggests that slight modifications of the model are needed, and that the ellipticity measurements could help to better constrain the shallow crustal structure of the region. Predictions for the Po Plain are larger than the observations by a factor of four or more and transition from retrograde to prograde Rayleigh wave motion at the surface for periods of T ~ 10–13 s is predicted for seismic stations in the plain. Analysis of corresponding real data indicates a possible detection of teleseismic prograde particle motion, but the weak teleseismic earthquake signals are mixed with ambient noise signals at the predicted, short, transition periods. Detection of the period of polarity inversion from the joint analysis of earthquake and ambient noise ellipticity measurements may provide further, stringent, constraints on the structure of sedimentary basins.

Description: Kawakatsu et al. and Kawakatsu introduced a new fifth parameter, , to describe transverse isotropy (TI). Considering that characterizes the incidence angle dependence of body wave phase velocities for TI models, its relevance for body wave seismology is obvious. Here, we derive expressions for partial derivatives (sensitivity kernels) of surface wave phase velocity and normal mode eigenfrequency for the new set of five parameters. The partial derivative for is about twice as large as that for the conventional , indicating that should be more readily resolved. While partial derivatives for S velocities are not so changed, those for P velocities are significantly modified; the sensitivity for anisotropic P velocities is greatly reduced. In contrary to the suggestion by Dziewonski & Anderson and Anderson & Dziewonski, there is not much control on the anisotropic P velocities. On the other hand, the significance of for long-period seismology has become clear.

Description: With the deployment of extensive seismic arrays, systematic and efficient parameter and uncertainty estimation is of increasing importance and can provide reliable, regional models for crustal and upper-mantle structure. We present an efficient Bayesian method for the joint inversion of surface-wave dispersion and receiver-function data that combines trans-dimensional (trans-D) model selection in an optimization phase with subsequent rigorous parameter uncertainty estimation. Parameter and uncertainty estimation depend strongly on the chosen parametrization such that meaningful regional comparison requires quantitative model selection that can be carried out efficiently at several sites. While significant progress has been made for model selection (e.g. trans-D inference) at individual sites, the lack of efficiency can prohibit application to large data volumes or cause questionable results due to lack of convergence. Studies that address large numbers of data sets have mostly ignored model selection in favour of more efficient/simple estimation techniques (i.e. focusing on uncertainty estimation but employing ad-hoc model choices). Our approach consists of a two-phase inversion that combines trans-D optimization to select the most probable parametrization with subsequent Bayesian sampling for uncertainty estimation given that parametrization. The trans-D optimization is implemented here by replacing the likelihood function with the Bayesian information criterion (BIC). The BIC provides constraints on model complexity that facilitate the search for an optimal parametrization. Parallel tempering (PT) is applied as an optimization algorithm. After optimization, the optimal model choice is identified by the minimum BIC value from all PT chains. Uncertainty estimation is then carried out in fixed dimension. Data errors are estimated as part of the inference problem by a combination of empirical and hierarchical estimation. Data covariance matrices are estimated from data residuals (the difference between prediction and observation) and periodically updated. In addition, a scaling factor for the covariance matrix magnitude is estimated as part of the inversion. The inversion is applied to both simulated and observed data that consist of phase- and group-velocity dispersion curves (Rayleigh wave), and receiver functions. The simulation results show that model complexity and important features are well estimated by the fixed dimensional posterior probability density. Observed data for stations in different tectonic regions of the southern Korean Peninsula are considered. The results are consistent with published results, but important features are better constrained than in previous regularized inversions and are more consistent across the stations. For example, resolution of crustal and Moho interfaces, and absolute values and gradients of velocities in lower crust and upper mantle are better constrained.

Description: NW Iran is a region of active deformation in the Eurasia–Arabia collision zone. This high strain field has caused intensive faulting accompanied by several major ( M 〉 6.5) earthquakes as it is evident from historical records. Whereas seismic data (i.e. instrumental and historical catalogues) are either short, or inaccurate and inhomogeneous, physics-based long-term simulations are beneficial to better assess seismic hazard. In this study, a deterministic seismicity model, which consists of major active faults, is first constructed, and used to generate a synthetic catalogue of large-magnitude ( M 〉 5.5) earthquakes. The frequency–magnitude distribution of the synthetic earthquake catalogue, which is based on the physical characteristic and slip rate of the mapped faults, is consistent with the empirical distribution evaluated using record of instrumental and historical events. The obtained results are also in accordance with palaeoseismic studies and other independent kinematic deformation models of the Iranian Plateau. Using the synthetic catalogue, characteristic magnitude for all 16 active faults in the study area is determined. Magnitude and epicentre of these earthquakes are comparable with the historical records. Large earthquake recurrence times and their variations are evaluated, either for an individual fault or for the region as a whole. Goodness-of-fitness tests revealed that recurrence times can be well described by the Weibull distribution. Time-dependent conditional probabilities for large earthquakes in the study area are also estimated for different time intervals. The resulting synthetic catalogue can be utilized as a useful data set for hazard and risk assessment instead of short, incomplete and inhomogeneous available catalogues.

Description: On 2012 April 11, a great strike-slip earthquake (moment magnitude of M w 8.6) occurred off the west coast of northern Sumatra area followed by an M w 8.2 aftershock 2 hr later. Different geophysical data and methods have been used to investigate the mechanism, faulting, seismic radiation and slip propagation of this event, but frequency-dependent features of its rupture process have not been discussed much. In this study, we use a compressive sensing method based on sparsity inversion in the frequency domain to study the frequency-dependent seismic radiation and rupture process of this event. Our results indicate a very complex rupture process concerning at least three different rupture stages on multiple subfaults with nearly conjugate geometries. The main shock has triggered seismicity on a series of ridge-perpendicular or ridge-parallel conjugate strike-slip faults around the Nighty East Ridge. Obvious frequency-dependent rupture process has been presented and discussed. Combining results from slip inversion based on the finite-fault model, we observe that in the beginning stage of the rupture lower frequency radiation appears to originate from the areas with large slip, while the high-frequency radiation is located at the boundary of large-slip region or rupture front. Some radiation probably originates from the repeating slip on the main faults or triggered events on some nearby faults in the rupture area. The complex frequency-dependent seismic radiation patterns observed in this study provide important information for future investigation of rupture physics of this complex strike-slip event.

Description: Following the M w 6.3 L'Aquila Earthquake of 2009 April 6, the Colle di Roio village, central Italy, suffered severe building damages. The village is located on top of an elongated carbonate ridge characterized by a complex subsurface structure, a condition prone to seismic amplification due to topographic and stratigraphic effects. We address the role of the subsurface structure and topography in the ground-motion amplification observed at the ridge top. To characterize the subsurface structure of the ridge we performed geological investigations and ambient vibration measurements in single-station as well as 2-D-array configuration. Geological investigations pointed out that the ridge top is characterized by the presence of fractured rock material as a consequence of its anticlinal fold structure. Horizontal-to-vertical spectral ratio (HVSR) processing of ambient vibration records showed a broad peak in the HVSR functions in the frequency range 4–6 Hz and 2-D-array data demonstrated that locally the subsurface structure at the ridge top cannot be considered homogeneous. In summer 2009, we further installed one accelerometric station on the ridge top to experimentally evaluate the site amplification. By means of HVSR analysis of a sample of 18 weak-motion records (H/V), we found that ground-motion amplification occurs in a narrow frequency range centred around 4 Hz with mean ratio amplitude of 6. We also analysed the dependence of seismic amplification on the azimuth by calculating H/V ratios for horizontal components rotated into a range of azimuths. This analysis showed that the higher level of horizontal amplification occurs in the direction perpendicular to the ridge trending direction. With the aim of evaluating the contribution of the topography and the local subsurface structure on the observed seismic amplification, we performed 2-D finite-difference modelling of wave propagation through the ridge, adopting both homogeneous and heterogeneous models. We were able to reproduce the frequency-dependent amplification only by using an heterogeneous 2-D model with a low-velocity layer (shear wave velocity of 400 m s –1 and maximum thickness of 20 m) at the ridge top. From the comparison between observations and numerical results, we suggest that the amplification at the Colle di Roio ridge top is due to the interplay between local shear wave velocity structure and topography, rather than to topography alone.

Description: The interaction between seismic waves and slopes is an important topic to provide reliable scenarios for earthquake-(re)triggered landslides. The physical properties of seismic waves as well as slope topography and geology can significantly modify the local seismic response, influencing landslide triggering. A novel approach is here applied to two case studies in Andalusia (southern Spain) for computing the expected earthquake-induced displacements of existing landslide masses. Towards this aim, dynamic stress–strain numerical modelling was carried out using a selection of seismic signals characterized by different spectral content and energy. In situ geophysical measurements, consisting of noise records and temporary seismometric arrays, were carried out to control the numerical outputs in terms of local seismic response. The results consist of relationships between the characteristic period, T m , of the seismic signals and the characteristic periods of the landslide masses, related to the thickness ( T s ) and length ( T l ), respectively. These relationships show that the larger the horizontal dimension (i.e. length of landslide) of a landslide is, the more effective the contribution (to the resulting coseismic displacement) of the long-period seismic waves is, as the maximum displacements are expected for a low T m at each energy level of the input. On the other hand, when the local seismic response mainly depends on stratigraphy (i.e. landslide thickness), the maximum expected displacements occur close to the resonance period of the landslide, except for high-energy seismic inputs.

Description: The aim of this paper is to assess the accuracy of the long-period (50–250 s) surface waves extracted from cross-correlation functions (CCF) of ambient noise. First, we compare waveforms of Empirical Green's functions (EGF) converted from CCF with their synthetics, and also compare seismograms from a ground truth earthquake with their synthetics, through numerical simulations using a common 3-D model. We then quantify the accuracy of EGFs by comparing two sets of time-shifts between the observed waveforms and the synthetics: one set for the ground truth earthquake and the other set for EGFs. Second, we compare Rayleigh wave phase velocity dispersion measurements from ambient noise and those from earthquake data in both global and regional studies. Through these comparisons, we conclude that both the dispersion curves and waveforms from noise data are consistent with their counterparts from earthquake data in the long-period band. The long-period surface waves from ambient noise are as accurate as those from earthquake data, and can be included in both global and regional ambient noise tomography and provide complementary data to constrain the lithospheric and asthenospheric structures.

Description: We investigate the cause of anomalously long duration (〉100 s), long-period ( T  = 13–14 s) phase packets following the Rayleigh wave, which are observed widely in eastern Japan during outer-rise earthquakes. As this phase is not seen in records of interplate earthquakes, the late long-period phase is assumed to develop as a result of propagation across the Japan Trench. Using 3-D finite-difference simulations of seismic wave propagation with a detailed structural model that includes seafloor topography and water layers, we demonstrate that the late phase could be generated by a fundamental-mode oceanic Rayleigh wave propagating in the deep Pacific (at a depth of 6 km), which converts to a fundamental-mode continental Rayleigh wave as it crosses the Japan Trench. Our models suggest that the conversion is caused by the change in bathymetry. Based on analysis of surface wave dispersion, we confirm that the first arrival of the Rayleigh wave on land can develop from the first higher mode oceanic Rayleigh wave. The successive arrivals of the two types of Rayleigh waves on land produce ground motions with very long durations after outer-rise earthquakes.

Description: A new method for finding solutions to ordinary differential equation boundary value problems is introduced, in which Sobolev gradient steepest descent is used to determine eigenfunctions and eigenvalues simultaneously in an iterative scheme. The technique is then applied to the 1-D Love wave problem. The algorithm has several advantages when computing dispersion curves. It avoids the problem of mode skipping, and can handle arbitrary Earth structure profiles in depth. For a given frequency range, computation times scale approximately as the square root of the number of frequencies, and the computation of dispersion curves can be implemented in a fully parallel manner over the modes involved. The steepest descent solutions are within a fraction of a per cent of the analytic solutions for the first 25 modes for a two-layer model. Since all corresponding eigenfunctions are computed along with the dispersion curves, the impact on group and phase velocity of the displacement behaviour with depth is thoroughly examined. The dispersion curves are used to compute synthetic Love wave seismograms that include many higher order modes. An example includes addition of attenuation to a model with a low-velocity zone, with values as low as Q = 20. Finally, a confirming comparison is made with a layer matrix method on the upper 700 km of a whole Earth model.

Description: The aspherical structure of the Earth is described in terms of lateral heterogeneity and anisotropy of the P - and S -wave velocities, density heterogeneity, ellipticity and rotation of the Earth and undulation of the discontinuity interfaces of the seismic wave velocities. Its structure significantly influences the normal mode spectra of the Earth's free oscillation in the form of cross-coupling between toroidal and spheroidal multiplets and self-coupling between the singlets forming them. Thus, the aspherical structure must be conversely estimated from the free oscillation spectra influenced by the cross-coupling and self-coupling. In the present study, we improve a spectral fitting inversion algorithm which was developed in a previous study to retrieve the global structures of the isotropic and anisotropic velocities of the P and S waves from the free oscillation spectra. The main improvement is that the geographical distribution of the intensity of the S -wave azimuthal anisotropy is represented by a nonlinear combination of structure coefficients for the anisotropic velocity structure, whereas in the previous study it was expanded into a generalized spherical harmonic series. Consequently, the improved inversion algorithm reduces the number of unknown parameters that must be determined compared to the previous inversion algorithm and employs a one-step inversion method by which the structure coefficients for the isotropic and anisotropic velocities are directly estimated from the fee oscillation spectra. The applicability of the improved inversion is examined by several numerical experiments using synthetic spectral data, which are produced by supposing a variety of isotropic and anisotropic velocity structures, earthquake source parameters and station-event pairs. Furthermore, the robustness of the inversion algorithm is investigated with respect to the back-ground noise contaminating the spectral data as well as truncating the series expansions by finite terms to represent the three-dimensional velocity structures. As a result, it is shown that the improved inversion can estimate not only the isotropic and anisotropic velocity structures but also the depth extent of the anisotropic regions in the Earth. In particular, the cross-coupling modes are essential to correctly estimate the isotropic and anisotropic velocity structures from the normal mode spectra. In addition, we argue that the effect of the seismic anisotropy is not negligible when estimating only the isotropic velocity structure from the spheroidal mode spectra.

Description: Fault zone head waves (FZHWs) are observed along major strike-slip faults and can provide high-resolution imaging of fault interface properties at seismogenic depth. In this paper, we present a new method to automatically detect FZHWs and pick direct P waves secondary arrivals (DWSAs). The algorithm identifies FZHWs by computing the amplitude ratios between the potential FZHWs and DSWAs. The polarities, polarizations and characteristic periods of FZHWs and DSWAs are then used to refine the picks or evaluate the pick quality. We apply the method to the Parkfield section of the San Andreas Fault where FZHWs have been identified before by manual picks. We compare results from automatically and manually picked arrivals and find general agreement between them. The obtained velocity contrast at Parkfield is generally 5–10 per cent near Middle Mountain while it decreases below 5 per cent near Gold Hill. We also find many FZHWs recorded by the stations within 1 km of the background seismicity (i.e. the Southwest Fracture Zone) that have not been reported before. These FZHWs could be generated within a relatively wide low velocity zone sandwiched between the fast Salinian block on the southwest side and the slow Franciscan Mélange on the northeast side. Station FROB on the southwest (fast) side also recorded a small portion of weak precursory signals before sharp P waves. However, the polarities of weak signals are consistent with the right-lateral strike-slip mechanisms, suggesting that they are unlikely genuine FZHW signals.

Description: Significant effort has been devoted over the last two decades to the development of various seismic velocity models for the region of southern California, United States. These models are mostly used in forward wave propagation simulation studies, but also as base models for tomographic and source inversions. Two of these models, the community velocity models CVM-S and CVM-H, are among the most commonly used for this region. This includes two alternative variations to the original models, the recently released CVM-S4.26 which incorporates results from a sequence of tomographic inversions into CVM-S, and the user-controlled option of CVM-H to replace the near-surface profiles with a V S 30 -based geotechnical model. Although either one of these models is regarded as acceptable by the modeling community, it is known that they have differences in their representation of the crustal structure and sedimentary deposits in the region, and thus can lead to different results in forward and inverse problems. In this paper, we evaluate the accuracy of these models when used to predict the ground motion in the greater Los Angeles region by means of an assessment of a collection of simulations of recent events. In total, we consider 30 moderate-magnitude earthquakes (3.5 〈 M w 〈 5.5) between 1998 and 2014, and compare synthetics with data recorded by seismic networks during these events. The simulations are done using a finite-element parallel code, with numerical models that satisfy a maximum frequency of 1 Hz and a minimum shear wave velocity of 200 m s –1 . The comparisons between data and synthetics are ranked quantitatively by means of a goodness-of-fit (GOF) criteria. We analyse the regional distribution of the GOF results for all events and all models, and draw conclusions from the results and how these correlate to the models. We find that, in light of our comparisons, the model CVM-S4.26 consistently yields better results.

Description: A multidomain indirect boundary element method (IBEM) is proposed to study the wave scattering of plane SH waves by complex local site in a layered half-space. The new method, using both the full-space and layered half-space Green's functions as its fundamental solutions can also be regarded as a coupled method of the full-space IBEM and half-space IBEM. First, the whole model is decomposed into independent closed regions and an opened layered half-space region with all of the irregular interfaces; then, fictitious uniformly distributed loads are applied separately on the boundaries of each region, and scattered fields of the closed regions and the opened layered half-space region are constructed by calculating the full-space and layered half-space Green's functions, respectively; finally, all of the regions are assembled to establish the linear algebraic system that arises from discretization. The densities of the distributed loads are determined directly by solving the algebraic system. The accuracy and capability of the new approach are verified extensively by comparing its results with those of published approaches for a class of hills, valleys and embedded inclusions. And the capability of the new method is further displayed when it is used to investigate a hill-triple layered valley-hill coupled topography in a multilayered half-space. All of the numerical calculations presented in this paper demonstrate that the new method is very suitable for solving multidomain coupled multilayered wave scattering problems with the merits of high accuracy and representing the scattered fields in different kinds of regions more reasonably and flexibly.

Description: Seismic interferometry and beam-forming techniques were applied to the ambient noise recorded during January 2014 at the ‘13 BB star’ array composed of thirteen seismic stations located in northern Poland. The circular and symmetric geometry of the array allowed the evaluation of the azimuths of noise sources and the velocities of recovered surface waves with a good reliability. After having pre-processed the raw records of the ambient noise in time- and frequency-domain, we studied the associated power spectral density to identify the frequency bands suitable for the recovery of the surface waves. Then the cross-correlation was performed between all the station pairs of the array to retrieve the Green's function, from which the velocity range of the surface waves can be determined. Making use of that analysis, the direction of the noise wavefield was linked to the maximum amplitude of the beam-power, estimated by the mixing in the frequency-domain of all the corresponding noise records. The results were related day by day to the mean wind velocity around Europe at 10 m above ground level obtained from global surveys carried out during the same month. Significant correlation between the direction of maximum beam-power associated to the ambient noise recorded at ‘13 BB star’ and the average wind velocity was found.

Description: We introduce a method for mapping the noise-source distribution of microseisms which uses information from the full length of covariograms (cross-correlations). We derive a forward calculation based on the plane-wave assumption in 2-D, to formulate an iterative, linearized inversion of covariogram envelopes in the time domain. The forward calculation involves bandpass filtering of the covariograms. The inversion exploits the well-known feature of noise cross-correlation, that is, an anomaly in the noise field that is oblique to the interstation direction appears as cross-correlation amplitude at a smaller time lag than the in-line, surface wave arrival. Therefore, the inversion extracts more information from the covariograms than that contained at the expected surface wave arrival, and this allows us to work with few stations to find the propagation directions of incoming energy. The inversion is naturally applied to data that retain physical units that are not amplitude normalized in any way. By dividing a network into groups of stations, we can constrain the source location by triangulation. We demonstrate results of the method with synthetic data and one year (2012) of data from the Swedish National Seismic Network and also look at the seasonal variation of source distribution around Scandinavia. After preprocessing and cross-correlation, the stations are divided into five groups of 9–12 stations. We invert the envelopes of each group in eight period ranges between 2 and 25 s. Results show that the noise sources at short periods (less than 12 s) lie predominantly in the North Atlantic Ocean and the Barents Sea, and at longer periods the energy appears to have a broader distribution. The strongly anisotropic source distribution in this area is estimated to cause significant biases of velocity measurements compared to the level of heterogeneity in the region. The amplitude of the primary microseisms varies little over the year, but secondary microseisms are much weaker in summer than in winter. Furthermore, the peak period of the secondary microseisms shifts from 5–6 s in winter to 4–5 s during the summer.

Description: We address the problem of obtaining a reliable seismic image without prior knowledge of the source wavelet, especially from data that contain strong surface-related multiples. Conventional reverse-time migration requires prior knowledge of the source wavelet, which is either technically or computationally challenging to accurately determine; inaccurate estimates of the source wavelet can result in seriously degraded reverse-time migrated images, and therefore wrong geological interpretations. To solve this problem, we present a ‘wavelet-free’ imaging procedure that simultaneously inverts for the source wavelet and the seismic image, by tightly integrating source estimation into a fast least-squares imaging framework, namely compressive imaging , given a reasonably accurate background velocity model. However, this joint inversion problem is difficult to solve as it is plagued with local minima and the ambiguity with respect to amplitude scalings because of the multiplicative, and therefore nonlinear, appearance of the source wavelet in the otherwise linear formalism. We have found a way to solve this nonlinear joint-inversion problem using a technique called variable projection, and a way to overcome the scaling ambiguity by including surface-related multiples in our imaging procedure following recent developments in surface-related multiple prediction by sparse inversion. As a result, we obtain without prior knowledge of the source wavelet high-resolution seismic images, comparable in quality to images obtained assuming the true source wavelet is known. By leveraging the computationally efficient compressive-imaging methodology, these results are obtained at affordable computational costs compared with conventional processing work flows that include surface-related multiple removal and reverse-time migration.

Description: To investigate the descend of subducted slabs we search for and analyse seismic arrivals that reflected off the surface of the slab. In order to distinguish between such arrivals and other seismic phases, we search for waves that reach a seismic array with a backazimuth deviating from the theoretical backazimuth of the earthquake. Source–receiver combinations are chosen in a way that their great circle paths do not intersect the slab region, hence the direct arrivals can serve as reference. We focus on the North and Northwest Pacific region by using earthquakes from Japan, the Philippines and the Hindu Kush area recorded at North American networks (e.g. USArray, Alaska and Canada). Using seismic array techniques for analysing the data and record information on slowness, backazimuth and traveltime of the observed out-of-plane arrivals we use these measurements to trace the wave back through a 1-D velocity model to its scattering/reflection location. We find a number of out-of-plane reflections. Assuming only single scattering, most out-of-plane signals have to travel as P -to- P phases and only a few as S -to- P phases, due to the length of the seismograms we processed. The located reflection points present a view of the 3-D structures within the mantle. In the upper mantle and the transition zone they correlate well with the edges of fast velocity regions in tomographic images. We also find reflection points in the mid- and lower mantle and their locations generally agree with fast velocities mapped by seismic tomography models suggesting that in the subduction regions we map, slabs enter the lower mantle. To validate our approach, we calculate and process synthetic seismograms for 3-D wave field propagation through a model containing a slab-like heterogeneity. We show, that depending on the source–receiver geometry relative to the reflection plane, it is indeed possible to observe and back-trace out-of-plane signals.

Description: We present a theoretical analysis of the error of synthetic seismograms computed by higher-order finite-element methods (ho-FEMs). We show the existence of a previously unrecognized type of error due to degenerate coupling between waves with the same frequency but different wavenumbers. These results are confirmed by simple numerical experiments using the spectral element method as an example of ho-FEMs. Errors of the type found by this study may occur generally in applications of ho-FEMs.

Description: We present a new methodology for detection and space–time location of seismic sources based on multiscale, frequency-selective coherence of the wave field recorded by dense large-scale seismic networks and local antennas. The method is designed to enhance coherence of the signal statistical features across the array of sensors and consists of three steps: signal processing, space–time imaging, and detection and location. The first step provides, for each station, a simplified representation of seismic signal by extracting multiscale non-stationary statistical characteristics, through multiband higher-order statistics or envelopes. This signal processing scheme is designed to account for a priori unknown transients, potentially associated with a variety of sources (e.g. earthquakes, tremors), and to prepare data for a better performance in posterior steps. Following space–time imaging is carried through 3-D spatial mapping and summation of station-pair time-delay estimate functions. This step produces time-series of 3-D spatial images representing the likelihood that each pixel makes part of a source. Detection and location is performed in the final step by extracting the local maxima from the 3-D spatial images. We demonstrate the efficiency of the method in detecting and locating seismic sources associated with low signal-to-noise ratio on an example of the aftershock earthquake records from local stations of International Maule Aftershock Deployment in Central Chile. The performance and potential of the method to detect, locate and characterize the energy release associated with possibly mixed seismic radiation from earthquakes and low-frequency tectonic tremors is further tested on continuous data from southwestern Japan.

Description: Information about seismic anisotropy is embedded in the variation of the amplitude of the Ps pulses as a function of the azimuth, on both the Radial and the Transverse components of teleseismic receiver functions (RF). We develop a semi-automatic method to constrain the presence and the depth of anisotropic layers beneath a single seismic broad-band station. An algorithm is specifically designed to avoid trial and error methods and subjective crustal parametrizations in RF inversions, providing a suitable tool for large-size data set analysis. The algorithm couples together information extracted from a 1-D V S profile and from a harmonic decomposition analysis of the RF data set. This information is used to determine the number of anisotropic layers and their approximate position at depth, which, in turn, can be used to, for example, narrow the search boundaries for layer thickness and S -wave velocity in a subsequent parameter space search. Here, the output of the algorithm is used to invert an RF data set by means of the Neighbourhood Algorithm (NA). To test our methodology, we apply the algorithm to both synthetic and observed data. We make use of synthetic RF with correlated Gaussian noise to investigate the resolution power for multiple and thin (1–3 km) anisotropic layers in the crust. The algorithm successfully identifies the number and position of anisotropic layers at depth prior the NA inversion step. In the NA inversion, strength of anisotropy and orientation of the symmetry axis are correctly retrieved. Then, the method is applied to field measurement from station BUDO in the Tibetan Plateau. Two consecutive layers of anisotropy are automatically identified with our method in the first 25–30 km of the crust. The data are then inverted with the retrieved parametrization. The direction of the anisotropic axis in the uppermost layer correlates well with the orientation of the major planar structure in the area. The deeper anisotropic layer is associated with an older phase of crustal deformation. Our results are compared with previous anisotropic RF studies at the same station, showing strong similarities.

Description: We report global phase-velocity maps of fundamental-mode Rayleigh waves at periods between 25 and 100 s based on finite-frequency tomography. Rayleigh wave dispersion measurements are made using a multitaper technique for both minor-arc and major-arc wave trains. The global phase-velocity maps confirm many features associated with surface tectonics including the ocean-continent dichotomy and the signature of lithospheric cooling in oceanic plates. In addition, the high-resolution phase-velocity maps reveal a major change in the distribution of small-scale anomalies in the Pacific at different wave periods. We calculate the global average of Rayleigh wave phase velocity in major tectonic regions and show that large discrepancies exist between our model and global crustal and mantle models: (1) In oceanic regions, short-period (〈~40 s) Rayleigh waves are faster than calculations based on models CRUST2.0 and S40RTS. The discrepancies could be explained by a thinner crust or faster wave speeds in the crust or upper mantle. The implementation of model CRUST1.0 significantly improves the agreement, with phase-velocity discrepancies less than 0.5 per cent on average. (2) In Archean cratons, Rayleigh wave phase velocities in our model are faster than calculations based on model S40RTS at periods longer than ~40 s; and the global average in orogenic belts is ~1–2 per cent slower than CRUST1.0 at periods shorter than ~50 s.

Description: Travelling-wave Green tensor has been widely used in calculations of synthetic seismograms and finite-frequency sensitivities of surface waves. The classic travelling-wave decomposition is based on a far-field approximation and may not be valid when applied to construct sensitivity kernels in regions close to the receiver. In this paper, we calculate synthetic seismograms and finite-frequency sensitivity kernels of Rayleigh waves based on travelling-wave representation of Green tensor that fully accounts for near-field effects. We show that far-field approximation is adequate for synthetic seismograms when the source–receiver epicentral distance is greater than the dominant wavelength. Errors in Rayleigh-wave sensitivity kernels introduced by far-field approximation are in general negligible for single-station measurements except for in a small region around the station, and the errors are more significant in sensitivity kernels for interstation measurements. In addition, interstation measurements are strongly sensitive to structures outside the region between the two stations, even for two stations along the same great circle path from the seismic source.

Description: The surface wave full ray theory (FRT) is an efficient tool to calculate synthetic waveforms of surface waves. It combines the concept of local modes with exact ray tracing as a function of frequency, providing a more complete description of surface wave propagation than the widely used great circle approximation (GCA). The purpose of this study is to evaluate the ability of the FRT approach to model teleseismic long-period surface waveforms ( T ~ 45–150 s) in the context of current 3-D Earth models to empirically assess its validity domain and its scope for future studies in seismic tomography. To achieve this goal, we compute vertical and horizontal component fundamental mode synthetic Rayleigh waveforms using the FRT, which are compared with calculations using the highly accurate spectral element method. We use 13 global earth models including 3-D crustal and mantle structure, which are derived by successively varying the strength and lengthscale of heterogeneity in current tomographic models. For completeness, GCA waveforms are also compared with the spectral element method. We find that the FRT accurately predicts the phase and amplitude of long-period Rayleigh waves ( T ~ 45–150 s) for almost all the models considered, with errors in the modelling of the phase (amplitude) of Rayleigh waves being smaller than 5 per cent (10 per cent) in most cases. The largest errors in phase and amplitude are observed for T ~ 45 s and for the three roughest earth models considered that exhibit shear wave anomalies of up to ~20 per cent, which is much larger than in current global tomographic models. In addition, we find that overall the GCA does not predict Rayleigh wave amplitudes well, except for the longest wave periods ( T ~ 150 s) and the smoothest models considered. Although the GCA accurately predicts Rayleigh wave phase for current earth models such as S20RTS and S40RTS, FRT's phase errors are smaller, notably for the shortest wave periods considered ( T ~ 45 s and T ~ 60 s). This suggests that the FRT approach is a useful means to build the next generation of elastic and anelastic surface wave tomography models. Finally, we observe a clear correlation between the FRT amplitude and phase errors and the roughness of the models. This allows us to quantify the limits of validity of the FRT in terms of model roughness thresholds, which can serve as useful guides in future seismic tomographic studies.

Description: Iranian Plateau reflects several different tectonic styles of collision, and large-scale strike-slip faults. We calculate a high-resolution 2-D M L shear velocity map for the Iranian Plateau to detect lateral crustal thickness changes associated with different tectonic boundaries. The M L velocity is very sensitive to strong lateral variations of crustal thickness and varies between the velocity of Lg and Sn phases. Our data set consists of 65 795 M L amplitude velocity measurements from 2531 precisely relocated events recorded by Iranian networks in the period 1996–2014. Using a constrained least-squares inversion scheme, we inverted the M L velocities for a 2-D shear velocity map of Iran. Our results show that the Zagros and South Caspian Basin (SCB) have shear wave velocities close to the Sn phase, and are thus Lg-blocking regions. High velocities in the High Zagros and the Simply Folded Belt imply significant crustal undulations within these zones. We note that in the central and south Zagros, the velocity border between the Zagros and central Iran is not coincident with the Zagros suture line that marks underthrusting of the Arabian plate beneath central Iran. The low plains of Gilan and Gorgan to the south of the Caspian Sea show high shear velocities similar to the SCB, implying that they are either underlain by an oceanic type crust or a transitional crust with a strong lateral crustal thickness gradient. The Lut block is an Lg-passing block implying that it is not surrounded by any sudden crustal thickness changes along its borders with central Iran. In the Alborz, NW Iran, Kopeh-Dagh, Binalud and most of the central Iran, low shear velocity near the Lg velocity is attributed to smooth or minor Moho undulations within these regions.

Description: We investigate the problem of finite-difference approximations of the velocity–stress formulation of the equation of motion and constitutive law on the staggered grid (SG) and collocated grid (CG). For approximating the first spatial and temporal derivatives, we use three approaches: Taylor expansion (TE), dispersion-relation preserving (DRP), and combined TE-DRP. The TE and DRP approaches represent two fundamental extremes. We derive useful formulae for DRP and TE-DRP approximations. We compare accuracy of the numerical wavenumbers and numerical frequencies of the basic TE, DRP and TE-DRP approximations. Based on the developed approximations, we construct and numerically investigate 14 basic TE, DRP and TE-DRP finite-difference schemes on SG and CG. We find that (1) the TE second-order in time, TE fourth-order in space, 2-point in time, 4-point in space SG scheme (that is the standard (2,4) VS SG scheme, say TE-2-4-2-4-SG) is the best scheme (of the 14 investigated) for large fractions of the maximum possible time step, or, in other words, in a homogeneous medium; (2) the TE second-order in time, combined TE-DRP second-order in space, 2-point in time, 4-point in space SG scheme (say TE-DRP-2-2-2-4-SG) is the best scheme for small fractions of the maximum possible time step, or, in other words, in models with large velocity contrasts if uniform spatial grid spacing and time step are used. The practical conclusion is that in computer codes based on standard TE-2-4-2-4-SG, it is enough to redefine the values of the approximation coefficients by those of TE-DRP-2-2-2-4-SG for increasing accuracy of modelling in models with large velocity contrast between rock and sediments.

Description: The multiple-taper correlation (MTC) algorithm for the estimation of teleseismic receiver functions (RFs) has desirable statistical properties. This paper presents several adaptations to the MTC algorithm that exploit its frequency-domain uncertainty estimates to generate stable RFs that include moveout corrections for deeper interfaces. Narrow-band frequency averaging implicit in spectral cross-correlation restricts the MTC-based RF estimates to resolve Ps converted phases only at short delay times, appropriate to the upper 100 km of Earth's lithosphere. The Ps conversions from deeper interfaces can be reconstructed by the MTC algorithm in two ways. Event cross-correlation computes a cross-correlation of single-taper spectrum estimates for a cluster of events rather than for a set of eigenspectrum estimates of a single P coda. To extend the reach of the algorithm, pre-stack moveout corrections in the frequency domain preserves the formal uncertainties of the RF estimates, which are used to weight RF stacks. Moving-window migration retains the multiple-taper approach, but cross-correlates the P -polarized motion with time-delayed SH and SV motion to focus on a Ps phase of interest. The frequency-domain uncertainties of bin-averaged RFs do not translate directly into the time domain. A jackknife over data records in each bin stack offers uncertainty estimates in the time domain while preserving uncertainty weighting in the frequency-domain RF stack.

Description: Amplification can occur in a graben as a result of strong earthquake-induced ground motion. Thus, in seismic hazard and seismic site response studies, it is of the utmost importance to determine the geometry of the bedrock depth. The main objectives of this study were to determine the bedrock depth and map the depth-to-bedrock ratio for use in land use planning in regard to the mitigation of earthquake hazards in the Eskisehir Basin. The fundamental resonance frequencies ( f r ) of 318 investigation sites in the Eskisehir Basin were determined through case studies, and the 2-D S -wave velocity structure down to the bedrock depth was explored. Single-station microtremor data were collected from the 318 sites, as well as microtremor array data from nine sites, seismic reflection data from six sites, deep-drilling log data from three sites and shallow drilling log data from ten sites in the Eskisehir Graben. The fundamental resonance frequencies of the Eskisehir Basin sites were obtained from the microtremor data using the horizontal-to vertical (H/V) spectral ratio (HVSR) method. The phase velocities of the Rayleigh waves were estimated from the microtremor data using the spatial autocorrelation (SPAC) method. The fundamental resonance frequency range at the deepest point of the Eskisehir Basin was found to be 0.23–0.35 Hz. Based on the microtremor array measurements and the 2-D S -wave velocity profiles obtained using the SPAC method, a bedrock level with an average velocity of 1300 m s –1 was accepted as the bedrock depth limit in the region. The log data from a deep borehole and a seismic reflection cross-section of the basement rocks of the Eskisehir Basin were obtained and permitted a comparison of bedrock levels. Tests carried out using a multichannel walk-away technique permitted a seismic reflection cross-section to be obtained up to a depth of 1500–2000 m using an explosive energy source. The relationship between the fundamental resonance frequency in the Eskisehir Basin and the results of deep drilling, shallow drilling, shear wave velocity measurement and sedimentary cover depth measurement obtained from the seismic reflection section was expressed in the form of a nonlinear regression equation. An empirical relationship between f r , the thickness of sediments and the bedrock depth is suggested for use in future microzonation studies of sites in the region. The results revealed a maximum basin depth of 1000 m, located in the northeast of the Eskisehir Basin, and the SPAC and HVSR results indicated that within the study area the basin is characterized by a thin local sedimentary cover with low shear wave velocity overlying stiff materials, resulting in a sharp velocity contrast. The thicknesses of the old Quaternary and Tertiary fluvial sediments within the basin serve as the primary data sources in seismic hazard and seismic site response studies, and these results add to the body of available seismic hazard data contributing to a seismic microzonation of the Eskisehir Graben in advance of the severe earthquakes expected in the Anatolian Region.

Description: We document space-dependent clustering properties of earthquakes with m ≥ 4 in the 1975–2015 worldwide seismic catalogue of the Northern California Earthquake Data Center. Earthquake clusters are identified using a nearest-neighbour distance in time–space–magnitude domain. Multiple cluster characteristics are compared with the heat flow level and type of deformation defined by parameters of the strain rate tensor. The analysis suggests that the dominant type of seismicity clusters in a region depends strongly on the heat flow, while the deformation style and intensity play a secondary role. The results show that there are two dominant types of global clustering: burst-like clusters that represent brittle fracture in relatively cold lithosphere (e.g. shallow events in subduction zones) and swarm-like clusters that represent brittle–ductile deformation in relatively hot lithosphere (e.g. mid-oceanic ridges). The global results are consistent with theoretical expectations and previous analyses of earthquake clustering in southern California based on higher quality catalogues. The observed region-specific deviations from average universal description of seismicity provide important constraints on the physics governing earthquakes and can be used to improve local seismic hazard assessments.

Description: Surface wave studies in the 1960s provided the first indication that the upper mantle was radially anisotropic. Resolving the anisotropic structure is important because it may yield information on deformation and flow patterns in the upper mantle. The existing radially anisotropic models are in poor agreement. Rayleigh waves have been studied extensively and recent models show general agreement. Less work has focused on Love waves and the models that do exist are less well-constrained than are Rayleigh wave models, suggesting it is the Love wave models that are responsible for the poor agreement in the radially anisotropic structure of the upper mantle. We have adapted the waveform inversion procedure of Debayle & Ricard to extract propagation information for the fundamental mode and up to the fifth overtone from Love waveforms in the 50–250 s period range. We have tomographically inverted these results for a mantle horizontal shear wave-speed model (β h ( z )) to transition zone depths. We include azimuthal anisotropy (2 and 4 terms) in the tomography, but in this paper we discuss only the isotropic β h ( z ) structure. The data set is significantly larger, almost 500 000 Love waveforms, than previously published Love wave data sets and provides ~17 000 000 constraints on the upper-mantle β h ( z ) structure. Sensitivity and resolution tests show that the horizontal resolution of the model is on the order of 800–1000 km to transition zone depths. The high wave-speed roots beneath the oldest parts of the continents appear to extend deeper for β h ( z ) than for β v ( z ) as in previous β h ( z ) models, but the resolution tests indicate that at least parts of these features could be artefacts. The low wave speeds beneath the mid-ocean ridges fade by ~150 km depth except for the upper mantle beneath the East Pacific Rise which remains slow to ~250 km depth. The resolution tests suggest that the low wave speeds at deeper depths beneath the East Pacific Rise are not solely due to vertical smearing of shallow, low wave speeds. Four prominent, low wave-speed features occur at transition zone depths—one aligned along the East African Rift, one centred south of the Indian peninsula, one located south of New Zealand and one in the south Pacific Ocean coinciding with the location of the South Pacific Superswell. The low wave-speed features south of New Zealand and south of the Indian peninsula correspond spatially with the two largest negative geoid lows on Earth.

Description: We assess the potential of additional rotational ground motions to increase the resolution of the full seismic moment tensor and its centroid depth during waveform inversion. For this purpose, we set up a test case of a shallow, medium-sized strike-slip source. In two scenarios, one based on theoretical station distribution and the other based on real station distribution, we compare the results based on inversion of translational ground motion data only and based on both, translational and rotational ground motion data. The inversion is done with a Bayesian approach to overcome the drawbacks of deterministic approaches and provide a comprehensive quantification of uncertainties. Our results indicate that the resolution of the moment tensor can be increased drastically by incorporating rotational ground motion data. Especially, the usually problematic components M xz and M yz as well as all components containing spatial derivatives with depth benefit most. Also, the resolution of the centroid depth is much better.

Description: We present a single station method for the determination of Rayleigh wave ellipticity, or Rayleigh wave horizontal to vertical amplitude ratio (H/V) using Frequency Dependent Polarization Analysis (FDPA). This procedure uses singular value decomposition of 3-by-3 spectral covariance matrices over 1-hr time windows to determine properties of the ambient seismic noise field such as particle motion and dominant wave-type. In FPDA, if the noise is mostly dominated by a primary singular value and the phase difference is roughly 90° between the major horizontal axis and the vertical axis of the corresponding singular vector, we infer that Rayleigh waves are dominant and measure an H/V ratio for that hour and frequency bin. We perform this analysis for all available data from the Earthscope Transportable Array between 2004 and 2014. We compare the observed Rayleigh wave H/V ratios with those previously measured by multicomponent, multistation noise cross-correlation (NCC), as well as classical noise spectrum H/V ratio analysis (NSHV). At 8 s the results from all three methods agree, suggesting that the ambient seismic noise field is Rayleigh wave dominated. Between 10 and 30 s, while the general pattern agrees well, the results from FDPA and NSHV are persistently slightly higher (~2 per cent) and significantly higher (〉20 per cent), respectively, than results from the array-based NCC. This is likely caused by contamination from other wave types (i.e. Love waves, body waves, and tilt noise) in the single station methods, but it could also reflect a small, persistent error in NCC. Additionally, we find that the single station method has difficulty retrieving robust Rayleigh wave H/V ratios within major sedimentary basins, such as the Williston Basin and Mississippi Embayment, where the noise field is likely dominated by reverberating Love waves and tilt noise.

Description: In this study, we use ambient seismic noise recorded at selected broad-band USArray Earthscope Transportable Array (TA) stations to obtain effective reflection seismograms using noise autocorrelations. In order to best retrieve the body-wave component of the Green's function beneath a station from ambient seismic noise, a number of processing steps are used, including temporal sign-bit normalization, spectral whitening and bandpass filtering. Hourly autocorrelations are stacked for different time periods including one day, one month and one year. On the final stack, different amplitude gain functions are applied, including automatic gain control (AGC), to equalize the correlation amplitudes. The robustness of the resulting ambient noise autocorrelations is first tested on a TA station in Nevada where we are able to identify arrivals similar to those found in an earlier study. We then investigated noise autocorrelations applied to several USArray TA stations in the central U.S., and the results were then compared with reflectivity synthetics for an average crustal model based on CRUST 1.0 where an AGC was used to enhance the later arrivals. Different stacking periods are also investigated in order to find stable correlation stacks.