ALBERT

Description: Current earthquake early warning (EEW) systems lack the ability to appropriately handle multiple concurrent earthquakes, which led to many false alarms during the 2011 Tohoku earthquake sequence in Japan. This paper uses a Bayesian probabilistic approach to handle multiple concurrent events for EEW. We implement the theory using a two-step algorithm. First, an efficient approximate Bayesian model class selection scheme is used to estimate the number of concurrent events. Then, the Rao-Blackwellized Importance Sampling method with a sequential proposal probability density function is used to estimate the earthquake parameters, that is hypocentre location, origin time, magnitude and local seismic intensity. A real data example based on 2 months data (2011 March 9–April 30) around the time of the 2011 M 9 Tohoku earthquake is studied to verify the proposed algorithm. Our algorithm results in over 90 per cent reduction in the number of incorrect warnings compared to the existing EEW system operating in Japan.

Description: Surface nuclear magnetic resonance (surface-NMR) is a promising technique for exploring shallow subsurface aquifer structures. Surface-NMR can be applied in environments that are characterized as a 1-D layered Earth. The technique utilizes a single loop and is referred to as magnetic resonance sounding. The technique referred to as magnetic resonance tomography (MRT) allows complex 2-D aquifer structures to be explored. Currently, MRT requires multiple loops and a roll along measurement scheme, which causes long survey time. We propose a loop layout using an elongated transmitter and an in-loop receiver arrays (ETRA) to conduct a 2-D survey with just one measurement. We present a comprehensive comparison between the new layout and the common approaches based on sensitivity and resolution analyses and show synthetic and field data. The results show that ETRA generates subsurface images at sufficient resolution with significantly lower survey times than other loop layouts.

Description: Measurements of ground deformation can be used to identify and interpret geophysical processes occurring at volcanoes. Most studies rely on a single geodetic technique, or fit a geophysical model to the results of multiple geodetic techniques. Here we present a methodology that combines GPS, Total Station measurements and InSAR into a single reference frame to produce an integrated 3-D geodetic velocity surface without any prior geophysical assumptions. The methodology consists of five steps: design of the network, acquisition and processing of the data, spatial integration of the measurements, time series computation and finally the integration of spatial and temporal measurements. The most significant improvements of this method are (1) the reduction of the required field time, (2) the unambiguous detection of outliers, (3) an increased measurement accuracy and (4) the construction of a 3-D geodetic velocity field. We apply this methodology to ongoing motion on Arenal's western flank. Integration of multiple measurement techniques at Arenal volcano revealed a deformation field that is more complex than that described by individual geodetic techniques, yet remains consistent with previous studies. This approach can be applied to volcano monitoring worldwide and has the potential to be extended to incorporate other geodetic techniques and to study transient deformation.

Description: Fractured rocks are known to exhibit seismic anisotropy and shear wave splitting (SWS). SWS is commonly used for fractured rock characterization and has been shown to be sensitive to fluid type. The presence of partial liquid/gas saturation is also known to affect the elastic properties of rocks. The combined effect of both fractures and partial liquid/gas saturation is still unknown. Using synthetic, silica-cemented sandstones with aligned penny-shaped voids, we conducted laboratory ultrasonic experiments to investigate the effect fractures aligned at an oblique angle to wave propagation would have on SWS under partial liquid/gas saturation conditions. The result for the fractured rock shows a saturation dependence which can be explained by combining a fractured rock model and a partial saturation model. At high to full water saturation values, SWS decreases as a result of the fluid bulk modulus effect on the quasi-shear wave. This bulk modulus effect is frequency dependent as a result of wave-induced fluid flow mechanisms, which would in turn lead to frequency dependent SWS. This result suggests the possible use of SWS for discriminating between full liquid saturation and partial liquid/gas saturation.

Description: The effect of network density and geometric distribution on kinematic non-linear source inversion is investigated by inverting synthetic ground motions from a buried strike-slip fault ( M w 6.5), that have been generated by dynamic spontaneous rupture modelling. For the inversion, we use a physics-based regularized Yoffe function as slip velocity function. We test three different cases of station network geometry: (i) single station, varying azimuth and epicentral distance; (ii) multistation circular configurations, that is stations at similar distances from the fault, and regularly spaced around the fault; (iii) irregular multistation configurations using different numbers of stations. Our results show: (1) single station tests suggest that it may be possible to obtain a relatively good source model even using a single station. The best source model using a single station is obtained with stations at which amplitude ratios between three components are not large. We infer that both azimuthal angle and source-to-station distance play an important role in the design of optimal seismic network for source inversion. (2) Multistation tests show that the quality of the inverted source systematically correlates neither with the number of stations, nor with waveform misfit. (3) Waveform misfit has a direct correlation with the number of stations, resulting in overfitting the observed data without any systematic improvement of the source. It suggests that the best source model is not necessarily derived from the model with minimum waveform misfit. (4) A seismic network with a small number of well-spaced stations around the fault may be sufficient to obtain acceptable source inversion.

Description: Fermat's interferometric principle is used to compute interior transmission traveltimes pq from exterior transmission traveltimes sp and sq . Here, the exterior traveltimes are computed for sources s on a boundary B that encloses a volume V of interior points p and q . Once the exterior traveltimes are computed, no further ray tracing is needed to calculate the interior times pq . Therefore this interferometric approach can be more efficient than explicitly computing interior traveltimes pq by ray tracing. Moreover, the memory requirement of the traveltimes is reduced by one dimension, because the boundary B is of one fewer dimension than the volume V . An application of this approach is demonstrated with interbed multiple (IM) elimination. Here, the IMs in the observed data are predicted from the migration image and are subsequently removed by adaptive subtraction. This prediction is enabled by the knowledge of interior transmission traveltimes pq computed according to Fermat's interferometric principle. We denote this principle as the ‘traveltime holographic principle’, by analogy with the holographic principle in cosmology where information in a volume is encoded on the region's boundary.

Description: Earthquake hypocentre locations are crucial in many domains of application (academic and industrial) as seismic event location maps are commonly used to delineate faults or fractures. The interpretation of these maps depends on location accuracy and on the reliability of the associated uncertainties. The largest contribution to location and uncertainty errors is due to the fact that the velocity model errors are usually not correctly taken into account. We propose a new Bayesian formulation that integrates properly the knowledge on the velocity model into the formulation of the probabilistic earthquake location. In this work, the velocity model uncertainties are first estimated with a Bayesian tomography of active shot data. We implement a sampling Monte Carlo type algorithm to generate velocity models distributed according to the posterior distribution. In a second step, we propagate the velocity model uncertainties to the seismic event location in a probabilistic framework. This enables to obtain more reliable hypocentre locations as well as their associated uncertainties accounting for picking and velocity model uncertainties. We illustrate the tomography results and the gain in accuracy of earthquake location for two synthetic examples and one real data case study in the context of induced microseismicity.

Description: An intraplate slow earthquake was detected in northernmost Hokkaido, Japan, by a dense network of the global navigation satellite system. Transient abnormal acceleration of 〈12 mm was observed during the period 2012 July to 2013 January (~5.5 months) at several sites. The spatial displacement distribution suggests that a localized tectonic event caused localized deformation. Estimated fault parameter indicates very shallow-dip reverse faulting in the uppermost crust, with a total seismic moment of 1.75E + 17 N m ( M w 5.4). This fault geometry is probably consistent with detachment structure indicated by geological studies. A simultaneous earthquake swarm with the maximum magnitude M 4.1 suggests a possibility that the slow slip triggered the seismic activity for unknown reasons. This slow earthquake is slower than its moment would indicate, with a duration–magnitude scaling relationship unlike either regular earthquakes or subduction slow slip events. This result indicates that even if the area is under different physical property from subduction zones, slow earthquake can occur by some causes. Slow earthquakes exist in remote regions away from subduction zones and might play an important role in strain release and tectonic activity.

Description: We studied ionospheric responses to the 2012 April 11 M w 8.6 North Sumatra earthquake using total electron content (TEC) measurements with the regional Global Navigation Satellite System network. This earthquake ruptured the oceanic lithosphere off the Indian Ocean coast of North Sumatra, and is known as the largest strike-slip earthquake ever recorded. Coseismic ionospheric disturbances (CIDs) with rapid TEC enhancement of a few TEC units propagated northward with a speed of acoustic waves (~1 km s –1 ). Resonant atmospheric oscillation with a frequency ~4 mHz have been found as monochromatic oscillation of TEC lasting for an hour after the main shock and the largest aftershock. We compared CID amplitudes of 21 earthquakes world-wide with moment magnitudes ( M w ) 6.6–9.2. They roughly obeyed a law such that CID amplitude increases by two orders of magnitude for the M w increase of three. The 2012 North Sumatra earthquakes slightly deviated negatively from the trend possibly reflecting their strike-slip mechanisms, that is small vertical crustal movements for their magnitudes.

Description: In autumn 2012, the new release 05 (RL05) of monthly geopotencial spherical harmonics Stokes coefficients (SC) from Gravity Recovery and Climate Experiment (GRACE) mission was published. This release reduces the noise in high degree and order SC, but they still need to be filtered. One of the most common filtering processing is the combination of decorrelation and Gaussian filters. Both of them are parameters dependent and must be tuned by the users. Previous studies have analyzed the parameters choice for the RL05 GRACE data for oceanic applications, and for RL04 data for global application. This study updates the latter for RL05 data extending the statistics analysis. The choice of the parameters of the decorrelation filter has been optimized to: (1) balance the noise reduction and the geophysical signal attenuation produced by the filtering process; (2) minimize the differences between GRACE and model-based data and (3) maximize the ratio of variability between continents and oceans. The Gaussian filter has been optimized following the latter criteria. Besides, an anisotropic filter, the fan filter, has been analyzed as an alternative to the Gauss filter, producing better statistics.

Description: We present high-resolution tomographic images in source areas of 26 large crustal earthquakes ( M 6.0–7.2) which occurred in Northeast Japan (Tohoku) during the past 120 yr from 1894 to 2014. Prominent low-velocity (low- V ) and high Poisson's ratio (high- ) anomalies are revealed in the crust and mantle wedge under the source areas. Beneath the volcanic front and backarc areas, the low- V and high- zones reflect arc-magma related high-temperature anomalies which are produced by joint effects of corner flow in the mantle wedge and fluids from dehydration of the subducting Pacific slab. The hot anomalies cause locally thinning and weakening of the brittle seismogenic layer above them. Low-frequency micro-earthquakes are observed in the lower crust and uppermost mantle in or around the low- V zones, which reflect ascending of arc magma and fluids from the mantle wedge to the crust. No volcano and magma exist in the forearc area due to low temperature there, hence the low- V zones in the forearc reflect fluids from the slab dehydration. The ascending fluids may have produced a ‘water wall’ in the mantle wedge and crust beneath the forearc area. When the water enters active faults in the crust, the fault-zone friction is reduced and so large earthquakes can be induced. These results indicate that the nucleation of a large earthquake is not entirely a mechanical process, but is closely associated with subduction dynamics and physical and chemical properties of rocks in the crust and upper mantle. In particular, arc magma and fluids play an important role in the seismogenesis.

Description: In spectral-induced polarization (SIP) studies of sites contaminated by organic hydrocarbons, it was shown that biodegradation by-products in general, and organic acids in particular, significantly alter the SIP signature of the subsurface. Still a systematic study that considers the effect of organic acid on the physicochemical and electrical (SIP) properties of the soil is missing. The goal of this work is to relate between the effect of organic acid on the physicochemical properties of the soil, and the soil electrical properties. To do so, we measured the temporal changes of the soil chemical (ion content) and electrical (low-frequency SIP) properties in response to influx of organic acid at different concentrations, gradually altering the soil pH. Our results show that organic acid reduces the soil pH, enhances mineral weathering and consequently reduces both the in-phase and quadrature conductivity. At the pH range where mineral weathering is most significant (pH 6–4.5) a negative linear relation between the soil pH and the soil formation factor was found, suggesting that mineral weathering changes the pore space geometry and hence affecting the in-phase electrical conductivity. In addition, we attribute the reduction in the quadrature conductivity to an exchange process between the natural cation adsorbed on the mineral surface and hydronium, and to changes in the width of the pore bottleneck that results from the mineral weathering. Overall, our results allow a better understanding of the SIP signature of soil undergoing acidification process in general and as biodegradation process in particular.

Description: Oil and gas reservoirs can cause anomalies in the energy and frequency of seismic signals. We can take advantage of these anomalies for hydrocarbon detection. Based on the Teager–Kaiser (TK) energy characteristics, a method that uses this energy in association with the Empirical Mode Decomposition (EMD) method is proposed for hydrocarbon detection in carbonate reservoirs. The EMD method, which can decompose the original seismic signals into a finite number of monocomponent intrinsic mode functions (IMFs) in the temporal domain, is used for multiband filtering. A TK energy separation algorithm is used to estimate the instantaneous frequency and amplitude of the selected IMFs from the EMD method. The proposed method can generate a joint time–frequency representation that can reflect the energy tracking of the seismic signals. The instantaneous spectrums produced by the EMD/TK method have the capability to detect hydrocarbon. The model results to the gas field, which located in the eastern Ordos Basin, China, show that the EMD/TK method can be adopted and they detect the gas-bearing reservoir efficiently. Application of the EMD/TK method in hydrocarbon detection in a gas field located in the Eastern Ordos Basin shows its effectiveness. The EMD/TK method can be used as a new analysis tool to determine the instantaneous spectral properties of a reservoir to detect hydrocarbon.

Description: A novel method for the effective identification of bedrock subsurface elevation from electrical resistivity tomography images is described. Identifying subsurface boundaries in the topographic data can be difficult due to smoothness constraints used in inversion, so a statistical population-based approach is used that extends previous work in calculating isoresistivity surfaces. The analysis framework involves a procedure for guiding a clustering approach based on the fuzzy c -means algorithm. An approximation of resistivity distributions, found using kernel density estimation, was utilized as a means of guiding the cluster centroids used to classify data. A fuzzy method was chosen over hard clustering due to uncertainty in hard edges in the topography data, and a measure of clustering uncertainty was identified based on the reciprocal of cluster membership. The algorithm was validated using a direct comparison of known observed bedrock depths at two 3-D survey sites, using real-time GPS information of exposed bedrock by quarrying on one site, and borehole logs at the other. Results show similarly accurate detection as a leading isosurface estimation method, and the proposed algorithm requires significantly less user input and prior site knowledge. Furthermore, the method is effectively dimension-independent and will scale to data of increased spatial dimensions without a significant effect on the runtime. A discussion on the results by automated versus supervised analysis is also presented.

Description: Significant deformations, followed by increased seismicity detected since 2011 July at El Hierro, Canary Islands, Spain, prompted the deployment of additional monitoring equipment. The climax of this unrest was a submarine eruption first detected on 2011 October 10, and located at about 2 km SW of La Restinga, southernmost village of El Hierro Island. The eruption ceased on 2012 March 5, after the volcanic tremor signals persistently weakened through 2012 February. However, the seismic activity did not end with the eruption, as several other seismic crises followed. The seismic episodes presented a characteristic pattern: over a few days the number and magnitude of seismic event increased persistently, culminating in seismic events severe enough to be felt all over the island. Those crises occurred in 2011 November, 2012 June and September, 2012 December to 2013 January and in 2013 March–April. In all cases the seismic unrest was preceded by significant deformations measured on the island's surface that continued during the whole episode. Analysis of the available GPS and seismic data suggests that several magma displacement processes occurred at depth from the beginning of the unrest. The first main magma movement or ‘injection’ culminated with the 2011 October submarine eruption. A model combining the geometry of the magma injection process and the variations in seismic energy release has allowed successful forecasting of the new-vent opening.

Description: Recent studies show that terrestrial and space-based observations of gravity agree over Europe. In this paper, we compare time-series of terrestrial gravity (including the contribution due to surface displacement) as measured by superconducting gravimeters (SGs), space-based observations from Gravity Recovery and Climate Experiment (GRACE) and predicted changes in gravity derived from two global hydrological models at 10 SG stations in central Europe. Despite the fact that all observations and models observe a maximum in the same season due to water storage changes, there is little agreement between the SG time-series even when they are separated by distances smaller than the spatial resolution of GRACE. We also demonstrate that GRACE and the SG observations and the water storage models do not display significant correlation at seasonal periods nor at interannual periods. These findings are consistent with the fact that the SGs are sensitive primarily to mass changes in the few hundred metres surrounding the station.

Description: Standard processing of seismic events for reporting in bulletins is usually done one-at-a time. State-of-the-art relative event methods, often involving cross correlation, are increasingly used and have improved estimates of event parameters for event detection, location and magnitude. This is because relative event techniques can simultaneously reduce measurement error and effects of model error. We show how cross correlation can be used to assign relative magnitudes for neighbouring seismic events distributed over a large region in east Asia and quantify to what extent the uncertainty in these values increases as waveform similarity breaks down. We find that cross correlation works well for magnitude comparison of two events when it is expected that they generate very similar signals even if these may be almost buried in large amounts of noise. This may be the case when investigating repeating earthquakes or nuclear explosions within a few kilometres of each other. Cross correlation is the optimal detector in these cases assuming noise is white and Gaussian, and also provides the least-squares solution for the relative amplitudes. However, when the waveform similarity of the underlying signals breaks down, due to interevent separation distance, source time function differences or focal mechanism differences, these assumptions are no longer valid and a bias is introduced into the relative magnitude measurement. This bias due to degradation of waveform similarity is modelled here with synthetics and an analytic expression for it is derived based on three terms—the cross-correlation coefficient (CC), and the signal-to-noise ratio (SNR) of the larger and smaller events. The analytic expression is a good match to the observed bias in the data. If the equation for relative magnitude is rewritten to correct for the bias due to the CC, a new equation results which is simply the log of the ratio of the L2 norms. The bias due to SNRs is still present because the observed waveforms inevitably contain both signal and noise. However, this bias is predicted to be minimal for typical detection thresholds. Making measurements of the ratio of the L2 norms is shown to remove the bias due to degradation of waveform similarity for real data. The scatter of these cross-correlation measurements of relative magnitude is much less than those obtained by differencing magnitudes in a traditional catalogue. Of 14 025 events in and near China, 34 per cent had over an order of magnitude reduction in the median standard deviation (0.0342 magnitude units) as compared to the estimated scatter in the catalogue (0.3454 magnitude units). And 78 per cent of the events show a factor 3 improvement or better in the precision of relative event size measured as the ratio of the L2 norms as compared to the precision of the catalogue for relative magnitudes. These results suggest that the ratio of the L2 norms is an appropriate measure of relative magnitudes for general seismicity of a monitoring region, when there is significant waveform dissimilarity for neighbouring events. This measure maintains a higher degree of measurement precision as compared to the catalogue.

Description: Seismic anisotropy in and around subducting Tonga slab (latitude ~20°S) is investigated by using three component broad-band seismograms of deep earthquakes at Tonga ( h 〉 550 km) recorded at the F-net, Japan. In the study area, the slab becomes stagnant when approaching the upper- and lower-mantle boundary, and the mantle transition zone and the shallowest lower mantle have been claimed to be anisotropic both the backarc side and the ocean side.We analyse shear wave splitting of teleseismic direct S waves from the deep earthquakes, and investigate a slightly different part of the Tonga subduction zone from previous studies. We find that anisotropy of an observable degree exists neither in the slab near the bottom of the upper mantle (below 600 km) nor in the lower mantle beneath the foci. The shear wave splitting lag time ( t ) attributable to the anisotropy inside the slab around the foci is less than 0.15 s, and the corresponding maximum degree of anisotropy is 0.9 per cent. The result is consistent with recent mineralogical studies which indicate that ringwoodite does not acquire significant preferred orientation of crystal lattice due to the deformation near the bottom of the upper mantle. The maximum t by large-scale anisotropy in the lower mantle traversed by the S waves from Tonga to Japan does not exceed 0.05 s, suggesting the absence of significant shear deformation near the top of the lower mantle.

Description: Full waveform inversion of ground-penetrating radar data is an emerging technique for the quantitative, high-resolution imaging of the near subsurface. Here, we present a 2-D frequency-domain full waveform inversion for the simultaneous reconstruction of the dielectric permittivity and of the electrical conductivity. The inverse problem is solved with a quasi-Newton optimization scheme, where the influence of the Hessian is approximated by the L-BFGS-B algorithm. This formulation can be considered to be fully multiparameter since it enables to update permittivity and conductivity values within the same descent step, provided we define scales of measurement through a reference permittivity, a reference conductivity, and an additional scaling factor. Numerical experiments on a benchmark from the literature demonstrate that the inversion is very sensitive to the parameter scaling, despite the consideration of the approximated Hessian that should correct for parameter dimensionalities. A proper scaling should respect the natural sensitivity of the misfit function and give priority to the parameter that has the most impact on the data (the permittivity, in our case). We also investigate the behaviour of the inversion with respect to frequency sampling, considering the selected frequencies either simultaneously or sequentially. As the relative imprint of permittivity and conductivity in the data varies with frequency, the simultaneous reconstruction of both parameters takes a significant benefit from broad frequency bandwidth data, so that simultaneous or cumulative strategies should be favoured. We illustrate our scaling approach with a realistic synthetic example for the imaging of a complex subsurface from on-ground multioffset data. Considering data acquired only from the ground surface increases the ill-posedness of the inverse problem and leads to a strong indetermination of the less-constrained conductivity parameters. A Tikhonov regularization can prevent the creation of high-wavenumber artifacts in the conductivity model that compensate for erroneous low-wavenumber structures, thus enabling to select model solutions. We propose a workflow for multiparameter imaging involving both parameter scaling and regularization. Optimal combinations of scaling factors and regularization weights can be identified by seeking regularization levels that exhibit a clear minimum of final data misfit with respect to parameter scaling. We confirm this workflow by inverting noise-contaminated synthetic data. In a surface-to-surface acquisition configuration, we have been able to reconstruct an accurate permittivity structure and a smooth version of the conductivity distribution, based entirely on the analysis of the data misfit with respect to parameter scaling, for different regularization levels.

Description: The aim of this paper is to determine the strain rate tensor (SRT) for the Iranian region. In this study, (1) we apply a method of computation of the SRT never used for the Iranian area and (2) we use a new GPS velocity field obtained from several previously published velocity fields. First, the method is described and tested on a synthetic case, which mimics the real Iranian case. The synthetic tests confirm that the method allows us to both retrieve high gradients of the strain rate field and reduce the effect of an erroneous velocity vector. Second, the method is applied to a real data set covering the Arabia-Eurasia collision zone in Iran. We particularly focus on the Zagros–Makran transition zone, the Central Iran region and the northernmost part of the Arabia-Eurasia collision zone (NW Iran–Caucasus–East Turkey). Whereas the main characteristics of the obtained SRT are consistent with known tectonic features, important new results are found in the Central Iran, with the strike-slip style along the Anar and Deshir faults, and the Zagros–Makran transition zone, with a north–south variation of the SRT along the Zendan–Minab–Palami fault system. We link these results to recent active tectonic studies.

Description: The Qaidam Basin is an ideal archive to study long-term climate and erosion histories at the NE Tibetan Plateau. We present a magnetostratigraphic study of the 723 m deep drill-core SG-1b of lacustrine sediments at the Jianshan anticline in the western Qaidam Basin. The polarity sequence shows 18 normal and 19 reverse polarity zones which can be readily correlated with chrons C1n-C3Br of the Geomagnetic Polarity Time Scale 2004 (GPTS 2004), dating the core at about 7.3–1.6 Ma. The resulting mean sediment accumulation rate (SAR) between polarity boundaries ranges from 6.5 to 30.4 cm ka –1 . High SARs occur within the intervals of 〉7.3–6.0, 5.2–4.2 and 3.6–2.6 Ma indicating three episodic phases of higher erosion. From the derived variation of SARs and previous results, we conclude that growth strata at the Jianshan anticline started to develop at ~1.6 Ma by limb rotation. All this we relate to pulse tectonic uplift of the NE Tibetan Plateau and fault-propagation-folding in the Qaidam Basin.

Description: When the results of geophysical models are compared with data, the uncertainties of the model are typically disregarded. This paper proposes a method for defining the uncertainty of a geophysical model based on a numerical procedure that estimates the empirical auto- and cross-covariances of model-estimated quantities. These empirical values are then fitted by proper covariance functions and used to compute the covariance matrix associated with the model predictions. The method is tested using a geophysical, spherical, thin-sheet finite element model of the Mediterranean region. Using a 2 analysis, the model's estimated horizontal velocities are compared with the velocities estimated from permanent GPS stations while taking into account the model uncertainty through its covariance structure and the covariance of the GPS estimates. The results indicate that including the estimated model covariance in the testing procedure leads to lower observed 2 values and might help a sharper identification of the best-fitting geophysical models.

Description: Optimizing an experimental design is a compromise between maximizing information we get about the target and limiting the cost of the experiment, providing a wide range of constraints. We present a statistical algorithm for experiment design that combines the use of linearized inverse theory and stochastic optimization technique. Linearized inverse theory is used to quantify the quality of one given experiment design while genetic algorithm (GA) enables us to examine a wide range of possible surveys. The particularity of our algorithm is the use of the multi-objective GA NSGA II that searches designs that fit several objective functions (OFs) simultaneously. This ability of NSGA II is helping us in defining an experiment design that focuses on a specified target area. We present a test of our algorithm using a 1-D electrical subsurface structure. The model we use represents a simple but realistic scenario in the context of CO 2 sequestration that motivates this study. Our first synthetic test using a single OF shows that a limited number of well-distributed observations from a chosen design have the potential to resolve the given model. This synthetic test also points out the importance of a well-chosen OF, depending on our target. In order to improve these results, we show how the combination of two OFs using a multi-objective GA enables us to determine an experimental design that maximizes information about the reservoir layer. Finally, we present several tests of our statistical algorithm in more challenging environments by exploring the influence of noise, specific site characteristics or its potential for reservoir monitoring.

Description: Gravity and height changes, reflecting magma accumulation in subsurface chambers, are evaluated using finite element models in order to resolve controversial relationships observed in some volcanic areas. When significant gravity changes occur without any significant deformation, or vice versa, it is often difficult, if not impossible, to explain the observations using the popular Mogi model. Here, we explore whether these discrepancies can be explained by magma compressibility and source geometry effects. Compression of resident magma and expansion of the chamber wall act concurrently to accommodate newly added magma. Gravity-height ratios are found to mainly depend on: (i) geometry of the sources, which control the volume expansion of the chamber, (ii) magma compressibility, which affects the contraction of the magma resident in the chamber, and (iii) depth of the sources. Our numerical results show that, when magma compressibility and non-spherical sources are taken into account, significant gravity variations can, indeed, be successfully reconciled with negligible height changes. This may be the case at Etna volcano, where gravity changes (about 40 μGal) without any significant deformation (below 5 cm) were observed during the 1994–1995 inflation period. The numerical results point to the accumulation of a 1.4 10 10 kg mass into an elongated source simulating a shallow storage region supplying the summit craters.

Description: Understanding the mode of deformation around the Eastern Himalayan Syntaxis (EHS) is crucial for models of Tibetan Plateau evolution. To constrain rotations and translation east of the EHS, we present new palaeomagnetic data from meta-basalt layers—possibly sills—of the Baoshan Block in western Yunnan, southeastern Tibet Plateau; the meta-basalts were either emplaced or metamorphosed to greenschist facies at ~30 Ma ( 40 Ar– 39 Ar whole rock ages). A detailed rock magnetic study—in combination with reflected and transmitted light microscopy—identifies magnetite and Ti-rich titanomagnetite as the main magnetic carriers. Using alternating field demagnetization, we separated well-clustered characteristic remanent magnetizations. A fold test indicates a syn-folding remanence acquisition, supporting an Oligocene magnetization due to evidence for regional Eocene–Oligocene shortening. The tilt corrected overall site mean direction at 30 per cent of unfolding yields declination/inclination = 042.2°/47.0°, corresponding to a clockwise rotation of 35.1° ± 12.7° with respect to stable Eurasia. This denotes an average rotation rate of 1.17 ± 0.42° Myr –1 since ~30 Ma, ranging at the lower limit of the present-day, GPS-derived rotation rates. We explain the clockwise rotation by tectonic escape of the Baoshan and the Lanping-Simao blocks along the Ailao Shan shear zone. With the onset of shearing along the Chong Shan and Gaoligong Shan shear zones, the Baoshan Block continued its southeastward escape, decoupled from the Lanping-Simao Block between these two major shear zones. Crustal flow could explain rotation and southward escape after ~20 Ma.

Description: We use body-waveform modelling to constrain the source parameters of earthquakes occurring globally in oceanic lithosphere beneath the subduction zone outer rise and outer trench slope. These data are then used to map the stress state in the lithosphere of the downgoing plate as it bends into the subduction zone. Our results provide new constraints on the faulting of oceanic lithosphere at the outer rise, which is important for understanding the transmission of plate-driving forces through the subduction system. In all cases, shallow normal-faulting earthquakes are observed at the top of the plate, and are separated in depth from any deeper thrust-faulting earthquakes. No temporal variation associated with large thrust-faulting earthquakes on the subduction interface is seen in the depth extent of each type of faulting at the outer rise. The transition depth from trench-normal extension to compression is found to vary in agreement with models in which deformation is driven by the combination of in-plane stresses and bending stresses, resulting principally from slab pull. Combining the seismologically derived constraints on the thickness of the elastic core of the plate with estimates of the plate curvature, we place upper bounds on the strength of the lithosphere at the outer rise, which is required to be 300 MPa for a constant yield stress model, or governed by an effective coefficient of friction of 0.3.

Description: Geodetic observations across most of the major strike-slip faults show an interseismic strain rate, which presents a sharp localization of elastic shear strain in the fault vicinity (20–60 km). The screw dislocation model of Savage & Burford is commonly used to fit these geodetic data and to retrieve the far field velocity and the locking depth. This model is very popular because of its inherent simplicity to derive fault slip rates, and mostly because it predicts locking depths, which are of the same order of magnitude as the base of the seismogenic zone (5–20 km). A first issue with the screw dislocation model is that localization is paradoxically introduced by imposing a step function in the velocity field at a depth where the crust is otherwise recognized to behave following a viscous rheology. A second issue with this model is that it is not consistent with the rheological model of the crust that is valid for both post-seismic (1–10 yr), interseismic (several thousand years) and long-term geodynamic (several thousand to several million years) timescales. Here we use numerical models to study how alternative and more geologically realistic boundary conditions and rheological structures can lead to the localization of elastic strain at the Earth's surface during the interseismic period. We find that simple elastic models resembling the Savage & Burford model but driven by far field plate velocity are inefficient at localizing strain unless this driving velocity is transmitted by a rigid indentor. We also find that models including a weak viscous heterogeneity beneath the fault zone are able to produce appropriate localization of the deformation near the fault. This alternative class of model is shown to be pertinent in regards to the boundary conditions and geological observations along exhumed ductile strike-slip shear zone.

Description: In this work, we present a study of the coseismic and post-seismic crustal deformation associated to the M w 6.3, 2009 April 6 L'Aquila earthquake from the analysis of GPS displacement time-series. We use a principal component decomposition-based inversion method to study the space- and time-dependent evolution of slip on faults without any a priori assumption on the model used to characterize the temporal evolution of crustal deformation. The method adopted allows us to account for the initial post-seismic deformation in estimating the coseismic displacements, in a consistent manner for the whole GPS network. We use elastic dislocation theory and a least-squares procedure to invert for the slip distribution on the mainshock fault (Paganica fault) and a second fault (Campotosto fault), where a M w 5.2 aftershock occurred on April 9. The geometries for these faults are obtained from a singular value decomposition of precisely relocated aftershocks. We find that the use of complex fault geometries is not justified by the GPS observations available. An inversion that accounts for post-seismic slip to occur on both the Paganica and Campotosto faults provides a better fit to the GPS time-series observations, than using only the Paganica fault segment, at a 95 per cent confidence level. Within our resolution, afterslip regions do not migrate over time and are localized on fault patches that are approximately complementary to those of coseismic slip. We find that the position of some relevant afterslip patches is different if the inversion is performed assuming a fixed rake or not. We estimate the parameter a – b of rate- and state-dependent friction on those fault regions accommodating afterslip that are robustly characterized in our inversions. We find values of the order of 10 –3 , which is near the transition from potentially unstable to nominally stable friction. These results are in agreement with laboratory measurements performed on typical rocks of the L'Aquila region.

Description: The weighted mean temperature ( T m ) is a critical parameter in research of GPS meteorology. Recently, we have established two global T m models, GTm-I and GTm-, using ground-based radiosonde data. These two models, which directly use the location and the day of year to calculate T m , consider the annual cycle only. In this study, we further take into account the variation characteristics of T m in semi-annual and diurnal periodicity and estimate the initial phase of each cycle. A more accurate global empirical T m model called GTm-III is constructed using high-precision Global Geodetic Observing System (GGOS) Atmosphere T m grid data, and then validated using a new set of T m grid data, radiosonde data and Constellation Observation System for Meteorology, Ionosphere, and Climate (COSMIC) radio occultation data. Results indicate that compared to the existing global T m models, GTm-III provides T m estimates of high accuracy on a global scale and the predicted values are close to the true values in every moment of every day.

Description: Increased displacement rates have been observed following manylarge earthquakes and magmatic events. Although an order of magnitude smaller than the displacements associated with the main event, the post-seismic or post-rifting deformation may continue for years to decades after the initial earthquake or dyke intrusion. Due to the rare occurrence of subaerial rifting events, there are very few observations to constrain models of post-rifting deformation. In 2005 September, a 60-km-long dyke was intruded along the Dabbahu segment of the Nubia-Arabia Plate boundary (Afar, Ethiopia), marking the beginning of an ongoing rifting episode. Continued activity has been monitored using satellite radar interferometry and data from global positioning system instruments deployed around the rift in response to the initial intrusion. Using multiple satellite passes, we are able to separate the rift perpendicular and vertical displacement fields around the Dabbahu segment. Rift perpendicular and vertical rates of up to 180 and 240 mm yr –1 , respectively. Here, we show that models of viscoelastic relaxation alone are insufficient to reproduce the observed deformation field and that a large portion of the observed signal is related to the movement of magma within the rift segment. Our models suggest upper mantle viscosities of 10 18–19 Pa s overlain by an elastic crust of between 15 and 30 km. To fit the observations, inflation and deflation of magma chambers in the centre of the rift and to the south east of the rift axis is required at rates of ~0.13 and –0.08 km 3  yr –1 .

Description: Convection in planetary cores can generate fluid flow and magnetic fields, and a number of sophisticated codes exist to simulate the dynamic behaviour of such systems. We report on the first community activity to compare numerical results of computer codes designed to calculate fluid flow within a whole sphere. The flows are incompressible and rapidly rotating and the forcing of the flow is either due to thermal convection or due to moving boundaries. All problems defined have solutions that allow easy comparison, since they are either steady, slowly drifting or perfectly periodic. The first two benchmarks are defined based on uniform internal heating within the sphere under the Boussinesq approximation with boundary conditions that are uniform in temperature and stress-free for the flow. Benchmark 1 is purely hydrodynamic, and has a drifting solution. Benchmark 2 is a magnetohydrodynamic benchmark that can generate oscillatory, purely periodic, flows and magnetic fields. In contrast, Benchmark 3 is a hydrodynamic rotating bubble benchmark using no slip boundary conditions that has a stationary solution. Results from a variety of types of code are reported, including codes that are fully spectral (based on spherical harmonic expansions in angular coordinates and polynomial expansions in radius), mixed spectral and finite difference, finite volume, finite element and also a mixed Fourier–finite element code. There is good agreement between codes. It is found that in Benchmarks 1 and 2, the approximation of a whole sphere problem by a domain that is a spherical shell (a sphere possessing an inner core) does not represent an adequate approximation to the system, since the results differ from whole sphere results.

Description: Large igneous provinces (LIPs) lie approximately above the margins of the African and Pacific large low shear velocity provinces (LLSVPs) in the deep mantle. This spatial correlation has been used to argue that plumes are preferentially generated at the margins of LLSVPs. We perform a series of Monte Carlo–based statistical tests to assess the uniqueness of this conclusion. These tests indicate that (1) the reconstructed locations of LIPs are significantly correlated with both slower-than-average shear wave velocity regions, which contain LLSVPs, and the margins of these structures; and (2) these correlations cannot be statistically distinguished. That is, given current constraints, if plumes were generated randomly throughout regions of slower-than-average shear wave velocity in the deep mantle, then statistical tests are expected to show a significant correlation between the locations of LIPs and the margins of LLSVPs. We therefore conclude that it is premature to argue that the margins of LLSVPs represent preferred zones of plume generation. This conclusion is reinforced in our analysis by a demonstration that the expected mean distance of a set of points randomly placed in slower-than-average shear wave velocity regions is consistent with the observed mean distance between LIPs and the margins of LLSVPs. Finally, we also test the correlation between the reconstructed locations of LIPs and the horizontal gradient in deep mantle shear velocity perturbations. We find, given the uncertainty implied by different tomography models, that there is no statistically significant correlation and that being in a slow region (i.e. in the region of LLSVPs) is a stronger geographic requirement for plume generation than being at a specific (high) gradient.

Description: Zagros is a relatively young and active fold-thrust belt, which has formed due to convergence between the Eurasian and Arabian plates. Magnetotelluric (MT) soundings along a transect were carried out to determine the crustal structure in the collision zone of the two Palaeocontinents. MT data were analysed and modelled using 2-D inversion schemes. The models show clear conductive and resistive domains along the MT profile consistent to a great extent with documented tectonic features and surface geology. The models obtained from the joint inversion of transverse electric and transverse magnetic modes as well as the inversion of the determinant data show similar features along the profile. The new MT results reveal that the transition between two continents at the surface coincides with the western boundary of Sanandaj-Sirjan Zone (SSZ) at the Main Zagros Thrust (MZT). Along the profile towards northeast the conductors at top indicate massive Neogene sediments of the central domain (CD) while the very thick, shallow-located, resistive body (5–25 km thick and 100 km long) beneath is unlikely to be of oceanic affinity, but continental. Another main feature along the profile is the main resistive and conductive parts of the Arabian Plate, which coincide with the tectonic events of High Zagros Fault and Mountain Front Fault. Two highly conductive thick zones are recognized at the southwest part and in the middle of the profile apparently extending to a depth of about 50 km, possibly related to a downward smearing effect due to the presence of thick sedimentary columns in the upper crust. Along the profile, conductive features are recognized at the metamorphic SSZ and Urumieh-Dokhtar Magmatic Assemblage units as well as at CD. Below site 31 along the surface trace of the MZT, the transition between the two continents is distinguished by a complex sequence of conductive and resistive zones both varying laterally as well as vertically. The main difference between the two domains is that the Eurasian Plate seems to be more resistive than the Arabian Plate, although some part of the difference can be related to the thick sequence of conductive sedimentary rocks on the Arabian Plate.

Description: Global Positioning System (GPS) has been proven to be an effective tool to retrieve high-precision displacement for the natural hazard monitoring. The network positioning and Precise Point Positioning (PPP) are the two basic approaches for its data solution, but the former one can only get a relative displacement within the local reference frame and requires a complex and continuously linked infrastructure, and the latter one with a long convergence time to obtain the absolute displacements within the global reference frame. To overcome these drawbacks, this paper proposed a method of fast determining the displacement by PPP velocity estimation (PPPVE). The key of the approach is that the velocity vector parameters are not correlated with other unknown parameters, such as ambiguities and atmosphere, so they can be fast and accurately estimated and integrated into displacements. The validation shows that the displacement can be provided with a precision of 1–2 cm in 1 min by PPPVE. In additional, the Kalman smoothing estimation can be used to improve the PPP solution.

Description: Rifting occurs as episodes of active deformation in individual rift segments of the Northern Volcanic Zone (NVZ) in Iceland. Here, we simulate deformation around the Krafla central volcano and rift system in the NVZ using a 3-D numerical model in order to explain synthetic aperture radar data acquired by the ERS and Envisat satellite missions between 1993 and 2008. The deformation is non-linear in time over the observed interval. The observed deformation can be explained by a combination of three processes, including: (i) secular plate spreading between the North American and Eurasian plates at a rate of 18.2 mm yr –1 , (ii) viscoelastic relaxation following the Krafla Fires rifting episode between 1975 and 1984 and (iii) inflation/deflation of shallow magma chambers beneath the Theistareykir and Krafla central volcanoes. We minimize the misfit between the observed and modelled values of the range change gradient, averaged over all samples, using a simulated annealing algorithm that uses a first-order Taylor series to approximate the fitting function. The calibration parameters include the locking depth of the plate boundary and the rheological properties of the lower crust and mantle. The 68-per cent confidence intervals for the parameters in the solution that best fits the data are: (i) a locking depth of 8.0 to 9.5 km, (ii) a viscosity of 19 to 49 EPa.s (1 EPa.s =10 18 Pa.s) in the lower crust at depths between 8 and 24 km and (iii) a viscosity of 5 to 9 EPa.s in the upper mantle below 24 km.

Description: In this paper, we use synthetic data sets from a profile to demonstrate the importance of aligning the 3-D mesh and data coordinate system with the dominant geo-electrical strike direction in 3-D inverse modelling. The resistivity model investigated consists of a regional, elongated 2-D conductive structure at 45° to the profile. We compare the results of full impedance tensor inversion with the results from inversion of only off-diagonal components of the magnetotelluric impedance tensor. The 3-D inversion result obtained with the complete tensor elements yields the subsurface model closest to the original model, whereas the result of inverting only off-diagonal components is the poor imaging of the continuity of the conductive 2-D body. However, the conductor can be correctly recovered using only the off-diagonal components if the model mesh and the data are aligned with quasi-2-D geo-electrical strike.

Description: Alfvén waves are a set of transverse waves that propagate in an electrically conducting fluid in the presence of an ambient magnetic field. Studies of such waves in the Earth's interior are important because they can be used to make inferences about the structure and physical properties of the core that would otherwise remain inaccessible. We produce 1-D forward models of cylindrical torsional Alfvén waves in the Earth's core, also known as torsional oscillations, and study their evolution in a full sphere and an equatorially symmetric spherical shell. Here, we find that travelling torsional waves undergo significant geometric dispersion that increases with successive reflections from the boundaries such that an initial wave pulse becomes unidentifiable within three transits of the core. Low amplitude wakes trail behind sharply defined pulses during propagation, a phenomenon that we interpret using the failure of Huygens’ principle in even dimensions. We investigate the relationship between geometric dispersion and wavelength, concluding that long-wavelength features are more dispersive than short-wavelength features. This result is particularly important because torsional waves that have been inferred in the Earth's core from secular variation are relatively long wavelength, and are therefore likely to undergo significant dispersion within the core. When stress-free boundary conditions on angular velocity are applied, waves are reflected at the equator of the core–mantle boundary with the same sign as the incident wave. Waves that pass through the rotation axis undergo a pseudo-reflection and display a more complicated behaviour due to a phase shift. In an equatorially symmetric shell, we identify a weak reflection at the tangent cylinder due to geometric effects.

Description: Magnetic experiments, a Monte Carlo simulation and transmission electron microscopy observations combine to confirm variable chemical phase separation during quench and annealing of metastable ferri-ilmenite compositions, caused by inhomogeneous Fe-Ti ordering and anti-ordering. Separation begins near interfaces between growing ordered and anti-ordered domains, the latter becoming progressively enriched in ilmenite component, moving the Ti-impoverished hematite component into Fe-enriched diffusion waves near the interfaces. Even when disordered regions are eliminated, Fe-enriched waves persist and enlarge on anti-phase boundaries between growing and shrinking ordered and anti-ordered domains. Magnetic results and conceptual models show that magnetic ordering with falling T initiates in the Fe-enriched wave crests. Although representing only a tiny fraction of material, identified at highest T s on a field-cooling curve, they control the ‘pre-destiny’ of progressive magnetization at lower T . They can provide a positive magnetic moment in a minority of ordered ferrimagnetic material, which, by exchange coupling, then creates a self-reversed negative moment in the remaining majority. Four T s or T ranges are recognized on typical field-cooling curves: T PD is the T range of ‘pre-destination’; T C is the predominant Curie T where major positive magnetization increases sharply; T MAX is where magnetization reaches a positive maximum, beyond which it is outweighed by self-reversed magnetization and T ZM is the T where total magnetization passes zero. Disposition of these T s on cooling curves indicate the fine structure of self-reversed thermoremanent magnetization. These results confirm much earlier suspicions that the ‘ x -phase’ responsible for self-reversed magnetization resides in Fe-enriched phase boundaries.

Description: Improving the accuracy of the marine gravity field requires both improved altimeter range precision and dense track coverage. After a hiatus of more than 15 yr, a wealth of suitable data is now available from the CryoSat-2, Envisat and Jason-1 satellites. The range precision of these data is significantly improved with respect to the conventional techniques used in operational oceanography by retracking the altimeter waveforms using an algorithm that is optimized for the recovery of the short-wavelength geodetic signal. We caution that this new approach, which provides optimal range precision, may introduce large-scale errors that would be unacceptable for other applications. In addition, CryoSat-2 has a new synthetic aperture radar (SAR) mode that should result in higher range precision. For this new mode we derived a simple, but approximate, analytic model for the shape of the SAR waveform that could be used in an iterative least-squares algorithm for estimating range. For the conventional waveforms, we demonstrate that a two-step retracking algorithm that was originally designed for data from prior missions (ERS-1 and Geosat) also improves precision on all three of the new satellites by about a factor of 1.5. The improved range precision and dense coverage from CryoSat-2, Envisat and Jason-1 should lead to a significant increase in the accuracy of the marine gravity field.

Description: Numerical techniques, as such as finite element method, allow for the inclusion of features, such as topography and/or mechanical heterogeneities, for the interpretation of volcanic deformation. However, models based on these numerical techniques usually are not suitable to be included in non-linear estimations of source parameters based on explorative optimization schemes because they require a calculation of the numerical approach for every evaluation of the misfit function. We present a procedure for finite element (FE) models that can be combined with explorative inversion schemes. The methodology is based on including a body force term representing an infinitesimal source in the model formulation that is responsible for pressure (volume) changes in the medium. This provides significant savings in both the time required for mesh generation and actual computational time of the numerical approach. Furthermore, we develop an inversion algorithm to estimate those parameters that characterize the changes in location and pressure (volume) of deformation sources. Both provide FE inversions in a single step, avoiding remeshing and assembly of the linear system of algebraic equations that define the numerical approach and/or the automatic mesh generation. After providing the theoretical basis for the model, the numerical approach and the algorithm for the inversions, we test the methodology using a synthetic example in a stratovolcano. Our results suggest that the FE inversion methodology can be considered suitable for efficiently save time in quantitative interpretations of volcano deformation.

Description: Probabilistic inversion methods based on Markov chain Monte Carlo (MCMC) simulation are well suited to quantify parameter and model uncertainty of nonlinear inverse problems. Yet, application of such methods to CPU-intensive forward models can be a daunting task, particularly if the parameter space is high dimensional. Here, we present a 2-D pixel-based MCMC inversion of plane-wave electromagnetic (EM) data. Using synthetic data, we investigate how model parameter uncertainty depends on model structure constraints using different norms of the likelihood function and the model constraints, and study the added benefits of joint inversion of EM and electrical resistivity tomography (ERT) data. Our results demonstrate that model structure constraints are necessary to stabilize the MCMC inversion results of a highly discretized model. These constraints decrease model parameter uncertainty and facilitate model interpretation. A drawback is that these constraints may lead to posterior distributions that do not fully include the true underlying model, because some of its features exhibit a low sensitivity to the EM data, and hence are difficult to resolve. This problem can be partly mitigated if the plane-wave EM data is augmented with ERT observations. The hierarchical Bayesian inverse formulation introduced and used herein is able to successfully recover the probabilistic properties of the measurement data errors and a model regularization weight. Application of the proposed inversion methodology to field data from an aquifer demonstrates that the posterior mean model realization is very similar to that derived from a deterministic inversion with similar model constraints.

Description: A robust magnetotelluric (MT) inversion algorithm has been developed on the basis of quantile-quantile (q-q) plotting with confidence band and statistical modelling of inversion residuals for the MT response function (apparent resistivity and phase). Once outliers in the inversion residuals are detected in the q-q plot with the confidence band and the statistical modelling with the Akaike information criterion, they are excluded from the inversion data set and a subsequent inversion is implemented with the culled data set. The exclusion of outliers and the subsequent inversion is repeated until the q-q plot is substantially linear within the confidence band, outliers predicted by the statistical modelling are unchanged from the prior inversion, and the misfit statistic is unchanged at a target level. The robust inversion algorithm was applied to synthetic data generated from a simple 2-D model and observational data from a 2-D transect in southern Africa. Outliers in the synthetic data, which come from extreme values added to the synthetic responses, produced spurious features in inversion models, but were detected by the robust algorithm and excluded to retrieve the true model. An application of the robust inversion algorithm to the field data demonstrates that the method is useful for data clean-up of outliers, which could include model as well as data inconsistency (for example, inability to fit a 2-D model to a 3-D data set), during inversion and for objectively obtaining a robust and optimal model. The present statistical method is available irrespective of the dimensionality of target structures (hence 2-D and 3-D structures) and of isotropy or anisotropy, and can operate as an external process to any inversion algorithm without modifications to the inversion program.

Description: A general approach for constructing numerical equivalents of time-reversal mirrors is introduced. These numerical mirrors can be used to regenerate an original wavefield locally within a confined volume of arbitrary shape. Though time-reversal mirrors were originally designed to reproduce a time-reversed version of an original wavefield, the proposed method is independent of the time direction and can be used to regenerate a wavefield going either forward in time or backward in time. Applications to computational seismology and tomographic imaging of such local wavefield reconstructions are discussed. The key idea of the method is to directly express the source terms constituting the time-reversal mirror by introducing a spatial window function into the wave equation. The method is usable with any numerical method based on the discrete form of the wave equation, for example, with finite difference (FD) methods and with finite/spectral elements methods. The obtained mirrors are perfect in the sense that no additional error is introduced into the reconstructed wavefields apart from rounding errors that are inherent in floating-point computations. They are fully transparent as they do not interact with waves that are not part of the original wavefield and are permeable to these. We establish a link between some hybrid methods introduced in seismology, such as wave-injection, and the proposed time-reversal mirrors. Numerical examples based on FD and spectral elements methods in the acoustic, the elastic and the visco-elastic cases are presented. They demonstrate the accuracy of the method and illustrate some possible applications. An alternative implementation of the time-reversal mirrors based on the discretization of the surface integrals in the representation theorem is also introduced. Though it is out of the scope of the paper, the proposed method also apply to numerical schemes for modelling of other types of waves such as electro-magnetic waves.

Description: Differentiation between natural and induced seismicity is crucial for the ability to safely and soundly carry out various underground experiments and operations. This paper defines an objective tool for one of the criteria used to discriminate between natural and induced seismicity. The qualitative correlation between earthquake rates and the injected volume has been an established tool for investigating the possibility of induced, or triggered, seismicity. We derive mathematically, and verify using numerical examples, that the definition of normalized cross-correlation (NCC) between positive random functions exhibits high values with a limit equal to one, if these functions (such as earthquake rates and injection volumes) have a large mean and low standard deviation. In such a case, the high NCC values do not necessarily imply temporal relationship between the phenomena. Instead of positive-value time histories, the functions with their running mean subtracted should be used for cross-correlation. The NCC of such functions (called here NCCEP) may be close to zero, or may oscillate between positive and negative values in cases where seismicity is not related to injection. We apply this method for case studies of seismicity in Colorado, the United Kingdom, Switzerland and south-central Oklahoma, and show that NCCEP reliably determines induced seismicity. Finally, we introduce a geomechanical model explaining the positive cross-correlation observed in the induced seismicity data sets.

Description: Superconducting gravimeter (SG) records after the 2004 Sumatra M w 9.0, 2010 Maule M w 8.8 and 2011 Tohoku M w 9.1 earthquakes are selected to observe the singlets of six spheroidal normal mode multiplets below 1 mHz ( 0 S 2 , 2 S 1 , 0 S 3 , 0 S 4 , 1 S 2 and 3 S 1 ). To clearly observe their spectral splitting, the ensemble empirical mode decomposition (EEMD) is applied to SG records as a dyadic filter bank. Comparisons of the product spectra obtained from the direct residual gravity records and those obtained after using EEMD clearly confirm the validity of EEMD. After using EEMD, all of the singlets of the six multiplets, particularly those of 0 S 4 and 1 S 2 , are completely observed with high signal-to-noise ratio, whereas some of singlets could not be well resolved without the application of EEMD. This study demonstrates that EEMD may be important in the detection of the splitting of some weak and low-frequency seismic modes. The relevant observation results may improve the Earth's density models.

Description: Using a short-baseline hydrophone array, persistent volcanoacoustic sources are identified within the ambient noise field of the Lau Basin during the period between 2009 January and 2010 April. The submarine volcano West Mata and adjacent volcanic terrains, including the northern Matas and Volcano O, are the most active acoustic sources during the 15-month period of observation. Other areas of long-term activity include the Niua hydrothermal field, the volcanic islands of Hunga Ha'apai, Founalei, Niuatoputapu and Niuafo'ou, two seamounts located along the southern Tofua Arc and at least three unknown sites within the northern Lau Basin. Following the great Samoan earthquake on 2009 September 29, seven of the volcanoacoustic sources identified exhibit increases in the rate of acoustic detection. These changes persist over timescales of days-to-months and are observed up to 900 km from the earthquake hypocentre. At least one of the volcanoacoustic sources that did not respond to the 2009 Samoan earthquake exhibits an increase in detection rate following the great M w 8.8 Chile earthquake that occurred at a distance of ~9500 km on 2010 February 27. These observations suggest that great earthquakes may have undocumented impacts on Earth's vast submarine volcanic systems, potentially increasing the short-term flux of magma and volcanic gas into the overlying ocean.

Description: We perform a detailed synthetic study on the resolution of non-double-couple (non-DC) components in the seismic moment tensors from short-period data observed at regional networks designed typically for monitoring aftershock sequences of large earthquakes. In addition, we test two different inversion approaches—a linear full moment tensor inversion and a non-linear moment tensor inversion constrained to a shear-tensile source model. The inversions are applied to synthetic first-motion P - and S -wave amplitudes, which mimic seismic observations of aftershocks of the 1999 M w  = 7.4 Izmit earthquake in northwestern Turkey adopting a shear-tensile source model. To analyse the resolution capability for the obtained non-DC components inverted, we contaminate synthetic amplitudes with random noise and incorporate realistic uncertainties in the velocity model as well as in the hypocentre locations. We find that the constrained moment tensor inversion yields significantly smaller errors in the non-DC components than the full moment tensor inversion. In particular, the errors in the compensated linear vector dipole (CLVD) component are reduced if the constrained inversion is applied. Furthermore, we show that including the S -wave amplitudes in addition to P -wave amplitudes into the inversion helps to obtain reliable non-DC components. For the studied station configurations, the resolution remains limited due to the lack of stations with epicentral distances less than 15 km. Assuming realistic noise in waveform data and uncertainties in the velocity model, the errors in the non-DC components are as high as ±15 per cent for the isotropic and CLVD components, respectively, thus being non-negligible in most applications. However, the orientation of P - and T -axes is well determined even when errors in the modelling procedure are high.

Description: Laboratory experiments using laser-based ultrasonic techniques can be used to simulate seismic surveys on highly controlled small-scale physical models of the subsurface. Most of the time, such models consist in assemblies of homogeneous and consolidated materials. To enable the physical modelling of unconsolidated, heterogeneous and porous media, the use of granular materials is suggested here. We describe a simple technique to build a two-layer physical model characterized by lateral variations, strong property contrasts and velocity gradients. We use this model to address the efficiency of an innovative surface-wave processing technique developed to retrieve 2-D structures from a limited number of receivers. A step by step inversion procedure of the extracted dispersion curves yields accurate results so that the 2-D structure of the physical model is satisfactorily reconstructed. The velocity gradients within each layer are accurately retrieved as well, confirming current theoretical and experimental studies regarding guided surface acoustic modes in unconsolidated granular media.

Description: On 2010 March 11, a sequence of large, shallow continental crust earthquakes shook central Chile. Two normal faulting events with magnitudes around M w 7.0 and M w 6.9 occurred just 15 min apart, located near the town of Pichilemu. These kinds of large intraplate, inland crustal earthquakes are rare above the Chilean subduction zone, and it is important to better understand their relationship with the 2010 February 27, M w 8.8, Maule earthquake, which ruptured the adjacent megathrust plate boundary. We present a broad seismological analysis of these earthquakes by using both teleseismic and regional data. We compute seismic moment tensors for both events via a W-phase inversion, and test sensitivities to various inversion parameters in order to assess the stability of the solutions. The first event, at 14 hr 39 min GMT, is well constrained, displaying a fault plane with strike of N145°E, and a preferred dip angle of 55°SW, consistent with the trend of aftershock locations and other published results. Teleseismic finite-fault inversions for this event show a large slip zone along the southern part of the fault, correlating well with the reported spatial density of aftershocks. The second earthquake (14 hr 55 min GMT) appears to have ruptured a fault branching southward from the previous ruptured fault, within the hanging wall of the first event. Modelling seismograms at regional to teleseismic distances ( 〉 10°) is quite challenging because the observed seismic wave fields of both events overlap, increasing apparent complexity for the second earthquake. We perform both point- and extended-source inversions at regional and teleseismic distances, assessing model sensitivities resulting from variations in fault orientation, dimension, and hypocentre location. Results show that the focal mechanism for the second event features a steeper dip angle and a strike rotated slightly clockwise with respect to the previous event. This kind of geological fault configuration, with secondary rupture in the hanging wall of a large normal fault, is commonly observed in extensional geological regimes. We propose that both earthquakes form part of a typical normal fault diverging splay, where the secondary fault connects to the main fault at depth. To ascertain more information on the spatial and temporal details of slip for both events, we gathered near-fault seismological and geodetic data. Through forward modelling of near-fault synthetic seismograms we build a kinematic k –2 earthquake source model with spatially distributed slip on the fault that, to first-order, explains both coseismic static displacement GPS vectors and short-period seismometer observations at the closest sites. As expected, the results for the first event agree with the focal mechanism derived from teleseismic modelling, with a magnitude M w 6.97. Similarly, near-fault modelling for the second event suggests rupture along a normal fault, M w 6.90, characterized by a steeper dip angle (dip = 74°) and a strike clockwise rotated (strike = 155°) with respect to the previous event.

Description: The Maule earthquake (2010 February 27, M w 8.8, Chile) broke the subduction megathrust along a previously locked segment. Based on an international aftershock deployment, catalogues of precisely located aftershocks have become available. Using 23 well-located aftershocks, we calibrate the classic teleseismic backprojection procedure to map the high-frequency seismic radiation emitted during the earthquake. The calibration corrects traveltimes in a standard earth model both with a static term specific to each station, and a ‘dynamic’ term specific to each combination of grid point and station. The second term has been interpolated over the whole slipping area by kriging, and is about an order of magnitude smaller than the static term. This procedure ensures that the teleseismic images of rupture development are properly located with respect to aftershocks recorded with local networks and does not depend on accurate hypocentre location of the main shock. We track a bilateral rupture propagation lasting ~160 s, with its dominant branch rupturing northeastwards at about 3 km s –1 . The area of maximum energy emission is offset from the maximum coseismic slip but matches the zone where most plate interface aftershocks occur. Along dip, energy is preferentially released from two disconnected interface belts, and a distinct jump from the shallower belt to the deeper one is visible after about 20 s from the onset. However, both belts keep on being active until the end of the rupture. These belts approximately match the position of the interface aftershocks, which are split into two clusters of events at different depths, thus suggesting the existence of a repeated transition from stick-slip to creeping frictional regime.

Description: Secular motion of Earth's rotation pole results in large-scale secular deformation of Earth. Here, we investigate the magnitude of the deformation that has resulted from the rapid motion of the rotation pole to the east since ~2005. We show that geodetic (GNSS, DORIS, VLBI and SLR) estimates of vertical velocity since ~2005 have been biased by up to ±0.38 mm yr –1 relative to the longer-term deformation pattern. The largest signals occur within regions that include the U.S. Pacific Coast, Europe and South Pacific islands where geodetic measurements provide essential measurements of tide-gauge vertical movement and important constraints on models of glacial isostatic adjustment. Consequently, geodetic vertical velocities based on recent data should not be interpreted as being identical to centennial or longer term vertical land movement. Since 2010 the effect is further amplified by the overprediction of the IERS polar motion model relative to the ongoing secular change in pole position—during this time geodetic vertical velocities based on the IERS pole tide model are not just biased relative to the long-term rates but also from actual post-2010 Earth deformation. For geophysical or reference frame studies seeking geodetic vertical velocities that are representative of decadal timescales, where interannual variation is considered noise, the correction for this non-linear effect is straightforward, requiring an elastic computation using a reference rate of polar motion that is linear over the timescales of interest.

Description: Broad-band magnetic susceptibility (MS) measurement, a novel magnetic method capable of quantifying a narrow grain size distribution (GSD) of superparamagnetic (SP) particles by measuring low-field MS at a number of frequency steps spanning four orders of magnitude, has been tested in a loess/palaeosol section at Luochuan in the Chinese Loess Plateau. The studied succession consists of sequences from the latest palaeosol unit (S0) to the upper part of the loess unit (L2), spanning the last glacial–interglacial cycle. Reconstructed GSDs consist of volume fractions on the order of 10 –24 m 3 , and the mean GSDs are modal but with distinctive skewness among the loess, the weakly developed palaeosols (weak palaeosols), and the mature palaeosols. This indicates that the mean volume of SP particles in this loess/palaeosol sequence tends to increase during the transition from loess -〉 weak palaeosol -〉 palaeosol, an indication of grain growth as pedogenesis progresses. Total frequency dependence, or TFD(per cent), the difference between 130 at the lowest (130 Hz) and 500k at the highest (500 kHz) frequencies normalized to 130 , is judged to be a more suitable index than previous frequency dependence parameters for the concentration of SP particles. TFD(per cent) has a strong correlation with 130 , showing a continuous ‘growth curve’ with the rate of increase being highest for the loess, moderate for the weak palaeosols, and saturated for the palaeosols. The characteristic curve suggests that smaller SP particles are preferentially formed in the earlier stage of pedogenesis rather than the later phase when even larger particles are formed in mature palaeosols. These results demonstrate that the broad-band MS measurement method will be useful for the quantitative assessment of magnetic nanoparticles in soils and sediments.

Description: Environmental magnetism can provide simple and fast methods to semi-quantify the contamination levels by heavy metals in urban areas from the relationships between the magnetic properties and the heavy metal concentrations. The aim of this study is to explore the interpretation of selected magnetic parameters as proxies for the evaluation of heavy metal contamination in urban dusts. Dust samples were collected from different districts of Nanjing, Southeast China. The magnetic properties of dusts were analysed and heavy metal (Fe, Al, As, Ba, Cd, Cr, Co, Cu, Hg, Mn, Mo, Ni, Pb, Sr, V and Zn) concentrations were measured. Magnetic minerals, mainly from anthropogenic activities, are dominated by coarse stable single domain-grained ferrimagnetic minerals. The – T curves indicated that the main magnetic minerals are magnetite and hematite. Magnetic properties exhibited similar spatial distribution and significant positive correlations with the concentrations of most metals and with the pollution load index. Magnetic properties had a strong link with the abundance of heavy metals derived from industrial emission and traffic activities, whereas there was a poor correlation with the concentrations of metals from multi-sources such as commercial/domestic sources and natural processes. Besides the low-frequency magnetic susceptibility ( LF ) and the saturation isothermal remanent magnetization, other magnetic parameters can be considered as efficient indicators in the assessment of heavy metal contamination. The linkage between the magnetic properties and heavy metal concentrations in street dusts depends on the sources and nature of their magnetic fraction. Environmental magnetism analyses may provide simple and rapid methods for the assessment of heavy metal contamination in urban street dusts.

Description: Geophysical and geological data suggest that Tibetan middle crust is a partially molten, mechanically weak layer, but it is debated whether this low-viscosity layer is present beneath the entire plateau, what its properties are, how it deforms, and what role it has played in the plateau's evolution. Broad-band seismic surface waves yield resolution in the entire depth range of the Tibetan crust and can be used to constrain its shear-wave velocity structure (indicative of crustal composition, temperature and partial melting) and radial anisotropy (indicative of the patterns of deformation). We measured Love- and Rayleigh-wave phase-velocity curves in broad period ranges (up to 7–200 s) for a few tens of pairs and groups of stations across Tibet, combining, in each case, hundreds of interstation measurements, made with cross-correlation and waveform-inversion methods. Shear-velocity profiles were then determined by extensive series of non-linear inversions of the data, designed to constrain the depth-dependent ranges of isotropic-average shear speeds and radial anisotropy. Shear wave speeds within the Tibetan middle crust are anomalously low and, also, show strong lateral variations across the plateau. The lowest mid-crustal shear speeds are found in the north and west of the plateau (~3.1–3.2 km s –1 ), within a pronounced low-velocity zone. In southeastern Tibet, crustal shear wave speeds increase gradually towards southeast, whereas in the north, the change across the Kunlun Fault is relatively sharp. The lateral variations of shear speeds within the crust are indicative of those in temperature. A mid-crustal temperature of 800 °C, reported previously, can account for the low shear velocities across Lhasa. In the north, the temperature is higher and exceeds the solidus, resulting in partial melting that we estimate at 3–6 per cent. Strong radial anisotropy is required by the data in western-central Tibet (〉5 per cent) but not in northeastern Tibet. The amplitude of radial anisotropy in the crust does not correlate with isotropic-average shear speed (and, by inference, with crustal rock viscosity) or with surface elevation. Instead, radial anisotropy is related to the deformation pattern and is the strongest in regions experiencing extension (crustal flattening), as noted previously. The growth of Tibet by the addition of Indian crustal rocks into its crust from the south is reflected in the higher crustal seismic velocities (and, thus, lower temperatures) in the southern compared to northern parts of the plateau (more recently added rocks having had less time to undergo radioactive heating within the thickened Tibetan crust). Gravity-driven flattening—the basic cause of extension and normal faulting in the southern, western and central Tibet—is evidenced by pervasive radial anisotropy in the middle crust beneath the regions undergoing extension; the overall eastward flow of the crust is directed by the boundaries and motions of the lithospheric blocks that surround Tibet.

Description: For a period of about 1 yr between the summers of 2010 and 2011, 25 broad-band seismographs were deployed in a roughly linear array across the eastern end of the Qaidam basin and the Qilian Shan in the northeastern Tibetan plateau. This region is probably the most suitable place to study the ongoing convergence interaction between the high Tibetan plateau and the main Asian continental plate. Low-frequency P receiver function analysis of the data provides an image of the crust and mantle down to 700 km depth. In addition to the Moho at 45–65 km depth beneath the profile, the 410 and 660 km discontinuities bounding the mantle transition zone can be identified at 400–410 and 650–660 km depths, respectively. A possible increase in temperature in the upper mantle thought to exist beneath the northern part of the high Tibetan plateau is thus confined to this part of the plateau and lower upper-mantle temperatures similar to those beneath southern Tibet occur beneath the Qaidam basin and Qilian Shan. When higher frequencies are included in the P receiver function analysis, a positive Ps converter dipping down to the south from 70–75 km depth at 37.9°N to about 110 km depth at 36°N is imaged. As this feature is only seen in high-frequency images and not in the low-frequency image, it is modelled as the positive Ps conversion from the base of an approximately 5-km-thick anisotropic layer at the top of the Asian mantle lithosphere which is currently subducting. This south-dipping converter continues to the south on the INDEPTH IV profile. S receiver function analysis completes the image of the structure below the Qilian Shan profile with the identification of the lithosphere–asthenosphere boundary (LAB). The LAB of the Asian Plate is identified for a reference slowness of 6.4 s deg –1 at 12–14 s (105–125 km depth) between 38 and 41°N below the northern part of the S receiver function profile. To the south it increases in depth such that it is at about 19 s (170 km depth) between 34 and 35°N at the southern end of the profile. The LAB of the Asian Plate occurs at similar depths on the INDEPTH IV profile at the latitudes where the INDEPTH IV and Qilian Shan profiles overlap. As on the INDEPTH IV profile to the south, between 34 and 35°N at the southern end of the Qilian Shan profile there is evidence from the S receiver functions for the LAB of a separate Tibetan Plate.

Description: The aseismic Murray ridge (MR) is a continuation of the Owen fracture zone which marks the boundary between the Indian and Arabian plates. Due to large variation in morphology and structure within this NE–SW trending ridge in the Arabian Sea a large variation of the bathymetry from few hundred metres to about 4000 m is seen. Observed seismicity on the ridge system is predominantly strike-slip. Tsunamis generated in the Makran subduction zone (MSZ) will propagate through the MR system due to its proximity. As the tsunami speed depends on the depth of the ocean, bathymetry plays a vital role on tsunami propagation. In this paper, the effect of tsunami propagation through the MR system is carried out with the existing bathymetry and comparing the results by removing the bathymetry. To study this phenomenon the 1945 Makran tsunami and worst possible tsunamigenic earthquakes form eastern and western MSZ are considered. The directivity of tsunami propagation with the ridge system is seen to change after crossing the MR towards the southeast direction for tsunamis generated in the eastern MSZ. For tsunami generated in the western MSZ there is no change in its directivity and is almost same as without the ridge with propagation being towards the open sea. Hence the MR not only affects the amplitude of the tsunamis but also the directionality and the arrival times.

Description: The effects of Earth's layered structure, gravity and curvature on coseismic deformation are systematically quantified for all fundamental point sources and some finite-fault sources, respectively. The point-source simulations show that the layering effect (about ≤25 per cent) is significantly higher than the gravity effect (about ≤11 per cent) and the curvature effect (about ≤5 per cent). A case study on the 2011 M w 9.0 Tohoku earthquake is made to quantify the uncertainties of the dislocation models of large earthquakes, in which the three different effects are neglected. Finally, it is investigated how geodetic finite-fault slip inversions are affected by neglecting the layering effect.

Description: Recent work on the Marchenko equation has shown that the scalar 3-D Green's function for a virtual source in the subsurface can be retrieved from the single-sided reflection response at the surface and an estimate of the direct arrival. Here, we discuss the first steps towards extending this result to multicomponent data. After introducing a unified multicomponent 3-D Green's function representation, we analyse its 1-D version for elastodynamic waves in more detail. It follows that the main additional requirement is that the multicomponent direct arrival, needed to initiate the iterative solution of the Marchenko equation, includes the forward-scattered field. Under this and other conditions, the multicomponent Green's function can be retrieved from single-sided reflection data, and this is demonstrated with a 1-D numerical example.

Description: Multichannel analysis of surface waves (MASW) method analyses high-frequency surface waves to determine shear ( S )-wave velocities of near-surface materials, which are usually unconsolidated and possess higher Poisson's ratios. One of key steps using the MASW method to obtain the near-surface S -wave velocities is to pick correct phase velocities in dispersive images. A high-frequency seismic survey conducted over near-surface materials with a higher Poisson's ratio will often result in data that contains non-geometric wave, which will raise an additional energy in the dispersion image. Failure to identify it may result in misidentification. In this paper, we have presented a description about leaky surface wave and the influence caused by the existence of leaky waves in a high-frequency seismic record. We first introduce leaky wave and non-geometric wave. Next, we use two synthetic tests to demonstrate that non-geometric wave is leaky wave and show the properties about leaky surface wave by eigenfunctions using Chen's algorithm. We show that misidentification may occur in picking the dispersion curves of normal Rayleigh wave modes because the leaky-wave energy normally connects energy of fundamental and/or higher modes. Meanwhile, we use a real-world example to demonstrate the influence of leaky wave. We also propose that muting and filtering should been applied to raw seismic records prior to generating dispersive images to prevent misidentifying leaky surface waves as modal surface waves by a real-world example. Finally, we use a three-layer model with a low-velocity half-space to illustrate that leaky surface waves appear on condition that the phase velocities are higher than maximum S -wave velocity of the earth model when solving the Rayleigh equation.

Description: We present in this paper an alternative approach to modelling and inversion of surface nuclear magnetic resonance (sNMR) data that has numerous attractive characteristics. By considering the forward and inverse problem in the frequency-domain (FD) it is possible to heavily truncate the band-limited data set with no appreciable loss of information. Furthermore, it introduces a natural method of signal demodulating and sidesteps the rotating frame transformation, which can be difficult to apply to large volumes. By instead remaining in the laboratory frame, phase shifts due to off-resonance transmission and static variation in the Larmor frequencies become simpler to account for. For these reasons, we present the first practical scheme for comprehensively inverting the complex sNMR data set. The use of complex data allows for improved depth and decay constant resolution. However, it requires the knowledge of subsurface electrical conductivity, and by corollary the ability to accurately model those effects in the sNMR record. Finally, the complex FD signal is sensitive to dephasing of the sNMR record due to static magnetic field inhomogeneity. We illustrate the ability to detect dephasing in synthetic and field data and demonstrate an effective inversion which is able to solve for static dephasing effects. This results in better estimation of meaningful decay constants as well as reduced error in extrapolated initial signal amplitude. Both of these parameters are critical for using sNMR in applications requiring estimation of hydraulic parameters.

Description: Recent studies have shown that many early aftershocks are missing in the regular earthquake catalogues immediately following a moderate to large earthquake. An effective way to recover those missing aftershocks is to use existing events as templates to scan through continuous waveforms to detect additional events, also known as the matched filter technique. Here, we use this technique to detect missing early aftershocks within 3 d after the 2010 M L 6.4 Jiashian earthquake in southern Taiwan. Using waveforms of 875 aftershocks listed in the Central Weather Bureau catalogue as templates, we have detected 2307 early aftershocks in the first 3 d following the main shock. The aftershock rate generally follows an Omori's type decay with a logarithmic slope p of ~0.8. In addition, aftershocks migrated along the northwestward direction with a velocity of ~6–7 km decade –1 and ~3–4 km decade –1 in the opposite direction in depth of 0–20 km. Regions with different seismic properties in the vicinity of the main shock hypocentre may affect the locations and migration patterns of the early aftershocks. We also found a spatial anticorrelation between the main shock slip and aftershock hypocentres in first half hour, although aftershocks at later times filled in regions with moderate amount of slips during the main shock.

Description: This note points out that the statement in ‘A new moment-tensor decomposition for seismic events in anisotropic media’ by Chapman and Leaney that ‘this interpretation is always possible for any moment tensor whatever the anisotropy’ needs qualification. The statement is true for isotropic and weakly anisotropic media, but the decomposition is not possible for all anisotropies, and for some anisotropies and moment tensors is possible but non-unique. However, it is possible for most expected anisotropies. Conditions when the statement is invalid are discussed, and examples of realistic anisotropy are given for which the statement is valid.

Description: The paper in question by Van Camp and co-authors [MVC] challenges previous work showing that ground gravity data arising from hydrology can provide a consistent signal for the comparison with satellite gravity data. The data sets used are similar to those used previously, that is, the gravity field as measured by the GRACE satellites versus ground-based data from superconducting gravimeters (SGs) over the same continental area, in this case Central Europe. One of the main impediments in this paper is the presentation that is frequently confusing and misleading as to what the data analysis really shows, for example, the irregular treatment of annual components that are first subtracted then reappear in the analysis. More importantly, we disagree on specific points. Two calculations are included in our comment to illustrate where we believe that the processing in [MVC] paper is deficient. The first deals with their erroneous treatment of the global hydrology using a truncated spherical harmonic approach which explains almost a factor 2 error in their computation of the loading. The second shows the effect of making the wrong assumption in the GRACE/hydrology/surface gravity comparison by inverting the whole of the hydrology loading for underground stations. We also challenge their claims that empirical orthogonal function techniques cannot be done in the presence of periodic components, and that SG data cannot be corrected for comparisons with GRACE data. The main conclusion of their paper, that there is little coherence between ground gravity stations and this invalidates GRACE comparisons, is therefore questionable. There is nothing in [MVC] that contradicts any of the previous papers that have shown clearly a strong relation between seasonal signals obtained from both ground gravity and GRACE satellite data.

Description: We present a 3-D model of P and S velocities beneath El Hierro Island, constructed using the traveltime data of more than 13 000 local earthquakes recorded by the Instituto Geográfico Nacional (IGN, Spain) in the period from 2011 July to 2012 September. The velocity models were performed using the LOTOS code for iterative passive source tomography. The results of inversion were thoroughly verified using different resolution and robustness tests. The results reveal that the majority of the onshore area of El Hierro is associated with a high-velocity anomaly observed down to 10–12-km depth. This anomaly is interpreted as the accumulation of solid igneous rocks erupted during the last 1 Myr and intrusive magmatic bodies. Below this high-velocity pattern, we observe a low-velocity anomaly, interpreted as a batch of magma coming from the mantle located beneath El Hierro. The boundary between the low- and high-velocity anomalies is marked by a prominent seismicity cluster, thought to represent anomalous stresses due to the interaction of the batch of magma with crust material. The areas of recent eruptions, Orchilla and La Restinga, are associated with low-velocity anomalies surrounding the main high-velocity block. These eruptions took place around the island where the crust is much weaker than the onshore area and where the melted material cannot penetrate. These results put constraints on the geological model that could explain the origin of the volcanism in oceanic islands, such as in the Canaries, which is not yet clearly understood.

Description: The b value of the Gutenberg–Richter (GR) distribution is estimated as a function of a threshold magnitude m th and it is found to depend on m th for magnitudes larger than the completeness magnitude m c . We identify a magnitude interval [ m c , m m ] where b is a decreasing function of m th followed by a regime of increasing b for large magnitudes. This is a common feature of experimental catalogues for different geographic areas. The increase at large m th is explained in terms of an upper magnitude cut-off in experimental catalogues due to finite size effects. We develop a rigorous mathematical framework to relate the decrease of b in the intermediate regime to the functional form of the distribution of the b values. We propose two hypotheses: The first is that the spatial and temporal variability of b leads to a b distribution peaked around its average value. The second is that main shocks and aftershocks are distributed according to the GR law with different b values, leading to a bimodal distribution of b . Simulated Epidemic Type Aftershock Sequences catalogues, generated according to this hypothesis, exhibit the same magnitude distribution of experimental ones. In alternative, we cannot exclude the b dependence on m caused by magnitudes not homogeneously evaluated in a seismic catalogue. In the latter scenario our results provide the correction terms to the estimated magnitudes.

Description: The influence of changes in surface ice-mass redistribution and associated viscoelastic response of the Earth, known as glacial isostatic adjustment (GIA), on the Earth's rotational dynamics has long been known. Equally important is the effect of the changes in the rotational dynamics on the viscoelastic deformation of the Earth. This signal, known as the rotational feedback, or more precisely, the rotational feedback on the sea level equation, has been mathematically described by the sea level equation extended for the term that is proportional to perturbation in the centrifugal potential and the second-degree tidal Love number. The perturbation in the centrifugal force due to changes in the Earth's rotational dynamics enters not only into the sea level equation, but also into the conservation law of linear momentum such that the internal viscoelastic force, the perturbation in the gravitational force and the perturbation in the centrifugal force are in balance. Adding the centrifugal-force perturbation to the linear-momentum balance creates an additional rotational feedback on the viscoelastic deformations of the Earth. We term this feedback mechanism, which is studied in this paper, as the rotational feedback on the linear-momentum balance. We extend both the time-domain method for modelling the GIA response of laterally heterogeneous earth models developed by Martinec and the traditional Laplace-domain method for modelling the GIA-induced rotational response to surface loading by considering the rotational feedback on linear-momentum balance. The correctness of the mathematical extensions of the methods is validated numerically by comparing the polar-motion response to the GIA process and the rotationally induced degree 2 and order 1 spherical harmonic component of the surface vertical displacement and gravity field. We present the difference between the case where the rotational feedback on linear-momentum balance is considered against that where it is not. Numerical simulations show that the resulting difference in radial displacement and sea level change between these situations since the Last Glacial Maximum reaches values of ±25 and ±1.8 m, respectively. Furthermore, the surface deformation pattern is modified by up to 10 per cent in areas of former or ongoing glaciation, but by up to 50 per cent at the bottom of the southern Indian ocean. This also results in the movement of coastlines during the last deglaciation to differ between the two cases due to the difference in the ocean loading, which is seen for instance in the area around Hudson Bay, Canada and along the Chinese, Australian or Argentinian coastlines.

Description: Passive interferometry technology is based on the relation between the reflection and the transmission responses of the subsurface. The transmission response can be received at surface in the presence of the ambient noise source in the subsurface with the cross-correlation (CC) or multidimensional deconvolution methods. We investigate the feasibility of electromagnetic (EM) wave passive interferometry with CC method. We design a 2-D finite-difference time domain (FDTD) algorithm to simulate the long-duration ground penetrating radar (GPR) measurements with random distribution of passive EM sources. The noise sources have random duration time, waveform and spatial distribution. We test the FDTD GPR passive interferometry code with above source characteristics and apply the method to light non-aqueous phase liquid (LNAPL) monitoring. Based on the model simulation data, by using common midpoint velocity analysis and normal move out correction to process the interferometry retrieve record, we can accurately obtain the dynamic changing characteristics of the target's permittivity. The LNAPL dynamic leakage model can be imaged as well. The synthetic results demonstrate that the GPR passive interferometry is feasible in subsurface LNAPL monitoring. Our work provides a foundation for a passive interferometry field application using GPR.

Description: Lower and upper bounds for present deformation rates across faults in central California between the San Andreas Fault and Pacific coast are estimated from a new Global Positioning System (GPS) velocity field for central, western California in light of geodetic evidence presented in a companion paper for slow, but significant deformation within the Pacific Plate between young seafloor in the eastern Pacific and older seafloor elsewhere on the plate. Transects of the GPS velocity field across the San Andreas Fault between Parkfield and San Juan Buatista, where fault slip is dominated by creep and the velocity field thus reveals the off-fault deformation, show that GPS sites in westernmost California move approximately parallel to the fault at an average rate of 3.4 ± 0.4 mm yr –1 relative to the older interior of the Pacific Plate, but only 1.8 ± 0.6 mm yr –1 if the Pacific Plate frame of reference is corrected for deformation within the plate. Modelled interseismic elastic deformation from the weakly coupled creeping segment of the San Andreas Fault is an order-of-magnitude too small to explain the southeastward motions of coastal sites in western California. Similarly, models that maximize residual viscoelastic deformation from the 1857 Fort Tejon and 1906 San Francisco earthquakes mismatch both the rates and directions of GPS site motions in central California relative to the Pacific Plate. Neither thus explains the site motions southwest of the San Andreas fault, indicating that the site motions measure deformation across faults and folds outboard of the San Andreas Fault. The non-zero site velocities thus constitute strong evidence for active folding and faulting outboard from the creeping segment of the San Andreas Fault and suggest limits of 0–2 mm yr –1 for the Rinconada Fault slip rate and 1.8 ± 0.6 to 3.4 ± 0.4 mm yr –1 for the slip rates integrated across near-coastal faults such as the Hosgri, San Gregorio and San Simeon faults.

Description: Continental China consists of a complex amalgamation of geotectonic units and has experienced strong and widespread tectonic deformation since the Mesozoic. To understand its geological structure better, we conducted a systematic receiver function analysis using a total of 83 509 teleseismic traces in the time period of 2009–2010 recorded by 798 broad-band stations, among which 749 stations are permanent digital seismic stations from China Earthquake Networks Center and 49 stations were temporarily deployed in northern Central Tibet. A standard H – stacking method is employed to determine Moho depth and V p / V s ratio underneath each station from teleseismic receiver function analysis. The obtained Moho depth variations are generally consistent with those determined from various deep seismic soundings profiles. We combine our results with those from previous receiver functions studies to produce a high-resolution map of Moho depth and V p / V s variation for continental China. Compared to previous studies, the new study concerns many more stations and the resulting Moho depth map has much higher lateral resolution, especially in the eastern China. Overall, the Moho depth variation has a remarkable correlation with major tectonic units in continental China. For example, across the well-known gravity gradient line in east China, there is a clear shift in Moho depths. In general, the map of V p / V s ratio shows relatively high anomalies underneath Tibetan Plateau, along the gravity gradient line, and under several volcanoes.

Description: 3-D electrical resistivity surveys and inversion models are required to accurately resolve structures in areas with very complex geology where 2-D models might suffer from artefacts. Many 3-D surveys use a grid where the number of electrodes along one direction ( x ) is much greater than in the perpendicular direction ( y ). Frequently, due to limitations in the number of independent electrodes in the multi-electrode system, the surveys use a roll-along system with a small number of parallel survey lines aligned along the x -direction. The ‘Compare R’ array optimization method previously used for 2-D surveys is adapted for such 3-D surveys. Offset versions of the inline arrays used in 2-D surveys are included in the number of possible arrays (the comprehensive data set) to improve the sensitivity to structures in between the lines. The array geometric factor and its relative error are used to filter out potentially unstable arrays in the construction of the comprehensive data set. Comparisons of the conventional (consisting of dipole-dipole and Wenner–Schlumberger arrays) and optimized arrays are made using a synthetic model and experimental measurements in a tank. The tests show that structures located between the lines are better resolved with the optimized arrays. The optimized arrays also have significantly better depth resolution compared to the conventional arrays.

Description: In this study, information carried by the ambient seismic field is exploited to extract impulse response functions between two seismic stations using one as a ‘virtual’ source. Interferometry by deconvolution method is used and validated by comparing the extracted ambient noise impulse response waveforms with records of moderate magnitude earthquakes (from M w 4 to 5.8) that occurred close to the virtual source station in Japan. As the information is only available at low frequencies (less than 0.25 Hz), the ambient seismic field approach is coupled to a non-stationary stochastic model to simulate time domain accelerograms up to 50 Hz. This coupling allows the predicted ground motion to have both the deterministic part at low frequencies coming from the source and the crust structure and the high-frequency random contribution from the seismic waves scattering. The resulting combined accelerograms for an M w 5.8 event show a good agreement with observed ground motions from a real earthquake.

Description: We apply a reversible-jump Markov chain Monte Carlo method to sample the Bayesian posterior model probability density function of 2-D seafloor resistivity as constrained by marine controlled source electromagnetic data. This density function of earth models conveys information on which parts of the model space are illuminated by the data. Whereas conventional gradient-based inversion approaches require subjective regularization choices to stabilize this highly non-linear and non-unique inverse problem and provide only a single solution with no model uncertainty information, the method we use entirely avoids model regularization. The result of our approach is an ensemble of models that can be visualized and queried to provide meaningful information about the sensitivity of the data to the subsurface, and the level of resolution of model parameters. We represent models in 2-D using a Voronoi cell parametrization. To make the 2-D problem practical, we use a source–receiver common midpoint approximation with 1-D forward modelling. Our algorithm is transdimensional and self-parametrizing where the number of resistivity cells within a 2-D depth section is variable, as are their positions and geometries. Two synthetic studies demonstrate the algorithm's use in the appraisal of a thin, segmented, resistive reservoir which makes for a challenging exploration target. As a demonstration example, we apply our method to survey data collected over the Scarborough gas field on the Northwest Australian shelf.

Description: Compressional and shear wave seismic measurements were performed in an old railway tunnel and in galleries excavated in a 250-m-thick Toarcian claystone formation in the Tournemire experimental station (France). Three component (3C) geophones and three orthogonal orientations of the vibroseismic force source were used. Additionally, vertical seismic profiling (VSP) measurements were recorded with a 3C borehole geophone, a hydrophone and a microphone in a 159 m deep borehole (ID180) in the tunnel. The seismic data show that Toarcian claystone has strong transverse isotropy (TI) with a vertical symmetry axis. The qP , SH and qSV wave propagation velocities in horizontal directions—the plane of isotropy of the TI medium—are measured as 3550, 1850 and 1290 m s –1 , respectively. The zero-offset VSP reveals that only one shear wave propagates in the vertical (depth) direction and the P - and S -wave velocities are 3100 and 1375 m s –1 , respectively. Four elastic moduli of the TI medium are determined from the seismic velocities and from the bulk density of 2.53 g cm –3 : c 11 = 31.9 GPa, c 33 = 24.3 GPa, c 44 = 4.5 GPa and c 66 = 8.7 GPa. A walkaway VSP with the borehole geophone at 50 m depth in borehole ID180 and shot points in the galleries leads to oblique seismic ray paths which allow us to determine the fifth elastic modulus of the TI medium to c 13 = 16 GPa. The tube wave recorded by a hydrophone in the water filled lower part of the borehole propagates with 1350 m s –1 , which confirms the estimate of the elastic constant c 66 . The analysis of body wave and surface wave data from a seismic experiment in Galerie Est shows reflections from several fracture zones in the gallery floor. The thickness of the excavation damaged zone (EDZ) in the floor of Galerie Est is estimated to 0.7 m.

Description: By relocation of shallow and intermediate depth earthquakes and joint evaluation of already available and properly estimated waveform inversion focal mechanisms we investigate the location and shallow kinematics of the residual subducting slab in the Calabrian Arc region, that is the only, apparently still active segment of the old subduction front of the western Mediterranean. In agreement with high P -wave velocity anomaly found at intermediate depths by previous local earthquake tomography, our shallow-to-intermediate earthquake hypocentre distribution shows that the Ionian subducting slab is still in-depth continuous only in a small internal segment of the Arc, while detachment or break-off processes have already developed elsewhere along the Arc. At the same time, the space distribution and the waveform inversion focal mechanisms of the earthquakes occurring at shallow depth (〈70 km) do not evidence Subduction Transform Edge Propagator (STEP) fault activity at the edges of the descending slab. In particular, no trace is found of dip-slip faulting along near vertical planes parallel to the slab edges, that is no seismic evidence is available of vertical motion between the subducting segment of the plate and the adjacent portion of it. Also, the seismicity distribution and mechanisms found at crustal depths in the study region do not match properly with the expected scenario of relative motion at the lateral borders of the overriding plate. Our earthquake locations and mechanisms together with GPS information taken from the literature highlight a residual, laterally very short subducting slab showing quasi-nil velocity of trench retreat and no present-day STEP activity, still capable however of causing strong normal-faulting earthquakes in the trench area through its gravity-induced shallow deformation in a weak-coupling scenario.

Description: We combine new, well-determined GPS velocities from Clarion, Guadalupe and Socorro islands on young seafloor in the eastern Pacific basin with newly estimated velocities for 26 GPS sites from older seafloor in the central, western and southern parts of the Pacific Plate to test for deformation within the interior of the Pacific Plate and estimate the viscosity of the asthenosphere below the plate. Relative to a Pacific Plate reference frame defined from the velocities of the 26 GPS sites in other areas of the Pacific Plate, GPS sites on Clarion and Guadalupe islands in the eastern Pacific move 1.2 ± 0.6 mm yr –1 (1) towards S09°W ± 38° and 1.9 ± 0.3 mm yr –1 towards S19°E ± 10°, respectively. The two velocities, which are consistent within their 95 per cent uncertainties, both differ significantly from Pacific Plate motion. Transient volcanic deformation related to a 1993–1996 eruption of the Socorro Island shield volcano renders our GPS velocity from that island unreliable for the tectonic analysis although its motion is also southward like those of Clarion and Guadalupe islands. We test but reject the possibilities that drift of Earth's origin in ITRF2008 or unmodelled elastic offsets due to large-magnitude earthquakes around the Pacific rim since 1993 can be invoked to explain the apparent slow southward motions of Clarion and Guadalupe islands. Similarly, corrections to the Pacific Plate GPS velocity field for possible viscoelastic deformation triggered by large-magnitude earthquakes since 1950 also fail to explain the southward motions of the two islands. Viscoelastic models with prescribed asthenospheric viscosities lower than 1  x 10 19 Pa s instead introduce statistically significant inconsistencies into the Pacific Plate velocity field, suggesting that the viscosity of the asthenosphere below the plate is higher than 1  x 10 19 Pa s. Elastic deformation from locked Pacific–North America Plate boundary faults is also too small to explain the southward motions of the two islands. Horizontal thermal contraction of the plate interior may explain the motion observed at Clarion and Guadalupe islands, as might long-term tectonic deformation of the plate interior.

Description: Rayleigh wave group velocity dispersion measurements from local and regional earthquakes are used to interpret the lithospheric structure in the Gulf of California region. We compute group velocity maps for Rayleigh waves from 10 to 150 s using earthquakes recorded by broad-band stations of the Network of Autonomously Recording Seismographs in Baja California and Mexico mainland, UNM in Mexico, BOR, DPP and GOR in southern California and TUC in Arizona. The study area is gridded in 120 longitude cells by 180 latitude cells, with an equal spacing of 10 10 km. Assuming that each gridpoint is laterally homogeneous, for each period the tomographic maps are inverted to produce a 3-D lithospheric shear wave velocity model for the region. Near the Gulf of California rift axis, we found three prominent low shear wave velocity regions, which are associated with mantle upwelling near the Cerro Prieto volcanic field, the Ballenas Transform Fault and the East Pacific Rise. Upwelling of the mantle at lithospheric and asthenospheric depths characterizes most of the Gulf. This more detailed finding is new when compared to previous surface wave studies in the region. A low-velocity zone in northcentral Baja at ~28ºN which extends east–south–eastwards is interpreted as an asthenospheric window. In addition, we also identify a well-defined high-velocity zone in the upper mantle beneath central-western Baja California, which correlates with the previously interpreted location of the stalled Guadalupe and Magdalena microplates. We interpret locations of the fossil slab and slab window in light of the distribution of unique post-subduction volcanic rocks in the Gulf of California and Baja California. We also observe a high-velocity anomaly at 50-km depth extending down to ~130 km near the southwestern Baja coastline and beneath Baja, which may represent another remnant of the Farallon slab.

Description: We show that superresolution imaging in the far-field region of the sources and receivers is theoretically and practically possible if migration of resonant multiples is employed. A resonant multiple is one that bounces back and forth between two scattering points; it can also be the multiple between two smoothly varying interfaces as long as the reflection wave paths partially overlap and reflect from the same Fresnel zone. For a source with frequency f , compared to a one-way trip, N round trips in propagating between two scatterers increase the effective frequency by 2 N   f and decrease the effective wavelength by /(2 N ). Thus, multiples can, in principle, be used as high-frequency probes to estimate detailed properties of layers. Tests with both synthetic and field data validate this claim. Improved resolution by multiple imaging is not only feasible for crustal reflections, but might be applicable to mantle and core reverberations recorded by earthquake seismologists.

Description: We have investigated basic properties of transition remanent magnetization of natural magnetite in granite samples collected from the Minnesota River Valley, North America. Transition remanence was imparted during cooling and/or warming through the Verwey transition around 120 K. Depending on magnetic field conditions during cooling and warming, three types of transition remanences have been categorized: (1) TrRM, acquired during a cycle of field cooling and field warming; (2) TrWRM, acquired during zero-field cooling and field warming and (3) TrCRM, imparted during field cooling and zero-field warming. These remanences fulfil basic laws of remanent magnetization: (1) directions of the transition remanences are parallel to direction of the applied field (the law of parallelism), (2) intensities of the transition remanences are proportional to the applied field intensity (the law of proportionality) and (3) sum of the partial transition remanences is equal to the total transitional remanence, that is, TrRM = TrWRM + TrCRM. In addition, the ratio of TrRM to TRM LTD (the demagnetized component of thermoremanent magnetization by low-temperature demagnetization) shows a nearly constant value of ~0.34. This relationship might reflect differences in equilibrium magnetic domain state at low and high temperature.

Description: The 127 station NorthEast China Extended SeiSmic Array (NECESSArray) provides large quantities of high quality seismic data in northeast China that allow us to resolve lateral variations of Lg Q or crustal attenuation at 1 Hz ( Q o ) to 2.0° or greater. Using the reverse two-station/event method with 11 642 Lg path-amplitudes from 78 crustal earthquakes, we obtain a 2-D tomographic image of Lg Q o with values ranging from ~50 to 1400. A high degree of detail in the lateral variation of Lg attenuation is revealed in our tomographic image. High Q o regions are found in the Great Xing'an, Lesser Xing'an and Songen-Zhangguangcai Ranges. Low Q o regions are observed in the Songliao, Sanjiang and Erlian Basins. The lowest Q o is found near the Wudalianchi volcanic field and other Quaternary volcanic fields, the southern Songliao Basin, the western edge of the Erlian Basin and the Sanjiang Basin. Low Q o values are measured for paths that cross sedimentary basins with thick, unconsolidated sediments. Most of the high Lg attenuation in the Songliao Basin correlates reasonably well with low crustal Rayleigh wave phase velocity anomalies. The highest attenuating regions also correlate well with regions of Holocene volcanism.

Description: In this study, a new method for computing the sensitivity of the glacial isostatic adjustment (GIA) forward solution with respect to the Earth's mantle viscosity, the so-called the forward sensitivity method (FSM), and a method for computing the gradient of data misfit with respect to viscosity parameters, the so-called adjoint-state method (ASM), are presented. These advanced formal methods complement each other in the inverse modelling of GIA-related observations. When solving this inverse problem, the first step is to calculate the forward sensitivities by the FSM and use them to fix the model parameters that do not affect the forward model solution, as well as identifying and removing redundant parts of the inferred viscosity structure. Once the viscosity model is optimized in view of the forward sensitivities, the minimization of the data misfit with respect to the viscosity parameters can be carried out by a gradient technique which makes use of the ASM. The aim is this paper is to derive the FSM and ASM in the forms that are closely associated with the forward solver of GIA developed by Martinec. Since this method is based on a continuous form of the forward model equations, which are then discretized by spectral and finite elements, we first derive the continuous forms of the FSM and ASM and then discretize them by the spectral and finite elements used in the discretization of the forward model equations. The advantage of this approach is that all three methods (forward, FSM and ASM) have the same matrix of equations and use the same methodology for the implementation of the time evolution of stresses. The only difference between the forward method and the FSM and ASM is that the different numerical differencing schemes for the time evolution of the Maxwell and generalized Maxwell viscous stresses are applied in the respective methods. However, it requires only a little extra computational time for carrying out the FSM and ASM numerically. An straightforward approach to compute the gradient of the data misfit is the brute-force method, whereby the partial derivatives of the misfit with respect to model parameters are approximated by the centred difference of two forward model runs. Although the brute-force method is useful for computing the gradient of the data misfit with respect to a small number of model parameters, it becomes expensive for a viscosity model with a large number of parameters. The ASM offers an efficient alternative for computing the gradient of the misfit since the computational time of the ASM is independent of the number of viscosity parameters. The ASM is thus highly efficient for calculating the gradient of the misfit for models with large numbers of parameters. However, the forward-model solution for each time step must be stored, hence the memory demands scale linearly with the number of time steps. This is the main drawback of the ASM.

Description: The Ricker wavelet has been widely used in the analysis of seismic data, as its asymmetrical amplitude spectrum can represent the attenuation feature of seismic wave propagation through viscoelastic homogeneous media. However, the frequency band of the Ricker wavelet is not analytically determined yet. The determination of the frequency band leads to an inverse exponential equation. To solve this equation analytically a special function, the Lambert W function, is needed. The latter provides a closed and elegant expression of the frequency band of the Ricker wavelet, which is a sample application of the Lambert W function in geophysics and there have been other applications in various scientific and engineering fields in the past decade. Moreover, the Lambert W function is a variation of the Ricker wavelet amplitude spectrum. Since the Ricker wavelet is the second derivative of a Gaussian function and its spectrum is a single-valued smooth curve, numerical evaluation of the Lambert W function can be implemented by a stable interpolation procedure, followed by a recursive computation for high precision.

Description: The availability of accurate and reliable age control is a crucial prerequisite for palaeoclimate studies, particularly when different archives are to be compared. Here we present a detailed depth-to-time transformation for the lacustrine sediments of a ~940-m-long drill core (SG-1) from the western Qaidam Basin (NE Tibetan Plateau). To establish a more precise age model than the one previously available, which was based solely on magnetostratigraphic dating using polarity boundaries as tie points, we applied time-series analysis on magnetic susceptibility ( ) variation. The data are available in high resolution and are considered to be closely linked to orbital forcing. Since the sediment accumulation rate (SAR) varies strongly throughout the succession of core SG-1, conventional cyclostratigraphy by bandpass filtering cannot be applied. We present two alternative age models based on spectral characteristics and orbital tuning. The first age model (TPspec) is based on the assumption that changes in SAR occurred when the frequency spectra revealed obviously different characteristics in the spectral pattern. For the second age model (SARA), SAR was adjusted every 2 m by comparing observed with expected orbital cycles in accordance with the age of magnetic reversals. This age model appears more robust and shows the most convincing spectral results in the frequency and wavelet power spectrum of . According to the SARA age model, SAR varies between 14 and 73 cm kyr –1 , and the bottom of SG-1 has an age of 2.69 Ma. Our results show that orbital tuning can be successfully applied for sequences with strongly variable SAR. The age model SARA can be used for a more detailed analysis of the existing multiproxy data set in terms of palaeoclimate evolution. The most prominent feature of spectra using the SARA age model is the identification of the middle Pleistocene transition.

Description: Clay minerals and Al-substituted haematite (Al-hm) usually coexist in soils and sediments. However, effects of clay minerals on Al-hm during thermal magnetic measurements in argon environment have not been well studied. In order to quantify such effects, a series of Al-hm samples were synthesized, and were then mixed with clay minerals (illite, chlorite, kaolinite and Ca-montmorillonite). The temperature dependence of magnetic susceptibility curves in an argon environment showed that Al-substituted magnetite was produced during the thermal treatment via the reduction of Al-hm by the clay mineral, which leads to a significant magnetic enhancement of the thermal products. In addition, the reductive capacity varies among different types of clay minerals, that is, illite 〉 chlorite 〉 kaolinite 〉 Ca-montmorillonite. Furthermore, the iron content in the clay minerals and Al content of Al-hm are two predominant factors controlling the reduced haematite content. The iron is released from the clay minerals and provides the reducing agent, while Al decreases the crystallinity of haematite and thus facilitates the chemical reaction. Therefore, the thermal magnetic measurements can be used to quantify the Al content of Al-hm in natural samples. Our study provides significant information for palaeomagnetism and environmental magnetism studies, such as thermal magnetic analysis and palaeomagnetic intensity reconstruction using ancient pottery and kilns.

Description: We report a palaeomagnetic investigation of the last full geomagnetic field reversal, the Matuyama-Brunhes (M-B) transition, as preserved in a continuous sequence of exposed lacustrine sediments in the Apennines of Central Italy. The palaeomagnetic record provides the most direct evidence for the tempo of transitional field behaviour yet obtained for the M-B transition. 40 Ar/ 39 Ar dating of tephra layers bracketing the M-B transition provides high-accuracy age constraints and indicates a mean sediment accumulation rate of about 0.2 mm yr –1 during the transition. Two relative palaeointensity (RPI) minima are present in the M-B transition. During the terminus of the upper RPI minimum, a directional change of about 180 ° occurred at an extremely fast rate, estimated to be less than 2 ° per year, with no intermediate virtual geomagnetic poles (VGPs) documented during the transit from the southern to northern hemisphere. Thus, the entry into the Brunhes Normal Chron as represented by the palaeomagnetic directions and VGPs developed in a time interval comparable to the duration of an average human life, which is an order of magnitude more rapid than suggested by current models. The reported investigation therefore provides high-resolution integrated palaeomagnetic and radioisotopic data that document the fine details of the anatomy and tempo of the M-B transition in Central Italy that in turn are crucial for a better understanding of Earth's magnetic field, and for the development of more sophisticated models that are able to describe its global structure and behaviour.

Description: We develop a frequency-based approach to earthquake slip inversion that requires no prior information on the rupture velocity or slip-rate functions. Because the inversion is linear and is performed separately at each frequency, it is computationally efficient and suited to imaging the finest resolvable spatial details of rupture. We demonstrate the approach on synthetic seismograms based on the Source Inversion Validation Exercise 1 (SIV1) of a crustal M w 6.6 strike-slip earthquake recorded locally. A robust inversion approach is obtained by applying a combination of damping, smoothing and forcing zero slip at the edge of the fault model. This approach achieves reasonable data fits, overall agreement to the SIV1 model, including slip-rate functions of each subfault, from which its total slip, slip time history and rupture velocity can be extracted. We demonstrate the method's robustness by exploring the effects of noise, random timing errors, and fault geometry errors. The worst effects on the inversion are seen from errors in the assumed fault geometry.

Description: Passive seismic methods using earthquakes can be applied for extracting body waves and obtaining information of subsurface structure. In this study, we retrieve direct and reflected plane waves by applying seismic interferometry to the recorded ground motion from a cluster of regional earthquakes. We apply upgoing/downgoing P/S wavefield decomposition, time windowing, and multidimensional deconvolution to improve the quality of the extraction of reflected waves with seismic interferometry. The wavefield separation and seismic interferometry based on multidimensional deconvolution allow us to reconstruct PP, PS, SP and SS reflected waves without unwanted crosstalk between P and S waves. From earthquake data, we obtain PP, PS and SS reflected plane waves that reflect off the same reflector, and estimate P - and S -wave velocities.

Description: During megathrust earthquakes, great ruptures are accompanied by large scale mass redistribution inside the solid Earth and by ocean mass redistribution due to bathymetry changes. These large scale mass displacements can be detected using the monthly gravity maps of the GRACE satellite mission. In recent years it has become increasingly common to use the long wavelength changes in the Earth's gravity field observed by GRACE to infer seismic source properties for large megathrust earthquakes. An important advantage of space gravimetry is that it is independent from the availability of land for its measurements. This is relevant for observation of megathrust earthquakes, which occur mostly offshore, such as the $M_{\text{w}} \sim 9$ 2004 Sumatra–Andaman, 2010 Maule (Chile) and 2011 Tohoku-Oki (Japan) events. In Broerse et al. , we examined the effect of the presence of an ocean above the rupture on long wavelength gravity changes and showed it to be of the first order. Here we revisit the implementation of an ocean layer through the sea level equation and compare the results with approximated methods that have been used in the literature. One of the simplifications usually lies in the assumption of a globally uniform ocean layer. We show that especially in the case of the 2010 Maule earthquake, due to the closeness of the South American continent, the uniform ocean assumption is not valid and causes errors up to 57 per cent for modelled peak geoid height changes (expressed at a spherical harmonic truncation degree of 40). In addition, we show that when a large amount of slip occurs close to the trench, horizontal motions of the ocean floor play a mayor role in the ocean contribution to gravity changes. Using a slip model of the 2011 Tohoku-Oki earthquake that places the majority of slip close to the surface, the peak value in geoid height change increases by 50 per cent due to horizontal ocean floor motion. Furthermore, we test the influence of the maximum spherical harmonic degree at which the sea level equation is performed for sea level changes occurring along coastlines, which shows to be important for relative sea level changes occurring along the shore. Finally, we demonstrate that ocean floor loading, self-gravitation of water and conservation of water mass are of second order importance for coseismic gravity changes. When GRACE observations are used to determine earthquake parameters such as seismic moment or source depth, the uniform ocean layer method introduces large biases, depending on the location of the rupture with respect to the continent. The same holds for interpreting shallow slip when horizontal motions are not properly accounted for in the ocean contribution. In both cases the depth at which slip occurs will be underestimated.

Description: The representation of viscoelastic media in the time domain becomes more challenging with greater bandwidth of the propagating waves and number of travelled wavelengths. With the continuously increasing computational power, more extreme parameter regimes become accessible, which requires the reassessment and improvement of the standard ‘memory variable’ methods to implement attenuation in time-domain seismic wave-propagation methods. In this paper, we propose a method to minimize the error in the wavefield for a fixed complexity of the anelastic medium. This method consists of defining an appropriate misfit criterion based on a first-order analysis of how errors in the discretized medium propagate into errors in the wavefield and a simulated annealing optimization scheme to find the globally optimal parametrization. Furthermore, we derive an analytical time-stepping scheme for the memory variables that encode the strain history of the medium. Then we develop the coarse grained memory variable approach for the spectral element method (SEM) and benchmark it using the 2.5-D code AxiSEM for global body waves up to 1 Hz. Showing very good agreement with a reference solution, it also leads to a speedup of a factor of 5 in the anelastic part of the code (factor 2 in total) in this 2.5-D approach. A factor of 15 (3 in total) can be expected for the 3-D case compared to conventional implementations.

Description: In the western Pacific, high-frequency seismic energy is carried to very great distances from the source. The Po and So phases with observed seismic velocities characteristic of the mantle lithosphere have complex and elongated waveforms that are well explained by a model of stochastic heterogeneity. However, in the eastern part of the Pacific Basin equivalent paths show muted Po and weak, or missing, So . Once established, it is hard to eliminate such guided Po and So energy in the mantle lithosphere by purely structural effects. Even sharp changes in lithospheric thickness or complex transitions at fracture zones only weaken the mantle ducted wave trains, but Po and So remain distinct. In contrast, the effect of attenuation is much more severe and can lead to suppression of the So phase to below the noise level after passage of a few hundred kilometres. The differing characteristics of Po and So across the Pacific can therefore be related directly to the thermal state via the enhanced attenuation in hotter regions, such as the spreading ridges and backarc regions.

Description: The anisotropy of magnetic susceptibility (AMS) of soft rocks was measured in order to distinguish between the effect of remote and local strain fields, determine the size of the related inelastic damage zone and resolve the fault-plane solutions of past earthquakes. The AMS fabrics were explored next to late Pleistocene syndepositional normal faults (total displacement up to ~3.5 m) that cross soft lacustrine rocks within the seismically active Dead Sea basin. ‘Deposition fabrics’ prevail meters away from the fault planes and are characterized by scattered maximum and intermediate principal AMS axes. ‘Deformation fabrics’ are detected up to tens of centimetres from the fault planes and are characterized by well-grouped AMS axes, in which one of the principal axes is parallel to the strike of the nearby fault. Variations in the AMS fabrics and magnetic lineations define the size of the inelastic damage zone around the faults. The results demonstrate that the deformation-driven magnetic fabrics and the associated inelastic damage zones are compatible with coseismic dynamic faulting and the effects of the local strain field during earthquakes. Most of the AMS fabrics show a conspicuous similarity to that of the fault-plane solutions, i.e., the principal AMS axes and instantaneous strain ellipsoids are coaxial. These results suggest a novel application of the AMS method for defining the shape and size of the damage zones surrounding dynamic faults and determining the full tensor of the local strain field.

Description: The Tonga-Kermadec forearc is deforming in response to on-going subduction of the Pacific Plate beneath the Indo-Australian Plate. Previous research has focussed on the structural development of the forearc where large bathymetric features such as the Hikurangi Plateau and Louisville Ridge seamount chain are being subducted. Consequently, knowledge of the ‘background’ forearc in regions of normal plate convergence is limited. We report on an ~250-km-long multichannel seismic reflection profile that was shot perpendicular to the Tonga-Kermadec trench at ~28°S to determine the lateral and temporal variations in the structure, stratigraphy and deformation of the Kermadec forearc resulting solely from Pacific Plate subduction. Interpretation of the seismic profile, in conjunction with regional swath bathymetry data, shows that the Pacific Plate exhibits horst and graben structures that accommodate bending-induced extensional stresses, generated as the trenchward dip of the crust increases. Trench infill is also much thicker than expected at 1 km which, we propose, results from increased sediment flux into and along the trench. Pervasive normal faulting of the mid-trench slope most likely accommodates the majority of the observed forearc extension in response to basal subduction erosion, and a structural high is located between the mid- and upper-trench slopes. We interpret this high as representing a dense and most likely structurally robust region of crust lying beneath this region. Sediment of the upper-trench slope documents depositional hiatuses and on-going uplift of the arc. Strong along-arc currents appear to erode the Kermadec volcanic arc and distribute this sediment to the surrounding basins, while currents over the forearc redistribute deposits as sediment waves. Minor uplift of the transitional Kermadec forearc, observed just to the north of the profile, appears to relate to an underlying structural trend as well as subduction of the Louisville Ridge seamount chain 250 km to the north. Relative uplift of the Kermadec arc is observed from changes in the tilt of upper-trench slope deposits and extensional faulting of the basement immediately surrounding the Louisville Ridge.

Description: Very low frequency (VLF) electromagnetic waves that penetrate conductive magma-filled dykes generate secondary fields on the surface that can be used to invert for dyke properties. The model used for the interpretation calculates currents induced in a conductive strip by an inducing field that decays exponentially with depth due to the conductivity of the surrounding medium. The differential equations are integrated to give an inhomogeneous Fredholm equation of the second kind with a kernel consisting of a modified Bessel function of the second kind. Numerical methods are typically used to solve for the induced currents in the strip. In this paper, we apply a modified Galerkin–Chebyshev method, which involves separating the kernel into source and field spectra and integrating the source terms to obtain a matrix equation for the unknown coefficients. The incident wave is expressed as a Chebyshev series. The modified Bessel function is separated into a logarithmic singularity and a non-singular remainder, both of which are expanded in complex Chebyshev polynomials. The Chebyshev coefficients for the remainder are evaluated using a fast Fourier transform, while the logarithmic term and incident field have analytic series. The deconvolution then involves a matrix inversion. The results depend on the ratio of strip-size to skin-depth. For infinite skin-depth and a singular conductivity distribution given by $\tau _0 a/\sqrt{a^2 - z^2 }$ (where 0 is the conductance, a is the half-length and z the distance from the centre), Parker gives an analytic solution. We present a similar analytic series solution for the finite skin-depth case, where the size to skin depth ratio is small. Results are presented for different ratios of size to skin depth that can be compared with numerical solutions. We compare full-space and half-space solutions. A fit of the model to VLF data taken above a magma filled dykes in Hawaii and Mt Etna demonstrates that while properties such as depth to top, conductivity ratio, tilt and dip can be determined, the depth to bottom is indeterminate due to the exponential decay of the inducing field.

Description: The radially anisotropic shear velocity structure of the Earth's mantle provides a critical window on the interior dynamics of the planet, with isotropic variations that are interpreted in terms of thermal and compositional heterogeneity and anisotropy in terms of flow. While significant progress has been made in the more than 30 yr since the advent of global seismic tomography, many open questions remain regarding the dual roles of temperature and composition in shaping mantle convection, as well as interactions between different dominant scales of convective phenomena. We believe that advanced seismic imaging techniques, such as waveform inversion using accurate numerical simulations of the seismic wavefield, represent a clear path forwards towards addressing these open questions through application to whole-mantle imaging. To this end, we employ a ‘hybrid’ waveform-inversion approach, which combines the accuracy and generality of the spectral finite element method (SEM) for forward modelling of the global wavefield, with non-linear asymptotic coupling theory for efficient inverse modelling. The resulting whole-mantle model (SEMUCB-WM1) builds on the earlier successful application of these techniques for global modelling at upper mantle and transition-zone depths (≤800 km) which delivered the models SEMum and SEMum2. Indeed, SEMUCB-WM1 is the first whole-mantle model derived from fully numerical SEM-based forward modelling. Here, we detail the technical aspects of the development of our whole-mantle model, as well as provide a broad discussion of isotropic and radially anisotropic model structure. We also include an extensive discussion of model uncertainties, specifically focused on assessing our results at transition-zone and lower-mantle depths.

Description: Anchor Bay and surrounding regions are located on the northwest coast of the island of Malta, Central Mediterranean. The area is characterized by a coastal cliff environment having an outcropping layer of hard coralline limestone (UCL) resting on a thick (up to 50 m) layer of clays and marls (Blue Clay, BC). This configuration gives rise to coastal instability effects, in particular lateral spreading phenomena and rock falls. Previous and ongoing studies have identified both lateral spreading rates and vertical motions of several millimetres per year. The area is an interesting natural laboratory as coastal detachment processes in a number of different stages can be identified and are easily accessible. We investigate the site dynamic characteristics of this study area by recording ambient noise time-series at more than 30 points, over an area of 0.07 km 2 , using a portable three-component seismograph. The time-series are processed to give both horizontal-to-vertical spectral ratio graphs (H/V) as well as frequency-dependent polarisation analysis. The H/V graphs illustrate and quantify aspects of site resonance effects due both to underlying geology as well as to mechanical resonance of partly or wholly detached blocks. The polarization diagrams indicate the degree of linearity and predominant directions of vibrational effects. H/V curves closer to the cliff edge show complex responses at higher frequencies, characteristic of the dynamic behaviour of individual detached blocks. Particle motion associated with the higher frequencies shows strongly directional polarization and a high degree of linearity at well-defined frequencies, indicative of normal-mode vibration. The stable plateau areas, on the other hand, show simple, single-peak H/V curves representative of the underlying stratification and no predominant polarization direction. These results, which will be compared with those from other experiments in the area, have important implications for the understanding of ongoing processes in geologically active and unstable coastal environments.

Description: Correlated data errors (noise) are common in magnetotelluric (MT) data, but MT inversions typically neglect error correlations without investigating the impact of this simplification on inversion results. This paper examines effects of neglecting frequency- and spatially correlated noise on MT inversion, based on a nonlinear Bayesian formulation which quantifies the uncertainties of inversion results in terms of marginal posterior probability densities and credibility intervals. To do so, data with frequency- and spatially correlated noise of differing degrees are generated for several layered (1-D) synthetic cases. Bayesian MT inversions are carried out for these data sets with and without accounting for error correlation (i.e. applying full and diagonal covariance matrices, respectively, in the inversion), and the results are compared. For cases with noise that is strongly correlated over frequency or space, parameter uncertainties estimated using the diagonal-covariance simplification (neglecting error correlations) are found to often be significantly underestimated compared with results computed using the full covariance matrix.

Description: We focus on the relation between seismic and total postseismic afterslip following the Maule M w 8.8 earthquake on 2010 February 27 in central Chile. First, we calculate the cumulative slip released by aftershock seismicity. We do this by summing up the aftershock regions and slip estimated from scaling relations. Comparing the cumulative seismic slip with afterslip models we show that seismic slip of individual aftershocks exceeds locally the inverted afterslip model from geodetic constraints. As the afterslip model implicitly contains the displacements from the aftershocks, this reflects the tendency of afterslip models to smear out the actual slip pattern. However, it also suggests that locally slip for a number of the larger aftershocks exceeds the aseismic slip in spite of the fact that the total equivalent moment of the afterslip exceeds the cumulative moment of aftershocks by a large factor. This effect, seen weakly for the Maule 2010 and also for the Tohoku 2011 earthquake, can be explained by taking into account the uncertainties of the seismicity and afterslip models. In spite of uncertainties, the hypocentral region of the Nias 2005 earthquake is suggested to release a large fraction of moment almost purely seismically. Therefore, these aftershocks are not driven solely by the afterslip but instead their slip areas have probably been stressed by interseismic loading and the mainshock rupture. In a second step, we divide the megathrust of the Maule 2010 rupture into discrete cells and count the number of aftershocks that occur within 50 km of the centre of each cell as a function of time. We then compare this number to a time-dependent afterslip model by defining the ‘afterslip to aftershock ratio’ (ASAR) for each cell as the slope of the best fitting line when the afterslip at time t is plotted against aftershock count. Although we find a linear relation between afterslip and aftershocks for most cells, there is significant variability in ASAR in both the downdip and along-strike directions of the megathrust. We compare the spatial distribution of ASAR with the spatial distribution of seismic coupling, coseismic slip and Bouguer gravity anomaly, and in each case we find no significant correlation.

Description: Following earlier studies, we present forward and inverse simulations of heat and fluid transport of the upper crust using a local 3-D model of the Kola area. We provide best estimates for palaeotemperatures and permeabilities, their errors and their dependencies. Our results allow discriminating between the two mentioned processes to a certain extent, partly resolving the non-uniqueness of the problem. We find clear indications for a significant contribution of advective heat transport, which, in turn, imply only slightly lower ground surface temperatures during the last glacial maximum relative to the present value. These findings are consistent with the general background knowledge of (i) the fracture zones and the corresponding fluid movements in the bedrock and (ii) the glacial history of the Kola area.

Description: We analysed the conversion problem between teleseismic magnitudes ( M s and m b ) provided by the Seismological Bulletin of the International Seismological Centre and moment magnitudes ( M w ) provided by online moment tensor (MT) catalogues using the chi-square general orthogonal regression method (CSQ) that, differently from the ordinary least-square regression method (OLS), accounts for the measurement errors of both the predictor and response variables. To account for the non-linearity of the relationships, we used two types of curvilinear models: (i) the exponential model (EXP), recently proposed by the authors of the Global Catalogue sponsored by the Global Earthquake Model (GEM) Foundation and (ii) a connected bilinear (CBL) model, similar to that proposed by Ekström & Dziewonski, where two different linear trends at low and high magnitudes are connected by an arc of circle that preserves the continuity of the function and of its first derivative at the connecting points. For M s , we found that the regression curves computed for a global data set (GBL) are likely to be biased by the incompleteness of global MT catalogues for M w 〈5.0–5.5. In fact, the GBL curves deviate significantly from a similar regression curve computed for a Euro-Mediterranean data set (MED) integrated with the data provided by two regional MT catalogues including many more events with M w 〈 5.0–5.5. The GLB regression curves overestimate the M w proxies computed from M s up to 0.5 magnitude units. Hence for computing M w proxies at the global scale of M s ≤ 5.5, we suggest to adopt the coefficients obtained from the MED regression. The analysis of the frequency–magnitude relationship of the resulting M w proxy catalogues confirms the validity of this choice as the behaviour of b­ -value as a function of cut-off magnitude of the GBL data set is much more stable using such approach. The incompleteness of M w 's provided from MT global catalogues also affects the m b GBL data set but in this case the use of the CSQ regression method, in place of the OLS, mitigates the bias and then, at low magnitudes, the EXP regression curve computed from the more complete MED data set almost coincides with that computed from the GBL data set. Our results also indicate that the slope at low magnitudes of the M w – M s relationship is substantially consistent with the hypothesized theoretical value of 2/3 for M s 〈 5.0 while the slope of the M w – m b relationship at high magnitudes probably reaches the theoretically expected value of 2 only in the proximity of the upper limit of m b determinations in our data set ( m b  = 7.2).