Showing papers in "Geophysical Journal International in 2010"

Geologically current plate motions

Abstract: SUMMARY We describe best-fitting angular velocities and MORVEL, a new closure-enforced set of angular velocities for the geologically current motions of 25 tectonic plates that collectively occupy 97 per cent of Earth's surface. Seafloor spreading rates and fault azimuths are used to determine the motions of 19 plates bordered by mid-ocean ridges, including all the major plates. Six smaller plates with little or no connection to the mid-ocean ridges are linked to MORVEL with GPS station velocities and azimuthal data. By design, almost no kinematic information is exchanged between the geologically determined and geodetically constrained subsets of the global circuit—MORVEL thus averages motion over geological intervals for all the major plates. Plate geometry changes relative to NUVEL-1A include the incorporation of Nubia, Lwandle and Somalia plates for the former Africa plate, Capricorn, Australia and Macquarie plates for the former Australia plate, and Sur and South America plates for the former South America plate. MORVEL also includes Amur, Philippine Sea, Sundaland and Yangtze plates, making it more useful than NUVEL-1A for studies of deformation in Asia and the western Pacific. Seafloor spreading rates are estimated over the past 0.78 Myr for intermediate and fast spreading centres and since 3.16 Ma for slow and ultraslow spreading centres. Rates are adjusted downward by 0.6–2.6 mm yr−1 to compensate for the several kilometre width of magnetic reversal zones. Nearly all the NUVEL-1A angular velocities differ significantly from the MORVEL angular velocities. The many new data, revised plate geometries, and correction for outward displacement thus significantly modify our knowledge of geologically current plate motions. MORVEL indicates significantly slower 0.78-Myr-average motion across the Nazca–Antarctic and Nazca–Pacific boundaries than does NUVEL-1A, consistent with a progressive slowdown in the eastward component of Nazca plate motion since 3.16 Ma. It also indicates that motions across the Caribbean–North America and Caribbean–South America plate boundaries are twice as fast as given by NUVEL-1A. Summed, least-squares differences between angular velocities estimated from GPS and those for MORVEL, NUVEL-1 and NUVEL-1A are, respectively, 260 per cent larger for NUVEL-1 and 50 per cent larger for NUVEL-1A than for MORVEL, suggesting that MORVEL more accurately describes historically current plate motions. Significant differences between geological and GPS estimates of Nazca plate motion and Arabia–Eurasia and India–Eurasia motion are reduced but not eliminated when using MORVEL instead of NUVEL-1A, possibly indicating that changes have occurred in those plate motions since 3.16 Ma. The MORVEL and GPS estimates of Pacific–North America plate motion in western North America differ by only 2.6 ± 1.7 mm yr−1, ≈25 per cent smaller than for NUVEL-1A. The remaining difference for this plate pair, assuming there are no unrecognized systematic errors and no measurable change in Pacific–North America motion over the past 1–3 Myr, indicates deformation of one or more plates in the global circuit. Tests for closure of six three-plate circuits indicate that two, Pacific–Cocos–Nazca and Sur–Nubia–Antarctic, fail closure, with respective linear velocities of non-closure of 14 ± 5 and 3 ± 1 mm yr−1 (95 per cent confidence limits) at their triple junctions. We conclude that the rigid plate approximation continues to be tremendously useful, but—absent any unrecognized systematic errors—the plates deform measurably, possibly by thermal contraction and wide plate boundaries with deformation rates near or beneath the level of noise in plate kinematic data.

International Geomagnetic Reference Field: the eleventh generation

Abstract: The eleventh generation of the International Geomagnetic Reference Field (IGRF) was adopted in December 2009 by the International Association of Geomagnetism and Aeronomy Working Group V-MOD. It updates the previous IGRF generation with a definitive main field model for epoch 2005.0, a main field model for epoch 2010.0, and a linear predictive secular variation model for 2010.0–2015.0. In this note the equations defining the IGRF model are provided along with the spherical harmonic coefficients for the eleventh generation. Maps of the magnetic declination, inclination and total intensity for epoch 2010.0 and their predicted rates of change for 2010.0–2015.0 are presented. The recent evolution of the South Atlantic Anomaly and magnetic pole positions are also examined.

Hydrolytic Weakening of Quartz and Other Silicates

Abstract: Summary Rocks at high pressure and temperature become weak when exposed to an aqueous environment. Four different weakening processes have been found: (1) Increase in pore pressure reduces the effective confining pressure, thus negating the strengthening effect of confining pressure (Terzaghi). (2) Penetration of water into the intergranular boundaries reduces the cohesive strength. This is believed to have occurred in the dehydration weakening found in gypsum (Heard & Rubey) and in serpentinite (Raleigh & Paterson). (3) Water promotes recrystallization, greatly reducing the strength at low strain-rates at temperatures above the critical temperature for recrystallization. (4) Water permeates most silicate crystals so far tested and causes hydrolytic weakening of the component crystals of rocks. The most revealing experiments on hydrolytic weakening have been those on single crystals of synthetic quartz with varying water content (Griggs & Blacic). The water content inferred from 3 micron infra-red absorption varies from 0.0015 to 0.13 wt% (100–9000 H/106 Si) in the nine crystals tested. All crystals became weak at a critical temperature dependent on the water content, varying from 380°C at 0.13% to 1070°C at 0.0015%. The weakening process is thermally activated, is reversible, and is rate dependent. The mechanism of deformation is intracrystalline glide. MLaren has shown by transmission electron microscopy that deformed weakened synthetic quartz contains dislocation densities similar to those found in strong dry quartz deformed in the ductile regime. Annealing of water-weakened quartz removes most of the dislocations and causes the formation of bubbles. These are believed to be water bubbles formed by the reaction Si(OH)4→ SiO2+2H2O, the water migrating through the lattice. The Frank-Griggs hypothesis of hydrolytic weakening is that easy slip occurs only when Si-O-Si bridges adjacent to a dislocation are hydrolyzed by the migration of water. No Si-O bonds need be broken, for dislocation motion can occur by hydrogen bond exchange. All the evidence in hand suggests that hydrolytic weakening is general in silicates containing very small amounts of water. This phenomenon is expected to be important in the mantle, leading to a rate-dependent strength.

Determination of permeability from spectral induced polarization in granular media

Abstract: SUMMARY The surface conductivity of porous rocks has two contributions: the first is associated with the diffuse layer coating the grains and is frequency-independent as long as the diffuse layer is above a percolation threshold. The second contribution is associated with the Stern layer of weakly sorbed counterions on the mineral surface and is frequency-dependent if the Stern layer is discontinuous at the scale of the representative elementary volume. In the frequency range 1 mHz–100 Hz, this second contribution is also associated with the main polarization mechanism observed by the spectral induced polarization method in granular media (neglecting the contribution of other polarization processes like those associated with redox processes and membrane polarization). At the macroscale, we connect the Stern layer contribution to the complex conductivity and to the expectation of the probability distribution of the inverse of the grain size. This is done by performing a convolution between the probability distribution of the inverse of the grain size and the surface conductivity response obtained when all the grains have the same size. Surface conductivity at the macroscopic scale is also connected to an effective pore size used to characterize permeability. From these relationships, a new equation is derived connecting this effective pore size, the electrical formation factor, and the expected value of the probability distribution for the inverse of the grain size, which is in turn related to the distribution of the relaxation times. These new relationships are consistent with various formula derived in the literature in the limit where the grain size distribution is given by the delta function or a log normal distribution and agree fairly well with various experimental data showing also some limitations of the induced polarization method to infer permeability. One of these limitations is the difficulty to detect the relaxation, in the phase, associated with the smaller grains, as this polarization may be hidden by the Maxwell–Wagner polarization at relatively high frequencies (>100 Hz). Also, cemented aggregates of grains can behave as coarser grains.

Seismic tomography of the southern California crust based on spectral‐element and adjoint methods

Abstract: We iteratively improve a 3-D tomographic model of the southern California crust using numerical simulations of seismic wave propagation based on a spectral-element method (SEM) in combination with an adjoint method. The initial 3-D model is provided by the Southern California Earthquake Center. The data set comprises three-component seismic waveforms (i.e. both body and surface waves), filtered over the period range 2–30 s, from 143 local earthquakes recorded by a network of 203 stations. Time windows for measurements are automatically selected by the FLEXWIN algorithm. The misfit function in the tomographic inversion is based on frequency-dependent multitaper traveltime differences. The gradient of the misfit function and related finite-frequency sensitivity kernels for each earthquake are computed using an adjoint technique. The kernels are combined using a source subspace projection method to compute a model update at each iteration of a gradient-based minimization algorithm. The inversion involved 16 iterations, which required 6800 wavefield simulations. The new crustal model, m_(16), is described in terms of independent shear (V_S) and bulk-sound (V_B) wave speed variations. It exhibits strong heterogeneity, including local changes of ±30 per cent with respect to the initial 3-D model. The model reveals several features that relate to geological observations, such as sedimentary basins, exhumed batholiths, and contrasting lithologies across faults. The quality of the new model is validated by quantifying waveform misfits of full-length seismograms from 91 earthquakes that were not used in the tomographic inversion. The new model provides more accurate synthetic seismograms that will benefit seismic hazard assessment.

The angular velocities of the plates and the velocity of Earth's centre from space geodesy

Abstract: SUMMARY Using space geodetic observations from four techniques (GPS, VLBI, SLR and DORIS), we simultaneously estimate the angular velocities of 11 major plates and the velocity of Earth's centre. We call this set of relative plate angular velocities GEODVEL (for GEODesy VELocity). Plate angular velocities depend on the estimate of the velocity of Earth's centre and on the assignment of sites to plates. Most geodetic estimates of the angular velocities of the plates are determined assuming that Earth's centre is fixed in an International Terrestrial Reference Frame (ITRF), and are therefore subject to errors in the estimate of the velocity of Earth's centre. In ITRF2005 and ITRF2000, Earth's centre is the centre of mass of Earth, oceans and atmosphere (CM); the velocity of CM is estimated by SLR observation of LAGEOS's orbit. Herein we define Earth's centre to be the centre of mass of solid Earth (CE); we determine the velocity of CE by assuming that the portions of plate interiors not near the late Pleistocene ice sheets move laterally as if they were part of a rigid spherical cap. The GEODVEL estimate of the velocity of CE is likely nearer the true velocity of CM than are the ITRF2005 and ITRF2000 estimates because (1) no phenomena can sustain a significant velocity between CM and CE, (2) the plates are indeed nearly rigid (aside from vertical motion) and (3) the velocity of CM differs between ITRF2005 and ITRF2000 by an unacceptably large speed of 1.8 mm yr−1. The velocity of Earth's centre in GEODVEL lies between that of ITRF2000 and that of ITRF2005, with the distance from ITRF2005 being about twice that from ITRF2000. Because the GEODVEL estimates of uncertainties in plate angular velocities account for uncertainty in the velocity of Earth's centre, they are more realistic than prior estimates of uncertainties. GEODVEL differs significantly from all prior global sets of relative plate angular velocities determined from space geodesy. For example, the 95 per cent confidence limits for the angular velocities of GEODVEL exclude those of REVEL (Sella et al.) for 34 of the 36 plate pairs that can be formed between any two of the nine plates with the best-constrained motion. The median angular velocity vector difference between GEODVEL and REVEL is 0.028° Myr−1, which is up to 3.1 mm yr−1 on Earth's surface. GEODVEL differs the least from the geodetic angular velocities that Altamimi et al. determine from ITRF2005. GEODVEL's 95 per cent confidence limits exclude 11 of 36 angular velocities of Altamimi et al., and the median difference is 0.015° Myr−1. GEODVEL differs significantly from nearly all relative plate angular velocities averaged over the past few million years, including those of NUVEL-1A. The difference of GEODVEL from updated 3.2 Myr angular velocities is statistically significant for all but two of 36 angular velocities with a median difference of 0.063° Myr−1. Across spreading centres, eight have slowed down while only two have sped up. We conclude that plate angular velocities over the past few decades differ significantly from the corresponding angular velocity averaged over the past 3.2 Myr.

Regression analysis of MCS intensity and ground motion parameters in Italy and its application in ShakeMap

Abstract: SUMMARY In Italy, the Mercalli–Cancani–Sieberg (MCS) is the intensity scale in use to describe the level of earthquake ground shaking, and its subsequent effects on communities and on the built environment. This scale differs to some extent from the Mercalli Modified scale in use in other countries and adopted as standard within the USGS-ShakeMap procedure to predict intensities from observed instrumental data. We have assembled a new PGM/MCS-intensity data set from the Italian database of macroseismic information, DBMI04, and the Italian accelerometric database, ITACA. We have determined new regression relations between intensities and PGM parameters (acceleration and velocity). Since both PGM parameters and intensities suffer of consistent uncertainties we have used the orthogonal distance regression technique. The new relations are and Tests designed to assess the robustness of the estimated coefficients have shown that single-line parametrizations for the regression are sufficient to model the data within the model uncertainties. The relations have been inserted in the Italian implementation of the USGS-ShakeMap to determine intensity maps from instrumental data and to determine PGM maps from the sole intensity values. Comparisons carried out for earthquakes where both kinds of data are available have shown the general effectiveness of the relations.

A correlation-based misfit criterion for wave-equation traveltime tomography

Abstract: Wave-equation traveltime tomography tries to obtain a subsurface velocity model from seismic data, either passive or active, that explains their traveltimes. A key step is the extraction of traveltime differences, or relative phase shifts, between observed and modelled finite-frequency waveforms. A standard approach involves a correlation of the observed and measured waveforms. When the amplitude spectra of the waveforms are identical, the maximum of the correlation is indicative of the relative phase shift. When the amplitude spectra are not identical, however, this argument is no longer valid. We propose an alternative criterion to measure the relative phase shift. This misfit criterion is a weighted norm of the correlation and is less sensitive to differences in the amplitude spectra. For practical application it is important to use a sensitivity kernel that is consistent with the way the misfit is measured. We derive this sensitivity kernel and show how it differs from the standard banana–doughnut sensitivity kernel. We illustrate the approach on a cross-well data set.

Palaeolatitude and age of the Indo–Asia collision: palaeomagnetic constraints

Abstract: SUMMARY Ongoing controversies on the timing and kinematics of the Indo–Asia collision can be solved by palaeomagnetically determined palaeolatitudes of terranes bounding the Indo–Asia suture zone. We show here, based on new palaeomagnetic data from the Linzizong volcanic rocks (54–47 Ma) near the city of Lhasa, that the latitude of the southern margin of Asia was 22.8 ± 4.2°N when these rocks were deposited. This result, combined with revised palaeomagnetic results from the northernmost sedimentary units of Greater India and with apparent polar wander paths of India and Eurasia, palaeomagnetically constrain the collision to have occurred at 46 ± 8 Ma (95 per cent confidence interval). These palaeomagnetic results are consistent with tomographic anomalies at 15–25°N that are interpreted to locate the Tethyan oceanic slab that detached following collision, and with independent 56–46 Ma collision age estimates inferred from the timing of slowing down of India, high pressure metamorphism, the end of marine sedimentation and the first occurrence of suture zone and arc detritus on the Greater Indian margin. When compared with apparent polar wander paths of India and Eurasia, the ∼46 Ma onset of collision at 22.8 ± 4.2°N implies 2900 ± 600 km subsequent latitudinal convergence between India and Asia divided into 1100 ± 500 km within Asia and 1800 ± 700 km within India.

Accelerating a three-dimensional finite-difference wave propagation code using GPU graphics cards

Abstract: We accelerate a three-dimensional finite-difference in the time domain (FDTD) wave propagation code by a factor between about 20 and 60 compared to a serial implementation using Graphics Processing Unit (GPU) computing on NVIDIA graphics cards with the CUDA programming language. We describe the implementation of the code in CUDA to simulate the propagation of seismic waves in a heterogeneous elastic medium. We also implement Convolution Perfectly Matched Layers (CPMLs) on the graphics cards to efficiently absorb outgoing waves on the fictitious edges of the grid. We show that the code that runs on a graphics card gives the expected results by comparing our results to those obtained by running the same simulation on a classical processor core. The methodology that we present can be used for Maxwell's equations as well because their form is similar to that of the seismic wave equation written in velocity vector and stress tensor.

A unified continuum representation of post-seismic relaxation mechanisms: semi-analytic models of afterslip, poroelastic rebound and viscoelastic flow

Abstract: We present a unified continuum mechanics representation of the mechanisms believed to be commonly involved in post-seismic transients such as viscoelasticity, fault creep and poroelasticity. The time-dependent relaxation that follows an earthquake, or any other static stress perturbation, is considered in a framework of a generalized viscoelastoplastic rheology whereby some inelastic strain relaxes a physical quantity in the material. The relaxed quantity is the deviatoric stress in case of viscoelastic relaxation, the shear stress in case of creep on a fault plane and the trace of the stress tensor in case of poroelastic rebound. In this framework, the instantaneous velocity field satisfies the linear inhomogeneous Navier's equation with sources parametrized as equivalent body forces and surface tractions. We evaluate the velocity field using the Fourier-domain Green's function for an elastic half-space with surface buoyancy boundary condition. The accuracy of the proposed method is demonstrated by comparisons with finite-element simulations of viscoelastic relaxation following strike-slip and dip-slip ruptures for linear and power-law rheologies. We also present comparisons with analytic solutions for afterslip driven by coseismic stress changes. Finally, we demonstrate that the proposed method can be used to model time-dependent poroelastic rebound by adopting a viscoelastic rheology with bulk viscosity and work hardening. The proposed method allows one to model post-seismic transients that involve multiple mechanisms (afterslip, poroelastic rebound, ductile flow) with an account for the effects of gravity, non-linear rheologies and arbitrary spatial variations in inelastic properties of rocks (e.g. the effective viscosity, rate-and-state frictional parameters and poroelastic properties).

A threshold-based earthquake early warning using dense accelerometer networks

Abstract: Most earthquake early warning systems (EEWS) developed so far are conceived as either ‘regional’ (network-based) or ‘on-site’ (stand-alone) systems. The recent implementation of nationwide, high dynamic range, dense accelerometer arrays makes now available, potentially in real time, unsaturated waveforms of moderate-to-large magnitude earthquakes recorded at very short epicentral distances (<10–20 km). This would allow for a drastic increase of the early warning lead-time, for example, the time between the alert notification and the arrival time of potentially destructive waves at a given target site. By analysing strong motion data from modern accelerograph networks in Japan, Taiwan and Italy, we propose an integrated regional/on-site early warning method, which can be used in the very first seconds after a moderate-to-large earthquake to map the most probable damaged zones. The method is based on the real-time measurement of the period (τ_c) and peak displacement (Pd) parameters at stations located at increasing distances from the earthquake epicentre. The recorded values of early warning parameters are compared to threshold values, which are set for a minimum magnitude 6 and instrumental intensity VII, according to the empirical regression analyses of strong motion data. At each recording site the alert level is assigned based on a decisional table with four alert levels defined upon critical values of the parameters Pd and τ_c, which are set according to the error bounds estimated on the derived prediction equations. Given a real time, evolutionary estimation of earthquake location from first P arrivals, the method furnishes an estimation of the extent of potential damage zone as inferred from continuously updated averages of the period parameter and from mapping of the alert levels determined at the near-source accelerometer stations. The off-line application of the method to strong motion records of the M_w 6.3, 2009 Central Italy earthquake shows a very consistent match between the rapidly predicted (within a few seconds from the first recorded P wave) and observed damage zone, the latter being mapped from detailed macroseismic surveys a few days after the event. The proposed approach is suitable for Italy, where, during the last two decades, a dense network of wide dynamic-range accelerometer arrays has been deployed by the Department of Civil Protection (DPC), the Istituto Nazionale di Geofisica e Vulcanologia (INGV) and other regional research agencies.

Ambient vibration analysis of an unstable mountain slope

Abstract: SUMMARY A field experiment with small aperture seismic arrays was performed on the unstable rock slope above the village of Randa in the southern Swiss Alps. The aim of this experiment was to constrain the seismic response of a potential future rockslide using ambient vibration recordings. Weak seismic events were identified on the recordings and site-to-reference spectral ratios were calculated using a reference site located on the stable part of the slope. Spectral ratios of up to 30 were observed at sites located within the unstable portion of the slope. A strong variation of spectral ratios with azimuth indicates a directional site effect. Neither amplification nor directionality were observed at sites located in the stable part of the slope. Furthermore, time-frequency polarization analysis of the ambient noise was performed to provide robust estimates of frequency dependent directions of the maximum polarization. It was found that the unstable part of the slope vibrates within a narrow range of directions (130 ± 10°) for the frequency range centred around 5 Hz. The polarization directions estimated from ambient seismic vibrations are in good agreement with the deformation directions obtained by geodetic and in situ measurements. No directionality of ambient vibrations was observed at sites within the stable part of the slope.

Extension on the Tibetan plateau: recent normal faulting measured by InSAR and body wave seismology

Abstract: We use InSAR and body wave modelling to determine the faulting parameters for a series of five Mw 5.9–7.1 normal faulting earthquakes that occurred during 2008, including the March 20 Yutian earthquake (Mw 7.1), one of the largest normal faulting events to have occurred recently on the continents. We also study three earlier normal faulting earthquakes that occurred in southern Tibet between 1992 and 2005. Coseismic deformation for each of these eight events is measured with ascending and descending interferograms from ENVISAT, ERS and ALOS SAR data. Elastic dislocation modelling of the line-of-sight InSAR displacements and body wave seismological modelling of P and SH waves are used to estimate fault parameters and are found to be in good agreement for all the events studied. The use of InSAR to measure deformation allows a relatively precise determination of the fault location in addition to resolving the focal plane ambiguity. Only five of the eight events are associated with a clear surface topographic feature, suggesting that an underestimation of the amount of extension would result from using the surface expressions of normal faulting alone. The observations, in all cases, are consistent with slip on planar surfaces, with dips in the range 40–50°, that penetrates the uppermost crust to a depth of 10–15 km. We find no evidence for active low-angle (dip less than 30°) normal faulting. The contribution of the normal faulting to overall extension estimated by summing seismic moments over earthquakes for the past 43 yr is 3–4 mm yr−1, or 15–20 per cent of the rates of extension measured across the plateau using GPS. 85 per cent of the moment release in normal faulting over the past 43 yr has occurred in regions whose surface height exceeds 5 km. This observation adds weight to the suggestions that the widespread normal faulting on the plateau is the result of variations in the gravitational potential energy of the lithosphere.

Seismic reflection character of the Hikurangi subduction interface, New Zealand, in the region of repeated Gisborne slow slip events

Abstract: SUMMARY We use seismic reflection data to map the geometry and character of the subduction interface in the Gisborne area of the Hikurangi subduction margin, New Zealand, which experiences repeated shallow (<15 km) slow slip events The reflection character and geometry in this area is highly variable, which we interpret to be related to the subduction of seamounts and underthrust sediments Three zones of high-amplitude interface reflectivity (HRZ-1, 2 and 3) are interpreted to be the result of fluid-rich sediments that have been entrained with subducting seamounts The interface above the HRZ zones is shallower than the surrounding areas by 2–4 km, due to the warping of the interface to accommodate seamount subduction These zones of high-amplitude reflectivity and shallower interface geometry correlate broadly with locations of recorded slow slip events from 2002 to 2008 We hypothesize that effective stresses on the interface may be lower along the northeast margin in areas of high-amplitude reflectivity due to; (1) the enhanced underthrusting of fluid-rich sediment, (2) reduced overburden stresses where the interface has been warped to shallower depths to accommodate seamount subduction and (3) potential fluid flow concentration effects leading to overpressure along these shallower interface corrugations From our observations we propose localized reductions in effective stress caused by interface structural relief may be a potential factor in promoting shallow slow slip events

Receiver function inversion by trans‐dimensional Monte Carlo sampling

Abstract: SUMMARY A key question in the analysis of an inverse problem is the quantification of the non-uniqueness of the solution. Non-uniqueness arises when properties of an earth model can be varied without significantly worsening the fit to observed data. In most geophysical inverse problems, subsurface properties are parameterized using a fixed number of unknowns, and non-uniqueness has been tackled with a Bayesian approach by determining a posterior probability distribution in the parameter space that combines ‘a priori’ information with information contained in the observed data. However, less consideration has been given to the question whether the data themselves can constrain the model complexity, that is the number of unknowns needed to fit the observations. Answering this question requires solving a trans-dimensional inverse problem, where the number of unknowns is an unknown itself. Recently, the Bayesian approach to parameter estimation has been extended to quantify the posterior probability of the model complexity (the number of model parameters) with a quantity called ‘evidence’. The evidence can be hard to estimate in a non-linear problem; a practical solution is to use a Monte Carlo sampling algorithm that samples models with different number of unknowns in proportion to their posterior probability. This study presents a method to solve in trans-dimensional fashion the non-linear inverse problem of inferring 1-D subsurface elastic properties from teleseismic receiver function data. The Earth parameterization consists of a variable number of horizontal layers, where little is assumed a priori about the elastic properties, the number of layers, and and their thicknesses. We developed a reversible jump Markov Chain Monte Carlo algorithm that draws samples from the posterior distribution of Earth models. The solution of the inverse problem is a posterior probability distribution of the number of layers, their thicknesses and the elastic properties as a function of depth. These posterior distributions quantify completely the non-uniqueness of the solution. We illustrate the algorithm by inverting synthetic and field measurements, and the results show that the data constrain the model complexity. In the synthetic example, the main features of the subsurface properties are recovered in the posterior probability distribution. The inversion results for actual measurements show a crustal structure that agrees with previous studies in both crustal thickness and presence of intracrustal low-velocity layers.

Geophysical muon imaging: feasibility and limits

Abstract: SUMMARY We study the possibility of muon radiography as a tool to investigate space and time changes in the internal density distribution inside geological structures. Previous work has shown the practical applicability of this method. Nevertheless, quantitative information on factors which impose limitations on it are still sorely lacking in the literature. We discuss the main issues that can influence the final result of a geophysical imaging experiment. In particular, with the view of optimizing the signal-to-noise ratio, we address issues concerning (i) the energy spectrum for muons arriving at different zenith angles, (ii) the muon propagation model through matter and (iii) the characteristics of the muon detector (telescope) that we have designed to perform experiments of muon radiography against the harsh environment usually encountered in the active zone of a volcano. We thus identify factors that can induce either static or dynamic effects and that should be taken into account. We also define a feasibility eq. (32) relating the geometrical characteristics of the telescope and the duration of the experiment to the expected density resolution, in turn a function of the geometrical characteristics of the target structure. This relation is especially important to define the applicability domain of muon radiography and it is utilized to test the suitability of the method to investigate the density distribution inside some candidate target structures.

Earthquake mechanisms and active tectonics of the Hellenic subduction zone

Abstract: SUMMARY We use improved focal mechanisms and centroid depth estimates of earthquakes, combined with GPS velocities, to examine the tectonics of the Hellenic subduction zone, and in particular the processes occurring at both ends of the Hellenic Arc. Nubia-Aegean convergence is accommodated by shallowly dipping thrust-faulting along the subduction–zone interface, as well as by steeper splay faults in the overriding material. From a comparison of observed and expected seismic moment release over the last 100 yr, combined with existing knowledge of the longer-term documented historical record, we confirm earlier suggestions that most (80 per cent) of this convergence is accommodated aseismically, that is, that the subduction zone is uncoupled. This conclusion is robust, even allowing for rare very large earthquakes on splay faults, such as that of AD 365, and also allowing for the contribution of small earthquakes. The downgoing Nubian plate deforms by arc-parallel contraction at all depths, from 200 km seaward of Crete to at least 100 km within the subducting slab. Extensional (T) axes of earthquakes are aligned downdip within the descending slab suggesting that, even if the aseismic prolongation of the slab has reached the 670 km mantle discontinuity, it does not transmit stresses to shallower depths. Shallow thrust-faulting earthquakes on the subduction interface show a divergence of slip vectors round the arc, and GPS measurements show that this is accommodated mainly by E–W extension on normal faults in the overriding Aegean material. The eastern end of the subduction zone, south of Rhodes, displays distributed deformation in the overriding material, including a mixture of strike-slip and splay-thrust faulting, and probably involves rotations about a vertical axes. Here slip on the interface itself is by thrust faulting with slip vectors oblique to the arc but parallel to the overall Nubia-Aegean convergence: there is no evidence for slip-partitioning in the traditional sense. In the west, the subduction zone terminates in a distributed zone of parallel NE–SW strike-slip faults, of which the most prominent is the Kefalonia Transform Fault (KTF). A flexural gravity anomaly confirms that the deep bathymetric escarpment of the KTF is a lateral ramp, formed as the Ionian islands are emplaced SW onto the Apulian lithosphere, and enhanced by minor thrust faulting with slip vectors perpendicular to the scarp. Distributed parallel strike-slip faults both SW and NE of mainland central Greece terminate in E–W graben in central Greece, which accommodate the overall NE–SW shear by clockwise block rotation. Central Greece therefore acts as a relay zone between the strike-slip faulting of the NE Aegean and the Ionian Islands–western Peloponnese.

Detecting seasonal variations in seismic velocities within Los Angeles basin from correlations of ambient seismic noise

Abstract: SUMMARY We analyse 3 yr of continuous seismic records from broad-band stations of the Caltech Regional Seismic Network (CI) in vicinity of the Los Angeles basin. Using correlations of ambient seismic noise, relative velocity variations in the order of 0.1 per cent can be measured between all interstation pairs. It is the first time that such an extensive study between 861 interstation pairs over such a large area has been carried out. We perform these measurements using the ‘stretching’ technique, assuming that one of the two waveforms is merely a stretched version of the other. Obviously this assumption is always violated and the two waveforms are generally decorrelated due to temporal changes in the Earth crust, due to different sources or simply because the cross-correlations are not fully converged. We investigate the stability of these measurements by repeating each measurement over various time-windows of equal length. On average between all interstation pairs in the Los Angeles basin a seasonal signal in the relative velocity variation is observed, with peaks and troughs during winter and summer time, respectively. Generally the observed signal decreases with increasing interstation distance and relative traveltime perturbations can only be measured up to an interstation distance of 60 km. Furthermore, the traveltime perturbations do not depend on azimuth of station pairs, suggesting that they are not related to seasonal variations of the noise sources. Performing a simple regionalization by laterally averaging measurements over a subset of stations we found the sedimentary basin showing the most consistent signal and conclude that the observed seasonality might be induced either by changes in the ground-water aquifer or thermo-elastic strain variations that persist down to a depth of 15–22 km.

Elastic moduli evolution and accompanying stress changes with increasing crack damage: implications for stress changes around fault zones and volcanoes during deformation

Abstract: SUMMARY The elastic moduli of rock in areas susceptible to crack damage, such as within fault zones or volcanic edifices, can be subject to large modifications. Knowledge of how elastic moduli may vary in such situations is important for both the reliable modelling of volcano deformation and stability and for linear and non-linear elastic crack models for earthquake rupture. Furthermore, it has previously been shown that changes in elastic moduli can induce changes in the stress field surrounding faults. Here we report both uniaxial experimental measurements of changes in elastic moduli during increasing-amplitude cyclic stressing experiments on a range of different rock types (basalts, sandstones and granite) and the results of modelled stress modifications. The trend in elastic moduli evolution with increasing damage was remarkably similar for each rock type, with the exception of essentially crack-free intrusive basalt that exhibited very minor changes. In general, Young’s modulus decreased by between 11 and 32 per cent and Poisson’s ratio increased by between 72 and 600 per cent over the total sequence of loading cycles. These changes are attributed to an increasing level of anisotropic crack damage within the samples. Our results also show that acoustic emission (AE) output during any loading cycle only commenced when new crack damage was generated. This corresponded to the level of stress where AE ceased during the unloading portion of the previous cycle. Using the multilayer elastic model of Faulkner et al. we demonstrate that the damage-induced changes in elastic moduli also result in significant decreases in differential stress, increases in mean stress and rotation of the applied greatest principal stress relative to the orientation of the mechanical layering. The similar trend in the evolution of the elastic moduli of all the rocks tested suggests that stress modification in the damage zone of faults might take the same form, regardless of the lithology through which the fault runs. These observations are discussed in terms of their applicability to both fault zones and deformation at volcanoes.

An hp-adaptive discontinuous Galerkin finite-element method for 3-D elastic wave modelling

Abstract: We present a discontinuous Galerkin finite-element method (DG-FEM) formulation with Convolutional Perfectly Matched Layer (CPML) absorbing boundary condition for 3-D elastic seismic wave modelling. This method makes use of unstructured tetrahedral meshes locally refined according to the medium properties (h-adaptivity), and of approximation orders that can change from one element to another according to an adequate criterion (p-adaptivity). These two features allow us to significantly reduce the computational cost of the simulations. Moreover, we have designed an efficient CPML absorbing boundary condition, both in terms of absorption and computational cost, by combining approximation orders in the numerical domain. A quadratic interpolation is typically used in the medium to obtain the required accuracy, while lower approximation orders are used in the CPMLs to reduce the total computational cost and to obtain a well-balanced workload over the processors. While the efficiency of DG-FEMs have been largely demonstrated for high approximation orders, we favour the use of low approximation orders as they are more appropriate to the applications we are interested in. In particular, we address the issues of seismic modelling and seismic imaging in cases of complex geological structures that require a fine discretization of the medium. We illustrate the efficiency of our approach within the framework of the EUROSEISTEST verification and validation project, which is designed to compare high-frequency (up to 4 Hz) numerical predictions of ground motion in the Volvi basin (Greece). Through the tetrahedral meshing, we have achieved fine discretization of the basin, which appears to be a sine qua non condition for accurate computation of surface waves diffracted at the basin edges. We compare our results with predictions computed with the spectral element method (SEM), and demonstrate that our method yields the same level of accuracy with computation times of the same order of magnitude.

Automated determination of P-phase arrival times at regional and local distances using higher order statistics

Abstract: P>We present an algorithm for automatic P-phase arrival time determination for local and regional seismic events based on higher order statistics (HOS). Using skewness or kurtosis a characteristic function is determined to which a new iterative picking algorithm is applied. For P-phase identification we apply the Akaike Information Criterion to the characteristic function, while for a precise determination of the P-phase arrival time a pragmatic picking algorithm is applied to a recalculated characteristic function. In addition, an automatic quality estimate is obtained, based on the slope and the signal-to-noise ratio, both calculated from the characteristic function. To get rid of erroneous picks, a Jackknife procedure and an envelope function analysis is used. The algorithm is applied to a large data set with very heterogeneous qualities of P-onsets acquired by a temporary, regional seismic network of the EGELADOS-project in the southern Aegean. The reliability and robustness of the proposed algorithm is tested by comparing more than 3000 manually derived P readings, serving as reference picks, with the corresponding automatically estimated P-wave arrival times. We find an average deviation from the reference picks of 0.26 +/- 0.64 s when using kurtosis and 0.38 +/- 0.75 s when using skewness. If automatically as excellent classified picks are considered only, the average difference from the reference picks is 0.07 +/- 0.31 s and 0.07 +/- 0.41 s, respectively. However, substantially more P-arrival times are determined when using kurtosis, indicating that the characteristic function derived from kurtosis estimation is to be preferred. Since the characteristic function is calculated recursively, the algorithm is very fast and hence suited for earthquake early warning purposes. Furthermore, a comparative study with automatically derived P-readings using Allen's and Baer & Kradolfer's picking algorithms applied to the same data set demonstrates better quantitative and qualitative performance of the HOS approach. This study shows, that precise automatic P-onset determination is feasible, even when using data sets with very heterogeneous signal-to-noise ratio.

GPS and gravity constraints on continental deformation in the Alborz mountain range, Iran

Abstract: A network of 54 survey GPS sites, 28 continuous GPS stations and three absolute gravity (AG) observation sites have been set up in the Alborz mountain range to quantify the present-day kinematics of the range. Our results allow us to accurately estimate the motion of the South Caspian block (SCB) for the first time, and indicate rotation of the SCB relative to Eurasia, accounting for the left lateral motion in the Alborz range. In light of these new results, it clearly appears that deformation rates vary along the range, the eastern part accommodating mainly left lateral strike slip (2 mm yr(-1) south of the range and 5 mm yr(-1) north of the range) with a very low range normal shortening rate on the Khazar thrust fault (similar to 2 mm yr(-1)), and the western part accommodating range normal shortening (similar to 6 mm yr-1) on the Khazar thrust fault with a left lateral component of similar to 2 mm yr(-1) north of the range and 1 mm yr(-1) south of the range. These present-day kinematics agree with geomorphologic estimated slip rates, but not the long-term deformation, corroborating the idea that the kinematics of the range have changed recently due to the change of SCB motion.;Modelling of the interseismic deformation suggests a deep locking depth on the central-western segment of the Khazar fault (similar to 30 km) in agreement with the Baladeh earthquake rupture and aftershock ranging between 10 and 30 km. Given this unusual deep locking depth and the 34 degrees dip of the thrust, a large part of the Alborz range is located above the seismically coupled part of the fault. Based on our AG measurements this part of the range seems to uplift at a rate of 1-5 mm yr(-1), in agreement with terrace uplift.

A Preliminary Investigation into the Relationship between Long-Period Seismic Noise and Local Fluctuations in the Atmospheric Pressure Field

Abstract: Summary The displacement response of an elastic half space to a plane pressure wave is examined in order to establish the conditions under which sources of this type can contribute significantly to the long-period seismic noise field. The study is restricted to pressure waves which propagate at velocities well below the seismic wave velocities characteristic of the half space. The numerical studies indicate that pressure waves with amplitudes of 100 pbar or more can contribute significantly to the long-period vertical background noise observed at the surface, provided that the detectors are located on sections of alluvial fill or poorly to moderately indurated sandstones and shales whose thicknesses are greater than about a kilometre. These same waves can also create significant tilt noise on long-period horizontal seismographs located at or near the surface, regardless of the rock type. The seismic disturbances created by pressure waves decay rapidly away from the surface. Therefore, it appears that it may be possible to eliminate the effects of atmospherically generated noise by placing the detectors at moderate depths.

Bending mechanics and mode selection in free subduction: a thin-sheet analysis

Abstract: SUMMARY To elucidate the dynamics of free (buoyancy-driven) subduction of oceanic lithosphere, I study a model in which a 2-D sheet of viscous fluid with thickness h and viscosity γη1 subducts in an infinitely deep ambient fluid with viscosity η1. Numerical solutions for the sheet's evolution are obtained using the boundary-element method (BEM), starting from an initial configuration comprising a short ‘protoslab’ attached to a longer horizontal ‘plate’ that is free to move laterally beneath an impermeable traction-free surface. Interpretation of the solutions using thin viscous sheet theory shows that the fundamental length scale controlling the subduction is the ‘bending length’lb, defined at each instant as the length of the portion of the sheet's midsurface where the rate of change of curvature is significant. Geophysically speaking, lb is the sum of the lengths of the slab and of the region seaward of the trench where flexural bulging occurs. The bending length in turn enters into the definition of the sheet's dimensionless ‘stiffness’S≡γ(h/lb)3, which controls whether the sinking speed of the slab is determined by the viscosity of the sheet itself (S≫ 1) or by that of the ambient fluid (S≤ 1). Motivated by laboratory observations of different modes of subduction (retreating versus advancing trench, folding versus no folding, etc.) in fluid layers with finite depth, I calculate numerically the dip θD of the slab's leading end as a function of γ and the normalized depth D/h to which it has penetrated. The contours of the function θD(γ, D/h) strongly resemble the intermode boundaries in the laboratory-based regime diagram of Schellart, supporting the hypothesis that the mode of subduction observed at long times in experiments is controlled by the dip of the slab's leading end when it reaches the bottom of the layer. In particular, the BEM solutions explain why trenches advance in the laboratory only when γ lies in an intermediate range, and why they retreat when γ is either smaller or larger than this. Application of the BEM model to Wu et al.'s compilation of the minimum curvature radii of subducted slabs suggests γ∈[140, 510] for the Earth. This is too small to permit the laboratory-type ‘trench advancing’ mode, in agreement with the lack of tomographic evidence for slabs that are ‘bent over backwards’.

Spectral induced polarization of partially saturated clay-rocks: a mechanistic approach

Abstract: We have developed a mechanistic model to interpret spectral induced polarization data of partially saturated clay-rocks. This model accounts for the polarization of the grains through an electrical double layer model with a polarization model of the inner part of the electrical double layer called the Stern layer. The polarization model accounts also for the Maxwell-Wagner polarization at frequencies higher than 100 Hz. The Maxwell-Wagner polarization is modelled by using a conductivity model modified to account for the presence of a non-wetting immiscible phase like air in the pore space. The resulting model is consistent with the first and second Archie's laws in the case where surface conductivity can be neglected. The volumetric charge density of the diffuse layer at saturation is divided by the saturation of the water phase to account for the partial water saturation of the porous material. The model comprises seven fundamental parameters: the formation factor, the second Archie's exponent, a critical water saturation level, the mean electrical potential of the pore space at saturation, the density of the counterions in the Stern layer, and at least two parameters describing the grain size distribution. Most of these parameters can be derived independently using alternative measurements and electrochemical models. Measurements were performed in the frequency range 10 mHz-45 kHz using five samples from the Callovo-Oxfordian formation in the eastern part of the Paris Basin, France. The model agrees fairly well with the experimental data at saturation and for partially saturated clay-rocks down to 1 Hz. Most of the seven physical parameters entering the model were independently evaluated.

Constraining models of postglacial rebound using space geodesy; a detailed assessment of model ICE-5G VM2 and its relatives

Abstract: SUMMARY Using global positioning system, very long baseline interferometry, satellite laser ranging and Doppler Orbitography and Radiopositioning Integrated by Satellite observations, including the Canadian Base Network and Fennoscandian BIFROST array, we constrain, in models of postglacial rebound, the thickness of the ice sheets as a function of position and time and the viscosity of the mantle as a function of depth. We test model ICE-5G VM2 T90 Rot, which well fits many hundred Holocene relative sea level histories in North America, Europe and worldwide. ICE-5G is the deglaciation history having more ice in western Canada than ICE-4G; VM2 is the mantle viscosity profile having a mean upper mantle viscosity of 0.5 × 1021 Pa s and a mean uppermost-lower mantle viscosity of 1.6 × 1021 Pa s; T90 is an elastic lithosphere thickness of 90 km; and Rot designates that the model includes (rotational feedback) Earth's response to the wander of the North Pole of Earth's spin axis towards Canada at a speed of ≈1° Myr−1. The vertical observations in North America show that, relative to ICE-5G, the Laurentide ice sheet at last glacial maximum (LGM) at ≈26 ka was (1) much thinner in southern Manitoba, (2) thinner near Yellowknife (Northwest Territories), (3) thicker in eastern and southern Quebec and (4) thicker along the northern British Columbia–Alberta border, or that ice was unloaded from these areas later (thicker) or earlier (thinner) than in ICE-5G. The data indicate that the western Laurentide ice sheet was intermediate in mass between ICE-5G and ICE-4G. The vertical observations and GRACE gravity data together suggest that the western Laurentide ice sheet was nearly as massive as that in ICE-5G but distributed more broadly across northwestern Canada. VM2 poorly fits the horizontal observations in North America, predicting places along the margins of the Laurentide ice sheet to be moving laterally away from the ice centre at 2 mm yr−1 in ICE-4G and 3 mm yr−1 in ICE-5G, in disagreement with the observation that the interior of the North American Plate is deforming more slowly than 1 mm yr−1. Substituting VM5a T60 for VM2 T90, that is, introducing into the lithosphere at its base a layer with a high viscosity of 10 × 1021 Pa s, greatly improves the fit of the horizontal observations in North America. ICE-4G VM5a T60 Rot predicts most of the North American Plate to be moving horizontally more slowly than ≈1 mm yr−1, in agreement with the data. ICE-5G VM5a T60 Rot well fits both the vertical and horizontal observations in Europe. The space geodetic data cannot distinguish between models with and without rotational feedback, in the vertical because the velocity of Earth’ centre is uncertain, and in the horizontal because the areas of the plate interiors having geodetic sites is not large enough to detect the small differences in the predictions of rotational feedback going across the plate interiors.

Advanced insights into magmatism and volcanism of the Mozambique Ridge and Mozambique Basin in the view of new potential field data

Abstract: SUMMARY A new plate tectonic model is presented which describes the emplacement of the Mozambique Ridge off southeast Africa as the result of long lasting volcanic activity (140–122 Ma) during the initial opening between Africa and Antarctica. Thus, an oceanic origin for the Mozambique Ridge is proposed. This model is based on a new and systematic high resolution magnetic anomaly data set acquired across the Mozambique Ridge and throughout the Mozambique Basin. Data from the Mozambique Basin allow the identification of Mesozoic magnetic anomalies from M0r to M26 (124.61–155.3 Ma) with previously unmatched accuracy. Small-scale fracture zones are recognized by offsets of magnetic anomalies in the westernmost part of the basin. Additionally, a bend in the major fracture zones ‘F’ and ‘E’ between M17r and M18n (142.84–144.04 Ma) and a recognized sinusoidal change in spreading direction with an amplitude of about 15° indicate that the basin experienced several small scale changes in spreading direction through time. A maximum change in spreading direction to almost 0° at around M11n (135.69 Ma) is followed by a short lived increase in spreading half rate from 23.5 km Ma−1 to about 27.5 km Ma−1 in the time frame from M10r to M9n (134.30–132.83 Ma). We propose that this is related to the initial opening of the South Atlantic Ocean represented by the onset of seafloor spreading between the Falkland Plateau and Africa in the conjugate Georgia and Natal basins. Across the Mozambique Ridge, high amplitude magnetic anomalies at major structural boundaries suggest that the different plateaus of the ridge were formed at different times. A simple 2-D gravity and magnetic model for the ridge supports the hypothesis of multiple volcanic episodes which formed the ridge though long lasting volcanic activity between about 140 and 122 Ma. Together with new and verified rotation parameters from the Mozambique Basin and its conjugate, the Riiser-Larsen Sea, Antarctica, a series of plate tectonic reconstructions are presented which demonstrate when and how the different parts of the ridge evolved through time.

Spatial distribution of scattering loss and intrinsic absorption of short‐period S waves in the lithosphere of Japan on the basis of the Multiple Lapse Time Window Analysis of Hi‐net data

Abstract: SUMMARY Propagation of short period S waves through the crust is strongly controlled by scattering loss and intrinsic absorption. Scattering loss controls the shape of the seismic envelope and the spatial distribution of the energy while intrinsic absorption produces an exponential decay as a function of time which is independent of the coordinates. Since both parameters have a different effect, it is possible to compute them by means of the Multiple Lapse Time Window Analysis (MLTWA). In this paper, high resolution maps of Japan for both seismic S-wave attenuation parameters are obtained by using the MLTWA with data provided by the Hi-net seismic network. The maps show strong regional variation of the parameters. The variations depend mainly on the tectonic setting of each region and volcanism mechanisms. In the scattering loss maps, in the 1–2 Hz map, there is a very different behaviour between the northeast (stronger scattering loss) and the southwest (much weaker scattering loss). In the northeast of Japan, areas with strong scattering loss correspond mainly to the location of the volcanic arcs. In Southwest Japan, Chugoku region and the west side of Kii peninsula show strong scattering levels for the 4–8 Hz and higher frequency bands. In Hokkaido, strong scattering loss and intrinsic absorption can be observed for the lower frequency bands. In the Sendai plain, we find strong scattering loss for all the frequency bands. Under all these areas low velocity anomalies have been detected in previous works. In Tohoku area the east side shows lower intrinsic attenuation than the west side. This characteristic seems to extend to the southern part of Hokkaido and to central Japan. In Kanto area, the values of intrinsic absorption are lower than in the surrounding areas; but in the region around the Itoigawa-Shizuoka tectonic line and in its western side strong absorption areas are observed. The behaviour of each area for intrinsic absorption is similar to the behaviour observed in the Q−1c maps. The MLTWA is based on the hypothesis of multiple isotropic scattering in a medium with the homogeneous distribution of scattering mean free path and intrinsic absorption. Although these hypotheses are very simple, the results are very informative on the characteristics of each region, and show the usefulness of studying the properties of scattered seismic waves in the understanding of the properties of the crust.

Retrieval of Moho‐reflected shear wave arrivals from ambient seismic noise

Abstract: Theoretical studies on ambient seismic noise (ASN) predict that complete Green's function between seismic stations can be retrieved from cross correlation. However, only fundamental mode surface waves emerge in most studies involving real data. Here we show that Moho-reflected body wave (SmS) and its multiples can be identified with ASN for station pairs near their critical distances in the short period band (1–5 s). We also show that an uneven distribution of noise sources, such as mining activity and wind–topography interaction, can cause surface wave precursors, which mask weaker body wave phases.