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

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

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

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

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

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

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

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

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

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

Description: 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: Relative to the gravitational potential energy of the Earth's monopole, the multipole energy has received far less attention. In this paper, we recapitulate the basic physics from first principles and derive the formulas for multipole energies in analogy to classical electrostatic theory. We focus on the zonal quadrupole energy associated with the Earth's oblateness, the dominant term in Earth's gravity field apart from the monopole. We find the gravitational energy E oblateness 10 –6 | E monopole | = +2.5 x 10 26 J. We examine the implications of E oblateness and its changes associated with long-term ‘secular’ decreases in the oblateness parameter J 2 . We find the rate of loss of E oblateness due to the Earth rounding induced by the present-day GIA is about –200 GW, an amount quite significant in the kinetic energy budget of the mantle heat engine that drives the plate tectonics that has been estimated to be ~1 TW. We also assert that the tidal braking and the global earthquake dislocations, both resulting in Earth rounding on long-term geological timescales, are accompanied with a secular decrease of E oblateness at nearly the same rate of several GW.

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: Large-scale chemical lateral heterogeneities are inferred in the Earth's lowermost mantle by seismological studies. We explore the model space of thermochemical convection that can maintain reservoirs of dense material for a long period of time, by using similar analysis in 3-D spherical geometry. In this study, we focus on the parameters thought to be important in controlling the stability and structure of primordial dense reservoirs in the lower mantle, including the chemical density contrast between the primordial dense material and the regular mantle material (buoyancy ratio), thermal and chemical viscosity contrasts, volume fraction of primordial dense material and the Clapeyron slope of the phase transition at 660 km depth. We find that most of the findings from the 3-D Cartesian study still apply to 3-D spherical cases after slight modifications. Varying buoyancy ratio leads to different flow patterns, from rapid upwelling to stable layering; and large thermal viscosity contrasts are required to generate long wavelength chemical structures in the lower mantle. Chemical viscosity contrasts in a reasonable range have a second-order role in modifying the stability of the dense anomalies. The volume fraction of the initial primordial dense material does not effect the results with large thermal viscosity contrasts, but has significant effects on calculations with intermediate and small thermal viscosity contrasts. The volume fraction of dense material at which the flow pattern changes from unstable to stable depends on buoyancy ratio and thermal viscosity contrast. An endothermic phase transition at 660 km depth acts as a ‘filter’ allowing cold slabs to penetrate while blocking most of the dense material from penetrating to the upper mantle.

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: Relative sea level curves contain coupled information about absolute sea level change and vertical lithospheric movement. Such curves may be constructed based on, for example tide gauge data for the most recent times and different types of geological data for ancient times. Correct account for vertical lithospheric movement is essential for estimation of reliable values of absolute sea level change from relative sea level data and vise versa. For modern times, estimates of vertical lithospheric movement may be constrained by data (e.g. GPS-based measurements), which are independent from the relative sea level data. Similar independent data do not exist for ancient times. The purpose of this study is to test two simple inversion approaches for simultaneous estimation of lithospheric uplift rates and absolute sea level change rates for ancient times in areas where a dense coverage of relative sea level data exists and well-constrained average lithospheric movement values are known from, for example glacial isostatic adjustment (GIA) models. The inversion approaches are tested and used for simultaneous estimation of lithospheric uplift rates and absolute sea level change rates in southwest Scandinavia from modern relative sea level data series that cover the period from 1900 to 2000. In both approaches, a priori information is required to solve the inverse problem. A priori information about the average vertical lithospheric movement in the area of interest is critical for the quality of the obtained results. The two tested inversion schemes result in estimated absolute sea level rise of ~1.2/1.3 mm yr –1 and vertical uplift rates ranging from approximately –1.4/–1.2 mm yr –1 (subsidence) to about 5.0/5.2 mm yr –1 if an a priori value of 1 mm yr –1 is used for the vertical lithospheric movement throughout the study area. In case the studied time interval is broken into two time intervals (before and after 1970), absolute sea level rise values of ~0.8/1.2 mm yr –1 (before 1970) and ~2.0 mm yr –1 (after 1970) are found. The uplift patterns resulting from the different inversions suggest that the lithospheric post-GIA response changes near the border between the Danish Basin and the Fennoscandian Shield. The obtained patterns of vertical lithospheric movement rates are comparable to results from other studies based on different and similar data types. Main differences between the inversion results and the results from other studies are caused by factors such as the simplifications included in the inversion approach, such as neglecting local sea level variation caused by the dominant wind patterns, and the a priori values chosen for the vertical uplift rates. The tests of the inversion schemes reveal that realistic values of absolute sea level rise and lithospheric uplift may be simultaneously estimated provided that reliable prior knowledge regarding the overall lithospheric uplift in the study area is available beforehand. In the presented parametrizations, only one absolute sea level change rate value is estimated for each studied time interval while several vertical movement rates are found, and the inverse estimate of absolute sea level change rate is practically insensitive with respect to the choice of a priori value of absolute sea level change, as long as the uncertainty assigned to this a priori value is kept sufficiently high.

Description: In 1356, a magnitude 6–7 earthquake occurred near Basel, in Switzerland. But recent compilations of GPS measurements reveal that measured horizontal deformation rates in northwestern continental Europe are smaller than error bars on the measurements, proving present tectonic activity, if any, is very small in this area. We propose to reconcile these apparently antinomic observations with a mechanical model of the lithosphere that takes into account the geometry of the lithosphere–asthenosphere boundary, assuming that the only loading mechanism is gravity. The lithosphere is considered to be an elastoplastic material satisfying a Von Mises plasticity criterion. The model, which is 400 km long, 360 km wide and 230 km thick, is centred near Belfort in eastern France, with its width oriented parallel to the N145°E direction. It also takes into account the real topography of both the ground surface and that of the Moho discontinuity. Not only does the model reproduce observed principal stress directions orientations, it also identifies a plastic zone that fits roughly the most seismically active domain of the region. Interestingly, a somewhat similar stress map may be produced by considering an elastic lithosphere and an ad-hoc horizontal ‘tectonic’ stress field. However, for the latter model, examination of the plasticity criterion suggests that plastic deformation should have taken place. It is concluded that the present-day stress field in this region is likely controlled by gravity and rheology, rather than by active Alpine tectonics.

Description: We have investigated variations in transition zone thickness under the Borborema Province of NE Brazil by migrating and stacking teleseismic P -wave receiver functions at 32 seismic stations in the region. The Borborema Province represents the western portion of a larger Neoproterozoic mobile belt that occupied much of northern Gondwana, where extensional processes in the Mesozoic lead to the formation of a number of intracontinental basins and ultimately continental breakup. Episodes of intraplate volcanism and uplift marked the evolution of the Province during the Cenozoic, but it is unclear whether those episodes originated from shallow or deep-seated magmatic sources. On one hand, the elliptical shape of the uplifted area, the stress pattern of the Cenozoic deformation and the time overlap between uplift and volcanism suggest doming from thermal activation due to a deep-seated mantle plume. On the other hand, geochronological dates of volcanic bodies in the Province are better understood if resulting from lithospheric erosion by a shallow, small-scale convection cell. Large temperature anomalies are expected to be associated with mantle upwellings, and constraints on the depth extent of the upwellings can be obtained from transition zone thickness. Thinning of the transition zone with respect to its nominal 250 km value is considered diagnostic for positive temperature anomalies, while thickening is considered diagnostic for negative anomalies. Our results show that transition zone thickness is normal, around 250 km, throughout the Province and suggest that thermal perturbations—if present—are confined to the upper mantle. We argue that our results are consistent with a local, shallow magmatic source for the Cenozoic intraplate volcanism of the Borborema Province, although other proposed scenarios—such as channeling of upwelling plume material along lithospheric thin spots—cannot be ruled out with our analysis.

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: Standard techniques for computed tomography imaging are not directly applicable to a carbonate rock because of the geometric complexity of its pore space. In this study, we first characterized the pore structure in Majella limestone with 30 per cent porosity. Microtomography data acquired on this rock was partitioned into three distinct domains: macropores, solid grains, and an intermediate domain made up of voxels of solid embedded with micropores below the resolution. A morphological analysis of the microtomography images shows that in Majella limestone both the solid and intermediate domains are interconnected in a manner similar to that reported previously in a less porous limestone. We however show that the macroporosity in Majella limestone is fundamentally different, in that it has a percolative backbone which may contribute significantly to its permeability. We then applied for the first time 3-D-volumetric digital image correlation (DIC) to characterize the mode of mechanical failure in this limestone. Samples were triaxially deformed over a wide range of confining pressures. Tomography imaging was performed on these samples before and after deformation. Inelastic compaction was observed at all tested pressures associated with both brittle and ductile behaviors. Our DIC analysis reveals the structure of compacting shear bands in Majella limestone deformed in the transitional regime. It also indicates an increase of geometric complexity with increasing confinement—from a planar shear band, to a curvilinear band, and ultimately to a diffuse multiplicity of bands, before shear localization is inhibited as the failure mode completes the transition to delocalized cataclastic flow.

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: Wide-angle reflection/refraction seismic profiles were recorded across the Cyprus Arc, the plate boundary between the African Plate and the Aegean–Anatolian microplate, from the Eratosthenes Seamount to the Hecataeus Rise immediately south of Cyprus. The resultant models were able to resolve detail of significant lateral velocity variations, though the deepest crust and Moho are not well resolved from the seismic data alone. Conclusions from the modelling suggest that (i) Eratosthenes Seamount consists of continental crust but exhibits a laterally variable velocity structure with a thicker middle crust and thinner lower crust to the northeast; (ii) the Hecataeus Rise has a thick sedimentary rock cover on an indeterminate crust (likely continental) and the crust is significantly thinner than Eratosthenes Seamount based on gravity modelling; (iii) high velocity basement blocks, coincident with highs in the magnetic field, occur in the deep water between Eratosthenes and Hecataeus, and are separated and bounded by deep low-velocity troughs and (iv) one of the high velocity blocks runs parallel to the Cyprus Arc, while the other two appear linked based on the magnetic data and run NW–SE, parallel to the margin of the Hecataeus Rise. The high velocity block beneath the edge of Eratosthenes Seamount is interpreted as an older magmatic intrusion while the linked high velocity blocks along Hecataeus Rise are interpreted as deformed remnant Tethyan oceanic crust or mafic intrusives from the NNW–SSE oriented transform margin marking the northern boundary of Eratosthenes Seamount. Eratosthenes Seamount, the northwestern limit of rifted continental crust from the Levant Margin, is part of a jagged rifted margin transected by transform faults on the northern edge of the lower African Plate that is being obliquely subducted under the Aegean–Anatolian upper plate. The thicker crust of Eratosthenes Seamount may be acting as an asperity on the subducting slab, locally locking up subduction of the Cyprus Arc on its northern margin, while deformed Tethyan oceanic crust remains trapped between its northeastern margin and the Hecataeus Rise.

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 present an up-to-date high resolution picture of the ongoing crustal deformation field of Italy, based on an extensive combination of permanent and non-permanent GPS observations carried out since 1994. In addition, we present an updated map of contemporary S Hmax orientations computed by a multidisciplinary data set of well-constrained stress indicators, including both published results and novel analyses. The comparison of stress and geodetic strain-rates directions reveals that both patterns are near-parallel over a large part of the investigated area, highlighting that crustal stress and surface deformation are driven by the same mechanism. The comparison of the azimuthal patterns of surface strain and mantle deformation shows a modest correlation on the Alps and a low correlation along the Apennines chain and the Calabro-Peloritan Arc. Along the Apennines chain, this feature suggests the occurrence of significant strain partitioning and crust–mantle mechanical decoupling. Along the Calabro-Peloritan Arc, the apparent low correlation reflects a different mantle–crust mechanism of deformation to the ongoing subduction and rollback of the Ionian slab. In addition, the superposition of regional/local effects related to second-order sources (crustal lateral density changes, strength contrasts), which at regional/local scale modulate the crustal stress/strain-rate pattern, cannot be ruled out.

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: A series of linear analysis was performed on the onset of thermal convection of highly compressible fluids, in order to deepen the fundamental insights into the mantle convection of massive super-Earths in the presence of strong adiabatic compression. We consider the temporal evolution (growth or decay) of an infinitesimal perturbation superimposed to a highly compressible fluid which is in a hydrostatic (motionless) and conductive state in a basally heated horizontal layer. As a model of pressure-dependence in material properties, we employed an exponential decrease in thermal expansivity α and exponential increase in (reference) density with depth. The linearized equations for conservation of mass, momentum and internal (thermal) energy are numerically solved for the critical Rayleigh number as well as the vertical profiles of eigenfunctions for infinitesimal perturbations. The above calculations are repeatedly carried out by systematically varying (i) the dissipation number (Di), (ii) the temperature at the top surface and (iii) the magnitude of pressure-dependence in α and . Our analysis demonstrated that the onset of thermal convection is strongly affected by the adiabatic compression, in response to the changes in the static stability of thermal stratification in the fluid layer. For sufficiently large Di where a thick sublayer of stable stratification develops in the layer, for example, the critical Rayleigh number explosively increases with Di, together with drastic decreases in the length scales of perturbations both in vertical and horizontal directions. In particular, for very large Di, a thick ‘stratosphere’ occurs in the fluid layer where the vertical motion is significantly suppressed, resulting in a shrink of the incipient convection in a thin sublayer of unstable thermal stratification. In addition, when Di exceeds a threshold value above which a thermal stratification becomes stable in the entire layer, no perturbation is allowed to grow with time regardless of the Rayleigh number and/or the horizontal wavelength. We also found that the effect of adiabatic compression becomes prominent for higher temperature at the top surface of the fluid layer. These findings may imply the crucial importance of adiabatic compression in understanding the dynamics and evolution of the mantles of massive super-Earths, particularly for those orbiting their parent stars very closely.

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: In this paper, we propose an approach to compute the coseismic Earth's volume change based on a spherical-Earth elastic dislocation theory. We present a general expression of the Earth's volume change for three typical dislocations: the shear, tensile and explosion sources. We conduct a case study for the 2004 Sumatra earthquake ( M w 9.3), the 2010 Chile earthquake ( M w 8.8), the 2011 Tohoku-Oki earthquake ( M w 9.0) and the 2013 Okhotsk Sea earthquake ( M w 8.3). The results show that mega-thrust earthquakes make the Earth expand and earthquakes along a normal fault make the Earth contract. We compare the volume changes computed for finite fault models and a point source of the 2011 Tohoku-Oki earthquake ( M w 9.0). The big difference of the results indicates that the coseismic changes in the Earth's volume (or the mean radius) are strongly dependent on the earthquakes’ focal mechanism, especially the depth and the dip angle. Then we estimate the cumulative volume changes by historical earthquakes ( M w ≥ 7.0) since 1960, and obtain an Earth mean radius expanding rate about 0.011 mm yr –1 .

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: We present the crustal resistivity structure of the Pamir and Southern Tian Shan orogenic belts at the northwestern promontory of the India–Asia collision zone. The magnetotelluric (MT) data were recorded along a roughly north–south trending, 350 km long corridor from the Pamir Plateau in southern Tajikistan across the Pamir frontal ranges, the Alai Valley and the southwestern Tian Shan to Osh in the Kyrgyz part of the Fergana Basin. In total, we measured at 178 sites, whereof 26 combine broad band and long period recordings. One of the most intriguing features of the 2-D and 3-D inversion results is a laterally extended zone of high electrical conductivity below the Pamir Plateau, with resistivities below 1 m, starting at a depth of ~10–15 km. The high conductivity can be explained with the presence of partially molten rocks at middle to lower crustal levels, possibly related to ongoing migmatization and/or middle/lower crustal flow underneath the Southern Pamir. This interpretation is consistent with a low velocity zone found from local earthquake tomography, relatively high v p / v s ratios, elevated surface heat flow, and thermomechanical modelling suggesting that melting temperatures are reached in the felsic middle crust. In the upper crust of the Pamir and Tian Shan, the Palaeozoic–Mesozoic suture zones appear as electrically conductive, whereas the compact metamorphic rocks of the Muskol-Shatput Dome of the Central Pamir are highly resistive. The intra-montane basin of the Alai Valley—sandwiched between the Pamir and Tian Shan—exhibits a generally conductive upper crust that bifurcates into two conductors at depth. One of them connects to the active Main Pamir Thrust, which is absorbing most of today's convergence between the Pamir and the Tian Shan. Several deeper zones of high conductivity in the middle and lower crust of Central and Northern Pamir likely record fluid release due to metamorphism associated with active continental subduction/delamination.

Description: A new model of the deglaciation history of Antarctica over the past 25 kyr has been developed, which we refer to herein as ICE-6G_C (VM5a). This revision of its predecessor ICE-5G (VM2) has been constrained to fit all available geological and geodetic observations, consisting of: (1) the present day uplift rates at 42 sites estimated from GPS measurements, (2) ice thickness change at 62 locations estimated from exposure-age dating, (3) Holocene relative sea level histories from 12 locations estimated on the basis of radiocarbon dating and (4) age of the onset of marine sedimentation at nine locations along the Antarctic shelf also estimated on the basis of 14 C dating. Our new model fits the totality of these data well. An additional nine GPS-determined site velocities are also estimated for locations known to be influenced by modern ice loss from the Pine Island Bay and Northern Antarctic Peninsula regions. At the 42 locations not influenced by modern ice loss, the quality of the fit of postglacial rebound model ICE-6G_C (VM5A) is characterized by a weighted root mean square residual of 0.9 mm yr –1 . The Southern Antarctic Peninsula is inferred to be rising at 2 mm yr –1 , requiring there to be less Holocene ice loss there than in the prior model ICE-5G (VM2). The East Antarctica coast is rising at approximately 1 mm yr –1 , requiring ice loss from this region to have been small since Last Glacial Maximum. The Ellsworth Mountains, at the base of the Antarctic Peninsula, are inferred to be rising at 5–8 mm yr –1 , indicating large ice loss from this area during deglaciation that is poorly sampled by geological data. Horizontal deformation of the Antarctic Plate is minor with two exceptions. First, O'Higgins, at the tip of the Antarctic Peninsula, is moving southeast at a significant 2 mm yr –1 relative to the Antarctic Plate. Secondly, the margins of the Ronne and Ross Ice Shelves are moving horizontally away from the shelf centres at an approximate rate of 0.8 mm yr –1 , in viscous response to the early Holocene unloading of ice from the current locations of the ice shelf centers. ICE-6G_C (VM5A) fits the horizontal observations well (wrms residual speed of 0.7 mm yr –1 ), there being no need to invoke any influence of lateral variation in mantle viscosity. ICE-6G_C (VM5A) differs in several respects from the recently published W12A model of Whitehouse et al. First, the upper-mantle viscosity in VM5a is 5 10 20 Pa s, half that in W12A. The VM5a profile, which is identical to that inferred on the basis of the Fennoscandian relaxation spectrum, North American relative sea level histories and Earth rotation constraints, when coupled with the revised ICE-6G_C deglaciation history, fits all of the available constraints. Secondly, the net contribution of Antarctica ice loss to global sea level rise is 13.6 m, 2/3 greater than the 8 m in W12A. Thirdly, ice loss occurs quickly from 12 to 5 ka, and the contribution to global sea level rise during Meltwater Pulse 1B (11.5 ka) is large (5 m), consistent with sedimentation constraints from cores from the Antarctica ice shelf. Fourthly, in ICE-6G_C there is no ice gain in the East Antarctica interior, as there is in W12A. Finally, the new model of Antarctic deglaciation reconciles the global constraint upon the global mass loss during deglaciation provided by the Barbados record of relative sea level history when coupled with the Northern Hemisphere counterpart of this new model.

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.