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: 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: 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: Several bioinformatics methods have been proposed for the detection and characterization of genomic structural variation (SV) from ultra high-throughput genome resequencing data. Recent surveys show that comprehensive detection of SV events of different types between an individual resequenced genome and a reference sequence is best achieved through the combination of methods based on different principles (split mapping, reassembly, read depth, insert size, etc.). The improvement of individual predictors is thus an important objective. In this study, we propose a new method that combines deviations from expected library insert sizes and additional information from local patterns of read mapping and uses supervised learning to predict the position and nature of structural variants. We show that our approach provides greatly increased sensitivity with respect to other tools based on paired end read mapping at no cost in specificity, and it makes reliable predictions of very short insertions and deletions in repetitive and low-complexity genomic contexts that can confound tools based on split mapping of reads.

Description: Exome sequencing strategy is promising for finding novel mutations of human monogenic disorders. However, pinpointing the casual mutation in a small number of samples is still a big challenge. Here, we propose a three-level filtration and prioritization framework to identify the casual mutation(s) in exome sequencing studies. This efficient and comprehensive framework successfully narrowed down whole exome variants to very small numbers of candidate variants in the proof-of-concept examples. The proposed framework, implemented in a user-friendly software package, named KGGSeq ( http://statgenpro.psychiatry.hku.hk/kggseq ), will play a very useful role in exome sequencing-based discovery of human Mendelian disease genes.

Description: Small RNAs (sRNAs) are a class of short (20–25 nt) non-coding RNAs that play important regulatory roles in gene expression. An essential first step in understanding their function is to confidently identify sRNA targets. In plants, several classes of sRNAs such as microRNAs (miRNAs) and trans-acting small interfering RNAs have been shown to bind with near-perfect complementarity to their messenger RNA (mRNA) targets, generally leading to cleavage of the mRNA. Recently, a high-throughput technique known as Parallel Analysis of RNA Ends (PARE) has made it possible to sequence mRNA cleavage products on a large-scale. Computational methods now exist to use these data to find targets of conserved and newly identified miRNAs. Due to speed limitations such methods rely on the user knowing which sRNA sequences are likely to target a transcript. By limiting the search to a tiny subset of sRNAs it is likely that many other sRNA/mRNA interactions will be missed. Here, we describe a new software tool called PAREsnip that allows users to search for potential targets of all sRNAs obtained from high-throughput sequencing experiments. By searching for targets of a complete ‘sRNAome’ we can facilitate large-scale identification of sRNA targets, allowing us to discover regulatory interaction networks.

Description: Gene fusions are common driver events in leukaemias and solid tumours; here we present FusionAnalyser, a tool dedicated to the identification of driver fusion rearrangements in human cancer through the analysis of paired-end high-throughput transcriptome sequencing data. We initially tested FusionAnalyser by using a set of in silico randomly generated sequencing data from 20 known human translocations occurring in cancer and subsequently using transcriptome data from three chronic and three acute myeloid leukaemia samples. in all the cases our tool was invariably able to detect the presence of the correct driver fusion event(s) with high specificity. In one of the acute myeloid leukaemia samples, FusionAnalyser identified a novel, cryptic, in-frame ETS2–ERG fusion. A fully event-driven graphical interface and a flexible filtering system allow complex analyses to be run in the absence of any a priori programming or scripting knowledge. Therefore, we propose FusionAnalyser as an efficient and robust graphical tool for the identification of functional rearrangements in the context of high-throughput transcriptome sequencing data.

Description: The use of a priori knowledge in the alignment of targeted sequencing data is investigated using computational experiments. Adapting a Needleman–Wunsch algorithm to incorporate the genomic position information from the targeted capture, we demonstrate that alignment can be done to just the target region of interest. When in addition use is made of direct string comparison, an improvement of up to a factor of 8 in alignment speed compared to the fastest conventional aligner (Bowtie) is obtained. This results in a total alignment time in targeted sequencing of around 7 min for aligning approximately 56 million captured reads. For conventional aligners such as Bowtie, BWA or MAQ, alignment to just the target region is not feasible as experiments show that this leads to an additional 88% SNP calls, the vast majority of which are false positives (~92%).

Description: We introduce Grinder ( http://sourceforge.net/projects/biogrinder/ ), an open-source bioinformatic tool to simulate amplicon and shotgun (genomic, metagenomic, transcriptomic and metatranscriptomic) datasets from reference sequences. This is the first tool to simulate amplicon datasets (e.g. 16S rRNA) widely used by microbial ecologists. Grinder can create sequence libraries with a specific community structure, α and β diversities and experimental biases (e.g. chimeras, gene copy number variation) for commonly used sequencing platforms. This versatility allows the creation of simple to complex read datasets necessary for hypothesis testing when developing bioinformatic software, benchmarking existing tools or designing sequence-based experiments. Grinder is particularly useful for simulating clinical or environmental microbial communities and complements the use of in vitro mock communities.

Description: The most crucial step in data processing from high-throughput sequencing applications is the accurate and sensitive alignment of the sequencing reads to reference genomes or transcriptomes. The accurate detection of insertions and deletions (indels) and errors introduced by the sequencing platform or by misreading of modified nucleotides is essential for the quantitative processing of the RNA-based sequencing (RNA-Seq) datasets and for the identification of genetic variations and modification patterns. We developed a new, fast and accurate algorithm for nucleic acid sequence analysis, FANSe, with adjustable mismatch allowance settings and ability to handle indels to accurately and quantitatively map millions of reads to small or large reference genomes. It is a seed-based algorithm which uses the whole read information for mapping and high sensitivity and low ambiguity are achieved by using short and non-overlapping reads. Furthermore, FANSe uses hotspot score to prioritize the processing of highly possible matches and implements modified Smith–Watermann refinement with reduced scoring matrix to accelerate the calculation without compromising its sensitivity. The FANSe algorithm stably processes datasets from various sequencing platforms, masked or unmasked and small or large genomes. It shows a remarkable coverage of low-abundance mRNAs which is important for quantitative processing of RNA-Seq datasets.

Description: Ultra-deep RNA sequencing has become a powerful approach for genome-wide analysis of pre-mRNA alternative splicing. We develop MATS (multivariate analysis of transcript splicing), a Bayesian statistical framework for flexible hypothesis testing of differential alternative splicing patterns on RNA-Seq data. MATS uses a multivariate uniform prior to model the between-sample correlation in exon splicing patterns, and a Markov chain Monte Carlo (MCMC) method coupled with a simulation-based adaptive sampling procedure to calculate the P -value and false discovery rate (FDR) of differential alternative splicing. Importantly, the MATS approach is applicable to almost any type of null hypotheses of interest, providing the flexibility to identify differential alternative splicing events that match a given user-defined pattern. We evaluated the performance of MATS using simulated and real RNA-Seq data sets. In the RNA-Seq analysis of alternative splicing events regulated by the epithelial-specific splicing factor ESRP1, we obtained a high RT–PCR validation rate of 86% for differential exon skipping events with a MATS FDR of 〈10%. Additionally, over the full list of RT–PCR tested exons, the MATS FDR estimates matched well with the experimental validation rate. Our results demonstrate that MATS is an effective and flexible approach for detecting differential alternative splicing from RNA-Seq data.

Description: Numerous algorithms have been developed to analyze ChIP-Seq data. However, the complexity of analyzing diverse patterns of ChIP-Seq signals, especially for epigenetic marks, still calls for the development of new algorithms and objective comparisons of existing methods. We developed Qeseq, an algorithm to detect regions of increased ChIP read density relative to background. Qeseq employs critical novel elements, such as iterative recalibration and neighbor joining of reads to identify enriched regions of any length. To objectively assess its performance relative to other 14 ChIP-Seq peak finders, we designed a novel protocol based on Validation Discriminant Analysis (VDA) to optimally select validation sites and generated two validation datasets, which are the most comprehensive to date for algorithmic benchmarking of key epigenetic marks. In addition, we systematically explored a total of 315 diverse parameter configurations from these algorithms and found that typically optimal parameters in one dataset do not generalize to other datasets. Nevertheless, default parameters show the most stable performance, suggesting that they should be used. This study also provides a reproducible and generalizable methodology for unbiased comparative analysis of high-throughput sequencing tools that can facilitate future algorithmic development.

Description: The informational content of RNA sequencing is currently far from being completely explored. Most of the analyses focus on processing tables of counts or finding isoform deconvolution via exon junctions. This article presents a comparison of several techniques that can be used to estimate differential expression of exons or small genomic regions of expression, based on their coverage function shapes. The problem is defined as finding the differentially expressed exons between two samples using local expression profile normalization and statistical measures to spot the differences between two profile shapes. Initial experiments have been done using synthetic data, and real data modified with synthetically created differential patterns. Then, 160 pipelines (5 types of generator x 4 normalizations x 8 difference measures) are compared. As a result, the best analysis pipelines are selected based on linearity of the differential expression estimation and the area under the ROC curve. These platform-independent techniques have been implemented in the Bioconductor package rnaSeqMap. They point out the exons with differential expression or internal splicing, even if the counts of reads may not show this. The areas of application include significant difference searches, splicing identification algorithms and finding suitable regions for QPCR primers.

Description: The rapid expansion in the quantity and quality of RNA-Seq data requires the development of sophisticated high-performance bioinformatics tools capable of rapidly transforming this data into meaningful information that is easily interpretable by biologists. Currently available analysis tools are often not easily installed by the general biologist and most of them lack inherent parallel processing capabilities widely recognized as an essential feature of next-generation bioinformatics tools. We present here a user-friendly and fully automated R NA- S eq a nalysis p ipeline (R-SAP) with built-in multi-threading capability to analyze and quantitate high-throughput RNA-Seq datasets. R-SAP follows a hierarchical decision making procedure to accurately characterize various classes of transcripts and achieves a near linear decrease in data processing time as a result of increased multi-threading. In addition, RNA expression level estimates obtained using R-SAP display high concordance with levels measured by microarrays.

Description: Deciphering the structure of gene regulatory networks across the tree of life remains one of the major challenges in postgenomic biology. We present a novel ChIP-seq workflow for the archaea using the model organism Halobacterium salinarum sp. NRC-1 and demonstrate its application for mapping the genome-wide binding sites of natively expressed transcription factors. This end-to-end pipeline is the first protocol for ChIP-seq in archaea, with methods and tools for each stage from gene tagging to data analysis and biological discovery. Genome-wide binding sites for transcription factors with many binding sites (TfbD) are identified with sensitivity, while retaining specificity in the identification the smaller regulons (bacteriorhodopsin-activator protein). Chromosomal tagging of target proteins with a compact epitope facilitates a standardized and cost-effective workflow that is compatible with high-throughput immunoprecipitation of natively expressed transcription factors. The Pique package, an open-source bioinformatics method, is presented for identification of binding events. Relative to ChIP-Chip and qPCR, this workflow offers a robust catalog of protein–DNA binding events with improved spatial resolution and significantly decreased cost. While this study focuses on the application of ChIP-seq in H. salinarum sp. NRC-1, our workflow can also be adapted for use in other archaea and bacteria with basic genetic tools.

Description: A flexible statistical framework is developed for the analysis of read counts from RNA-Seq gene expression studies. It provides the ability to analyse complex experiments involving multiple treatment conditions and blocking variables while still taking full account of biological variation. Biological variation between RNA samples is estimated separately from the technical variation associated with sequencing technologies. Novel empirical Bayes methods allow each gene to have its own specific variability, even when there are relatively few biological replicates from which to estimate such variability. The pipeline is implemented in the edgeR package of the Bioconductor project. A case study analysis of carcinoma data demonstrates the ability of generalized linear model methods (GLMs) to detect differential expression in a paired design, and even to detect tumour-specific expression changes. The case study demonstrates the need to allow for gene-specific variability, rather than assuming a common dispersion across genes or a fixed relationship between abundance and variability. Genewise dispersions de-prioritize genes with inconsistent results and allow the main analysis to focus on changes that are consistent between biological replicates. Parallel computational approaches are developed to make non-linear model fitting faster and more reliable, making the application of GLMs to genomic data more convenient and practical. Simulations demonstrate the ability of adjusted profile likelihood estimators to return accurate estimators of biological variability in complex situations. When variation is gene-specific, empirical Bayes estimators provide an advantageous compromise between the extremes of assuming common dispersion or separate genewise dispersion. The methods developed here can also be applied to count data arising from DNA-Seq applications, including ChIP-Seq for epigenetic marks and DNA methylation analyses.

Description: ChIP-seq is increasingly used to characterize transcription factor binding and chromatin marks at a genomic scale. Various tools are now available to extract binding motifs from peak data sets. However, most approaches are only available as command-line programs, or via a website but with size restrictions. We present peak-motifs , a computational pipeline that discovers motifs in peak sequences, compares them with databases, exports putative binding sites for visualization in the UCSC genome browser and generates an extensive report suited for both naive and expert users. It relies on time- and memory-efficient algorithms enabling the treatment of several thousand peaks within minutes. Regarding time efficiency, peak-motifs outperforms all comparable tools by several orders of magnitude. We demonstrate its accuracy by analyzing data sets ranging from 4000 to 1 28 000 peaks for 12 embryonic stem cell-specific transcription factors. In all cases, the program finds the expected motifs and returns additional motifs potentially bound by cofactors. We further apply peak-motifs to discover tissue-specific motifs in peak collections for the p300 transcriptional co-activator. To our knowledge, peak-motifs is the only tool that performs a complete motif analysis and offers a user-friendly web interface without any restriction on sequence size or number of peaks.

Description: Genetic disorders can be detected by prenatal diagnosis using Chorionic Villus Sampling, but the 1:100 chance to result in miscarriage restricts the use to fetuses that are suspected to have an aberration. Detection of trisomy 21 cases noninvasively is now possible owing to the upswing of next-generation sequencing (NGS) because a small percentage of fetal DNA is present in maternal plasma. However, detecting other trisomies and smaller aberrations can only be realized using high-coverage NGS, making it too expensive for routine practice. We present a method, WISECONDOR (WIthin-SamplE COpy Number aberration DetectOR), which detects small aberrations using low-coverage NGS. The increased detection resolution was achieved by comparing read counts within the tested sample of each genomic region with regions on other chromosomes that behave similarly in control samples. This within-sample comparison avoids the need to re-sequence control samples. WISECONDOR correctly identified all T13, T18 and T21 cases while coverages were as low as 0.15–1.66. No false positives were identified. Moreover, WISECONDOR also identified smaller aberrations, down to 20 Mb, such as del(13)(q12.3q14.3), +i(12)(p10) and i(18)(q10). This shows that prevalent fetal copy number aberrations can be detected accurately and affordably by shallow sequencing maternal plasma. WISECONDOR is available at bioinformatics.tudelft.nl/wisecondor.

Description: Alternative splicing is the main mechanism governing protein diversity. The recent developments in RNA-Seq technology have enabled the study of the global impact and regulation of this biological process. However, the lack of standardized protocols constitutes a major bottleneck in the analysis of alternative splicing. This is particularly important for the identification of exon–exon junctions, which is a critical step in any analysis workflow. Here we performed a systematic benchmarking of alignment tools to dissect the impact of design and method on the mapping, detection and quantification of splice junctions from multi-exon reads. Accordingly, we devised a novel pipeline based on TopHat2 combined with a splice junction detection algorithm, which we have named FineSplice. FineSplice allows effective elimination of spurious junction hits arising from artefactual alignments, achieving up to 99% precision in both real and simulated data sets and yielding superior F 1 scores under most tested conditions. The proposed strategy conjugates an efficient mapping solution with a semi-supervised anomaly detection scheme to filter out false positives and allows reliable estimation of expressed junctions from the alignment output. Ultimately this provides more accurate information to identify meaningful splicing patterns. FineSplice is freely available at https://sourceforge.net/p/finesplice/ .

Description: Combinatorial interactions among transcription factors (TFs) are critical for integrating diverse intrinsic and extrinsic signals, fine-tuning regulatory output and increasing the robustness and plasticity of regulatory systems. Current knowledge about combinatorial regulation is rather limited due to the lack of suitable experimental technologies and bioinformatics tools. The rapid accumulation of ChIP-Seq data has provided genome-wide occupancy maps for a large number of TFs and chromatin modification marks for identifying enhancers without knowing individual TF binding sites. Integration of the two data types has not been researched extensively, resulting in underused data and missed opportunities. We describe a novel method for discovering frequent combinatorial occupancy patterns by multiple TFs at enhancers. Our method is based on probabilistic item set mining and takes into account uncertainty in both types of ChIP-Seq data. By joint analysis of 108 TFs in four human cell types, we found that cell–type-specific interactions among TFs are abundant and that the majority of enhancers have flexible architecture. We show that several families of transposable elements disproportionally overlap with enhancers with combinatorial patterns, suggesting that these transposable element families play an important role in the evolution of combinatorial regulation.

Description: Recent advances in high-throughput sequencing (HTS) technologies and computing capacity have produced unprecedented amounts of genomic data that have unraveled the genetics of phenotypic variability in several species. However, operating and integrating current software tools for data analysis still require important investments in highly skilled personnel. Developing accurate, efficient and user-friendly software packages for HTS data analysis will lead to a more rapid discovery of genomic elements relevant to medical, agricultural and industrial applications. We therefore developed Next-Generation Sequencing Eclipse Plug-in (NGSEP), a new software tool for integrated, efficient and user-friendly detection of single nucleotide variants (SNVs), indels and copy number variants (CNVs). NGSEP includes modules for read alignment, sorting, merging, functional annotation of variants, filtering and quality statistics. Analysis of sequencing experiments in yeast, rice and human samples shows that NGSEP has superior accuracy and efficiency, compared with currently available packages for variants detection. We also show that only a comprehensive and accurate identification of repeat regions and CNVs allows researchers to properly separate SNVs from differences between copies of repeat elements. We expect that NGSEP will become a strong support tool to empower the analysis of sequencing data in a wide range of research projects on different species.

Description: Next-generation sequencing (NGS) technologies-based transcriptomic profiling method often called RNA-seq has been widely used to study global gene expression, alternative exon usage, new exon discovery, novel transcriptional isoforms and genomic sequence variations. However, this technique also poses many biological and informatics challenges to extracting meaningful biological information. The RNA-seq data analysis is built on the foundation of high quality initial genome localization and alignment information for RNA-seq sequences. Toward this goal, we have developed RNASEQR to accurately and effectively map millions of RNA-seq sequences. We have systematically compared RNASEQR with four of the most widely used tools using a simulated data set created from the Consensus CDS project and two experimental RNA-seq data sets generated from a human glioblastoma patient. Our results showed that RNASEQR yields more accurate estimates for gene expression, complete gene structures and new transcript isoforms, as well as more accurate detection of single nucleotide variants (SNVs). RNASEQR analyzes raw data from RNA-seq experiments effectively and outputs results in a manner that is compatible with a wide variety of specialized downstream analyses on desktop computers.

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: 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: 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: 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: 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: 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: 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: 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: 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: 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 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 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: 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: 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).

Description: In this study, we present a new method for simulating the 3-D dynamic rupture process occurring on a non-planar fault. The method is based on the curved-grid finite-difference method (CG-FDM) proposed by Zhang & Chen and Zhang et al. to simulate the propagation of seismic waves in media with arbitrary irregular surface topography. While keeping the advantages of conventional FDM, that is computational efficiency and easy implementation, the CG-FDM also is flexible in modelling the complex fault model by using general curvilinear grids, and thus is able to model the rupture dynamics of a fault with complex geometry, such as oblique dipping fault, non-planar fault, fault with step-over, fault branching, even if irregular topography exists. The accuracy and robustness of this new method have been validated by comparing with the previous results of Day et al. , and benchmarks for rupture dynamics simulations. Finally, two simulations of rupture dynamics with complex fault geometry, that is a non-planar fault and a fault rupturing a free surface with topography, are presented. A very interesting phenomenon was observed that topography can weaken the tendency for supershear transition to occur when rupture breaks out at a free surface. Undoubtedly, this new method provides an effective, at least an alternative, tool to simulate the rupture dynamics of a complex non-planar fault, and can be applied to model the rupture dynamics of a real earthquake with complex geometry.

Description: Long-term volcanic subsidence provides insight into intereruptive processes, which comprise the longest portion of the eruptive cycle. Ground-based geodetic surveys of Medicine Lake Volcano (MLV), northern CA, document subsidence at rates of ~–10 mm yr –1 between 1954 and 2004. The long observation period plus the duration and stable magnitude of this signal presents an ideal opportunity to study long-term volcanic deformation, but this first requires accurate knowledge of the geometry and magnitude of the source. Best-fitting analytical source models to past levelling and GPS data sets show conflicting source parameters—primarily the model depth. To overcome this, we combine multiple tracks of InSAR data, each with a different look angle, to improve upon the spatial resolution of ground-based measurements. We compare the results from InSAR to those of past geodetic studies, extending the geodetic record to 2011 and demonstrating that subsidence at MLV continues at ~–10 mm yr –1 . Using geophysical inversions, we obtain the best-fitting analytical source model—a sill located at 9–10 km depth beneath the caldera. This model geometry is similar to those of past studies, providing a good fit to the high spatial density of InSAR measurements, while accounting for the high ratio of vertical to horizontal deformation derived from InSAR and recorded by existing levelling and GPS data sets. We discuss possible causes of subsidence and show that this model supports the hypothesis that deformation at MLV is driven by tectonic extension, gravitational loading, plus a component of volume loss at depth, most likely due to cooling and crystallization within the intrusive complex that underlies the edifice. Past InSAR surveys at MLV, and throughout the Cascades, are of variable success due to dense vegetation, snow cover and atmospheric artefacts. In this study, we demonstrate how InSAR may be successfully used in this setting by applying a suite of multitemporal analysis methods that account for atmospheric and orbital noise sources. These methods include: a stacking strategy based upon the noise characteristics of each data set; pixelwise rate-map formation (-RATE) and persistent scatterer InSAR (StaMPS).

Description: We present a numerical method to compute 3-D elastic waves in fully anisotropic axisymmetric media. This method is based on a decomposition of the wave equation into a series of uncoupled 2-D equations for which the dependence of the wavefield on the azimuth can be solved analytically. Four independent equations up to quadrupole order appear as solutions for moment-tensor sources located on the symmetry axis while single forces can be accommodated by two separate solutions up to dipole order. This decomposition gives rise to an efficient solution of the 3-D wave equation in a 2-D axisymmetric medium. First, we prove the validity of the decomposition of the wavefield in the presence of general anisotropy. Then we use it to derive the reduced 2-D equations of motions and discretize them using the spectral element method. Finally, we benchmark the numerical implementation for global wave propagation at 1 Hz and consider inner core anisotropy as an application for high-frequency wave propagation in anisotropic media at frequencies up to 2 Hz.

Description: We have developed a wavelet-based double-difference (DD) seismic tomography method. Instead of solving for the velocity model itself, the new method inverts for its wavelet coefficients in the wavelet domain. This method takes advantage of the multiscale property of the wavelet representation and solves the model at different scales. A sparsity constraint is applied to the inversion system to make the set of wavelet coefficients of the velocity model sparse. This considers the fact that the background velocity variation is generally smooth and the inversion proceeds in a multiscale way with larger scale features resolved first and finer scale features resolved later, which naturally leads to the sparsity of the wavelet coefficients of the model. The method is both data- and model-adaptive because wavelet coefficients are non-zero in the regions where the model changes abruptly when they are well sampled by ray paths and the model is resolved from coarser to finer scales. An iteratively reweighted least squares procedure is adopted to solve the inversion system with the sparsity regularization. A synthetic test for an idealized fault zone model shows that the new method can better resolve the discontinuous boundaries of the fault zone and the velocity values are also better recovered compared to the original DD tomography method that uses the first-order Tikhonov regularization.

Description: We show that it is possible to capture the oscillatory ground motion induced by the Tohoku-Oki event for periods ranging from 3 to 100 s using precise point positioning. We find that the ground motions of the sedimentary basins of Japan were large (respectively 〉0.15 m s –1 and 〉0.15 m s –2 for velocity and acceleration) even for periods larger than 3 s. We compare geodetic observables with a ground motion prediction equation designed for Japan seismicity and find that the spectral acceleration is well estimated for periods larger than 3 s and distances ranging from 100 to 500 km. At last, through the analysis of the displacement attenuation plots, we show that the 2011 Tohoku-Oki event is likely composed of multiple rupture patches as suggested before by time-reversal inversions of seismic data.

Description: In the literature, the inverted coseismic slip models from seismological and geodetic data for the 2011 Tohoku-Oki earthquake portray significant discrepancies, in particular regarding the intensity and the distribution of the rupture near the trench. For a megathrust earthquake, it is difficult to discern the slip along the shallow part of the fault from the geodetic data, which are often acquired on land. In this paper, we discuss the uncertainties in the slip distribution inversion using the geodetic data for the 2011 Tohoku earthquake and the Fully Bayesian Inversion method. These uncertainties are due to the prior information regarding the boundary conditions at the edges of the fault, the dip subduction angle and the smoothing operator. Using continuous GPS data from the Japan Island, the results for the rigid and free boundary conditions show that they produce remarkably different slip distributions at shallow depths, with the latter producing a large slip exceeding 30 m near the surface. These results indicate that the smoothing operator (gradient or Laplacian schemes) does not severely affect the slip pattern. To better invert the coseismic slip, we then introduce the ocean bottom GPS (OB-GPS) data, which improve the resolution of the shallow part of the fault. We obtain a near-trench slip greater than 40 m that reaches the Earth's surface, regardless of which boundary condition is used. Additionally, we show that using a mean dip angle for the fault as derived from subduction models is adequate if the goal is to invert for the general features of the slip pattern of this megathrust event.

Description: A geotechnical and geophysical campaign was performed at sites located in the alluvial plain of the river of Beirut (Lebanon), which is characterized by a significant lateral and vertical geological variability, along with anthropogenic disturbances in the first metres. The method combination has allowed detecting a shallow conductive low velocity layer of varying depth and thickness, corresponding to a soft clay layer embedded in coarser formations. This layer was found to exert strong control on the experimental dispersion curves (estimated from both active and passive experiments) characterized by a continuous mode superposition at high frequency, associated with an increase in phase velocity. Vs profiles in boreholes turned out to be of prime importance for adequately defining the parametrization before inversion and for ensuring the reliability of the inversion dispersive estimates at low frequency. A major output of this study is also to show that this low velocity layer, along with the strong shear velocity contrast at its bottom, significantly contributes to the site seismic response, and could make it difficult to use the measured H / V peak frequency as a proxy for the soil thickness over bedrock.

Description: We present an analysis of the M w = 5.3 earthquake that occurred in the Southeast Indian Ridge on 2010 February 11 using USArray data. The epicentre of this event is antipodal to the USArray, providing us with an opportunity to observe in details the antipodal focusing of seismic waves in space and time. We compare the observed signals with synthetic seismograms computed for a spherically symmetric earth model (PREM). A beamforming analysis is performed over the different seismic phases detected at antipodal distances. Direct spatial snapshots of the signals and the beamforming results show that the focusing is well predicted for the first P -wave phases such as PKP or PP . However, converted phases ( SKSP , PPS ) show a deviation of the energy focusing to the south, likely caused by the Earth's heterogeneity. Focusing of multiple S -wave phases strongly deteriorates and is barely observable.

Description: Global seismic tomography suffers from uncertainties in earthquake parameters routinely published in seismic catalogues. In particular, errors in earthquake location and origin-time may lead to strong biases in measured body wave delay-times and significantly pollute tomographic models. Common ways of dealing with this issue are to incorporate source parameters as additional unknowns into the linear tomographic equations, or to seek combinations of data to minimize the influence of source mislocations. We propose an alternative, physically-based method to desensitize direct S -wave delay-times to errors in earthquake location and origin-time. Our approach takes advantage of the fact that mislocation delay-time biases depend to first order on the earthquake-receiver azimuth, and to second order on the epicentral distance. Therefore, for every earthquake, we compute S -wave differential delay-times between optimized receiver pairs, such that a large part of their mislocation delay-time biases cancels out (for example origin-time fully subtracts out), while the difference of their sensitivity kernels remains sensitive to the model parameters of interest. Considering realistic, randomly distributed source mislocation vectors, as well as various levels of data noise and different synthetic Earths, we demonstrate that mislocation-related model errors are highly reduced when inverting for such differential delay-times, compared to absolute ones. The reduction is particularly rewarding for imaging the upper-mantle and transition zone. We conclude that using optimized receiver pairs is a suitable, low cost alternative to get rid of errors on earthquake location and origin-time for teleseismic direct S -wave traveltimes. Moreover, it can partly remove unilateral rupture propagation effects in cross-correlation delay-times, since they are similar to mislocation effects.

Description: As the Earth's inhomogeneous and viscoelastic properties, seismic signal attenuation we are trying to mitigate is a long-standing problem facing with high-resolution techniques. For addressing such a problem in the fields of time–frequency transform, Gabor transform methods such as atom-window method (AWM) and molecular window method (MWM) have been reported recently. However, we observed that these methods might be much better if we partition the non-stationary seismic data into adaptive stationary segments based on the amplitude and frequency information of the seismic signal. In this study, we present a new method called amplitude-frequency partition (AFP) to implement this process in the time–frequency domain. Cases of a synthetic and field seismic data indicated that the AFP method could partition the non-stationary seismic data into stationary segments approximately, and significantly, a high-resolution result would be achieved by combining the AFP method with conventional spectral-whitening method, which could be considered superior to previous resolution-enhancement methods like time-variant spectral whitening method, the AWM and the MWM as well. This AFP method presented in this study would be an effective resolution-enhancement tool for the non-stationary seismic data in the fields of an adaptive time–frequency transform.

Description: Full waveform inversion (FWI) of 3-D data sets has recently been possible thanks to the development of high performance computing. However, FWI remains a computationally intensive task when high frequencies are injected in the inversion or more complex wave physics (viscoelastic) is accounted for. The highest computational cost results from the numerical solution of the wave equation for each seismic source. To reduce the computational burden, one well-known technique is to employ a random linear combination of the sources, rather that using each source independently. This technique, known as source encoding, has shown to successfully reduce the computational cost when applied to real data. Up to now, the inversion is normally carried out using gradient descent algorithms. With the idea of achieving a fast and robust frequency-domain FWI, we assess the performance of the random source encoding method when it is interfaced with second-order optimization methods (quasi-Newton l -BFGS, truncated Newton). Because of the additional seismic modelings required to compute the Newton descent direction, it is not clear beforehand if truncated Newton methods can indeed further reduce the computational cost compared to gradient algorithms. We design precise stopping criteria of iterations to fairly assess the computational cost and the speed-up provided by the source encoding method for each optimization method. We perform experiment on synthetic and real data sets. In both cases, we confirm that combining source encoding with second-order optimization methods reduces the computational cost compared to the case where source encoding is interfaced with gradient descent algorithms. For the synthetic data set, inspired from the geology of Gulf of Mexico, we show that the quasi-Newton l -BFGS algorithm requires the lowest computational cost. For the real data set application on the Valhall data, we show that the truncated Newton methods provide the most robust direction of descent.

Description: We use the Rayleigh integral (RI) as an approximation to the Helmholtz–Kirchoff integral to model infrasound generation and propagation from underground chemical explosions at distances of 250 m out to 5 km as part of the Source Physics Experiment (SPE). Using a sparse network of surface accelerometers installed above ground zero, we are able to accurately create synthetic acoustic waveforms and compare them to the observed data. Although the underground explosive sources were designed to be symmetric, the resulting seismic wave at the surface shows an asymmetric propagation pattern that is stronger to the northeast of the borehole. This asymmetric bias may be attributed to the subsurface geology and faulting of the area and is observed in the acoustic waveforms. We compare observed and modelled results from two of the underground SPE tests with a sensitivity study to evaluate the asymmetry observed in the data. This work shows that it is possible to model infrasound signals from underground explosive sources using the RI and that asymmetries observed in the data can be modelled with this technique.

Description: The Terre Adélie and George V Land (East Antarctica) represent key areas for understanding tectonic relationships between terranes forming the Neoarchean-Palaeoproterozoic Terre Adélie Craton (TAC) and the neighbouring lithospheric blocks, together with the nature of its boundary. This region that represents the eastern border of the TAC is limited on its eastern side by the Mertz shear zone (MSZ) separating more recent Palaeozoic units from the craton. The MSZ, that recorded dextral strike-slip movement at 1.7 and 1.5 Ga, is likely correlated with the Kalinjala or Coorong shear zone in South Australia, east of the Gawler Craton and may therefore represent a frozen lithospheric-scale structure. In order to investigate the lithospheric structure of the TAC and the MSZ, we deployed from 2009 October to 2011 October four temporary seismic stations, which sampled the various lithospheric units of the TAC and of the neighbouring Palaeozoic block, together with the MSZ. We used receiver function method to deduce Moho depths and seismic anisotropy technique to infer the upper mantle deformation. Results from receiver functions analysis reveal Moho at 40–44 km depth beneath the TAC, at 36 km under the MSZ and at 28 km beneath the eastern Palaeozoic domain. The MSZ therefore delimits two crustal blocks of different thicknesses with a vertical offset of the Moho of 12 km. Seismic anisotropy deduced from SKS splitting at stations on the TAC shows fast polarisation directions () trending E–W, that is, parallel to the continental margin, and delay times ( t ) ranging from 0.8 to 1.6 s. These results are similar to the splitting parameters observed at the permanent GEOSCOPE Dumont D'Urville station (DRV: 95°N, t 1.1 s) located in the Palaeoproterozoic domain of TAC. On the MSZ, the small number of good quality measurements limits the investigation of the deep signature of the shear zone. However, the station in the Palaeozoic domain shows trending N60°E, which is significantly different to the trending measurements from stations on the TAC, suggesting that the MSZ may also represent a major frontier between the Neoarchean-Palaeoproterozoic and Palaeozoic terranes.

Description: The seismic investigation of the lowermost mantle is in many places hampered by the lack of suitable source–receiver combinations that sample the D '' region and have to meet the requirements of a suitable epicentral distance range. The low velocity regions beneath the central Pacific and Atlantic Oceans in particular have been sampled in fewer places than circum Pacific regions. In this study, we use data from two recent ocean bottom seismometer (OBS) deployments for the Plume-Lithosphere Undersea Mantle Experiment (PLUME) around Hawaii to increase the coverage of the lower mantle with reflected P waves. Through stacking of the data we achieve significant reduction in noise levels. The most favourable epicentral distances to detect D '' reflections are around 70–79°. Most of our source–receiver combinations have distances less than that, thereby limiting the number of candidate observed reflections. Nevertheless, using array methods, we are able to test approximately 70 events for arrivals with slowness values and arrival times that would be consistent with a top-side reflection off a hypothetical D '' structure (PdP wave). Modelling these data with a 1-D reflectivity method, we identify a few places of detectable PdP waves, for which the velocity contrast in P - and S -wave velocity across the D '' reflector have to be relatively large (around 3–5 per cent increase and decrease, respectively) compared to other regions (e.g. beneath the Caribbean or Eurasia where the contrast is closer to 1–2 per cent). For larger distance ranges, smaller velocity contrasts are sufficient to cause observable reflections. This study shows that, despite the possible dominance of microseisms on OBS records, it is possible to use relatively short-period waves, with dominant periods as short as 3–7 s. Our findings suggest that, with future such deployments, OBS deployments will help to extend D '' studies to previously unmapped regions.

Description: In the context of elastic wave propagation, the non-periodic homogenization asymptotic method allows to find a smooth effective medium and equations that correspond to the wave propagation in a given complex elastic or acoustic medium down to a given minimum wavelength. By smoothing all discontinuities and fine scales of the original medium, the homogenization technique considerably reduces meshing difficulties as well as the numerical cost associated with the wave equation solver, while producing the same waveform as for the original medium (up to the wanted accuracy). We present here a variation of the original method, allowing to homogenize the difference, or residual, between an original medium and a reference medium. This makes it possible to, for example, homogenize some specific parts of a model or to leave unchanged a specific interface while homogenizing the rest of the model. We present two examples of applications, one implying a complex geological shallow structure and the other involving the combination of deterministic and stochastic elastic models.

Description: Converted phase (CP) elastic seismic signals are comparable in amplitude to the primary signals recorded at large offsets and have the potential to be used in seismic imaging and velocity analysis. We present an approach for CP elastic wave equation velocity analysis that does not use source information and is applicable to surface-seismic, microseismic, teleseismic and vertical seismic profile (VSP) studies. Our approach is based on the cross-correlation between reflected or transmitted PP and CP PS (and/or SS and CP PS) waves propagated backward in time, and is formulated as an optimization problem with a differential semblance criterion objective function for the simultaneous update of both P - and S -wave velocity models. The merit of this approach is that it is fully data-driven, uses full waveform information, and requires only one elastic backward propagation to form an image rather than the two (one forward and one backward) propagations needed for standard reverse-time migration. Moreover, as the method does not require forward propagation, it does not suffer from migration operator source aliasing when a small number of shots are used. We present a derivation of the method and test it with a synthetic model and field micro-seismic data.

Description: We examine P -wave velocity structure at the base of the mantle beneath the western Pacific, near the western edge of the Pacific Large-Low Shear Velocity Province (LLSVP), using high-quality seismograms provided by a large-scale mobile broad-band seismic observation in northeastern China (the NECESSArray project). Forward modelling using the reflectivity method is conducted to explain the variation of P -wave traveltimes as a function of epicentral distance near the core shadow zone. Additionally, PcP–P traveltimes are examined to enlarge the survey area. As a result, a rapid variation of P -wave velocity is detected at the base of the mantle. Regions of thin (20–50 km thick) and low velocity (–2 to –5 per cent) layers at the base of the mantle are intersected by an 80-km-thick region with a high velocity (+2 per cent). A slightly fast region exists at the northwest of the region with the thin low-velocity layer. These layers are typically separated by several hundred kilometres and would be difficult to explain by thermal effects alone. These observations suggest that very complicated thermochemical reactions occur near the edge of the Pacific LLSVP.

Description: Complications arise in the interpretation of gravity fields because of interference from systematic degradations, such as boundary blurring and distortion. The major sources of these degradations are the various systematic errors that inevitably occur during gravity field data acquisition, discretization and geophysical forward modelling. To address this problem, we evaluate deconvolution method that aim to detect the clear horizontal boundaries of anomalous sources by the suppression of systematic errors. A convolution-based multilayer projection model, based on the classical 3-D gravity field forward model, is innovatively derived to model the systematic error degradation. Our deconvolution algorithm is specifically designed based on this multilayer projection model, in which three types of systematic error are defined. The degradations of the different systematic errors are considered in the deconvolution algorithm. As the primary source of degradation, the convolution-based systematic error is the main object of the multilayer projection model. Both the random systematic error and the projection systematic error are shown to form an integral part of the multilayer projection model, and the mixed norm regularization method and the primal-dual optimization method are therefore employed to control these errors and stabilize the deconvolution solution. We herein analyse the parameter identification and convergence of the proposed algorithms, and synthetic and field data sets are both used to illustrate their effectiveness. Additional synthetic examples are specifically designed to analyse the effects of the projection systematic error, which is caused by the uncertainty associated with the estimation of the impulse response function.

Description: We successfully apply the semi-global inverse method of simulated annealing to determine the best-fitting 1-D anisotropy model for use in acoustic frequency domain waveform tomography. Our forward problem is based on a numerical solution of the frequency domain acoustic wave equation, and we minimize wavefield phase residuals through random perturbations to a 1-D vertically varying anisotropy profile. Both real and synthetic examples are presented in order to demonstrate and validate the approach. For the real data example, we processed and inverted a cross-borehole data set acquired by Vale Technology Development (Canada) Ltd. in the Eastern Deeps deposit, located in Voisey's Bay, Labrador, Canada. The inversion workflow comprises the full suite of acquisition, data processing, starting model building through traveltime tomography, simulated annealing and finally waveform tomography. Waveform tomography is a high resolution method that requires an accurate starting model. A cycle-skipping issue observed in our initial starting model was hypothesized to be due to an erroneous anisotropy model from traveltime tomography. This motivated the use of simulated annealing as a semi-global method for anisotropy estimation. We initially tested the simulated annealing approach on a synthetic data set based on the Voisey's Bay environment; these tests were successful and led to the application of the simulated annealing approach to the real data set. Similar behaviour was observed in the anisotropy models obtained through traveltime tomography in both the real and synthetic data sets, where simulated annealing produced an anisotropy model which solved the cycle-skipping issue. In the real data example, simulated annealing led to a final model that compares well with the velocities independently estimated from borehole logs. By comparing the calculated ray paths and wave paths, we attributed the failure of anisotropic traveltime tomography to the breakdown of the ray-theoretical approximation in the vicinity of strong velocity discontinuities.

Description: Seismic traveltimes and their spatial derivatives are the basis of many imaging methods such as pre-stack depth migration and tomography. A common approach to compute these quantities is to solve the eikonal equation with a finite-difference scheme. If many recently published algorithms for resolving the eikonal equation do now yield fairly accurate traveltimes for most applications, the spatial derivatives of traveltimes remain very approximate. To address this accuracy issue, we develop a new hybrid eikonal solver that combines a spherical approximation when close to the source and a plane wave approximation when far away. This algorithm reproduces properly the spherical behaviour of wave fronts in the vicinity of the source. We implement a combination of 16 local operators that enables us to handle velocity models with sharp vertical and horizontal velocity contrasts. We associate to these local operators a global fast sweeping method to take into account all possible directions of wave propagation. Our formulation allows us to introduce a variable grid spacing in all three directions of space. We demonstrate the efficiency of this algorithm in terms of computational time and the gain in accuracy of the computed traveltimes and their derivatives on several numerical examples.

Description: Our study compares the seismic properties between the flat and normal subduction regions in central Chile, to better understand the links between the slab geometry, surface deformation and the deeper structures. In comparison with previous studies, we show the most complete 3-D regional seismic tomography images for this region, in which we use (1) a larger seismic data set compiled from several short-term seismic catalogues, (2) a denser seismic array allowing a better resolution of the subduction zone from the trench to the backarc and into the upper ~30 km of the slab and (3) a starting 1-D background velocity model specifically calculated for this region and refined over the years. We assess and discuss our tomography results using regional seismic attenuation models and estimating rock types on the basis of pressure and temperature conditions computed from thermomechanical models. Our results show significant seismic differences between the flat and normal subduction zones. As expected, the faster seismic velocities and increased seismicity within the flat slab and overriding lithosphere are generally consistent with a cooler thermal state. Our results are also consistent with dehydration of the mantle above the subducted Juan Fernandez Ridge at the eastern tip of the flat slab segment, indicating that the latter retains some fluids during subduction. However, fluids in the upper portion of the flat slab segment are not seismically detected, since we report instead fast slab seismic velocities which contradict the argument of its buoyancy being the cause of horizontal subduction. The forearc region, above the flat slab, exhibits high Vs and very low Vp / Vs ratios, uncorrelated with typical rock compositions, increased density or reduced temperature; this feature is possibly linked with the aftershock effects of the M w7.1 1997 Punitaqui earthquake, the flat slab geometry and/or seismic anisotropy. At the surface, the seismic variations correlate with the geological terranes. The Andean crust is strongly reduced in seismic velocities along the La Ramada–Aconcagua deformation belt, suggesting structural damage. Slow seismic velocities along the Andean Moho match non-eclogitized hydrated rocks, consistent with a previous delamination event or a felsic composition, which in turn supports the extent of the Chilenia terrane at these depths. We confirm previous studies that suggest that the Cuyania terrane in the backarc region is mafic and contains an eclogitized lower crust below 50-km depth. We also hypothesize major Andean basement detachment faults (or shear zones) to extend towards the plate interface and canalize slab-derived fluids into the continental crust.

Description: The most common type of waves used for probing anisotropy of rocks in laboratory is the direct P wave. Information potential of the measured P -wave velocity, however, is limited. In rocks displaying weak triclinic anisotropy, the P -wave velocity depends just on 15 linear combinations of 21 elastic parameters, called the weak-anisotropy parameters. In strong triclinic anisotropy, the P -wave velocity depends on the whole set of 21 elastic parameters, but inversion for six of them is ill-conditioned and these parameters are retrieved with a low accuracy. Therefore, in order to retrieve the complete elastic tensor accurately, velocities of S waves must also be measured and inverted. For this purpose, we developed a lab facility which allows the P - and S -wave ultrasonic sounding of spherical rock samples in 132 directions distributed regularly over the sphere. The velocities are measured using a pair of P -wave sensors with the transmitter and receiver polarized along the radial direction and using two pairs of S -wave sensors with the transmitter and receiver polarized tangentially to the spherical sample in mutually perpendicular directions. We present inversion methods of phase and ray velocities for elastic parameters describing general triclinic anisotropy. We demonstrate on synthetic tests that the inversion becomes more robust and stable if the S -wave velocities are included. This applies even to the case when the velocity of the S waves is measured in a limited number of directions and with a significantly lower accuracy than that of the P wave. Finally, we analyse velocities measured on a rock sample from the Outokumpu deep drill hole, Finland. We present complete sets of elastic parameters of the sample including the error analysis for several levels of confining pressure ranging from 0.1 to 70 MPa.

Description: We use normal-mode splitting functions in addition to surface wave phase anomalies, body wave traveltimes and long-period waveforms to construct a 3-D model of anisotropic shear wave velocity in the Earth's mantle. Our modelling approach inverts for mantle velocity and anisotropy as well as transition-zone discontinuity topographies, and incorporates new crustal corrections for the splitting functions that are consistent with the non-linear corrections we employ for the waveforms. Our preferred anisotropic model, S362ANI+M, is an update to the earlier model S362ANI, which did not include normal-mode splitting functions in its derivation. The new model has stronger isotropic velocity anomalies in the transition zone and slightly smaller anomalies in the lowermost mantle, as compared with S362ANI. The differences in the mid- to lowermost mantle are primarily restricted to features in the Southern Hemisphere. We compare the isotropic part of S362ANI+M with other recent global tomographic models and show that the level of agreement is higher now than in the earlier generation of models, especially in the transition zone and the lower mantle. The anisotropic part of S362ANI+M is restricted to the upper 300 km in the mantle and is similar to S362ANI. When radial anisotropy is allowed throughout the mantle, large-scale anisotropic patterns are observed in the lowermost mantle with v SV  〉  v SH beneath Africa and South Pacific and v SH  〉  v SV beneath several circum-Pacific regions. The transition zone exhibits localized anisotropic anomalies of ~3 per cent v SH  〉  v SV beneath North America and the Northwest Pacific and ~2 per cent v SV  〉  v SH beneath South America. However, small improvements in fits to the data on adding anisotropy at depth leave the question open on whether large-scale radial anisotropy is required in the transition zone and in the lower mantle. We demonstrate the potential of mode-splitting data in reducing the trade-offs between isotropic velocity and anisotropy in the lowermost mantle for the even-degree variations. Spurious anisotropic variations in the mid-mantle are also suppressed with the addition of mode-splitting data.

Description: We describe a methodology for quantitatively characterizing the fractured nature of a hydrocarbon or geothermal reservoir from surface seismic data under a Bayesian inference framework. The method combines different kinds of measurements of fracture properties to find a best-fitting model while providing estimates of the uncertainty of model parameters. Fractures provide pathways for fluid flow in a reservoir, and hence, knowledge about a reservoir's fractured nature can be used to enhance production from the reservoir. The fracture properties of interest in this study (to be inferred) are fracture orientation and excess compliance, where each of these properties are assumed to vary spatially over a 2-D horizontal grid which is assumed to represent the top of a reservoir. The Bayesian framework in which the inference problem is cast has the key benefits of (1) utilization of a prior model that allows geological information to be incorporated, (2) providing a straightforward means of incorporating all measurements (across the 2-D spatial grid) into the estimates at each gridpoint, (3) allowing different types of measurements to be combined under a single inference procedure and (4) providing a measure of uncertainty in the estimates. The observed data are taken from a 2-D array of surface seismic receivers responding to an array of surface sources. Well understood features from the seismic traces are extracted and treated as the observed data, namely the P -wave reflection amplitude variation with acquisition azimuth and offset (amplitude versus azimuth data) and fracture transfer function (FTF) data. Amplitude versus azimuth data are known to be more sensitive to fracture properties when the fracture spacing is significantly smaller than the seismic wavelength, whereas FTF data are more sensitive to fracture properties when the fracture spacing is on the order of the seismic wavelength. Combining these two measurements has the benefit of allowing inferences to be made about fracture properties over a larger range of fracture spacing than otherwise attainable. Geophysical forward models for the measurements are used to arrive at likelihood models for the data. The prior distribution for the fracture variables is obtained by defining a Markov random field over the lateral 2-D grid where we wish to obtain fracture properties, where this method for defining the prior has the added benefit of allowing for non-stationarity in the resulting model covariance. The fracture variables are then inferred by application of loopy belief propagation to yield approximations for the posterior marginal distributions of the fracture properties, as well as the maximum a posteriori and Bayes least-squares (posterior mean) estimates of these properties. Verification of the inference procedure is performed using a synthetic data set, where the estimates obtained are shown to be at or near ground truth for the full range of fracture spacings for fracture orientation and at low fracture spacings for excess compliance estimates.

Description: In this study, new intensity-distance attenuation law is presented, using directly the intensity observations, rather than the subjective, and sometimes controversial, isoseismal lines. This intensity-distance attenuation law is the only one defined for Portugal mainland, which is expressed as a function of magnitude. We computed this attenuation law using the slope and the intercept of the logarithmic regression of 25 events, with magnitudes between 4.4 and 6.2. Using the Bakun and Wentworth method (1997), this new attenuation law allows performing better results in the earthquake epicentral position and magnitude estimations of the 1909 Benavente event than the Atkinston and Boore attenuation law (1997). This law also gives good results in the study of site effects, presenting good matches between intensity residuals and geological structures where site effects are expected.

Description: We infer 3-D localized shear velocity structure in the lowermost 400 km of the mantle at the western edge of the Pacific large low shear velocity province (LLSVP) by applying waveform inversion to transverse component body-wave waveforms from the F-net seismic array in Japan. Our data set consists of relatively long period (12.5–200 s) broad-band seismic waveforms of Tonga-Fiji deep focus and intermediate deep earthquakes. We conduct several tests to confirm the robustness of the inversion results. We find two low-velocity zones at the bottom of the target region, with a high-velocity zone in the middle, and a low-velocity zone above the high-velocity zone and contiguous with the two deeper low-velocity zones at a depth of 200–300 km above the core-mantle boundary (CMB). This supports the idea that the Pacific LLSVP may be an aggregation of small upwelling plumes rather than a single large thermochemical pile.

Description: Samples of shales from the Ordovician Bongabinni and Goldwyer source rock formations were recovered from the Canning Basin (Western Australia). Attenuation was experimentally measured on preserved plugs from these formations in the frequency range between 10 –2 and 10 2 Hz. Samples cored with different orientations with respect to the sedimentary bedding were prepared and tested in their native saturated state and after drying in the oven at 105 °C for 24 hr to assess the effect of fluids and of the sediment anisotropy on attenuation. To aid the interpretation of the experimental results, the clay-rich samples were characterized in terms of mineralogy, water content, porosity, permeability and microstructure. The two shales have significantly different quality factors; and this is seen to be dependent on both the saturation state of the samples and the propagation direction of the oscillatory signal. The attenuation coefficient for compression/extension parallel to bedding is less than that vertical to bedding in both the preserved and partially dehydrated situations. No frequency dependency is observed in the preserved samples within the range of frequencies explored in this study. On the other hand partially saturated samples show peaks in attenuation at around 40 Hz when the stress perturbation is transmitted normal to the macroscopic bedding. The interpretation of the attenuation measurements in terms of well-established theoretical models is discussed in view of the physical characteristics and microstructure of the tested rocks.

Description: We examine seismic attenuation in the northern South Island, New Zealand, where subduction transitions to oblique collision, and plate motion occurs along multiple crustal faults overlying the subducted slab. The 3-D inversions derive Q (1/attenuation) using path attenuation t * observations for 334 distributed earthquakes recorded on permanent and temporary stations, including both velocity and acceleration records. A 2.5 s window was used for P spectra, but for S spectra longer varied lengths were selected around the predicted S arrival, using the 5–95 per cent energy integral. The Q results highlight many aspects of the structure more clearly than previously derived seismic velocity models, including a high Q slab, low Q basins, moderately low Q active fault regions, and thick lithosphere. Qs tends to be greater than Qp , except in low Q shallow upper crust. The mantle above the slab does not exhibit low Q , unlike mantle to the north beneath the volcanic region of central North Island, and is inferred to be cool and stagnant with some vertical flux of slab dehydration fluid. In the brittle crust, low Q is imaged along those faults with most recent seismicity and may be related to distributed microfractures. In the ductile crust of the greywacke terranes, zones of low Q under the faults are attributed to localized ductile deformation with high strain-rate and grain size reduction, consistent with numerical models showing the development of enhanced strain-rate zones above the strong underlying slab. In contrast the Christchurch region has no ductile lower crust and instead has high Q indicative of strong mafic rocks at 12 km depth.

Description: Numerical simulation of repeated occurrences of slip events on a fault patch (asperity) is used to interpret the mechanism of irregular sequences of slip events. The fault is uniformly shear loaded at a constant rate, and the frictional stress acting on the fault is assumed to obey a rate- and state-dependent friction (RSF) law. A circular patch with velocity-weakening frictional property is embedded in the fault, which apart from this has velocity-strengthening frictional property. The numerical simulations are conducted using various characteristic slip distances L of the RSF law. For small values of L seismic slip events (earthquakes) repeatedly occur at regular intervals. With increasing L , the recurrence of slip events becomes more complex. A period doubled slip pattern, where seismic and aseismic slip events alternately occur, multiperiodic patterns and aperiodic patterns occur. At the same time, slip tends to become aseismic with increasing L . The distributions of shear stress on the fault just before slip events are variable because of variation in the residual stress of the preceding slip event and the shear stress generated by aseismic sliding during interseismic periods. These variations in shear stress cause the complex sequence of slip events seen here. An iteration map of the recurrence intervals of slip events for an aperiodic sequence of slip events is expressed by a simple curve, indicating that the timing of an event is predictable from the previous time interval, and the sequence of slip events exhibits deterministic chaos. To help interpret these results for a sequence of slip events on a velocity-weakening patch embedded in a velocity-strengthening region, a numerical simulation is conducted of slip on a velocity-weakening patch enclosed by a permanently locked region. In this case, no complex recurrence of slip events is observed. When L is less than a critical value, seismic slip events repeatedly occur at a constant interval. Stable sliding occurs when L is larger than the critical value. This result indicates that the complex slip behaviour seen for a velocity-weakening patch embedded in a velocity-strengthening region is caused by the interaction between the velocity-weakening and velocity-strengthening regions.

Description: We carry out a joint inversion of surface wave dispersion curves and teleseismic shear wave arrival times across the Tibetan collision zone, from just south of the Himalaya to the Qaidam Basin at the northeastern margin of the plateau, and from the surface to 600 km depth. The surface wave data consist of Rayleigh-wave group dispersion curves, mainly in the period range from 10 to 70 s, with a maximum of 2877 source–receiver pairs. The body wave data consist of more than 8000 S -wave arrival times recorded from 356 telesesmic events. The tomographic images show a ‘wedge’ of fast seismic velocities beneath central Tibet that starts underneath the Himalaya and reaches as far as the Bangong–Nujiang Suture (BNS). In our preferred interpretation, in central Tibet the Indian lithosphere underthrusts the plateau to approximately the BNS, and then subducts steeply. Further east, Indian lithosphere appears to be subducting at an angle of ~45°. We see fast seismic velocities under much of the plateau, as far as the BNS in central Tibet, and as far as the Xiangshuihe–Xiaojiang Fault in the east. At 150 km depth, the fast region is broken by an area ~300 km wide that stretches from the northern edge of central Tibet southeastwards as far as the Himalaya. We suggest that this gap, which has been observed previously by other investigators, represents the northernmost edge of the Indian lithosphere, and is a consequence of the steepening of the subduction zone from central to eastern Tibet. This also implies that the fast velocities in the northeast have a different origin, and are likely to be caused by lithospheric thickening or small-scale subduction of Asian lithosphere. Slow velocities observed to the south of the Qaidam suggest that the basin is not subducting. Finally, we interpret fast velocities below 400 km as subducted material from an earlier stage of the collision that has stalled in the transition zone. Its position to the south of the present subduction is likely to be due to the relative motion of India to the northeast.

Description: Probabilistic seismic hazard models (PSHM) are used for quantifying the seismic hazard at a site or a grid of sites. In this study, a methodology is proposed to compare the distribution of the expected number of sites with exceedance with the observed number considering an acceleration threshold at a set of recording sites. The method is applied to France and Turkey. The French accelerometric database is checked to produce a reliable accelerometric data set. In addition, we also used a synthetic data set inferred from an instrumental catalogue combined with a ground-motion prediction equation. The results show that the MEDD2002 and AFPS2006 PSH models overestimate the number of sites with exceedance for low acceleration levels (below 40 cm s –2 ) or short return periods (smaller than 50 yr for AFPS2006 and 475 yr for MEDD2002). For larger acceleration levels, there are few observations and none of the models is rejected. In Turkey, the SHARE hazard estimates can be tested against ground-motion levels of interest in earthquake engineering. As the completeness issue is crucial, the recorded data at each station is analysed to detect potential gaps in the recording. As most accelerometric stations are located on soil, accelerations at rock are estimated using a site-amplification model. Different minimum intersite distances and station configurations are considered. The observed numbers of sites with exceedance are well within the bounds of the predicted distribution for accelerations between 103 and 397 cm s –2 . For higher levels, both the observed number and the predicted percentile 2.5 are zero, and no conclusion can be drawn.

Description: This study presents a depth inversion of Scholte wave group and phase velocity maps obtained from cross-correlation of 6.5 hr of noise data from the Valhall Life of Field Seismic network. More than 2 600 000 vertical–vertical component cross-correlations are computed from the 2320 available sensors, turning each sensor into a virtual source emitting Scholte waves. We used a traditional straight-ray surface wave tomography to compute the group velocity map. The phase velocity maps have been computed using the Eikonal tomography method. The inversion of these maps in depth are done with the Neighbourhood Algorithm. To reduce the number of free parameters to invert, geological a priori information are used to propose a power-law 1-D velocity profile parametrization extended with a gaussian high-velocity layer where needed. These parametrizations allowed us to create a high-resolution 3-D S -wave model of the first 600 m of the Valhall subsurface and to precise the locations of geological structures at depth. These results would have important implication for shear wave statics and monitoring of seafloor subsidence due to oil extraction. The 3-D model could also be a good candidate for a starting model used in full-waveform inversions.

Description: Cross-correlation and cross-spectral time delays often exhibit strong outliers due to ambiguities or cycle jumps in the correlation function. Their number increases when signal-to-noise, signal similarity or spectral bandwidth decreases. Such outliers heavily determine the time-delay probability density function and the results of further computations (e.g. double-difference location and tomography) using these time delays. In the present research we expressed cross-correlation as a function of the squared difference between signal amplitudes and show that they are closely related. We used this difference as a cost function whose minimum is reached when signals are aligned. Ambiguities may be removed in this function by using a priori information. We propose using the traveltime difference as a priori time-delay information. By modelling the probability density function of the traveltime difference by a Cauchy distribution and the probability density function of the data (differences of seismic signal amplitudes) by a Laplace distribution we were able to find explicitly the time-delay a posteriori probability density function. The location of the maximum of this a posteriori probability density function is the maximum a p osteriori time-delay estimation for earthquake signals. Using this estimation to calculate time delays for earthquakes on the south flank of Kilauea statistically improved the cross-correlation time-delay estimation for these data and resulted in successful double-difference relocation for an increased number of earthquakes. This robust time-delay estimation improves the spatiotemporal resolution of seismicity rates in the south flank of Kilauea.

Description: Fracture shear-dilatancy is an essential process for enhancing the permeability of deep geothermal reservoirs, and is usually accompanied by the radiation of seismic waves. However, the hazard and risk perspective of induced seismicity research typically focuses only on the question of how to reduce the occurrence of induced earthquakes. Here we present a quantitative analysis of seismic hazard as a function of the two key factors defining an enhanced geothermal system: The permeability enhancement, and the size of the stimulated reservoir. Our model has two coupled components: (1) a pressure diffusion model and (2) a stochastic seismicity model. Permeability is increased in the source area of each induced earthquake depending on the amount of slip, which is determined by the magnitude. We show that the few largest earthquakes (i.e. 5–10 events with M ≥ 1.5) contribute more than half of the total reservoir stimulation. The results further indicate that planning and controlling of reservoir engineering operations may be compromised by the considerable variability of maximum observed magnitude, reservoir size, the Gutenberg–Richter b -value and Shapiro's seismogenic index (i.e. a measure of seismic reactivity of a reservoir) that arises from the intrinsic stochastic nature of induced seismicity. We also find that injection volume has a large impact on both reservoir size and seismic hazard. Injection rate and injection scheme have a negligible effect. The impact of site-specific parameters on seismicity and reservoir properties is greater than that of the injected volume. In particular, conditions that lead to high b -values—possibly a low differential stress level—have a high impact on seismic hazard, but also reduce the efficiency of the stimulation in terms of permeability enhancement. Under such conditions, target reservoir permeability can still be achieved without reaching an unacceptable level of seismic hazard, if either the initial reservoir permeability is high or if several fractures are stimulated. The proposed methodology is a first step towards including induced seismic hazard analysis into the design of reservoir stimulation in a quantitative and robust manner.

Description: Three large earthquakes occurred in the Ashikule stepover zone, NW Tibet in 2008 ( M w 7.1), 2012 ( M w 6.2) and 2014 ( M w 6.9). In this paper we use InSAR data to examine the 2008 event in detail, and place it in the context of both recent/possible future events and regional tectonic environment. We present InSAR data covering the coseismic and post-seismic phases of the 2008 event, and invert for slip distributions during both phases in order to examine the spatial relationship between earthquake rupture and afterslip. To account for the curved rupture trace as mapped in the field, we model the slip on a non-planar fault. A change in slip sense from left-lateral in the north to extensional in the south is consistent with the focal mechanism of aftershocks. We find that afterslip occurs around the perimeter of the coseismic slip, although in a patchy way, as seen in other normal faulting case studies. We calculate the Coulomb stress changes imparted by the 2008 earthquake on neighbouring faults, including stress loading on the subsequent rupture zones of the 2012 and 2014 earthquakes. We also compute the combined loading effect of the 2008, 2012 and 2014 earthquakes on the nearby Altyn Tagh Fault, which hasn't experienced a major earthquake since early last century. This study demonstrates that stepover zone failure can increase stresses on bounding strike-slip faults, which likely has relevance in analogous tectonic environments elsewhere in the world.

Description: Source–receiver interferometric imaging can be used to synthesize a subsurface acoustic or elastic image, consisting of a zero-time, zero-offset response (or Green's function) between a colocated pseudo-source and pseudo-receiver placed at each point in the subsurface image. However, if the imaging process does not properly account for multiple reflections, and enclosing boundaries of sources and receivers are not available, the image shows artefacts, poorly illuminated areas and distorted image amplitudes. Here we demonstrate with numerical examples that two-sided non-linear imaging provides the best elastic pure-mode ( PP and SS ) and converted-mode ( PS ) images, having higher resolution and more uniform illumination than those obtained from both one-sided linear imaging and from other intermediate steps of imaging (e.g. non-linear one-sided, linear two-sided). We also propose practical approaches to construct the additional fields required by two-sided non-linear imaging without the need for a detailed velocity model and receivers (and/or sources) in the subsurface. Moreover, when conversions are used for imaging, ‘true-amplitude’ images (here true-amplitude means properly retrieving amplitudes that represent the zero-time, zero-offset elastic response) should theoretically vanish because neither P -to- S or S -to- P conversions arise at zero-time and zero-offset. Applying a correction procedure that accounts for the polarity reversal in PS (or SP ) single-shot images helps with their structural interpretation but results in an unphysical estimate of the subsurface response and uninterpretable amplitudes. This suggests that there are advantages in exploiting pure-mode SS reflections/transmissions, in addition to converted waves only, because they require no polarity correction and the resulting image contains meaningful amplitudes that are proportional to the local shear-wave properties of the medium.

Description: We present a time-independent gridded earthquake rate forecast for the European region including Turkey. The spatial component of our model is based on kernel density estimation techniques, which we applied to both past earthquake locations and fault moment release on mapped crustal faults and subduction zone interfaces with assigned slip rates. Our forecast relies on the assumption that the locations of past seismicity is a good guide to future seismicity, and that future large-magnitude events occur more likely in the vicinity of known faults. We show that the optimal weighted sum of the corresponding two spatial densities depends on the magnitude range considered. The kernel bandwidths and density weighting function are optimized using retrospective likelihood-based forecast experiments. We computed earthquake activity rates ( a - and b -value) of the truncated Gutenberg–Richter distribution separately for crustal and subduction seismicity based on a maximum likelihood approach that considers the spatial and temporal completeness history of the catalogue. The final annual rate of our forecast is purely driven by the maximum likelihood fit of activity rates to the catalogue data, whereas its spatial component incorporates contributions from both earthquake and fault moment-rate densities. Our model constitutes one branch of the earthquake source model logic tree of the 2013 European seismic hazard model released by the EU-FP7 project ‘Seismic HAzard haRmonization in Europe’ (SHARE) and contributes to the assessment of epistemic uncertainties in earthquake activity rates. We performed retrospective and pseudo-prospective likelihood consistency tests to underline the reliability of our model and SHARE's area source model (ASM) using the testing algorithms applied in the collaboratory for the study of earthquake predictability (CSEP). We comparatively tested our model's forecasting skill against the ASM and find a statistically significant better performance for testing periods of 10–20 yr. The testing results suggest that our model is a viable candidate model to serve for long-term forecasting on timescales of years to decades for the European region.

Description: To estimate the slowness vector of infrasound waves propagating across an infrasound array, it is often considered that the sensors are located in a perfectly horizontal plane (2-D approximation). However, the arrays are not planar and differences between sensor altitudes cannot be neglected without introducing biases when estimating the parameters of interest (backazimuth and incidence). A closed-formula of this error depending on the geometry of the array and the wave incidence is presented. Since the unbiased 3-D estimation of the slowness parameters results in a significant alteration of the variance, a metric based on the mean square error is proposed. We found it useful to compare the 3-D to the 2-D processed results. We show the 3-D estimator does not always give the best estimates. We propose also a formulation of the boundary that allow us to choose between these two estimators depending on the situation (arrays, incidence and signal-to-noise ratios). Finally, we compare these two approaches (2-D/3-D) to all International Monitoring System arrays with synthetic data, and perform a comparative estimation of backazimuth for real data.

Description: Microseismic monitoring is an essential tool for the characterization of hydraulic fractures. Fast estimation of the parameters that define a microseismic event is relevant to understand and control fracture development. The amount of data contained in the microseismic records however, poses a challenge for fast continuous detection and evaluation of the microseismic source parameters. Work inspired by the emerging field of Compressive Sensing has showed that it is possible to evaluate source parameters in a compressed domain, thereby reducing processing time. This technique performs well in scenarios where the amplitudes of the signal are above the noise level, as is often the case in microseismic monitoring using downhole tools. This paper extends the idea of the compressed domain processing to scenarios of microseismic monitoring using surface arrays, where the signal amplitudes are commonly at the same level as, or below, the noise amplitudes. To achieve this, we resort to the use of an imaging operator, which has previously been found to produce better results in detection and location of microseismic events from surface arrays. The operator in our method is formed by full-waveform elastodynamic Green's functions that are band-limited by a source time function and represented in the frequency domain. Where full-waveform Green's functions are not available, ray tracing can also be used to compute the required Green's functions. Additionally, we introduce the concept of the compressed inverse, which derives directly from the compression of the migration operator using a random matrix. The described methodology reduces processing time at a cost of introducing distortions into the results. However, the amount of distortion can be managed by controlling the level of compression applied to the operator. Numerical experiments using synthetic and real data demonstrate the reductions in processing time that can be achieved and exemplify the process of selecting the compression rate that produces a tolerable amount of distortion into the results.

Description: The topography of the core–mantle boundary (CMB) is directly linked to the dynamics of both the mantle and the outer core, although it is poorly constrained and understood. Recent studies have produced topography models with mutual agreement up to degree 2. A broad-band waveform inversion strategy is introduced and applied here, with relatively low computational cost and based on a first-order Born approximation. Its performance is validated using synthetic waveforms calculated in theoretical earth models that include different topography patterns with varying lateral wavelengths, from 600 to 2500 km, and magnitudes (~10 km peak-to-peak). The source–receiver geometry focuses mainly on the P diff , PKP , PcP and ScS phases. The results show that PKP branches, PcP and ScS generally perform well and in a similar fashion, while P diff yields unsatisfactory results. We investigate also how 3-D mantle correction influences the output models, and find that despite the disturbance introduced, the models recovered do not appear to be biased, provided that the 3-D model is correct. Using cross-correlated traveltimes, we derive new topography models from both P and S waves. The static corrections used to remove the mantle effect are likely to affect the inversion, compromising the agreement between models derived from P and S data. By modelling traveltime residuals starting from sensitivity kernels, we show how the simultaneous use of volumetric and boundary kernels can reduce the bias coming from mantle structures. The joint inversion approach should be the only reliable method to invert for CMB topography using absolute cross-correlation traveltimes.

Description: We introduce the single-station cross-correlation (SC) technique of processing ambient seismic noise and compare its results with the established cross-correlation (CC) and autocorrelation (AC) techniques. While CC is the correlation of the signals of two seismic stations with each other and AC is the correlation of a signal with itself, SC is the correlation of two different components of a single three-component seismic sensor. The comparison of the three different correlation techniques shows that CCs give the best results at frequencies below 0.5 Hz and that SCs give the best results at higher frequencies. In all three processing techniques, ambient seismic noise is correlated in order to reconstruct the Green's function describing the wave propagation between the first and the second sensor. By relating the coda parts of the daily Green's functions with the long-term reference Green's functions, shear wave velocity changes are determined. Here, we apply this technique to the data of 20 seismic stations in the surroundings of the fault zone of the Iwate-Miyagi Nairiku earthquake ( M W  = 6.9), which occurred on 2008 June 13, UTC (2008 June 14, Japan Standard Time) in the northern part of the Japanese island Honshu. The data range from 2008 January to 2011 June and therefore include the Tohoku earthquake ( M W  = 9.0), which occurred on 2011 March 11, off the coast of northern Honshu. The data are analysed in five different frequency ranges between 0.125 and 4.0 Hz. The data show coseismic velocity changes for both earthquakes followed by a post-seismic velocity recovery. In general, the coseismic velocity changes increase with frequency. For the Iwate-Miyagi Nairiku earthquake, the strongest velocity changes occur close to the fault zone. Quickly recovering coseismic velocity changes can be separated from changes not recovering during the study period. For the Tohoku earthquake, the complete area is affected by coseismic velocity changes. A modelling of the depth of the coseismic velocity changes indicates that the Iwate-Miyagi Nairiku earthquake can be explained either by large shallow velocity changes or by small, but deep changes. For one station, the observations can only be explained by assuming deeper changes. For the Tohoku earthquake, the modelling shows that different parts of the study area are affected in different ways, some showing shallow changes, others deeper changes. Furthermore, seasonal velocity variations occur, which are compatible for the different stations above 0.5 Hz, with velocity maxima in autumn.

Description: We report in this paper an original analysis of microseismic events (MSEs) induced by an excavation operation in the clay environment of the Mont Terri underground rock laboratory. In order to identify the MSEs with confidence, we develop a restrictive but efficient multistep method for filtering the recorded events. We deduce the spatial distribution and processes associated with the excavation-induced damage from the spatial location and focal mechanisms of the MSEs. We observe an asymmetric geometry of the excavation damaged zone around the excavated gallery, without notable microseismic activity in the sandy facies sidewall, in contrast with the shaly facies sidewall where a first burst of events is recorded, followed by two smaller bursts: one locates ahead of the excavation front and is associated with a dominant double-couple component, suggesting bedding plane reworking, that is, shear fracture mode, and the MSEs of the other cluster inside the shaly sidewall of the gallery, with a dominant compensated linear vector dipole component, suggesting extensive cracking. We identify and discuss four major factors that seem to control the MSEs source mechanisms: lithology, geometry of the geological features, gallery orientation and direction of the main compressive stress.

Description: High-rate GNSS has attracted increasing attention and numerous applications in geohazard monitoring and early warning. In this paper, we investigate three current existing single-receiver approaches for real-time GNSS seismology, comparing their observation models for equivalence and assessing the impact of main error components. We propose some refinements to the variometric approach and especially consider compensating the geometry error component by using the accurate initial coordinates before the earthquake to eliminate the drift trend in the integrated coseismic displacements. After careful corrections of satellite ephemeris, ionospheric delay, tropospheric delay and geometry errors, the refined variometric approach and the temporal point positioning (TPP) method have equivalent mathematical model with the converged precise point positioning (PPP). We evaluated the precision of the variometric and TPP approaches with various error correction schemes and duration time using numerous data sets and demonstrated that few centimetres accuracy of coseismic displacements is achievable even for 20 min interval. We applied these single-receiver approaches to process 1 Hz GPS data collected from the Tohoku-Oki earthquake ( M w 9.0, 2011 March 11) in Japan to capture coseismic displacement, and further, inverted the obtained displacement fields for fault slip distribution and moment magnitude. Comparisons of the results obtained using the refined variometric approach and TPP, as well as the converged PPP, displayed very good consistence both in coseismic displacements within few centimetres and in the slip distribution patterns and moment magnitudes.

Description: Macroseismic intensities are the only available data for most historical earthquakes and often represent the unique source of information for crucial events in the definition of seismic hazard. In this paper, we attempt at getting insight into source characteristics by reproducing the observed intensity field. As a test case, we study the source of 1908 Messina Straits earthquake ( M W  = 7.1), by testing three distinct fault models deduced from the analysis of geodetic data. Starting from the static slip distribution, we develop kinematic source models for the investigated fault and compute full waveform synthetic seismograms in a 1-D structural model, also accounting for anelastic attenuation. Then, we convert both computed peak-ground acceleration (PGA) and peak-ground velocity (PGV) to macroseismic intensity at 100 selected sites, by means of specific empirical relations for the Italian region. By comparing the original data separately with PGA- and PGV-based intensity fields, we discriminate among the tested faults and determine the best values for the investigated kinematic parameters of the source. We also perform a misfit analysis for the best source model, in order to investigate the dependence of the results on the selected parametrization. The results of the analysis indicate that among the tested models, the one characterized by an east-dipping fault, with strike-oriented NS slightly rotated clockwise, better explains the observed macroseismic field of the 1908 Messina Straits earthquake. Besides, the fracture nucleated at the southern end of the fault and ruptured northward, producing considerable directivity effects. This is in agreement with the published results obtained from the investigation of the historical seismograms. We also determine realistic values for the rupture velocity and the rise-time. Our study confirms the great potential of the macroseismic data, demonstrating that they contain enough information to constrain important characteristics of the fault, which can be retrieved by using complex source models and computing complete wavefield. Moreover, we also show that the simultaneous comparison of both PGA- and PGV-based synthetic macroseismic fields with the original intensities provides tighter constraints for discriminating among different source models, with respect to what attainable from each of them.

Description: We propose to test if gravimetry can prove useful in discriminating different models of long-term deep crustal processes in the case of the Taiwan mountain belt. We discuss two existing tectonic models that differ in the deep processes proposed to sustain the long-term growth of the orogen. One model assumes underplating of the uppermost Eurasian crust with subduction of the deeper part of the crust into the mantle. The other one suggests the accretion of the whole Eurasian crust above crustal-scale ramps, the lower crust being accreted into the collisional orogen. We compute the temporal gravity changes caused only by long-term rock mass transfers at depth for each of them. We show that the underplating model implies a rate of gravity change of –6 x 10 –2 μGal yr –1 , a value that increases to 2 x 10 –2 μGal yr –1 if crustal subduction is neglected. If the accretion of the whole Eurasian crust occurs, a rate of 7 x 10 –2 μGal yr –1 is obtained. The two models tested differ both in signal amplitude and spatial distribution. The yearly gravity changes expected by long-term deep crustal mass processes in Taiwan are two orders of magnitude below the present-day uncertainty of land-based gravity measurements. Assuming that these annually averaged long-term gravity changes will linearly accumulate with ongoing mountain building, multidecadal time-series are needed to identify comparable rates of gravity change. However, as gravity is sensitive to any mass redistribution, effects of short-term processes such as seismicity and surface mass transfers (erosion, sedimentation, ground-water) may prevent from detecting any long-term deep signal. This study indicates that temporal gravity is not appropriate for deciphering the long-term deep crustal processes involved in the Taiwan mountain belt.

Description: Perfectly matched layer (PML) is an efficient absorbing technique for numerical wave simulations. Since it appeared, various improvements have been made. The complex frequency-shifted PML (CFS-PML) improves the absorbing performance for near-grazing incident waves and evanescent waves. The auxiliary differential equation (ADE) formulation of the PML provides a convenient unsplit-field PML implementation that can be directly used with high order time marching schemes. The multi-axial PML (MPML) stabilizes the PML on anisotropic media. However, these improvements were generally developed for Cartesian grids. In this paper, we extend the ADE CFS-PML to general curvilinear (non-orthogonal) grids for elastic wave modelling. Unlike the common implementations to absorb the waves in the computational space, we apply the damping along the perpendicular direction of the PML layer in the local Cartesian coordinates. Further, we relate the perpendicular and parallel components of the gradient operator in the local Cartesian coordinates to the derivatives in the curvilinear coordinates, to avoid mapping the wavefield to the local Cartesian coordinates. It is thus easy to be incorporated with numerical schemes on curvilinear grids. We derive the PML equations for the interior region and for the free surface separately because the free surface boundary condition modifies the elastic wave equations. We show that the elastic wave modelling on curvilinear grids exhibits anisotropic effects in the computational space, which may lead to unstable simulations. To stabilize the simulation, we adapt the MPML strategy to also absorb the wavefield along the two parallel directions of the PML. We illustrate the stability of this ADE CFS-MPML for finite-difference elastic wave simulations on curvilinear grids by two numerical experiments.

Description: Seismic monitoring can greatly benefit from imaging events with a low signal-to-noise ratio (SNR) as the number of the events with a low signal grows exponentially. One way to detect weaker events is improvement of a SNR by migration-type stacking of waveforms from multiple stations. We have developed a new method of location of seismic events that involves stacking of seismic phases and amplitudes along diffraction traveltime curves to suppress noise and detect seismic events with a SNR lower than that on individual receivers. The stacking includes polarity correction based on a simultaneous seismic moment tensor inversion and detection algorithm on the stack function. We applied this method to locate microseismicity induced by hydraulic fracturing. First we calibrated the velocity model by locating perforation shots at known locations. Then we processed 3 d of data from microseismic monitoring of shale stimulation and benchmarked migration-type locations of the largest events that were manually located. The detected and located events induced by hydraulic fracturing in this case study are mostly shear events forming narrow bands along the maximum horizontal stress direction approximately 100 m above the injection intervals. The proposed technique is fully automated and feasible for real-time seismic monitoring.

Description: A series of six seismic events that occurred in one panel of Rudna copper-ore mine in Poland is studied. Although the events had comparable magnitudes, from 3.0 to 3.7, their ground effects were very diverse. Comparing the accelerations observed at various locations with their estimates from ground motion prediction equation the events split into three distinct pairs. The first pair ground effects exceed considerably the estimates at most of observation points, the second pair effects are abnormally high at short epicentral distances, whereas the third pair effects are much less than the estimates at most of observation points. The similarities in ground effects correlate with the fault mechanisms similarities. The first two pairs’ events, whose ground effects were strong, exhibit thrust faulting and the third pair events, which caused unexpectedly low ground motion exhibit normal faulting mechanisms. The paired events have also similar apparent stress values. These stress values of the two events of very weak ground effects are distinctly lower than the values of other four events. All events demonstrate dominating non-double-couple components in the overall mechanisms. A kinematic source analysis indicates that these events have extended rather than point sources, and five of them have distinct directivity effects. A static stress transfer analysis signifies interrelations between these events. The rupture of all events started in areas of Coulomb failure function increase due to the cumulative effect of previous events. Linking results of the ground effects, source, rupture and interaction analyses tentative geodynamic conclusions are formulated. The untypical and diverse ground effects of the studied events result likely from the events’ complexity expressed by tensile source mechanisms, finite sources, directivity of ruptures and nearly horizontal rupture planes. The above features seem to be implied by a superposition of coseismic alterations of stress field and stress changes due to mining. One cannot, therefore, exclude the possibility of other cases of significant deviations from the expected ground motion amplitudes, due to specific geodynamics in another parts of the mine. An analysis like that done within this work can allow, however, foreseeing such extreme surface impacts.

Description: We present an epidemic type aftershock sequenc (ETAS) model where the offspring rates vary both spatially and temporally. This is achieved by distinguishing between those space–time volumes where the interpoint space and time distances are small, and those where they are considerably larger. We also question the nature of the background component in the ETAS model. Is it simply a temporal boundary correction ( t  = 0) or does it represent an additional tectonic process not described by the aftershock component? The form of these stochastic models should not be considered to be fixed. As we accumulate larger and better earthquake catalogues, GPS data, strain rates, etc., we have the ability to ask more complex questions about the nature of the process. By fitting modified models consistent with such questions, we should gain a better insight into the earthquake process. Hence, we consider a sequence of incrementally modified ETAS type models rather than ‘the’ ETAS model.

Description: The computation of quasi-static deformation for axisymmetric viscoelastic structures on a gravitating spherical earth is addressed using the spectral element method (SEM). A 2-D spectral element domain is defined with respect to spherical coordinates of radius and angular distance from a pole of symmetry, and 3-D viscoelastic structure is assumed to be azimuthally symmetric with respect to this pole. A point dislocation source that is periodic in azimuth is implemented with a truncated sequence of azimuthal order numbers. Viscoelasticity is limited to linear rheologies and is implemented with the correspondence principle in the Laplace transform domain. This leads to a series of decoupled 2-D problems which are solved with the SEM. Inverse Laplace transform of the independent 2-D solutions leads to the time-domain solution of the 3-D equations of quasi-static equilibrium imposed on a 2-D structure. The numerical procedure is verified through comparison with analytic solutions for finite faults embedded in a laterally homogeneous viscoelastic structure. This methodology is applicable to situations where the predominant structure varies in one horizontal direction, such as a structural contrast across (or parallel to) a long strike-slip fault.

Description: The calculation of first-order P -wave ray synthetic seismograms based on first-order ray tracing (FORT) and dynamic ray tracing (FODRT) for P -waves propagating in inhomogeneous, weakly anisotropic media is extended from smooth to layered media. All the basic formulae necessary to calculate the P -wave FORT and FODRT quantities inside layers and to transform them at the points of reflection/transmission are given. The proposed formulae are applicable in subcritical as well as overcritical regions. The accuracy of the results is tested by comparing the approximate (FORT) results with the results obtained from a standard ray tracer for anisotropic media. The tests indicate that, except for critical regions, where the ray theory provides incorrect results anyway, the accuracy of FORT and FODRT in layered media is comparable with the accuracy in smooth media.