ALBERT

Description: Elastic reverse time migration (RTM) can yield accurate subsurface information (e.g. PP and PS reflectivity) by imaging the multicomponent seismic data. However, the existing RTM methods are still insufficient to provide satisfactory results because of the finite recording aperture, limited bandwidth and imperfect illumination. Besides, the P - and S -wave separation and the polarity reversal correction are indispensable in conventional elastic RTM. Here, we propose an iterative elastic least-squares RTM (LSRTM) method, in which the imaging accuracy is improved gradually with iteration. We first use the Born approximation to formulate the elastic de-migration operator, and employ the Lagrange multiplier method to derive the adjoint equations and gradients with respect to reflectivity. Then, an efficient inversion workflow (only four forward computations needed in each iteration) is introduced to update the reflectivity. Synthetic and field data examples reveal that the proposed LSRTM method can obtain higher-quality images than the conventional elastic RTM. We also analyse the influence of model parametrizations and misfit functions in elastic LSRTM. We observe that Lamé parameters, velocity and impedance parametrizations have similar and plausible migration results when the structures of different models are correlated. For an uncorrelated subsurface model, velocity and impedance parametrizations produce fewer artefacts caused by parameter crosstalk than the Lamé coefficient parametrization. Correlation- and convolution-type misfit functions are effective when amplitude errors are involved and the source wavelet is unknown, respectively. Finally, we discuss the dependence of elastic LSRTM on migration velocities and its antinoise ability. Imaging results determine that the new elastic LSRTM method performs well as long as the low-frequency components of migration velocities are correct. The quality of images of elastic LSRTM degrades with increasing noise.

Description: Seismic noise measurements (ambient vibrations) have been increasingly used in rock slope stability assessment for both investigation and monitoring purposes. Recent studies made on gravitational hazard revealed significant spectral amplification at given frequencies and polarization of the wave-field in the direction of maximum rock slope displacement. Different properties (resonance frequencies, polarization and spectral ratio amplitudes) can be derived from the spectral analysis of the seismic noise to characterize unstable rock masses. The objective here is to identify the dynamic parameters that could be used to gain information on prone-to-fall rock columns’ geometry. To do so, the dynamic response of prone-to-fall columns to seismic noise has been studied on two different sites exhibiting cliff-like geometry. Dynamic parameters (main resonance frequency and spectral ratio amplitudes) that could characterize the column decoupling were extracted from seismic noise and their variations were studied taking into account the external environmental parameter fluctuations. Based on this analysis, a two-dimensional numerical model has been set up to assess the influence of the rear vertical fractures identified on both sites on the rock column motion response. Although a simple relation was found between spectral ratio amplitudes and the rock column slenderness, it turned out that the resonance frequency is more stable than the spectral ratio amplitudes to characterize this column decoupling, provided that the elastic properties of the column can be estimated. The study also revealed the effect of additional remote fractures on the dynamic parameters, which in turn could be used for detecting the presence of such discontinuities.

Description: In this study, we present a new synthesis of GPS velocities for tectonic deformation within the Tibetan Plateau and its surrounding areas, a combined data set of ~1854 GPS-derived horizontal velocity vectors. Assuming that crustal deformation is localized along major faults, a block modelling approach is employed to interpret the GPS velocity field. We construct a 30-element block model to describe present-day deformation in western China, with half of them located within the Tibetan Plateau, and the remainder located in its surrounding areas. We model the GPS velocities simultaneously for the effects of block rotations and elastic strain induced by the bounding faults. Our model yields a good fit to the GPS data with a mean residual of 1.08 mm a –1 compared to the mean uncertainty of 1.36 mm a –1 for each velocity component, indicating a good agreement between the predicted and observed velocities. The major strike-slip faults such as the Altyn Tagh, Xianshuihe, Kunlun and Haiyuan faults have relatively uniform slip rates in a range of 5–12 mm a –1 along most of their segments, and the estimated fault slip rates agree well with previous geologic and geodetic results. Blocks having significant residuals are located at the southern and southeastern Tibetan Plateau, suggesting complex tectonic settings and further refinement of accurate definition of block geometry in these regions.

Description: Surface-related multiples have been utilized in the reverse-time migration (RTM) procedures, and additional illumination for subsurface can be provided. Meanwhile, many cross-talks are generated from undesired interactions between forward- and backward-propagated seismic waves. In this paper, subsequent to analysing and categorizing these cross-talks, we propose RTM of first-order multiples to avoid most undesired interactions in RTM of all-order multiples, where only primaries are forward-propagated and crosscorrelated with the backward-propagated first-order multiples. With primaries and multiples separated during regular seismic data processing as the input data, first-order multiples can be obtained by a two-step scheme: (1) the dual-prediction of higher-order multiples; and (2) the adaptive subtraction of predicted higher-order multiples from all-order multiples within local offset-time windows. In numerical experiments, two synthetic and a marine field data sets are used, where different cross-talks generated by RTM of all-order multiples can be identified and the proposed RTM of first-order multiples can provide a very interpretable image with a few cross-talks.

Description: We investigated post-seismic velocity changes within the fault zone of the 2008 M 7.9 Wenchuan earthquake using coda wave data of repeating small earthquakes. We employed template matching and grid search methods to identify well-defined repeating earthquakes in order to minimize artefacts induced by variations in source location. We identified a total of 12 isolated patches in the fault zone that ruptured more than twice in a 1 yr period after the M 7.9 earthquake. We applied the coda wave interferometry technique to the waveform data of the 34 identified repeating earthquakes to estimate velocity changes between the first and subsequent events in each cluster. We found that major post-seismic velocity changes occurred in the southwestern part of the rupture area, where the main rupture was initiated and characterized by thrust motion, while the Beichuan area in the northeastern part of the rupture zone appears to experience very little post-seismic velocity changes.

Description: In this study, we demonstrate the feasibility of imaging broad-band (10–150 s) Rayleigh wave phase velocity maps on a continental scale using ambient noise tomography (ANT). We obtain broad-band Rayleigh waves from cross-correlations of ambient noise data between all station pairs of USArray and measure the dispersion curves from these cross-correlations at a period band of 10–150 s. The large-scale dense USArray enables us to obtain over 500 000 surface wave paths which cover the contiguous United States densely. Using these paths, we generate Rayleigh wave phase velocity maps at 10–150 s periods. Our phase velocity maps are similar to other reported phase velocity maps based on ambient noise data at short periods (〈50 s) and based on earthquake data at intermediate/long periods (50–90 s). This study extends ANT from short/intermediate periods (〈50 s) to long periods up to 150 s in a continental scale of the USA. These broad-band phase velocity maps from ANT can be used to construct 3-D lithospheric and asthenospheric velocity structures.

Description: In finite-difference (FD) method, numerical dispersion is the dominant factor influencing the accuracy of seismic modelling. Various optimized FD schemes for scalar wave modelling have been proposed to reduce grid dispersion, while the optimized time–space domain FD schemes for elastic wave modelling have not been fully investigated yet. In this paper, an optimized FD scheme with Equivalent Staggered Grid (ESG) for elastic modelling has been developed. We start from the constant P - and S -wave speed elastic wave equations and then deduce analytical plane wave solutions in the wavenumber domain with eigenvalue decomposition method. Based on the elastic plane wave solutions, three new time–space domain dispersion relations of ESG elastic modelling are obtained, which are represented by three equations corresponding to P -, S - and converted-wave terms in the elastic equations, respectively. By using these new relations, we can study the dispersion errors of different spatial FD terms independently. The dispersion analysis showed that different spatial FD terms have different errors. It is therefore suggested that different FD coefficients to be used to approximate the three spatial derivative terms. In addition, the relative dispersion error in L2-norm is minimized through optimizing FD coefficients using Newton's method. Synthetic examples have demonstrated that this new optimal FD schemes have superior accuracy for elastic wave modelling compared to Taylor-series expansion and optimized space domain FD schemes.

Description: Markov chain Monte-Carlo (McMC) sampling generates correlated random samples such that their distribution would converge to the true distribution only as the number of samples tends to infinity. In practice, McMC is found to be slow to converge, convergence is not guaranteed to be achieved in finite time, and detection of convergence requires the use of subjective criteria. Although McMC has been used for decades as the algorithm of choice for inference in complex probability distributions, there is a need to seek alternative approaches, particularly in high dimensional problems. Walker & Curtis ( 2014 ) developed a method for Bayesian inversion of 2-D spatial data using an exact sampling alternative to McMC which always draws independent samples of the target distribution. Their method thus obviates the need for convergence and removes the concomitant bias exhibited by finite sample sets. Their algorithm is nevertheless computationally intensive and requires large memory. We propose a more efficient method for Bayesian inversion of categorical variables, such as geological facies that requires no sampling at all. The method is based on a 2-D Hidden Markov Model (2D-HMM) over a grid of cells where observations represent localized data constraining each cell. The data in our example application are seismic attributes such as P - and S -wave impedances and rock density; our categorical variables are the hidden states and represent the geological rock types in each cell—facies of distinct subsets of lithology and fluid combinations such as shale, brine-sand and gas-sand. The observations at each location are assumed to be generated from a random function of the hidden state (facies) at that location, and to be distributed according to a certain probability distribution that is independent of hidden states at other locations – an assumption referred to as ‘localized likelihoods’. The hidden state (facies) at a location cannot be determined solely by the observation at that location as it also depends on prior information concerning the spatial distribution of other hidden states elsewhere. The prior information is included in the inversion in the form of a training image which represents a conceptual depiction of local geologies that might be expected, but other forms of prior information can be used in the method as desired. The method provides direct estimates of posterior marginal probability distributions over each variable, so these do not need to be estimated from samples such as in McMC. Nevertheless, in the case that samples are desired, these can be generated. On a 2-D test example the method is shown to outperform previous methods significantly, at a fraction of the computational cost. In many foreseeable applications there are therefore no serious impediments to extending the method to 3-D cases.

Description: Using an up-to-date global plate rotation model, applied to the endpoints of preserved major spreading ridge isochrons, we have calculated the explicitly reconstructable length-weighted mean global half-spreading rate (HSR), ridge length and area production as a function of time since the end of the Cretaceous Normal Superchron at 83.0 Ma. Our calculations integrate uncertainties in rotation parameters and chron boundary ages with the partial sampling uncertainties arising from progressive subduction of older oceanic lithosphere and its preserved spreading record. This record of directly reconstructable oceanic ridge production provides a well-constrained baseline that can be compared to reconstructions that include the largely unconstrained extrapolated histories of entirely subducted oceanic plates. The directly reconstructable global mean HSR has not varied by more than ±15 per cent about an average rate of 28.4 ± 4.6 mm a –1 since 83 Ma. No long-term secular trend is evident: a maximum global mean half-rate of 32 ± 6 mm a –1 occurred from 33.1 Ma to about 25.8 Ma, with minima of 26 ± 5 mm a –1 between about 56 and 40.2 Ma, and 24 ± 1 mm a –1 since 3.2 Ma. Only this most recent interval has a rate that differs significantly (at ±2) from the long-term mean. The global, reconstructable ridge length at 56 Ma decreases by less than 15 per cent relative to the modern ridge system; by 83 Ma it has decreased by 38 per cent. These relatively high preserved ridge fractions mean that the estimated uncertainty due to partial sampling stays roughly equivalent to the estimated rotation model uncertainties, allowing long-term spreading rate variations of 〉20 per cent since the Late Cretaceous to be ruled out. In contrast, prior to 83 Ma too little oceanic lithosphere is preserved to reliably reconstruct global spreading rates.

Description: A double-correlation method is introduced to locate tremor sources based on stacks of complex, doubly-correlated tremor records of multiple triplets of seismographs back projected to hypothetical source locations in a geographic grid. Peaks in the resulting stack of moduli are inferred source locations. The stack of the moduli is a robust measure of energy radiated from a point source or point sources even when the velocity information is imprecise. Application to real data shows how double correlation focuses the source mapping compared to the common single correlation approach. Synthetic tests demonstrate the robustness of the method and its resolution limitations which are controlled by the station geometry, the finite frequency of the signal, the quality of the used velocity information and noise level. Both random noise and signal or noise correlated at time shifts that are inconsistent with the assumed velocity structure can be effectively suppressed. Assuming a surface wave velocity, we can constrain the source location even if the surface wave component does not dominate. The method can also in principle be used with body waves in 3-D, although this requires more data and seismographs placed near the source for depth resolution.

Description: The possibility of applying one explicit finite-difference (FD) scheme to all interior grid points (points not lying on a grid border) no matter what their positions are with respect to the material interface is one of the key factors of the computational efficiency of the FD modelling. Smooth or discontinuous heterogeneity of the medium is accounted for only by values of the effective grid moduli and densities. Accuracy of modelling thus very much depends on how these effective grid parameters are evaluated. We present an orthorhombic representation of a heterogeneous medium for the FD modelling. We numerically demonstrate its superior accuracy. Compared to the harmonic-averaging representation the orthorhombic representation is more accurate mainly in the case of strong surface waves that are especially important in local surface sedimentary basins. The orthorhombic representation is applicable to modelling seismic wave propagation and earthquake motion in isotropic models with material interfaces and smooth heterogeneities using velocity–stress, displacement–stress and displacement FD schemes on staggered, partly staggered, Lebedev and collocated grids.

Description: We present two independent, automated methods for estimating the absolute horizontal misorientation of seismic sensors. We apply both methods to 44 free-fall ocean-bottom seismometers (OBSs) of the RHUM-RUM experiment ( http://www.rhum-rum.net/ ). The techniques measure the 3-D directions of particle motion of (1) P -waves and (2) Rayleigh waves of earthquake recordings. For P -waves, we used a principal component analysis to determine the directions of particle motions (polarizations) in multiple frequency passbands. We correct for polarization deviations due to seismic anisotropy and dipping discontinuities using a simple fit equation, which yields significantly more accurate OBS orientations. For Rayleigh waves, we evaluated the degree of elliptical polarization in the vertical plane in the time and frequency domain. The results obtained for the RHUM-RUM OBS stations differed, on average, by 3.1° and 3.7° between the methods, using circular mean and median statistics, which is within the methods’ estimate uncertainties. Using P -waves, we obtained orientation estimates for 31 ocean-bottom seismometers with an average uncertainty (95 per cent confidence interval) of 11° per station. For 7 of these OBS, data coverage was sufficient to correct polarization measurements for underlying seismic anisotropy and dipping discontinuities, improving their average orientation uncertainty from 11° to 6° per station. Using Rayleigh waves, we obtained misorientation estimates for 40 OBS, with an average uncertainty of 16° per station. The good agreement of results obtained using the two methods indicates that they should also be useful for detecting misorientations of terrestrial seismic stations.

Description: We investigate the elastic and anelastic structure of the lowermost mantle at the western edge of the Pacific large low shear velocity province (LLSVP) by inverting a collection of S and ScS waveforms. The transverse component data were obtained from F-net for 31 deep earthquakes beneath Tonga and Fiji, filtered between 12.5 and 200 s. We observe a regional variation of S and ScS arrival times and amplitude ratios, according to which we divide our region of interest into three subregions. For each of these subregions, we then perform 1-D (depth-dependent) waveform inversions simultaneously for radial profiles of shear wave velocity ( V S ) and seismic quality factor ( Q ). Models for all three subregions show low V S and low Q structures from 2000 km depth down to the core–mantle boundary. We further find that V S and Q in the central subregion, sampling the Caroline plume, are substantially lower than in the surrounding regions, whatever the depth. In the central subregion, V S -anomalies with respect to PREM (d V S ) and Q are about –2.5 per cent and 216 at a depth of 2850 km, and –0.6 per cent and 263 at a depth of 2000 km. By contrast, in the two other regions, d V S and Q are –2.2 per cent and 261 at a depth of 2850 km, and –0.3 per cent and 291 at a depth of 2000 km. At depths greater than ~2500 km, these differences may indicate lateral variations in temperature of ~100 K within the Pacific LLSVP. At shallower depths, they may be due to the temperature difference between the Caroline plume and its surroundings, and possibly to a small fraction of iron-rich material entrained by the plume.

Description: Earth's free oscillations excited by a mega-thrust earthquake were observed by a continent-scale array of groundwater monitoring sites for the first time. After the occurrence of the 2011 Tohoku M w 9.0 earthquake, water level records at 43 out of 216 wells in the China mainland revealed long-period free oscillation signals. In the time domain, these free oscillations exhibit globe circling Rayleigh surface waves. In some single wells, even the globe-circling Rayleigh wave R7 was visible, which travels three times around the Earth after the first arrival and appears about 10 hr after the earthquake occurrence in the present case. The spectral analysis shows that the principal oscillatory fluctuations seen in the water level records correspond to the spheroidal modes 0 S l ( l  = 2–31 for frequencies between 0.3 and 5.0 mHz) of the Earth's free oscillation. Especially at quiet sites, the spheroidal modes at very low frequencies (〈1.5 mHz) can be identified with high signal-to-noise ratios. Using signal enhancement methods (product spectrum over 43 wells), even the gravest modes of these oscillations can be detected. The results suggest that groundwater level arrays can be considered as a low-cost complementary tool to study the Earth's free oscillations excited by great earthquakes. Additionally, the site-specific aquifer response may provide further insight into local hydrogeological conditions.

Description: We performed numerical simulations of the 2011 deep-seated Akatani landslide in central Japan to understand the dynamic evolution of friction of the landslide. By comparing the forces obtained from numerical simulation to those resolved from seismic waveform inversion, the coefficient of the friction during sliding was investigated in the range of 0.1–0.4. The simulation assuming standard Coulomb friction shows that the forces obtained by the seismic waveform inversion are well explained using a constant friction of μ = 0.3. A small difference between the residuals of Coulomb simulation and a velocity-dependent simulation suggests that the coefficient of friction over the volume is well constrained as 0.3 most of time during sliding. It suggests the sudden loss of shearing resistance at the onset of sliding, that is, sudden drop of the initial coefficient of friction in our model, which accelerates the deep-seated landslide. Our numerical simulation calibrated by seismic data provides the evolution of dynamic friction with a reasonable resolution in time, which is difficult to obtain from a conventional runout simulation, or seismic waveform inversion alone.

Description: We explore thermal convection of a fluid with a temperature-dependent viscosity in a basally heated 3-D spherical shell using linear stability analyses and numerical experiments, while considering the application of our results to terrestrial planets. The inner to outer radius ratio of the shell f assumed in the linear stability analyses is in the range of 0.11–0.88. The critical Rayleigh number R c for the onset of thermal convection decreases by two orders of magnitude as f increases from 0.11 to 0.88, when the viscosity depends sensitively on the temperature, as is the case for real mantle materials. Numerical simulations carried out in the range of f  = 0.11–0.55 show that a thermal boundary layer (TBL) develops both along the surface and bottom boundaries to induce cold and hot plumes, respectively, when f is 0.33 or larger. However, for smaller f values, a TBL develops only on the bottom boundary. Convection occurs in the stagnant-lid regime where the root mean square velocity on the surface boundary is less than 1 per cent of its maximum at depth, when the ratio of the viscosity at the surface boundary to that at the bottom boundary exceeds a threshold that depends on f . The threshold decreases from 10 6.5 at f  = 0.11 to 10 4 at f  = 0.55. If the viscosity at the base of the convecting mantle is 10 20 –10 21  Pa s, the Rayleigh number exceeds R c for Mars, Venus and the Earth, but does not for the Moon and Mercury; convection is unlikely to occur in the latter planets unless the mantle viscosity is much lower than 10 20  Pa s and/or the mantle contains a strong internal heat source.

Description: We introduce a technique to compute exact anelastic sensitivity kernels in the time domain using parsimonious disk storage. The method is based on a reordering of the time loop of time-domain forward/adjoint wave propagation solvers combined with the use of a memory buffer. It avoids instabilities that occur when time-reversing dissipative wave propagation simulations. The total number of required time steps is unchanged compared to usual acoustic or elastic approaches. The cost is reduced by a factor of 4/3 compared to the case in which anelasticity is partially accounted for by accommodating the effects of physical dispersion. We validate our technique by performing a test in which we compare the K α sensitivity kernel to the exact kernel obtained by saving the entire forward calculation. This benchmark confirms that our approach is also exact. We illustrate the importance of including full attenuation in the calculation of sensitivity kernels by showing significant differences with physical-dispersion-only kernels.

Description: Subducting oceanic lithosphere is an example of a thin sheet-like object whose characteristic lateral dimension greatly exceeds its thickness. Here we exploit this property to derive a new hybrid boundary-integral/thin sheet (BITS) representation of subduction that combines in a single equation all the forces acting on the sheet: gravity, internal resistance to bending and stretching, and the tractions exerted by the ambient mantle. For simplicity, we limit ourselves to 2-D. We solve the BITS equations using a discrete Lagrangian approach in which the sheet is represented by a set of vertices connected by edges. Instantaneous solutions for the sinking speed of a slab attached to a trailing flat sheet obey a scaling law of the form V / V Stokes  = fct(St), where V Stokes is a characteristic Stokes sinking speed and St is the sheet's flexural stiffness. Time-dependent solutions for the evolution of the sheet's shape and thickness show that these are controlled by the viscosity ratio between the sheet and its surroundings. An important advantage of the BITS approach is the possibility of generalizing the sheet's rheology, either to a viscosity that varies along the sheet or to a non-Newtonian shear-thinning rheology.

Description: Real Earth media are not perfectly elastic. Instead, they attenuate propagating mechanical waves. This anelastic phenomenon in wave propagation can be modeled by a viscoelastic mechanical model consisting of several standard linear solids. Using this viscoelastic model, we approximate a constant Q over a frequency band of interest. We use a four-element viscoelastic model with a trade-off between accuracy and computational costs to incorporate Q into 2-D time-domain first-order velocity–stress wave equations. To improve the computational efficiency, we limit the Q in the model to a list of discrete values between 2 and 1000. The related stress and strain relaxation times that characterize the viscoelastic model are pre-calculated and stored in a database for use by the finite-difference calculation. A viscoelastic finite-difference scheme that is second order in time and fourth order in space is developed based on the MacCormack algorithm. The new method is validated by comparing the numerical result with analytical solutions that are calculated using the generalized reflection/transmission coefficient method. The synthetic seismograms exhibit greater than 95 per cent consistency in a two-layer viscoelastic model. The dispersion generated from the simulation is consistent with the Kolsky–Futterman dispersion relationship.

Description: After 83 yr, the great normal-faulting earthquake of 1933 March 2, which took place off the Japan Trench and produced a devastating tsunami on the Sanriku coast and damaging waves in Hawaii, remains the largest recorded normal-faulting earthquake. This study uses advanced methods to investigate this event using far-field seismological and tsunami data and complements a sister study by Uchida et al. which used exclusively arrival times at Japanese stations. Our relocation of the main shock (39.22°N, 144.45°E, with a poorly constrained depth of less than 40 km) places it in the outer trench slope, below a seafloor depth of ~6500 m, in a region of horst-and-graben structure, with fault scarps approximately parallel to the axis of the Japan Trench. Relocated aftershocks show a band of genuine shallow aftershocks parallel to the Japan Trench under the outer trench slope and a region of post-mainshock events landward of the trench axis that occur over roughly the same latitude range and are thought to be the result of stress transfer to the interplate thrust boundary following the normal-faulting rupture. Based on a combination of P -wave first motions and inversion of surface wave spectral amplitudes, we propose a normal-faulting focal mechanism ( = 200°, = 61° and = 271°) and a seismic moment M 0 = (7 ± 1) x 10 28 dyn cm ( M w = 8.5). A wide variety of data, including the distribution of isoseismals, the large magnitudes (up to 8.9) proposed by early investigators before the standardization of magnitude scales, estimates of energy-to-moment ratios and the tentative identification of a T wave at Pasadena (and possibly Riverside), clearly indicate that this seismic source was exceptionally rich in high-frequency wave energy, suggesting a large apparent stress and a sharp rise time, and consistent with the behaviour of many smaller shallow normal-faulting earthquakes. Hydrodynamic simulations based on a range of possible sources consistent with the above findings, including a compound rupture on two opposite-facing normal-faulting segments, are in satisfactory agreement with tsunami observations in Hawaii, where run-up reached 3 m, causing significant damage. This study emphasizes the need to include off-trench normal-faulting earthquake sources in global assessments of tsunami hazards emanating from the subduction of old and cold plates, whose total length of trenches exceed 20 000 km, even though only a handful of great such events are known with confidence in the instrumental record.

Description: To refine the 3-D seismic velocity model in the greater Parkfield, California region, a new data set including regular earthquakes, shots, quarry blasts and low-frequency earthquakes (LFEs) was assembled. Hundreds of traces of each LFE family at two temporary arrays were stacked with time–frequency domain phase weighted stacking method to improve signal-to-noise ratio. We extend our model resolution to lower crustal depth with LFE data. Our result images not only previously identified features but also low velocity zones (LVZs) in the area around the LFEs and the lower crust beneath the southern Rinconada Fault. The former LVZ is consistent with high fluid pressure that can account for several aspects of LFE behaviour. The latter LVZ is consistent with a high conductivity zone in magnetotelluric studies. A new Vs model was developed with S picks that were obtained with a new autopicker. At shallow depth, the low Vs areas underlie the strongest shaking areas in the 2004 Parkfield earthquake. We relocate LFE families and analyse the location uncertainties with the NonLinLoc and tomoDD codes. The two methods yield similar results.

Description: The aftershock distribution of the 1933 Sanriku-oki outer trench earthquake is estimated by using modern relocation methods and a newly developed velocity structure to examine the spatial extent of the source-fault and the possibility of a triggered interplate seismicity. In this study, we first examined the regional data quality of the 1933 earthquake based on smoked-paper records and then relocated the earthquakes by using the 3-D velocity structure and double-difference method. The improvements of hypocentre locations using these methods were confirmed by the examination of recent earthquakes that are accurately located based on ocean bottom seismometer data. The results show that the 1933 aftershocks occurred under both the outer- and inner-trench-slope regions. In the outer-trench-slope region, aftershocks are distributed in a ~280-km-long area and their depths are shallower than 50 km. Although we could not constrain the fault geometry from the hypocentre distribution, the depth distribution suggests the whole lithosphere is probably not under deviatoric tension at the time of the 1933 earthquake. The occurrence of aftershocks under the inner trench slope was also confirmed by an investigation of waveform frequency difference between outer and inner trench earthquakes as recorded at Mizusawa. The earthquakes under the inner trench slope were shallow (depth less double equals30 km) and the waveforms show a low-frequency character similar to the waveforms of recent, precisely located earthquakes in the same area. They are also located where recent activity of interplate thrust earthquakes is high. These suggest that the 1933 outer-trench-slope main shock triggered interplate earthquakes, which is an unusual case in the order of occurrence in contrast with the more common pairing of a large initial interplate shock with subsequent outer-slope earthquakes. The off-trench earthquakes are distributed about 80 km width in the trench perpendicular direction. This wide width cannot be explained from a single high-angle fault confined at a shallow depth (depth less double equals50 km). The upward motion of the 1933 tsunami waveform records observed at Sanriku coast also cannot be explained from a single high-angle west-dipping normal fault. If we consider additional fault, involvement of high-angle, east-dipping normal faults can better explain the tsunami first motion and triggering of the aftershock in a wide area under the outer trench slope. Therefore multiple off-trench normal faults may have activated during the 1933 earthquake. We also relocated recent (2001–2012) seismicity by the same method. The results show that the present seismicity in the outer-trench-slope region can be divided into several groups along the trench. Comparison of the 1933 rupture dimensions based on our aftershock relocations with the morphologies of fault scarps in the outer trench slope suggest that the rupture was limited to the region where fault scarps are largely trench parallel and cross cut the seafloor spreading fabric. These findings imply that bending geometry and structural segmentation of the incoming plate largely controls the spatial extent of the 1933 seismogenic faulting. In this shallow rupture model for this largest outer trench earthquake, triggered seismicity in the forearc and structural control of faulting represent an important deformation styles for off-trench and shallow megathrust zones.

Description: The Cadzow rank-reduction method can be effectively utilized in simultaneously denoising and reconstructing 5-D seismic data that depend on four spatial dimensions. The classic version of Cadzow rank-reduction method arranges the 4-D spatial data into a level-four block Hankel/Toeplitz matrix and then applies truncated singular value decomposition (TSVD) for rank reduction. When the observed data are extremely noisy, which is often the feature of real seismic data, traditional TSVD cannot be adequate for attenuating the noise and reconstructing the signals. The reconstructed data tend to contain a significant amount of residual noise using the traditional TSVD method, which can be explained by the fact that the reconstructed data space is a mixture of both signal subspace and noise subspace. In order to better decompose the block Hankel matrix into signal and noise components, we introduced a damping operator into the traditional TSVD formula, which we call the damped rank-reduction method. The damped rank-reduction method can obtain a perfect reconstruction performance even when the observed data have extremely low signal-to-noise ratio. The feasibility of the improved 5-D seismic data reconstruction method was validated via both 5-D synthetic and field data examples. We presented comprehensive analysis of the data examples and obtained valuable experience and guidelines in better utilizing the proposed method in practice. Since the proposed method is convenient to implement and can achieve immediate improvement, we suggest its wide application in the industry.

Description: An exact analytical solution is presented for the effective dynamic transverse shear modulus in a heterogeneous fluid-filled porous solid containing cylindrical inclusions. The complex and frequency-dependent properties of the dynamic shear modulus are caused by the physical mechanism of mesoscopic-scale wave-induced fluid flow whose scale is smaller than wavelength but larger than the size of pores. Our model consists of three phases: a long cylindrical inclusion, a cylindrical shell of poroelastic matrix material with different mechanical and/or hydraulic properties than the inclusion and an outer region of effective homogeneous medium of laterally infinite extent. The behavior of both the inclusion and the matrix is described by Biot's consolidation equations, whereas the surrounding effective medium which is used to describe the effective transverse shear properties of the inner poroelastic composite is assumed to be a viscoelastic solid whose complex transverse shear modulus needs to be determined. The determined effective transverse shear modulus is used to quantify the S -wave attenuation and velocity dispersion in heterogeneous fluid-filled poroelastic rocks. The calculation shows the relaxation frequency and relative position of various fluid saturation dispersion curves predicted by this study exhibit very good agreement with those of a previous 2-D finite-element simulation. For the double-porosity model (inclusions having a different solid frame than the matrix but the same pore fluid as the matrix) the effective shear modulus also exhibits a size-dependent characteristic that the relaxation frequency moves to lower frequencies by two orders of magnitude if the radius of the cylindrical poroelastic composite increases by one order of magnitude. For the patchy-saturation model (inclusions having the same solid frame as the matrix but with a different pore fluid from the matrix), the heterogeneity in pore fluid cannot cause any attenuation in the transverse shear modulus at all. A comparison with the case of spherical inclusions illustrates that the transverse shear modulus for the cylindrical inclusion exhibits more S -wave attenuation than spherical inclusions.

Description: Secondary microseism sources are pressure fluctuations close to the ocean surface. They generate acoustic P waves that propagate in water down to the ocean bottom where they are partly reflected and partly transmitted into the crust to continue their propagation through the Earth. We present the theory for computing the displacement power spectral density of secondary microseism P waves recorded by receivers in the far field. In the frequency domain, the P -wave displacement can be modeled as the product of (1) the pressure source, (2) the source site effect that accounts for the constructive interference of multiply reflected P waves in the ocean, (3) the propagation from the ocean bottom to the stations and (4) the receiver site effect. Secondary microseism P waves have weak amplitudes, but they can be investigated by beamforming analysis. We validate our approach by analysing the seismic signals generated by typhoon Ioke (2006) and recorded by the Southern California Seismic Network. Backprojecting the beam onto the ocean surface enables to follow the source motion. The observed beam centroid is in the vicinity of the pressure source derived from the ocean wave model WAVEWATCH III R . The pressure source is then used for modeling the beam and a good agreement is obtained between measured and modeled beam amplitude variation over time. This modeling approach can be used to invert P -wave noise data and retrieve the source intensity and lateral extent.

Description: In the context of the verification of the Comprehensive Nuclear-Test Ban Treaty in the marine environment, we present a new discriminant based on the empirical observation that hydroacoustic phases recorded at T -phase stations from explosive sources in the water column feature a systematic inverse dispersion, with lower frequencies traveling slower, which is absent from signals emanating from earthquake sources. This difference is present even in the case of the so-called ‘hotspot earthquakes’ occurring inside volcanic edifices featuring steep slopes leading to efficient seismic–acoustic conversions, which can lead to misidentification of such events as explosions when using more classical duration-amplitude discriminants. We propose an algorithm for the compensation of the effect of dispersion over the hydroacoustic path based on a correction to the spectral phase of the ground velocity recorded by the T -phase station, computed individually from the dispersion observed on each record. We show that the application of a standard amplitude-duration algorithm to the resulting compensated time-series satisfactorily identifies records from hotspot earthquakes as generated by dislocation sources, and present a full algorithm, lending itself to automation, for the discrimination of explosive and earthquake sources of hydroacoustic signals at T -phase stations. The only sources not readily identifiable consist of a handful of complex explosions which occurred in the 1970s, believed to involve the testing of advanced weaponry, and which should be independently identifiable through routine vetting by analysts. While we presently cannot provide a theoretical justification to the observation that only explosive sources generate dispersed T phases, we hint that this probably reflects a simpler, and more coherent distribution of acoustic energy among the various modes constituting the wave train, than in the case of dislocation sources embedded in the solid Earth.

Description: This work presents an innovative strategy to enhance the resolution of surface wave tomography obtained from ambient noise cross-correlation ( C 1 ) by bridging asynchronous seismic networks through the correlation of coda of correlations ( C 3 ). Rayleigh wave group dispersion curves show consistent results between synchronous and asynchronous stations. Rayleigh wave group traveltimes are inverted to construct velocity–period maps with unprecedented resolution for a region covering Mexico and the southern United States. The resulting period maps are then used to regionalize dispersion curves in order to obtain local 1-D shear velocity models ( V S ) of the crust and uppermost mantle in every cell of a grid of 0.4°. The 1-D structures are obtained by iteratively adding layers until reaching a given misfit, and a global tomography model is considered as an input for depths below 150 km. Finally, a high-resolution 3-D V S model is obtained from these inversions. The major structures observed in the 3-D model are in agreement with the tectonic-geodynamic features and with previous regional and local studies. It also offers new insights to understand the present and past tectonic evolution of the region.

Description: Archaeomagnetic field models cover longer timescales than historical models and may therefore resolve the motion of geomagnetic features on the core–mantle boundary (CMB) in a more meaningful statistical sense. Here we perform a detailed appraisal of archaeomagnetic field models to infer some aspects of the physics of the outer core. We characterize and compare the identification and tracking of reversed flux patches (RFPs) in order to assess the RFPs robustness. We find similar behaviour within a family of models but differences among different families, suggesting that modelling strategy is more influential than data set. Similarities involve recurrent positions of RFPs, but no preferred direction of motion is found. The tracking of normal flux patches shows similar qualitative behaviour confirming that RFPs identification and tracking is not strongly biased by their relative weakness. We also compare the tracking of RFPs with that of the historical field model gufm1 and with seismic anomalies of the lowermost mantle to explore the possibility that RFPs have preferred locations prescribed by lower mantle lateral heterogeneity. The archaeomagnetic field model that most resembles the historical field is interpreted in terms of core dynamics and core–mantle thermal interactions. This model exhibits correlation between RFPs and low seismic shear velocity in co-latitude and a shift in longitude. These results shed light on core processes, in particular we infer toroidal field lines with azimuthal orientation below the CMB and large fluid upwelling structures with a width of about 80° (Africa) and 110° (Pacific) at the top of the core. Finally, similar preferred locations of RFPs in the past 9 and 3 kyr of the same archaeomagnetic field model suggest that a 3 kyr period is sufficiently long to reliably detect mantle control on core dynamics. This allows estimating an upper bound of 220–310 km for the magnetic boundary layer thickness below the CMB.

Description: Determination of a response of the sea water column to teleseismic plane wave is important to suppress adverse effects of water reverberations in calculating receiver functions (RFs) using ocean-bottom seismometer (OBS) records. We present a novel non-linear waveform analysis method using the simulated annealing algorithm to determine such a water-layer response recorded by an OBS array. We then demonstrate its usefulness for the RF estimation through its application to synthetic and observed data. Synthetic experiments suggest that the water-layer response constrained in this way has a potential to improve RFs of OBS records drastically even in the high-frequency range (to 4 Hz). By applying it to data observed by the OBS array around the Kii Peninsula, southwestern Japan, we identified a low-velocity zone at the top of the subducting Philippine Sea plate. This zone may represent the incoming fluid-rich sediment layer that has been reported by active-source seismic survey.

Description: We develop an automated strategy for discriminating deep microseismic events from shallow ones on the basis of the waveforms recorded on a limited number of surface receivers. Machine-learning techniques are employed to explore the relationship between event hypocentres and seismic features of the recorded signals in time, frequency and time–frequency domains. We applied the technique to 440 microearthquakes –1.7 〈  M w  〈 1.29, induced by an underground cavern collapse in the Napoleonville Salt Dome in Bayou Corne, Louisiana. Forty different seismic attributes of whole seismograms including degree of polarization and spectral attributes were measured. A selected set of features was then used to train the system to discriminate between deep and shallow events based on the knowledge gained from existing patterns. The cross-validation test showed that events with depth shallower than 250 m can be discriminated from events with hypocentral depth between 1000 and 2000 m with 88 per cent and 90.7 per cent accuracy using logistic regression and artificial neural network models, respectively. Similar results were obtained using single station seismograms. The results show that the spectral features have the highest correlation to source depth. Spectral centroids and 2-D cross-correlations in the time–frequency domain are two new seismic features used in this study that showed to be promising measures for seismic event classification. The used machine-learning techniques have application for efficient automatic classification of low energy signals recorded at one or more seismic stations.

Description: We conduct a numerical experiment to investigate potential bias in measurements of S -wave splitting (apparent differences between the arrival times of SH and SV phases) for waves propagating close to the core–mantle boundary (CMB) in the D'' layer. The bias is defined as the discrepancy between shear wave splitting measured from finite frequency synthetic seismograms (‘apparent splitting’) and the splitting predicted by ray theory, which is a high-frequency approximation. For simple isotropic models, we find biases which are typically between 0.5 and 4 s, depending on the model, the Q structure and the dominant period of the synthetics. The bias increases for lower frequencies or lower Q values. The epicentral distance at which the bias starts depends on the frequency and the Q structure. We also compute synthetics for models based on mineral physics (using the elastic constants under lower-mantle pressure and temperature conditions, taking into account the phase transition from Mg-perovskite to Mg-post-perovskite) and geodynamics (the thermal boundary layer) and find that the depth of the positive velocity jump associated with the phase transition and the depth range over which the velocity decreases (due to temperature increases) in the thermal boundary layer significantly influence the wavefield in the lowermost mantle. For example, in cold regions beneath subduction zones, wavefields for SH and SV differ greatly due to the steep velocity decrease close to the CMB. For complex models, apparent splitting can also arise from the possibility that low amplitude direct phases might be overlooked, and larger amplitude later phases might instead incorrectly be picked as the direct arrival. Biases of the type investigated in this study combine with other sources of uncertainty for splitting in D'' (e.g. the correction for upper-mantle anisotropy and the difference between SH and SV ray paths) to make a precise evaluation of the anisotropy in D'' difficult.

Description: Geodetic, geologic and palaeomagnetic data reveal that Oregon (western USA) rotates clockwise at 0.3 to 1.0° Ma –1 (relative to North America) about an axis near the Idaho–Oregon–Washington border, while northeast Washington is relatively fixed. This rotation has been going on for at least 15 Ma. The Yakima fold and thrust belt (YFTB) forms the boundary between northern Oregon and central Washington where convergence of the clockwise-rotating Oregon block is apparently accommodated. North–south shortening across the YFTB has been thought to occur in a fan-like manner, increasing in rate to the west. We obtained high-accuracy, high-density geodetic GPS measurements in 2012–2014 that are used with earlier GPS measurements from the 1990s to characterize YFTB kinematics. The new results show that the deformation associated with the YFTB starts at the Blue Mountains Anticline in northern Oregon and extends north beyond the Frenchman Hills in Washington, past the epicentre of the 1872 M w 7.0 Entiat earthquake to 49°N. The north–south strain rate across the region is 2 to 3 x 10 –9 yr –1 between the volcanic arc and the eastern edge of the YFTB (241.0°E); east of there it drops to about 10 –9 yr –1 . At the eastern boundary of the YFTB, faults and earthquake activity are truncated by a north-trending, narrow zone of deformation that runs along the Pasco Basin and Moses Lake regions near 240.9°E. This zone, abutting the Department of Energy Hanford Nuclear Reservation, accommodates about 0.5 mm yr –1 of east to northeast shortening. A similar zone of N-trending transpression is seen along 239.9°E where there is a change in the strike of the Yakima folds. The modern deformation of the YFTB is about 600 km wide from south to north and internally may be controlled by pre-existing crustal structure.

Description: The Central Anatolian orogenic plateau is represented by young volcanism, rapid plateau uplift and distinctive (past and active) tectonic deformation. In this study, we consider observational data in terms of regional present-day geodynamics in the region. The residual topography of Central Anatolia was derived to define the regional isostatic conditions according to Airy isostasy and infer the potential role of ‘dynamic topography’. 2-D thermomechanical forward models for coupled mantle-lithosphere flow/deformation were conducted along an N–S directional profile through the region (e.g. northern/Pontides, interior and southern/Taurides). These models were based on seismic tomography data that provide estimates about the present-day mantle thermal structure beneath the Anatolian plate. We compare the modelling results with calculated residual topography and independent data sets of geological deformation, gravity and high surface heat flow/widespread geothermal activity. Model results suggest that there is ~1 km of mantle flow induced dynamic topography associated with the sublithospheric flow driven by the seismically inferred mantle structure. The uprising mantle may have also driven the asthenospheric source of volcanism in the north (e.g. Galatia volcanic province) and the Cappadocia volcanic province in the south while elevating the surface in the last 10 Myr. Our dynamic topography calculations emphasize the role of vertical forcing under other orogenic plateaux underlain by relatively thin crust and low-density asthenospheric mantle.

Description: Post-depositional reductive diagenesis usually results in partial or entire cleansing of the pristine palaeomagnetic signal, therefore, its intensity is important to be assessed for sediments that are in the purpose of retrieving palaeomagnetic information. Grain size, rock magnetic and geochemical studies on the entire core, along with scanning electron microscope observations and X-ray diffraction analyses for representative samples were carried out on a Holocene sediment core retrieved from the deep water part of Huguangyan maar lake (HGY), southeast China. The pristine magnetic mineral assemblage of the studied core is domianted by superparamagnetic (SP) and stable single domain titanomagnetite, and high coercivity minerals are not detectable. Based on down-core variations of the average grain size ( M Z ), total organic carbon (TOC), detrital elements (Al, Ti, Fe and Mn) and the concentration and mineralogy of magnetic minerals, the studied core could be divided into three subsections. The uppermost subsection is the least affected by diagenesis, with detrital titanomagnetite as the dominant magnetic mineral. This is owing to low TOC contents, but high detrital input generated by weak Asian summer monsoon intensity during the late Holocene. The intermediate subsection shows down-core progressively enhanced dissolution of detrital titanomagnetite, and concomitant formation of authigenic pyrite and siderite, which indicates down-core progressively enhanced diagenesis generated by down-core progressive increasing TOC content, but decreasing detrital input as the result of down-core progressively strengthened Asian summer monsoon intensity. The pristine magnetic mineral assemblage has been profoundly modified in the lowermost subsection. At certain positions of the lowermost subsection, detrital titanomagnetite has been even completely dissolved via diagenesis, giving place to authigenic pyrite and siderite. High TOC content, but low detrital input generated from strong Asian summer monsoon intensity during the early Holocene are accountable for intensive diagenesis in the lowermost subsection. Complete erasing of detrital magnetic input signal at certain positions of the lowermost subsection, and considerable formation of authigenic siderite indicate that palaeomagnetic records of the studied core have been significantly compromised. The studied core has relatively higher TOC content, lower detrital matter content, calmer sedimentary environments, and less DO available at its water–sediment interface than the cores retrieved at relatively shallower water depths, which all contribute to its relatively stronger diagenesis. Progressive thickening of the upper two subsections with increasing water depth is owing to progressive increase in sedimentation rate with increasing water depth, which is the key factor in determining the thickness of each diagenetic subsection of cores from HGY. It would be better that lake sediments for palaeomagnetic investigations collected at a water depth shallower than the depth of its thermocline.

Description: A velocity ( Vs ) and structure model is derived for the Los Angeles Basin, California based on ambient-noise surface wave and receiver-function analysis, using data from a low-cost, short-duration, dense broad-band survey (LASSIE) deployed across the basin. The shear wave velocities show lateral variations at the Compton-Los Alamitos and the Whittier Faults. The basement beneath the Puente Hills–San Gabriel Valley shows an unusually high velocity (~4.0 km s –1 ) and indicates the presence of schist. The structure of the model shows that the basin is a maximum of 8 km deep along the profile and that the Moho rises to a depth of 17 km under the basin. The basin has a stretch factor of 2.6 in the centre grading to 1.3 at the edges and is in approximate isostatic equilibrium.

Description: In many geological systems, inversion of density stratification sets in Rayleigh–Taylor (RT) instabilities, leading to an ascent of relatively low-density materials through the high-density overburden in the form of diapirs. These diapirs often originate from dipping low-density layers. This study aims to show how the initial tilt of such source layers can control the ascent behaviour of diapirs initiated by RT instabilities. Using two-layer viscous models we produced RT instabilities in physical experiments, and investigated the effects of source-layer tilts ( β ). Our experiments suggest that these diapirs ascend with contrasting lateral spreading rates in the up and down slope directions, resulting in their axi-asymmetric geometry. However, their heads retain a circular outline on the horizontal top surface, where the upwelling axis is located away from their geometric centre in the upslope direction. In this paper, we present a series of experimental models to demonstrate the spectrum of axi-symmetric to -asymmetric geometrical transitions with increasing β . Our experiments also reveal that when β is large (〉4°) the diapirs become unstable, resulting in a continuous migration of their stems in the upslope direction. Using the volume of fluid method we ran computational fluid dynamic (CFD) simulations to study the underlying hydrodynamics of axi-asymmetric diapiric growth. The CFD simulations show that β 〉 0° conditions develop stronger flow vortices on the downslope side of an ascending diapir, leading to a pressure difference between the up- and downslope flanks. Such a differential pressure causes the diapir head to spread at a faster rate in the tilt direction. An estimate of the asymmetric spreading rates is given as a function of β . Our present study provides a fundamental understanding of the hydrodynamic flow structure responsible for the asymmetric growth of RT instabilities on tilted source layers, as applicable to a wide range of large-scale geological settings, such as sedimentary basins and subduction zones.

Description: We present a catalogue of full seismic moment tensors for 63 events from Uturuncu volcano in Bolivia. The events were recorded during 2011–2012 in the PLUTONS seismic array of 24 broad-band stations. Most events had magnitudes between 0.5 and 2.0 and did not generate discernible surface waves; the largest event was M w 2.8. For each event we computed the misfit between observed and synthetic waveforms, and we used first-motion polarity measurements to reduce the number of possible solutions. Each moment tensor solution was obtained using a grid search over the 6-D space of moment tensors. For each event, we show the misfit function in eigenvalue space, represented by a lune. We identify three subsets of the catalogue: (1) six isotropic events, (2) five tensional crack events, and (3) a swarm of 14 events southeast of the volcanic centre that appear to be double couples. The occurrence of positively isotropic events is consistent with other published results from volcanic and geothermal regions. Several of these previous results, as well as our results, cannot be interpreted within the context of either an oblique opening crack or a crack-plus-double-couple model. Proper characterization of uncertainties for full moment tensors is critical for distinguishing among physical models of source processes.

Description: The 2-D acoustic wave equation is commonly solved numerically by finite-difference (FD) methods in which the accuracy of solution is significantly affected by the FD stencils. The commonly used cross stencil can reach either only second-order accuracy for space domain dispersion-relation-based FD method or (2 M )th-order accuracy along eight specific propagation directions for time–space domain dispersion-relation-based FD method, if the conventional (2 M )th-order spatial FD and second-order temporal FD are used to discretize the equation. One other newly developed rhombus stencil can reach arbitrary even-order accuracy. However, this stencil adds significantly to computational cost when the operator length is large. To achieve a balance between the solution accuracy and efficiency, we develop a new FD stencil to solve the 2-D acoustic wave equation. This stencil is a combination of the cross stencil and rhombus stencil. A cross stencil with an operator length parameter M is used to approximate the spatial partial derivatives while a rhombus stencil with an operator length parameter N together with the conventional second-order temporal FD is employed in approximating the temporal partial derivatives. Using this stencil, a new FD scheme is developed; we demonstrate that this scheme can reach (2 M )th-order accuracy in space and (2 N )th-order accuracy in time when spatial FD coefficients and temporal FD coefficients are derived from respective dispersion relation using Taylor-series expansion (TE) method. To further increase the accuracy, we derive the FD coefficients by employing the time–space domain dispersion relation of this FD scheme using TE. We also use least-squares (LS) optimization method to reduce dispersion at high wavenumbers. Dispersion analysis, stability analysis and modelling examples demonstrate that our new scheme has greater accuracy and better stability than conventional FD schemes, and thus can adopt large time steps. To reduce the extra computational cost resulting from adopting the new stencil, we apply the variable spatial operator length schemes. Adopting our new FD scheme, characterized by new stencil, LS-based optimization, variable operator lengths and larger time step, modelling efficiency is significantly improved.

Description: The geomechanical analysis of a highly compartmentalized reservoir is performed to simulate the seafloor subsidence due to gas production. The available observations over the hydrocarbon reservoir consist of bathymetric surveys carried out before and at the end of a 10-yr production life. The main goal is the calibration of the reservoir compressibility c M , that is, the main geomechanical parameter controlling the surface response. Two conceptual models are considered: in one (i) c M varies only with the depth and the vertical effective stress (heterogeneity due to lithostratigraphic variability); in another (ii) c M varies also in the horizontal plane, that is, it is spatially distributed within the reservoir stratigraphic units. The latter hypothesis accounts for a possible partitioning of the reservoir due to the presence of sealing faults and thrusts that suggests the idea of a block heterogeneous system with the number of reservoir blocks equal to the number of uncertain parameters. The method applied here relies on an ensemble-based data assimilation (DA) algorithm (i.e. the ensemble smoother, ES), which incorporates the information from the bathymetric measurements into the geomechanical model response to infer and reduce the uncertainty of the parameter c M . The outcome from conceptual model (i) indicates that DA is effective in reducing the c M uncertainty. However, the maximum settlement still remains underestimated, while the areal extent of the subsidence bowl is overestimated. We demonstrate that the selection of the heterogeneous conceptual model (ii) allows to reproduce much better the observations thus removing a clear bias of the model structure. DA allows significantly reducing the c M uncertainty in the five blocks (out of the seven) characterized by large volume and large pressure decline. Conversely, the assimilation of land displacements only partially constrains the prior c M uncertainty in the reservoir blocks marginally contributing to the cumulative seafloor subsidence, that is, blocks with low pressure.

Description: Near-surface geophysical imaging is often performed by generating surface waves, and estimating the subsurface properties through inversion, that is, iteratively matching experimentally observed dispersion curves with predicted curves from a layered half-space model of the subsurface. Key to the effectiveness of inversion is the efficiency and accuracy of computing the dispersion curves and their derivatives. This paper presents improved methodologies for both dispersion curve and derivative computation. First, it is shown that the dispersion curves can be computed more efficiently by combining an unconventional complex-length finite element method (CFEM) to model the finite depth layers, with perfectly matched discrete layers (PMDL) to model the unbounded half-space. Second, based on analytical derivatives for theoretical dispersion curves, an approximate derivative is derived for the so-called effective dispersion curve for realistic geophysical surface response data. The new derivative computation has a smoothing effect on the computation of derivatives, in comparison with traditional finite difference (FD) approach, and results in faster convergence. In addition, while the computational cost of FD differentiation is proportional to the number of model parameters, the new differentiation formula has a computational cost that is almost independent of the number of model parameters. At the end, as confirmed by synthetic and real-life imaging examples, the combination of CFEM + PMDL for dispersion calculation and the new differentiation formula results in more accurate estimates of the subsurface characteristics than the traditional methods, at a small fraction of computational effort.

Description: Powerful subduction zone earthquakes rupture thousands of square kilometres along continental margins but at certain locations earthquake rupture terminates. To date, detailed knowledge of the parameters that govern seismic rupture and aftershocks is still incomplete. On 2015 September 16, the M w 8.3 Illapel earthquake ruptured a 200 km long stretch of the Central Chilean subduction zone, triggering a tsunami and causing significant damage. Here, we analyse the temporal and spatial pattern of the coseismic rupture and aftershocks in relation to the tectonic setting in the earthquake area. Aftershocks cluster around the area of maximum coseismic slip, in particular in lateral and downdip direction. During the first 24 hr after the main shock, aftershocks migrated in both lateral directions with velocities of approximately 2.5 and 5 km hr –1 . At the southern rupture boundary, aftershocks cluster around individual subducted seamounts that are related to the downthrusting Juan Fernández Ridge. In the northern part of the rupture area, aftershocks separate into an upper cluster (above 25 km depth) and a lower cluster (below 35 km depth). This dual seismic–aseismic transition in downdip direction is also observed in the interseismic period suggesting that it may represent a persistent feature for the Central Chilean subduction zone.

Description: This work aims to explore the ongoing tectonic activity of structures in the outermost sector of the Northern Apennines, which represents the active leading edge of the thrust belt and is dominated by compressive deformation. We have applied the Persistent Scatterer Interferometry (PSI) technique to obtain new insights into the present-day deformation pattern of the frontal area of the Northern Apennine. PSI has proved to be effective in detecting surface deformation of wide regions involved in low tectonic movements. We used 34 Envisat images in descending geometry over the period of time between 2004 and 2010, performing about 300 interferometric pairs. The analysis of the velocity maps and of the PSI time-series has allowed to observe ground deformation over the sector of the Po Plain between Piacenza and Reggio Emilia. The time-series of permanent GPS stations located in the study area, validated the results of the PSI technique, showing a good correlation with the PS time-series. The PS analysis reveals the occurrence of a well-known subsidence area on the rear of the Ferrara arc, mostly connected to the exploitation of water resources. In some instances, the PS velocity pattern reveals ground uplift (with mean velocities ranging from 1 to 2.8 mm yr –1 ) above active thrust-related anticlines of the Emilia and Ferrara folds, and part of the Pede-Apennine margin. We hypothesize a correlation between the observed uplift deformation pattern and the growth of the thrust-related anticlines. As the uplift pattern corresponds to known geological features, it can be used to constrain the seismo-tectonic setting, and a working hypothesis may involve that the active Emilia and Ferrara thrust folds would be characterized by interseismic periods possibly dominated by aseismic creep.

Description: Least-squares inversion of seismic arrivals can provide remarkably detailed models of the Earth's subsurface. However, cycle skipping associated with these oscillatory arrivals is the main cause for local minima in the least-squares objective function. Therefore, it is often difficult for descent methods to converge to the solution without an accurate initial large-scale velocity estimate. The low frequencies in the arrivals, needed to update the large-scale components in the velocity model, are usually unreliable or absent. To overcome this difficulty, we propose a multi-objective inversion scheme that uses the conventional least-squares functional along with an auxiliary data-domain objective. As the auxiliary objective effectively replaces the seismic arrivals by bumps, we call it the bump functional. The bump functional minimization can be made far less sensitive to cycle skipping and can deal with multiple arrivals in the data. However, it can only be used as an auxiliary objective since it usually does not provide a unique model after minimization even when the regularized-least-squares functional has a unique global minimum and hence a unique solution. The role of the bump functional during the multi-objective inversion is to guide the optimization towards the global minimum by pulling the trapped solution out of the local minima associated with the least-squares functional whenever necessary. The computational complexity of the bump functional is equivalent to that of the least-squares functional. In this paper, we describe various characteristics of the bump functional using simple and illustrative numerical examples. We also demonstrate the effectiveness of the proposed multi-objective inversion scheme by considering more realistic examples. These include synthetic and field data from a cross-well experiment, surface-seismic synthetic data with reflections and synthetic data with refracted arrivals at long offsets.

Description: In Europe, common input data types for seismic hazard evaluation include earthquake catalogues, seismic zonation models and ground motion models, all with well-constrained epistemic uncertainties. In contrast, neotectonic deformation models and their related uncertainties are rarely considered in earthquake forecasting and seismic hazard studies. In this study, for the first time in Europe, we developed a seismic hazard model based exclusively on active fault and geodynamic deformation models. We applied it to the External Dinarides, a slow-deforming fold-and-thrust belt in the Central Mediterranean. The two deformation models furnish consistent long-term earthquake rates above the M w 4.7 threshold on a latitude/longitude grid with 0.2° spacing. Results suggest that the use of deformation models is a valid alternative to empirical-statistical approaches in earthquake forecasting in slow-deforming regions of Europe. Furthermore, we show that the variability of different deformation models has a comparable effect on the peak ground motion acceleration uncertainty as do the ground motion prediction equations.

Description: The rocks in the crust and the upper mantle of the Earth are believed to exhibit electrical anisotropy to some extent. It is beneficial to further understand and recognize the propagation of the electromagnetic waves in the Earth by investigating the magnetotelluric (which is one of the main geophysical techniques to probe the deep structures in the Earth) responses of the media with anisotropic conductivity structures. In this study, we examine the magnetotelluric fields over an idealized 2-D model consisting of two segments with axially anisotropic conductivity structures overlying a perfect conductor basement by a quasi-static analytic approach. The resulting analytic solution could not only contribute to the electromagnetic induction theory in the anisotropic Earth but also serve as at least an initial standard solution which could be used to validate the reliability and accuracy of the numerical algorithms developed for modelling the magnetotelluric responses of the 2-D media with much more general anisotropic conductivity.

Description: Owing to the increasing availability of computational resources, in recent years the probabilistic solution of non-linear, geophysical inverse problems by means of sampling methods has become increasingly feasible. Nevertheless, we still face situations in which a Monte Carlo approach is not practical. This is particularly true in cases where the evaluation of the forward problem is computationally intensive or where inversions have to be carried out repeatedly or in a timely manner, as in natural hazards monitoring tasks such as earthquake early warning. Here, we present an alternative to Monte Carlo sampling, in which inferences are entirely based on a set of prior samples—that is, samples that have been obtained independent of a particular observed datum. This has the advantage that the computationally expensive sampling stage becomes separated from the inversion stage, and the set of prior samples—once obtained—can be reused for repeated evaluations of the inverse mapping without additional computational effort. This property is useful if the problem is such that repeated inversions of independent data have to be carried out. We formulate the inverse problem in a Bayesian framework and present a practical way to make posterior inferences based on a set of prior samples. We compare the prior sampling based approach to a Markov Chain Monte Carlo approach that samples from the posterior probability distribution. We show results for both a toy example, and a realistic seismological source parameter estimation problem. We find that the posterior uncertainty estimates obtained based on prior sampling can be considered conservative estimates of the uncertainties obtained by directly sampling from the posterior distribution.

Description: The Afar Depression and its adjacent areas are underlain by an upper mantle marked by some of the world's largest negative velocity anomalies, which are frequently attributed to the thermal influences of a lower-mantle plume. In spite of numerous studies, however, the existence of a plume beneath the area remains enigmatic, partially due to inadequate quantities of broad-band seismic data and the limited vertical resolution at the mantle transition zone (MTZ) depth of the techniques employed by previous investigations. In this study, we use an unprecedented quantity (over 14 500) of P -to- S receiver functions (RFs) recorded by 139 stations from 12 networks to image the 410 and 660 km discontinuities and map the spatial variation of the thickness of the MTZ. Non-linear stacking of the RFs under a 1-D velocity model shows robust P -to- S conversions from both discontinuities, and their apparent depths indicate the presence of an upper-mantle low-velocity zone beneath the entire study area. The Afar Depression and the northern Main Ethiopian Rift are characterized by an apparent 40–60 km depression of both MTZ discontinuities and a normal MTZ thickness. The simplest and most probable interpretation of these observations is that the apparent depressions are solely caused by velocity perturbations in the upper mantle and not by deeper processes causing temperature or hydration anomalies within the MTZ. Thickening of the MTZ on the order of 15 km beneath the southern Arabian Plate, southern Red Sea and western Gulf of Aden, which comprise the southward extension of the Afro-Arabian Dome, could reflect long-term hydration of the MTZ. A 20 km thinning of the MTZ beneath the western Ethiopian Plateau is observed and interpreted as evidence for a possible mantle plume stem originating from the lower mantle.

Description: We propose a procedure for uncertainty quantification in Probabilistic Tsunami Hazard Analysis (PTHA), with a special emphasis on the uncertainty related to statistical modelling of the earthquake source in Seismic PTHA (SPTHA), and on the separate treatment of subduction and crustal earthquakes (treated as background seismicity). An event tree approach and ensemble modelling are used in spite of more classical approaches, such as the hazard integral and the logic tree. This procedure consists of four steps: (1) exploration of aleatory uncertainty through an event tree, with alternative implementations for exploring epistemic uncertainty; (2) numerical computation of tsunami generation and propagation up to a given offshore isobath; (3) (optional) site-specific quantification of inundation; (4) simultaneous quantification of aleatory and epistemic uncertainty through ensemble modelling. The proposed procedure is general and independent of the kind of tsunami source considered; however, we implement step 1, the event tree, specifically for SPTHA, focusing on seismic source uncertainty. To exemplify the procedure, we develop a case study considering seismic sources in the Ionian Sea (central-eastern Mediterranean Sea), using the coasts of Southern Italy as a target zone. The results show that an efficient and complete quantification of all the uncertainties is feasible even when treating a large number of potential sources and a large set of alternative model formulations. We also find that (i) treating separately subduction and background (crustal) earthquakes allows for optimal use of available information and for avoiding significant biases; (ii) both subduction interface and crustal faults contribute to the SPTHA, with different proportions that depend on source-target position and tsunami intensity; (iii) the proposed framework allows sensitivity and deaggregation analyses, demonstrating the applicability of the method for operational assessments.

Description: The expanding fleet of broad-band ocean-bottom seismograph (OBS) stations is facilitating the study of the structure and seismicity of oceanic plates at regional scales. For continental studies, an important tool to characterize continental crust and mantle structure is the analysis of teleseismic P receiver functions. In the oceans, however, receiver functions potentially suffer from several limiting factors that are unique to ocean sites and plate structures. In this study, we model receiver functions for a variety of oceanic lithospheric structures to investigate the possibilities and limitations of receiver functions using OBS data. Several potentially contaminating effects are examined, including pressure reverberations from the water column for various ocean-floor depths and the effects of a layer of low-velocity marine sediments. These modelling results indicate that receiver functions from OBS data are difficult to interpret in the presence of marine sediments, but shallow-water sites in subduction zone forearcs may be suitable for constraining various crustal elements around the locked megathrust fault. We propose using a complementary approach based on transfer function modelling combined with a grid search approach that bypasses receiver functions altogether and estimates model properties directly from minimally processed waveforms. Using real data examples from the Cascadia Initiative, we show how receiver and transfer functions can be used to infer seismic properties of the oceanic plate in both shallow (Cascadia forearc) and deep (Juan de Fuca Ridge) ocean settings.

Description: Several works have reported that haematite has non-linear initial susceptibility at room temperature, like pyrrhotite or titanomagnetite, but there is no explanation for the observed behaviours yet. This study sets out to determine which physical property (grain size, foreign cations content and domain walls displacements) controls the initial susceptibility. The performed measurements include microprobe analysis to determine magnetic phases different to haematite; initial susceptibility (300 K); hysteresis loops, SIRM and backfield curves at 77 and 300 K to calculate magnetic parameters and minor loops at 77 K, to analyse initial susceptibility and magnetization behaviours below Morin transition. The magnetic moment study at low temperature is completed with measurements of zero field cooled–field cooled and AC susceptibility in a range from 5 to 300 K. The minor loops show that the non-linearity of initial susceptibility is closely related to Barkhausen jumps. Because of initial magnetic susceptibility is controlled by domain structure it is difficult to establish a mathematical model to separate magnetic subfabrics in haematite-bearing rocks.

Description: A new implementation of indirect boundary element method allows simulating the elastic wave propagation in complex configurations made of embedded regions that are homogeneous with irregular boundaries or flat layers. In an older implementation, each layer of a flat layered region would have been treated as a separated homogeneous region without taking into account the flat boundary information. For both types of regions, the scattered field results from fictitious sources positioned along their boundaries. For the homogeneous regions, the fictitious sources emit as in a full-space and the wave field is given by analytical Green's functions. For flat layered regions, fictitious sources emit as in an unbounded flat layered region and the wave field is given by Green's functions obtained from the discrete wavenumber (DWN) method. The new implementation allows then reducing the length of the discretized boundaries but DWN Green's functions require much more computation time than the full-space Green's functions. Several optimization steps are then implemented and commented. Validations are presented for 2-D and 3-D problems. Higher efficiency is achieved in 3-D.

Description: An important real world application of doublet flow occurs in well design of both geothermal and hydrocarbon reservoirs. A guiding principle for fluid management of injection and extraction wells is that mass balance is commonly assumed between the injected and produced fluid. Because the doublets are considered closed loops, the injection fluid is assumed to eventually reach the producer well and all the produced fluid ideally comes from stream tubes connected to the injector of the well pair making up the doublet. We show that when an aquifer background flow occurs, doublets will rarely retain closed loops of fluid recirculation. When the far-field flow rate increases relative to the doublet's strength, the area occupied by the doublet will diminish and eventually vanishes. Alternatively, rather than using a single injector (source) and single producer (sink), a linear array of multiple injectors separated by some distance from a parallel array of producers can be used in geothermal energy projects as well as in waterflooding of hydrocarbon reservoirs. Fluid flow in such an arrangement of parallel source-sink arrays is shown to be macroscopically equivalent to that of a line doublet. Again, any far-field flow that is strong enough will breach through the line doublet, which then splits into two vortices. Apart from fundamental insight into elementary flow dynamics, our new results provide practical clues that may contribute to improve the planning and design of doublets and direct line drives commonly used for flow management of groundwater, geothermal and hydrocarbon reservoirs.

Description: Though well known for layer boundaries, the use of amplitude-versus-offset (AVO) variations for non-welded boundaries like fractures is not yet investigated. Depending on the seismic wavelength used, fractures can be regarded as thin, compliant zones in rocks, in different scales. We explore the potential of multiangle AVO inversion of P-P and P-S reflections from a fracture to estimate fracture properties. We conduct laboratory experiments to measure reflection responses of dry and wet fractures. The observed P-P reflections of the wet fracture and the fracture aperture are very well predicted by the non-welded interface model. We invert the angle-dependent P-P reflectivity of the fracture to estimate both normal and tangential fracture compliances. The estimated value of the normal compliance is accurate, and it is also possible to obtain the value of the non-zero tangential compliance. We find that supplementing the information of converted P-S reflections in the AVO inversion greatly improves the estimate of the tangential compliance. The calculated compliance ratio clearly shows the existence of fluid in the fracture. This finding can be crucial for new applications in a wide range of scale—from earthquake seismology, deep and shallow seismic exploration, to non-destructive material testing.

Description: The present study evaluates the capacity of the Boom Clay as a host rock for disposal purposes, more precisely its seismic characterization, which may assess its long-term performance to store radioactive wastes. Although the formation is relatively uniform and homogeneous, there are embedded thin layers of septaria (carbonates) that may affect the integrity of the Boom Clay. Therefore, it is essential to locate these geobodies. The seismic data to characterize the Boom Clay has been acquired at the Kruibeke test site. The inversion, which allowed us to obtain the anisotropy parameters and seismic velocities of the clay, is complemented with further information such as log and laboratory data. The attenuation properties have been estimated from equivalent formations (having similar composition and seismic velocities). The inversion yields quite consistent results although the symmetry of the medium is unusual but physically possible, since the anisotropy parameter is negative. According to a time-domain calculation of the energy velocity at four frequency bands up to 900 Hz, velocity increases with frequency, a behaviour described by the Zener model. Then, we use this model to describe anisotropy and anelasticity that are implemented into the equation of motion to compute synthetic seismograms in the space–time domain. The technique is based on memory variables and the Fourier pseudospectral method. We have computed reflection coefficients of the septaria thin layer. At normal incidence, the P -wave coefficient vanishes at specific thicknesses of the layer and there is no conversion to the S wave. For example, calculations at 600 Hz show that for thicknesses of 1 m the septarium can be detected more easily since the amplitudes are higher (nearly 0.8). Converted PS waves have a high amplitude at large offsets (between 30° and 80°) and can be useful to identify the target on this basis. Moreover, we have investigated the effect of septaria embedded in the Boom Clay with several simulations, by considering a lateral partial continuity of the calcareous thin inclusions. The simulations with layers of calcareous material show continuity of the reflections even when the percentage of carbonate within the layer is very small (5–15 per cent), while for low content of the calcareous material, isolated septaria boulders generate diffraction events. We have also simulated the stacked seismic section obtained from processing of the field data. The matching between the field and synthetic sections is acceptable.

Description: A new Matched Filtering Algorithm (MFA) is proposed for detecting and analysing microseismic events recorded by downhole monitoring of hydraulic fracturing. This method requires a set of well-located template (‘parent’) events, which are obtained using conventional microseismic processing and selected on the basis of high signal-to-noise (S/N) ratio and representative spatial distribution of the recorded microseismicity. Detection and extraction of ‘child’ events are based on stacked, multichannel cross-correlation of the continuous waveform data, using the parent events as reference signals. The location of a child event relative to its parent is determined using an automated process, by rotation of the multicomponent waveforms into the ray-centred co-ordinates of the parent and maximizing the energy of the stacked amplitude envelope within a search volume around the parent's hypocentre. After correction for geometrical spreading and attenuation, the relative magnitude of the child event is obtained automatically using the ratio of stacked envelope peak with respect to its parent. Since only a small number of parent events require interactive analysis such as picking P - and S -wave arrivals, the MFA approach offers the potential for significant reduction in effort for downhole microseismic processing. Our algorithm also facilitates the analysis of single-phase child events, that is, microseismic events for which only one of the S - or P -wave arrivals is evident due to unfavourable S/N conditions. A real-data example using microseismic monitoring data from four stages of an open-hole slickwater hydraulic fracture treatment in western Canada demonstrates that a sparse set of parents (in this case, 4.6 per cent of the originally located events) yields a significant (more than fourfold increase) in the number of located events compared with the original catalogue. Moreover, analysis of the new MFA catalogue suggests that this approach leads to more robust interpretation of the induced microseismicity and novel insights into dynamic rupture processes based on the average temporal (foreshock–aftershock) relationship of child events to parents.

Description: Azimuthal anisotropy is a powerful tool to reveal information about both the present structure and past evolution of the mantle. Anisotropic images of the upper mantle are usually obtained by analysing various types of seismic observables, such as surface wave dispersion curves or waveforms, SKS splitting data, or receiver functions. These different data types sample different volumes of the earth, they are sensitive to different length scales, and hence are associated with different levels of uncertainties. They are traditionally interpreted separately, and often result in incompatible models. We present a Bayesian inversion approach to jointly invert these different data types. Seismograms for SKS and P phases are directly inverted using a cross-convolution approach, thus avoiding intermediate processing steps, such as numerical deconvolution or computation of splitting parameters. Probabilistic 1-D profiles are obtained with a transdimensional Markov chain Monte Carlo scheme, in which the number of layers, as well as the presence or absence of anisotropy in each layer, are treated as unknown parameters. In this way, seismic anisotropy is only introduced if required by the data. The algorithm is used to resolve both isotropic and anisotropic layering down to a depth of 350 km beneath two seismic stations in North America in two different tectonic settings: the stable Canadian shield (station FFC) and the tectonically active southern Basin and Range Province (station TA-214A). In both cases, the lithosphere–asthenosphere boundary is clearly visible, and marked by a change in direction of the fast axis of anisotropy. Our study confirms that azimuthal anisotropy is a powerful tool for detecting layering in the upper mantle.

Description: We have obtained new results in the statistical analysis of global earthquake catalogues with special attention to the largest earthquakes, and we examined the statistical behaviour of earthquake rate variations. These results can serve as an input for updating our recent earthquake forecast, known as the ‘Global Earthquake Activity Rate 1’ model (GEAR1), which is based on past earthquakes and geodetic strain rates. The GEAR1 forecast is expressed as the rate density of all earthquakes above magnitude 5.8 within 70 km of sea level everywhere on earth at 0.1 x 0.1 degree resolution, and it is currently being tested by the Collaboratory for Study of Earthquake Predictability. The seismic component of the present model is based on a smoothed version of the Global Centroid Moment Tensor (GCMT) catalogue from 1977 through 2013. The tectonic component is based on the Global Strain Rate Map, a ‘General Earthquake Model’ (GEM) product. The forecast was optimized to fit the GCMT data from 2005 through 2012, but it also fit well the earthquake locations from 1918 to 1976 reported in the International Seismological Centre-Global Earthquake Model (ISC-GEM) global catalogue of instrumental and pre-instrumental magnitude determinations. We have improved the recent forecast by optimizing the treatment of larger magnitudes and including a longer duration (1918–2011) ISC-GEM catalogue of large earthquakes to estimate smoothed seismicity. We revised our estimates of upper magnitude limits, described as corner magnitudes, based on the massive earthquakes since 2004 and the seismic moment conservation principle. The new corner magnitude estimates are somewhat larger than but consistent with our previous estimates. For major subduction zones we find the best estimates of corner magnitude to be in the range 8.9 to 9.6 and consistent with a uniform average of 9.35. Statistical estimates tend to grow with time as larger earthquakes occur. However, by using the moment conservation principle that equates the seismic moment rate with the tectonic moment rate inferred from geodesy and geology, we obtain a consistent estimate of the corner moment largely independent of seismic history. These evaluations confirm the above-mentioned corner magnitude value. The new estimates of corner magnitudes are important both for the forecast part based on seismicity as well as the part based on geodetic strain rates. We examine rate variations as expressed by annual earthquake numbers. Earthquakes larger than magnitude 6.5 obey the Poisson distribution. For smaller events the negative-binomial distribution fits much better because it allows for earthquake clustering.

Description: In order to improve our understanding of hazardous underground cavities, the development and collapse of a ~200 m wide salt solution mining cavity was seismically monitored in the Lorraine basin in northeastern France. The microseismic events show a swarm-like behaviour, with clustering sequences lasting from seconds to days, and distinct spatiotemporal migration. Observed microseismic signals are interpreted as the result of detachment and block breakage processes occurring at the cavity roof. Body wave amplitude patterns indicated the presence of relatively stable source mechanisms, either associated with dip-slip and/or tensile faulting. Signal overlaps during swarm activity due to short interevent times, the high-frequency geophone recordings and the limited network station coverage often limit the application of classical source analysis techniques. To overcome these shortcomings, we investigated the source mechanisms through different procedures including modelling of observed and synthetic waveforms and amplitude spectra of some well-located events, as well as modelling of peak-to-peak amplitude ratios for the majority of the detected events. We extended the latter approach to infer the average source mechanism of many swarming events at once, using multiple events recorded at a single three component station. This methodology is applied here for the first time and represents a useful tool for source studies of seismic swarms and seismicity clusters. The results obtained with different methods are consistent and indicate that the source mechanisms for at least 50 per cent of the microseismic events are remarkably stable, with a predominant thrust faulting regime with faults similarly oriented, striking NW–SE and dipping around 35°–55°. This dominance of consistent source mechanisms might be related to the presence of a preferential direction of pre-existing crack or fault structures. As an interesting byproduct, we demonstrate, for the first time directly on seismic data, that the source radiation pattern significantly controls the detection capability of a seismic station and network.

Description: We show that higher modes are an important component of high-frequency Rayleigh waves in the cross-correlations over sedimentary basins. The particle motions provide a good test for distinguishing and separating the fundamental from the first higher mode, with the fundamental mode having retrograde and the first higher mode having prograde motion in the 1–10 s period of interest. The basement depth controls the cut-off period of the first higher mode, which coincides with a rapid increase (over period) in the particle-motion ellipticity or H / V ratio of the fundamental mode. The strong higher mode we observed is not only due to the low-velocity sedimentary layer but also due to the noise sources with significant radial component such as the basin edge scattering. It is important to correctly identify the mode order when inverting the dispersion curves because misidentifying the higher mode as fundamental will lead to an anomalous high V SV velocity.

Description: We introduce a ‘double-difference’ method for the inversion for seismic wave speed structure based on adjoint tomography. Differences between seismic observations and model predictions at individual stations may arise from factors other than structural heterogeneity, such as errors in the assumed source-time function, inaccurate timings and systematic uncertainties. To alleviate the corresponding non-uniqueness in the inverse problem, we construct differential measurements between stations, thereby reducing the influence of the source signature and systematic errors. We minimize the discrepancy between observations and simulations in terms of the differential measurements made on station pairs. We show how to implement the double-difference concept in adjoint tomography, both theoretically and practically. We compare the sensitivities of absolute and differential measurements. The former provide absolute information on structure along the ray paths between stations and sources, whereas the latter explain relative (and thus higher resolution) structural variations in areas close to the stations. Whereas in conventional tomography a measurement made on a single earthquake-station pair provides very limited structural information, in double-difference tomography one earthquake can actually resolve significant details of the structure. The double-difference methodology can be incorporated into the usual adjoint tomography workflow by simply pairing up all conventional measurements; the computational cost of the necessary adjoint simulations is largely unaffected. Rather than adding to the computational burden, the inversion of double-difference measurements merely modifies the construction of the adjoint sources for data assimilation.

Description: The creation of a magnitude-homogenized catalogue is often one of the most fundamental steps in seismic hazard analysis. The process of homogenizing multiple catalogues of earthquakes into a single unified catalogue typically requires careful appraisal of available bulletins, identification of common events within multiple bulletins and the development and application of empirical models to convert from each catalogue's native scale into the required target. The database of the International Seismological Center (ISC) provides the most exhaustive compilation of records from local bulletins, in addition to its reviewed global bulletin. New open-source tools are developed that can utilize this, or any other compiled database, to explore the relations between earthquake solutions provided by different recording networks, and to build and apply empirical models in order to harmonize magnitude scales for the purpose of creating magnitude-homogeneous earthquake catalogues. These tools are described and their application illustrated in two different contexts. The first is a simple application in the Sub-Saharan Africa region where the spatial coverage and magnitude scales for different local recording networks are compared, and their relation to global magnitude scales explored. In the second application the tools are used on a global scale for the purpose of creating an extended magnitude-homogeneous global earthquake catalogue. Several existing high-quality earthquake databases, such as the ISC-GEM and the ISC Reviewed Bulletins, are harmonized into moment magnitude to form a catalogue of more than 562 840 events. This extended catalogue, while not an appropriate substitute for a locally calibrated analysis, can help in studying global patterns in seismicity and hazard, and is therefore released with the accompanying software.

Description: This paper introduces a novel approach to constructing an effective pre-conditioner for finite-difference (FD) electromagnetic modelling in geophysical applications. This approach is based on introducing an FD contraction operator, similar to one developed for integral equation formulation of Maxwell's equation. The properties of the FD contraction operator were established using an FD analogue of the energy equality for the anomalous electromagnetic field. A new pre-conditioner uses a discrete Green's function of a 1-D layered background conductivity. We also developed the formulae for an estimation of the condition number of the system of FD equations pre-conditioned with the introduced FD contraction operator. Based on this estimation, we have established that the condition number is bounded by the maximum conductivity contrast between the background conductivity and actual conductivity. When there are both resistive and conductive anomalies relative to the background, the new pre-conditioner is advantageous over using the 1-D discrete Green's function directly. In our numerical experiments with both resistive and conductive anomalies, for a land geoelectrical model with 1:10 contrast, the method accelerates convergence of an iterative method (BiCGStab) by factors of 2–2.5, and in a marine example with 1:50 contrast, by a factor of 4.6, compared to direct use of the discrete 1-D Green's function as a pre-conditioner.

Description: A series of important geological events occurred in the Tibetan Plateau area during the Jurassic, such as the collision of the Lhasa and Qiangtang Terranes, the closure of the Meso-Tethyan Ocean, the opening of the Neo-Tethyan Ocean and the cessation of the mega-monsoon. The ~3000 m thick Jurassic sedimentary sequence in the Qiangtang Basin on the central Tibetan Plateau, which is called the Yanshiping (YSP) Group, recorded these geological events. However, the chronology of the sequence is surprisingly poorly constrained. Here, we perform a detailed palaeomagnetic analysis on the ~1060 m thick middle and upper portions of the YSP Group (the Xiali and Suowa Formations) in the YSP section of the eastern Qiangtang Basin. Three bivalve zones at stratigraphic intervals of ~40–140, 640–800 and 940–1040 m are identified, which yield a Bathonian–Callovian age for the Lower Xiali Fm., a Callovian–Oxfordian age for the Lower Suowa Fm. and an Oxfordian–Kimmeridgian age for the Upper Suowa Fm. A total of 544 oriented palaeomagnetic samples were collected from the section. By combining thermal and alternating field demagnetizations, clear characteristic remanent magnetization (ChRM) directions are isolated for most of the samples. The robust ChRM directions pass fold and reversals tests, which support the primary nature of the ChRMs and yield a palaeopole at 76.8°N/297.2°E (dp = 2.2°, dm = 3.7°). A total of 27 normal and 26 reversed polarity zones were successfully recorded in the section. Combined with fossil age constraints, results suggest that the section is plausibly composed of a Callovian-Early Kimmeridgian age sedimentary sequence.

Description: Reaching the global minimum of a waveform misfit function requires careful choices about the nonlinear optimization, preconditioning and regularization methods underlying an inversion. Because waveform inversion problems are susceptible to erratic convergence associated with strong nonlinearity, one or two test cases are not enough to reliably inform such decisions. We identify best practices, instead, using four seismic near-surface problems, one regional problem and two global problems. To make meaningful quantitative comparisons between methods, we carry out hundreds of inversions, varying one aspect of the implementation at a time. Comparing nonlinear optimization algorithms, we find that limited-memory BFGS provides computational savings over nonlinear conjugate gradient methods in a wide range of test cases. Comparing preconditioners, we show that a new diagonal scaling derived from the adjoint of the forward operator provides better performance than two conventional preconditioning schemes. Comparing regularization strategies, we find that projection, convolution, Tikhonov regularization and total variation regularization are effective in different contexts. Besides questions of one strategy or another, reliability and efficiency in waveform inversion depend on close numerical attention and care. Implementation details involving the line search and restart conditions have a strong effect on computational cost, regardless of the chosen nonlinear optimization algorithm.

Description: The conventional finite-difference time-domain (FDTD) method for elastic waves suffers from the staircasing error when applied to model a curved free surface because of its structured grid. In this work, an improved, stable and accurate 3-D FDTD method for elastic wave modelling on a curved free surface is developed based on the finite volume method and enlarged cell technique (ECT). To achieve a sufficiently accurate implementation, a finite volume scheme is applied to the curved free surface to remove the staircasing error; in the mean time, to achieve the same stability as the FDTD method without reducing the time step increment, the ECT is introduced to preserve the solution stability by enlarging small irregular cells into adjacent cells under the condition of conservation of force. This method is verified by several 3-D numerical examples. Results show that the method is stable at the Courant stability limit for a regular FDTD grid, and has much higher accuracy than the conventional FDTD method.

Description: Pressure solution creep (PSC) is an important elementary process in rock friction at high temperatures where solubilities of rock-forming minerals are significantly large. It significantly changes the frictional resistance and enhances time-dependent strengthening. A recent microphysical model for PSC-involved friction of clay–quartz mixtures, which can explain a transition between dilatant and non-dilatant deformation (d-nd transition), was modified here and implemented in dynamic earthquake sequence simulations. The original model resulted in essentially a kind of rate- and state-dependent friction (RSF) law, but assumed a constant friction coefficient for clay resulting in zero instantaneous rate dependency in the dilatant regime. In this study, an instantaneous rate dependency for the clay friction coefficient was introduced, consistent with experiments, resulting in a friction law suitable for earthquake sequence simulations. In addition, a term for time-dependent strengthening due to PSC was added which makes the friction law logarithmically rate-weakening in the dilatant regime. The width of the zone in which clasts overlap or, equivalently, the interface porosity involved in PSC plays a role as the state variable. Such a concrete physical meaning of the state variable is a great advantage in future modelling studies incorporating other physical processes such as hydraulic effects. Earthquake sequence simulations with different pore pressure distributions demonstrated that excess pore pressure at depth causes deeper rupture propagation with smaller slip per event and a shorter recurrence interval. The simulated ruptures were arrested a few kilometres below the point of pre-seismic peak stress at the d-nd transition and did not propagate spontaneously into the region of pre-seismic non-dilatant deformation. PSC weakens the fault against slow deformation and thus such a region cannot produce a dynamic stress drop. Dynamic rupture propagation further down to brittle-plastic transition, evidenced by geological observations, would require even smaller frictional resistance at coseismic slip rate, suggesting the importance of implementation of dynamic weakening activated at coseismic slip rates for more realistic simulation of earthquake sequences. The present models produced much smaller afterslip at deeper parts of arrested ruptures than those with logarithmic RSF laws because of a more significant rate-strengthening effect due to linearly viscous PSC. Detailed investigation of afterslip would give a clue to understand the deformation mechanism which controls shear resistance of the fault in a region of arrest of earthquake ruptures.

Description: Surface wave tomography routinely uses empirically scaled density model in the inversion of dispersion curves for shear wave speeds of the crust and uppermost mantle. An improperly selected empirical scaling relationship between density and shear wave speed can lead to unrealistic density models beneath certain tectonic formations such as sedimentary basins. Taking the Sichuan basin east to the Tibetan plateau as an example, we investigate the differences between density profiles calculated from four scaling methods and their effects on Rayleigh wave phase velocities. Analytical equations for 1-D layered models and adjoint tomography for 3-D models are used to examine the trade-off between density and S -wave velocity structures at different depth ranges. We demonstrate that shallow density structure can significantly influence phase velocities at short periods, and thereby affect the shear wave speed inversion from phase velocity data. In particular, a deviation of 25 per cent in the initial density model can introduce an error up to 5 per cent in the inverted shear velocity at middle and lower crustal depths. Therefore one must pay enough attention in choosing a proper velocity–density scaling relationship in constructing initial density model in Rayleigh wave inversion for crustal shear velocity structure.

Description: Application of the ambient noise surface wave tomography method (ANT) for determination of the upper-mantle structure requires data on long-periodic noise ( T 〉 40 s). The ANT technique implies that noise sources are distributed almost uniformly over the surface. This is practically true for short-periodic noise, however, it is not so in the case of long periods. In this paper we show that the main contribution to noise at long periods is caused by signals from earthquakes. In some cases, they may strongly distort noise cross-correlation. This leads to an incorrect determination of surface wave velocity dispersion curves. To minimize such a distortion we propose two means: (1) to use records of noise for the periods when there is no clustering of earthquakes, such as aftershocks of strong events; (2) to stack cross-correlation functions for a period of at least three years in order to achieve sufficient uniformity of earthquake locations. Validity of this approach is demonstrated by ANT results for Europe. Tomographic reconstruction of Rayleigh wave group velocities for 10–100 s measured along interstation paths was carried out in a central part of Western Europe where resolving power of the data was the highest. Locally averaged dispersion curves were inverted to vertical S -wave velocity sections in this area. The results correspond closely to known features of the structure of the region, namely: strong difference of the crust and upper-mantle structure at the opposite sides from the Tornquist–Teisseyre Line down to ~ 250 km, penetration of high-velocity material of East European Platform lithosphere under Carpathians, as well as penetration of low-velocity asthenospheric layer from the Carpathian region towards the northeast.

Description: Data on the evolution of Earth's magnetic field intensity are important for understanding the geodynamo and planetary evolution. However, the paleomagnetic record in rocks may be adversely affected by many physical processes, which must be taken into account when analysing the palaeointensity database. This is especially important in the light of an ongoing debate regarding core thermal conductivity values, and how these relate to the Precambrian geodynamo. Here, we demonstrate that several data sets in the Precambrian palaeointensity database overestimate the true paleofield strength due to the presence of non-ideal carriers of palaeointensity signals and/or viscous re-magnetizations. When the palaeointensity overestimates are removed, the Precambrian database does not indicate a robust change in geomagnetic field intensity during the Mesoproterozoic. These findings call into question the recent claim that the solid inner core formed in the Mesoproterozoic, hence constraining the thermal conductivity in the core to ‘moderate’ values. Instead, our analyses indicate that the presently available palaeointensity data are insufficient in number and quality to constrain the timing of solid inner core formation, or the outstanding problem of core thermal conductivity. Very young or very old inner core ages (and attendant high or low core thermal conductivity values) are consistent with the presently known history of Earth's field strength. More promising available data sets that reflect long-term core structure are geomagnetic reversal rate and field morphology. The latter suggests changes that may reflect differences in Archean to Proterozoic core stratification, whereas the former suggest an interval of geodynamo hyperactivity at ca. 550 Ma.

Description: Sensitivity analysis with synthetic models is widely used in seismic tomography as a means for assessing the spatial resolution of solutions produced by, in most cases, linear or iterative nonlinear inversion schemes. The most common type of synthetic reconstruction test is the so-called checkerboard resolution test in which the synthetic model comprises an alternating pattern of higher and lower wave speed (or some other seismic property such as attenuation) in 2-D or 3-D. Although originally introduced for application to large inverse problems for which formal resolution and covariance could not be computed, these tests have achieved popularity, even when resolution and covariance can be computed, by virtue of being simple to implement and providing rapid and intuitive insight into the reliability of the recovered model. However, checkerboard tests have a number of potential drawbacks, including (1) only providing indirect evidence of quantitative measures of reliability such as resolution and uncertainty, (2) giving a potentially misleading impression of the range of scale-lengths that can be resolved, and (3) not giving a true picture of the structural distortion or smearing that can be caused by the data coverage. The widespread use of synthetic reconstruction tests in seismic tomography is likely to continue for some time yet, so it is important to implement best practice where possible. The goal of this paper is to develop the underlying theory and carry out a series of numerical experiments in order to establish best practice and identify some common pitfalls. Based on our findings, we recommend (1) the use of a discrete spike test involving a sparse distribution of spikes, rather than the use of the conventional tightly spaced checkerboard; (2) using data coverage (e.g. ray-path geometry) inherited from the model constrained by the observations (i.e. the same forward operator or matrix), rather than the data coverage obtained by solving the forward problem through the synthetic model; (3) carrying out multiple tests using structures of different scale length; (4) taking special care with regard to what can be inferred when using synthetic structures that closely mimic what has been recovered in the observation-based model; (5) investigating the range of structural wavelengths that can be recovered using realistic levels of imposed data noise; and (6) where feasible, assessing the influence of model parametrization error, which arises from making a choice as to how structure is to be represented.

Description: In November 2010, intense seismic activity including 29 events with a magnitude above 5.0, started in the western part of the Gulf of Aden, where the structure of the oceanic spreading ridge is characterized by a series of N115°-trending slow-spreading segments set within an EW-trending rift. Using signals recorded by permanent and temporary networks in Djibouti and Yemen, we located 1122 earthquakes, with a magnitude ranging from 2.1 to 5.6 from 2010 November 1 to 2011 March 31. By looking in detail at the space–time distribution of the overall seismicity, and both the frequency and the moment tensor of large earthquakes, we re-examine the chronology of this episode. In addition, we also interpret the origin of the activity using high-resolution bathymetric data, as well as from observations of seafloor cable damage caused by high temperatures and lava flows. The analysis allows us to identify distinct active areas. First, we interpret that this episode is mainly related to a diking event along a specific ridge segment, located at E044°. In light of previous diking episodes in nearby subaerial rift segments, for which field constraints and both seismic and geodetic data exist, we interpret the space–time evolution of the seismicity of the first few days. Migration of earthquakes suggests initial magma ascent below the segment centre. This is followed by a southeastward dike propagation below the rift immediately followed by a northwestward dike propagation below the rift ending below the northern ridge wall. The cumulative seismic moment associated with this sequence reaches 9.1 x 10 17 Nm, and taking into account a very low seismic versus geodetic moment, we estimate a horizontal opening of ~0.58–2.9 m. The seismic activity that followed occurred through several bursts of earthquakes aligned along the segment axis, which are interpreted as short dike intrusions implying fast replenishment of the crustal magma reservoir feeding the dikes. Over the whole period, the opening is estimated to be ~1.76–8.8 m across the segment. A striking feature of this episode is that the seismicity remained confined within one individual segment, whereas the adjacent en-echelon segments were totally quiescent, suggesting that the magma supply system of one segment is disconnected from those of the neighbouring segments. Second, we identify activity induced by the first intrusion with epicentres aligned along an N035°E-trending, ~30 km long at the northwestern end of the active opening segment. This group encompasses more than seven earthquakes with magnitude larger than 5.0, and with strike-slip focal mechanisms consistent with the faults identified in the bathymetry and the structural pattern of the area. We propose that a transform fault is currently in formation which indicates an early stage of the ridge segmentation, at the locus of the trend change of the spreading ridge, which also corresponds to the boundary between a clear oceanic lithosphere and the zone of transform between continental and oceanic crust.

Description: We develop an approach for simulating acousto-elastic wave phenomena, including scattering from fluid–solid boundaries, where the solid is allowed to be anisotropic, with the discontinuous Galerkin method. We use a coupled first-order elastic strain-velocity, acoustic velocity–pressure formulation, and append penalty terms based on interior boundary continuity conditions to the numerical (central) flux so that the consistency condition holds for the discretized discontinuous Galerkin weak formulation. We incorporate the fluid–solid boundaries through these penalty terms and obtain a stable algorithm. Our approach avoids the diagonalization into polarized wave constituents such as in the approach based on solving elementwise Riemann problems.

Description: We challenge the perspective that seismicity could contribute to polar motion by arguing quantitatively that, in first approximation and on the average, interseismic deformations can compensate for it. This point is important because what we must simulate and observe in Earth Orientation Parameter time-series over intermediate timescales of decades or centuries is the residual polar motion resulting from the two opposing processes of coseismic and interseismic deformations. In this framework, we first simulate the polar motion caused by only coseismic deformations during the longest period available of instrumental seismicity, from 1900 to present, using both the CMT and ISC-GEM catalogues. The instrumental seismicity covering a little longer than one century does not represent yet the average seismicity that we should expect on the long term. Indeed, although the simulation shows a tendency to move the Earth rotation pole towards 133°E at the average rate of 16.5 mm yr –1 , this trend is still sensitive to individual megathrust earthquakes, particularly to the 1960 Chile and 1964 Alaska earthquakes. In order to further investigate this issue, we develop a global seismicity model (GSM) that is independent from any earthquake catalogue and that describes the average seismicity along plate boundaries on the long term by combining information about present-day plate kinematics with the Anderson theory of faulting, the seismic moment conservation principle and a few other assumptions. Within this framework, we obtain a secular polar motion of 8 mm yr –1 towards 112.5°E that is comparable with that estimated from 1900 to present using the earthquake catalogues, although smaller by a factor of 2 in amplitude and different by 20° in direction. Afterwards, in order to reconcile the idea of a secular polar motion caused by earthquakes with our simplest understanding of the seismic cycle, we adapt the GSM in order to account for interseismic deformations and we use it to quantify, for the first time ever, their contribution to polar motion. Taken together, coseismic and interseismic deformations make the rotation pole wander around the north pole with maximum polar excursions of about 1 m. In particular, the rotation pole moves towards about Newfoundland when the interseismic contribution dominates over the coseismic ones (i.e. during phases of low seismicity or, equivalently, when most of the fault system associated with plate boundaries is locked). When megathrust earthquakes occur, instead, the rotation pole is suddenly shifted in an almost opposite direction, towards about 133°E.

Description: It is well known that large earthquakes generally trigger aftershock sequences. However, the duration of those sequences is unclear due to the gradual power-law decay with time. The triggering time is assumed to be infinite in the epidemic type aftershock sequence (ETAS) model, a widely used statistical model to describe clustering phenomena in observed earthquake catalogues. This assumption leads to the constraint that the power-law exponent p of the Omori-Utsu decay has to be larger than one to avoid supercritical conditions with accelerating seismic activity on long timescales. In contrast, seismicity models based on rate- and state-dependent friction observed in laboratory experiments predict p ≤ 1 and a finite triggering time scaling inversely to the tectonic stressing rate. To investigate this conflict, we analyse an ETAS model with finite triggering times, which allow smaller values of p . We use synthetic earthquake sequences to show that the assumption of infinite triggering times can lead to a significant bias in the maximum likelihood estimates of the ETAS parameters. Furthermore, it is shown that the triggering time can be reasonably estimated using real earthquake catalogue data, although the uncertainties are large. The analysis of real earthquake catalogues indicates mainly finite triggering times in the order of 100 days to 10 years with a weak negative correlation to the background rate, in agreement with expectations of the rate- and state-friction model. The triggering time is not the same as the apparent duration, which is the time period in which aftershocks dominate the seismicity. The apparent duration is shown to be strongly dependent on the mainshock magnitude and the level of background activity. It can be much shorter than the triggering time. Finally, we perform forward simulations to estimate the effective forecasting period, which is the time period following a mainshock, in which ETAS simulations can improve rate estimates after the occurrence of a mainshock. We find that this effective forecasting period is only in the order of 100 days for moderate mainshocks and in the order of a few years for large events, even if the underlying triggering process lasts much longer.

Description: We present a new upper-mantle tomographic model derived solely from hum seismic data. Phase correlograms between station pairs are computed to extract phase-coherent signals. Correlograms are then stacked using the time–frequency phase-weighted stack method to build-up empirical Green's functions. Group velocities and uncertainties are measured in the wide period band of 30–250 s, following a resampling approach. Less data are required to extract reliable group velocities at short periods than at long periods, and 2 yr of data are necessary to measure reliable group velocities for the entire period band. Group velocities are first regionalized and then inverted versus depth using a simulated annealing method in which the number and shape of splines that describes the S -wave velocity model are variable. The new S -wave velocity tomographic model is well correlated with models derived from earthquakes in most areas, although in India, the Dharwar craton is shallower than in other published models.

Description: Accurate prediction of ground motion intensity and its variability is an important element in seismic hazard assessment. Simulation-based ground motion prediction has become popular in light of earthquake rupture and wave-propagation modelling methods and the rapid growth in computing power. Earthquake rupture modelling needs to be physics-based and also computationally efficient. It also requires the ability to quantify the variability of finite source models for future scenario events. A generalized pseudo-dynamic source model (e.g. earthquake source statistics model) for M w 6.5–7.0 has been constructed by analysing a number of dynamic rupture models based on 1-point and 2-point statistics of kinematic source parameters. Synthetic broad-band ground motions derived from the pseudo-dynamic source model are validated against empirical ground motion prediction equations (GMPEs). The pseudo-dynamic source model produces ground motion intensities mostly compatible with the empirical GMPEs in the broad frequency range. The perturbation analysis of 1-point and 2-point statistics helps to elucidate the effect of source statistics on ground motions in a systematic way. The constructed pseudo-dynamic source model may be used to generate a feasible range of rupture scenarios for future events and to simulate expected ground motions for seismic hazard assessment.

Description: We evaluate the bias in parameter estimates of the ETAS model. We show that when a simulated catalogue is magnitude-truncated there is considerable bias, whereas when it is not truncated there is no discernible bias. We also discuss two further implied assumptions in the ETAS and other self-exciting models. First, that the triggering boundary magnitude is equivalent to the catalogue completeness magnitude. Secondly, the assumption in the Gutenberg–Richter relationship that numbers of events increase exponentially as magnitude decreases. These two assumptions are confounded with the magnitude truncation effect. We discuss the effect of these problems on analyses of real earthquake catalogues.

Description: The frequency dependence of the quality factor Q has long been predicted by mathematical modelling and laboratory measurements; however, in situ evidence from seismic surveys is still lacking. We have conducted the cross-hole seismic surveys to investigate the near-surface seismic attenuation in the Daqing oilfield in northeastern China. The seismic waves were fired in a source hole of 40 m from the bottom to the surface at an interval of 1 m and were recorded in a receiver hole of 40 m by two geophones with one at the surface and the other one at the bottom. The direct waves were extracted to avoid the noise disturbance and the reflection interference, and the attenuations without the effects of the source signature and the receiver coupling were estimated by a method we proposed. The nonlinear attenuations were observed and fitted using the power-law-based Q . The reliability of Q estimate was verified by the high similarity between the real and the simulated attenuations. Therefore, the experiment we have conducted can be treated as a reliable evidence for the frequency dependence of near-surface  Q .

Description: The southeastern European cities of Sofia and Thessaloniki are explored as example site-specific scenarios by geographically zoning their individual localized seismic sources based on the highest probabilities of magnitude exceedance. This is with the aim of determining the major components contributing to each city's seismic hazard. Discrete contributions from the selected input earthquake catalogue are investigated to determine those areas that dominate each city's prevailing seismic hazard with respect to magnitude and source-to-site distance. This work is based on an earthquake catalogue developed and described in a previously published paper by the author and components of a magnitude probability density function. Binned magnitude and distance classes are defined using a joint magnitude–distance distribution. The prevailing seismicity to each city—as defined by a child data set extracted from the parent earthquake catalogue for each city considered—is divided into distinct constrained data bins of small discrete magnitude and source-to-site distance intervals. These are then used to describe seismic hazard in terms of uni-variate modal values; that is, M * and D * which are the modal magnitude and modal source-to-site distance in each city's local historical seismicity. This work highlights that Sofia's dominating seismic hazard—that is, the modal magnitudes possessing the highest probabilities of occurrence—is located in zones confined to two regions at 60–80 km and 170–180 km from this city, for magnitude intervals of 5.75–6.00 M w and 6.00–6.25 M w respectively. Similarly, Thessaloniki appears prone to highest levels of hazard over a wider epicentral distance interval, from 80 to 200 km in the moment magnitude range 6.00–6.25 M w .

Description: Red Clay underlying the loess-palaeosol sequences on the Chinese Loess Plateau is an eolian deposit. There is a controversy over whether magnetic susceptibility ( ) variations in Red Clay sequence can be used as an indicator of summer palaeomonsoon intensity. This study investigates the magnetic mineralogy, magnetic concentration and magnetic grain size distribution of Jiaxian Red Clay with multimagnetic methods. Our results indicate that the magnetic properties of Jiaxian Red Clay are similar to those of the Quaternary loess-palaeosol sequences, and ultrafine ferrimagnetic grains produced during pedogenesis are responsible for an increase in susceptibility, therefore the enhancement mechanism of Red Clay is similar to that of the overlying loess-palaeosol sequences. This paper explores variations in the Red Clay sequence through spatial and temporal analysis. The susceptibility variation of six sites along a NNE to SSW transect correlate to palaeoclimatic cycles, so can be used to trace the summer palaeomonsoon intensity from a spatial perspective. However, a simple loess-derived calibration function cannot be used to quantitative reconstruct the palaeomonsoon intensity variations thought time. An adjusted calibration function for palaeosols from Red Clay sequence needs to be developed, so that can be used to quantitative reconstruct palaeomonsoon intensity. Further study is necessary to develop such a transfer function.

Description: Palaeointensity experiments were carried out to a sample collection from two sections of basalt lava flow sequences of Pliocene age in north central Iceland (Chron C2An) to further refine the knowledge of the behaviour of the palaeomagnetic field. Selection of samples was mainly based on their stability of remanence to thermal demagnetization as well as good reversibility in variations of magnetic susceptibility and saturation magnetization with temperature, which would indicate the presence of magnetite as a product of deuteric oxidation of titanomagnetite. Among 167 lava flows from two sections, 44 flows were selected for the Königsberger–Thellier–Thellier experiment in vacuum. In spite of careful pre-selection of samples, an Arai plot with two linear segments, or a concave-up appearance, was often encountered during the experiments. This non-ideal behaviour was probably caused by an irreversible change in the domain state of the magnetic grains of the pseudo-single-domain (PSD) range. This is assumed because an ideal linear plot was obtained in the second run of the palaeointensity experiment in which a laboratory thermoremanence acquired after the final step of the first run was used as a natural remanence. This experiment was conducted on six selected samples, and no clear difference between the magnetic grains of the experimented and pristine sister samples was found by scanning electron microscope and hysteresis measurements, that is, no occurrence of notable chemical/mineralogical alteration, suggesting that no change in the grain size distribution had occurred. Hence, the two-segment Arai plot was not caused by the reversible multidomain/PSD effect in which the curvature of the Arai plot is dependent on the grain size. Considering that the irreversible change in domain state must have affected data points at not only high temperatures but also low temperatures, f v ≥ 0.5 was adopted as one of the acceptance criteria where f v is a vectorially defined fraction of the linear segment. A measure of curvature k ' was also used to check the linearity of the selected linear segment. It was avoided, however, to reject the result out of hand by the large curvature k of the entire data points because it might still include a linear segment with a large fraction. Combining with the results of Shaw's experiments, 52 palaeointensities were obtained out of 192 specimens, or 11 flow means were obtained out of the 44 lava flows. Most of the palaeointensities were from the upper part of the lava section (Chron C2An.1n) and ranged between 30 and 66 μT. Including two results from the bottom part of the lava section, the mean virtual dipole moment for 2.5–3.5 Ma is 6.3 ± 1.4  x  10 22 Am 2 ( N  = 11), which is ~19 per cent smaller than the present-day dipole moment.

Description: We collect two months of ambient noise data recorded by 35 broad-band seismic stations in a 9 x 11 km area (1–3 km station interval) near Karamay, China, and do cross-correlation of noise data between all station pairs. Array beamforming analysis of the ambient noise data shows that ambient noise sources are unevenly distributed and the most energetic ambient noise mainly comes from azimuths of 40°–70°. As a consequence of the strong directional noise sources, surface wave components of the cross-correlations at 1–5 Hz show clearly azimuthal dependence, and direct dispersion measurements from cross-correlations are strongly biased by the dominant noise energy. This bias renders that the dispersion measurements from cross-correlations do not accurately reflect the interstation velocities of surface waves propagating directly from one station to the other, that is, the cross-correlation functions do not retrieve empirical Green's functions accurately. To correct the bias caused by unevenly distributed noise sources, we adopt an iterative inversion procedure. The iterative inversion procedure, based on plane-wave modeling, includes three steps: (1) surface wave tomography, (2) estimation of ambient noise energy and biases and (3) phase velocities correction. First, we use synthesized data to test the efficiency and stability of the iterative procedure for both homogeneous and heterogeneous media. The testing results show that: (1) the amplitudes of phase velocity bias caused by directional noise sources are significant, reaching ~2 and ~10 per cent for homogeneous and heterogeneous media, respectively; (2) phase velocity bias can be corrected by the iterative inversion procedure and the convergence of inversion depends on the starting phase velocity map and the complexity of the media. By applying the iterative approach to the real data in Karamay, we further show that phase velocity maps converge after 10 iterations and the phase velocity maps obtained using corrected interstation dispersion measurements are more consistent with results from geology surveys than those based on uncorrected data. As ambient noise in high-frequency band (〉1 Hz) is mostly related to human activities or climate events, both of which have strong directivity, the iterative approach demonstrated here helps improve the accuracy and resolution of ANT in imaging shallow earth structures.

Description: While elasticity is a defining characteristic of the Earth's lithosphere, it is often ignored in numerical models of long-term tectonic processes in favour of a simpler viscoplastic description. Here we assess the consequences of this assumption on a well-studied geodynamic problem: the growth of normal faults at an extensional plate boundary. We conduct 2-D numerical simulations of extension in elastoplastic and viscoplastic layers using a finite difference, particle-in-cell numerical approach. Our models simulate a range of faulted layer thicknesses and extension rates, allowing us to quantify the role of elasticity on three key observables: fault-induced topography, fault rotation, and fault life span. In agreement with earlier studies, simulations carried out in elastoplastic layers produce rate-independent lithospheric flexure accompanied by rapid fault rotation and an inverse relationship between fault life span and faulted layer thickness. By contrast, models carried out with a viscoplastic lithosphere produce results that may qualitatively resemble the elastoplastic case, but depend strongly on the product of extension rate and layer viscosity U x L . When this product is high, fault growth initially generates little deformation of the footwall and hanging wall blocks, resulting in unrealistic, rigid block-offset in topography across the fault. This configuration progressively transitions into a regime where topographic decay associated with flexure is fully accommodated within the numerical domain. In addition, high U x L favours the sequential growth of multiple short-offset faults as opposed to a large-offset detachment. We interpret these results by comparing them to an analytical model for the fault-induced flexure of a thin viscous plate. The key to understanding the viscoplastic model results lies in the rate-dependence of the flexural wavelength of a viscous plate, and the strain rate dependence of the force increase associated with footwall and hanging wall bending. This behaviour produces unrealistic deformation patterns that can hinder the geological relevance of long-term rifting models that assume a viscoplastic rheology.

Description: We developed an improved method for the separation of intrinsic and scattering attenuation of seismic shear waves by envelope inversion called Qopen . The method optimizes the fit between Green's functions for the acoustic, isotropic radiative transfer theory and observed energy densities of earthquakes. The inversion allows the determination of scattering and intrinsic attenuation, site corrections and spectral source energies for the investigated frequency bands. Source displacement spectrum and the seismic moment of the analysed events can be estimated from the obtained spectral source energies. We report intrinsic and scattering attenuation coefficients of shear waves near three geothermal reservoirs in Germany for frequencies between 1 and 70 Hz. The geothermal reservoirs are located in Insheim, Landau (both Upper Rhine Graben) and Unterhaching (Molasse basin). We compare these three sedimentary sites to two sites located in crystalline rock with respect to scattering and intrinsic attenuation. The inverse quality factor for intrinsic attenuation is constant in sediments for frequencies smaller than $10\hspace{1.00006pt}\mathrm{Hz}$ and decreasing for higher frequencies. For crystalline rock, it is on a lower level and strictly monotonic decreasing with frequency. Intrinsic attenuation dominates scattering except for the Upper Rhine Graben, where scattering is dominant for frequencies below $10\hspace{1.00006pt}\mathrm{Hz}$ . Observed source displacement spectra show a high-frequency fall-off greater than or equal to 3.

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

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

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

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

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

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

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

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

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

Description: Understanding the contribution of biogenic magnetic particles into sedimentary assemblages is a current challenge in palaeomagnetism. It has been demonstrated recently that magnetic particles produced through biologically controlled mineralization processes, such as magnetosomes from magnetotactic bacteria, contribute to the recording of natural remanent magnetization in marine and lacustrian sediments. Contributions from other, biologically induced, mineralization types, which are known from multiple laboratory experiments to include magnetic minerals, remain largely unknown. Here, we report magnetic properties of iron minerals formed by a community of iron- and manganese-reducing bacteria isolated from a natural groundwater deposit during a 2 yr long incubation experiment. The main iron phases of the biomineralized mass are lepidocrocite, goethite and magnetite, each of which has environmental significance. Unlike the majority of the previous studies that reported superparamagnetic grain size, and thus no remanence carrying capacity of biologically induced magnetite, hysteresis and first-order reversal curves measurements in our study have not detected significant superparamagnetic contribution. The biomineralized mass, instead, contains a mixture of single-domain to pseudo-single-domain and multidomain magnetite particles that are capable of carrying a stable chemical remanent magnetization. Isothermal remanent magnetization acquisition parameters and first-order reversal curves signatures of the biomineralized samples deviate from previously proposed criteria for the distinction of extracellular (biologically induced) magnetic particles in mixtures. Given its potential significance as a carrier of natural remanent magnetization, environmental requirements, distribution in nature and the efficiency in the geomagnetic field recording by biologically induced mineralization need comprehensive investigation.

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

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

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

Description: In the context of the 2014 realization of the International Terrestrial Reference Frame, the International DORIS (Doppler Orbitography Radiopositioning Integrated by Satellite) Service (IDS) has delivered to the IERS a set of 1140 weekly SINEX files including station coordinates and Earth orientation parameters, covering the time period from 1993.0 to 2015.0. From this set of weekly SINEX files, the IDS combination centre estimated a cumulative DORIS position and velocity solution to obtain mean horizontal and vertical motion of 160 stations at 71 DORIS sites. The main objective of this study is to validate the velocities of the DORIS sites by comparison with external models or time-series. Horizontal velocities are compared with two recent global plate models (GEODVEL 2010 and NNR-MORVEL56). Prior to the comparisons, DORIS horizontal velocities were corrected for Global Isostatic Adjustment from the ICE-6G (VM5a) model. For more than half of the sites, the DORIS horizontal velocities differ from the global plate models by less than 2–3 mm yr –1 . For five of the sites (Arequipa, Dionysos/Gavdos, Manila and Santiago) with horizontal velocity differences with respect to these models larger than 10 mm yr –1 , comparisons with GNSS estimates show the veracity of the DORIS motions. Vertical motions from the DORIS cumulative solution are compared with the vertical velocities derived from the latest GPS cumulative solution over the time span 1995.0–2014.0 from the University of La Rochelle solution at 31 co-located DORIS-GPS sites. These two sets of vertical velocities show a correlation coefficient of 0.83. Vertical differences are larger than 2 mm yr –1 at 23 percent of the sites. At Thule, the disagreement is explained by fine-tuned DORIS discontinuities in line with the mass variations of outlet glaciers. Furthermore, the time evolution of the vertical time-series from the DORIS station in Thule show similar trends to the GRACE equivalent water height.

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

Description: We have used ensemble averages of satellite-derived free-air gravity anomaly data, together with inverse modelling techniques, to determine the effective elastic thickness, T e , of circum-Pacific subducting oceanic lithosphere and its relationship to plate age. Synthetic modelling tests show that T e can be recovered best using gravity anomaly, rather than bathymetry, data and profiles that are at least 750 km long. Inverse modelling based on a uniform T e elastic plate suggests that T e increases with age of the subducting oceanic lithosphere and is given approximately by the depth to the 390 ± 10 °C oceanic isotherm based on a cooling plate model. Misfits between the observed and calculated gravity anomalies are significantly improved if a mechanically weak zone is included between the trench axis and the outer rise. This weak zone is coincident with observations of bend-faulting and seismicity. Inverse modelling shows that T e landward of the outer rise is generally 40–65 per cent less than the T e seaward of the outer rise. Both landward and seaward T e increases with age of the lithosphere and are given by the depth to the 342–349 °C and 671–714 °C oceanic isotherm, respectively. A dependence of T e on age is consistent with models for the cooling of oceanic lithosphere as it moves away from a mid-ocean ridge and the temperature-dependent ductile creep of oceanic lithospheric minerals such as olivine. By comparing the observed T e to the predicted T e based on laboratory-derived yield strength envelopes and an assumption of elastic-perfectly plastic deformation, we have attempted to constrain the rheology of oceanic lithosphere. Regardless of the assumed friction coefficient, the dry-olivine low-temperature plasticity flow laws of Goetze, Evans & Goetze, Raterron et al . and Mei et al . all provide quite a good fit to the observed T e at circum-Pacific subduction zones. This result contrasts with the Hawaiian Islands, where these flow laws are generally too strong to fit the observations. The discrepancy in rheology within Pacific plate may be caused by differences in the timescale of loading and therefore the amount of viscoelastic stress relaxation that has occurred. Other possibilities include thermal rejuvenation and magma-assisted flexure at the Hawaiian Islands.