Showing papers in "Geophysical Journal International in 2008"

Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps

Abstract: SUMMARY We present the results of Rayleigh wave and Love wave phase velocity tomography in the western United States using ambient seismic noise observed at over 250 broad-band stations from the EarthScope/USArray Transportable Array and regional networks. All available threecomponent time-series for the 12-month span between 2005 November 1 and 2006 October 31 have been cross-correlated to yield estimated empirical Rayleigh and Love wave Green’s functions. The Love wave signals were observed with higher average signal-to-noise ratio (SNR) than Rayleigh wave signals and hence cannot be fully explained by the scattering of Rayleigh waves. Phase velocity dispersion curves for both Rayleigh and Love waves between 5 and 40 speriod were measured for each interstation path by applying frequency‐time analysis. The average uncertainty and systematic bias of the measurements are estimated using a method based on analysing thousands of nearly linearly aligned station-triplets. We find that empirical Green’s functions can be estimated accurately from the negative time derivative of the symmetric component ambient noise cross-correlation without explicit knowledge of the source distribution. The average traveltime uncertainty is less than 1 s at periods shorter than 24 s. We present Rayleigh and Love wave phase speed maps at periods of 8, 12, 16,and 20 s. The maps show clear correlations with major geological structures and qualitative agreement with previous results based on Rayleigh wave group speeds.

Surface wave array tomography in SE Tibet from ambient seismic noise and two-station analysis - II. Crustal and upper-mantle structure

Abstract: SUMMARY We determine the 3-D shear wave speed variations in the crust and upper mantle in the southeastern borderland of the Tibetan Plateau, SW China, with data from 25 temporary broad-band stations and one permanent station. Interstation Rayleigh wave (phase velocity) dispersion curves were obtained at periods from 10 to 50 s from empirical Green’s function (EGF) derived from (ambient noise) interferometry and from 20 to 150 s from traditional two-station (TS) analysis. Here, we use these measurements to construct phase velocity maps (from 10 to 150 s, using the average interstation dispersion from the EGF and TS methods between 20 and 50 s) and estimate from them (with the Neighbourhood Algorithm) the 3-D wave speed variations and their uncertainty. The crust structure, parametrized in three layers, can be well resolved with a horizontal resolution about of 100 km or less. Because of the possible effect of mechanically weak layers on regional deformation, of particular interest is the existence and geometry of low (shear) velocity layers (LVLs). In some regions prominent LVLs occur in the middle crust, in others they may appear in the lower crust. In some cases the lateral transition of shear wave speed coincides with major fault zones. The spatial variation in strength and depth of crustal LVLs suggests that the 3-D geometry of weak layers is complex and that unhindered crustal flow over large regions may not occur. Consideration of such complexity may be the key to a better understanding of relative block motion and patterns of seismicity.

Non-parametric seismic data recovery with curvelet frames

Abstract: SUMMARY Seismic data recovery from data with missing traces on otherwise regular acquisition grids forms a crucial step in the seismic processing flow. For instance, unsuccessful recovery leads to imaging artefacts and to erroneous predictions for the multiples, adversely affecting the performance of multiple elimination. A non-parametric transform-based recovery method is presented that exploits the compression of seismic data volumes by recently developed curvelet frames. The elements of this transform are multidimensional and directional and locally resemble wave fronts present in the data, which leads to a compressible representation for seismic data. This compression enables us to formulate a new curvelet-based seismic data recovery algorithm through sparsity-promoting inversion. The concept of sparsity-promoting inversion is in itself not new to geophysics. However, the recent insights from the field of ‘compressed sensing’ are new since they clearly identify the three main ingredients that go into a successful formulation of a recovery problem, namely a sparsifying transform, a sampling strategy that subdues coherent aliases and a sparsity-promoting program that recovers the largest entries of the curvelet-domain vector while explaining the measurements. These concepts are illustrated with a stylized experiment that stresses the importance of the degree of compression by the sparsifying transform. With these findings, a curvelet-based recovery algorithms is developed, which recovers seismic wavefields from seismic data volumes with large percentages of traces missing. During this construction, we benefit from the main three ingredients of compressive sampling, namely the curvelet compression of seismic data, the existence of a favourable sampling scheme and the formulation of a large-scale sparsity-promoting solver based on a cooling method. The recovery performs well on synthetic as well as real data by virtue of the sparsifying property of curvelets. Our results are applicable to other areas such as global seismology.

Waveform inversion in the Laplace domain

Abstract: SUMMARY For the last 30 yr, since Tarantola’s pioneering theoretical study of waveform inversion, geophysicists and applied mathematicians have utilized a waveform inversion to delineate the earth’s structures. However, successful applications of waveform inversion to real data are nominal. The failures are mainly caused by the high non-linearity of the waveform inversion and the real data containing insufficient low-frequency components. We propose a waveform inversion algorithm that is robust and is not sensitive to the initial model by exploiting the wavefield in the Laplace domain and the adjoint property of the wave equation. The wavefield in the Laplace domain is equivalent to the zero frequency component of the damped wavefield. Therefore, the inversion of Poisson’s equation in electrical prospecting can be viewed as a waveform inversion problem, exploiting the zero frequency component of an undamped wavefield. Since our inversion algorithm in the Laplace domain is strongly associated with the nature of DC inversion for Poisson’s equation in electrical prospecting, our algorithm can generate a velocity model that is equivalent to a long-wavelength velocity model (a smooth velocity model). Through numerical tests, we found that the Laplace-transformed wavefields for optimally low Laplace damping constants are critically needed in our algorithm as the low-frequency components are necessary in the waveform inversion of the frequency domain. We note that the objective function by l2 norm of the logarithmic wavefield in the Laplace domain behaves as if it has no local minimum points in low and high Laplace damping constants. Moreover, we note that the forward modelling in the Laplace domain could be accurately calculated by using a coarser grid of several hundreds of metres than that of the frequency domain. Numerical tests of the salt dome model with high-velocity contrast demonstrate the robustness of our algorithm to both synthetic and real data. To apply our algorithm to real data more successfully, we need to improve the accuracy of the Laplace transformation for the wavefields in the time domain. However, our waveform inversion in the Laplace domain can be successfully applied to real data since our algorithm has more advantages in forward modelling and inversion than those in the frequency domain.

Tectonic evolution of the North Sea basin: crustal stretching and subsidence

Abstract: Summary. The lithospheric stretching model for the formation of sedimentary basins was tested in the central North Sea by a combined study of crustal thinning and basement subsidence patterns. A profile of crustal structure was obtained by shooting a long-range seismic experiment across the Central Graben, the main axis of subsidence. A seabed array of 12 seismometers in the graben was used to record shots fired in a line 530 km long across the basin. The data collected during the experiment were interpreted by modelling synthetic seismograms from a laterally varying structure, and the final model showed substantial crustal thinning beneath the graben. Subsidence data from 19 exploration wells were analysed to obtain subsidence patterns in the central North Sea since Jurassic times. Changes in water depth were quantified using foraminiferal assemblages where possible, and observed basement subsidence paths were corrected for sediment loading, compaction and changes in water depth through time. The seismic model is shown to be compatible with the observed gravity field, and the small size of observed gravity anomalies is used to argue that the basin is in local isostatic equilibrium. Both crustal thinning and basement subsidence studies indicate about 70 km of stretching across the Central Graben during the mid-Jurassic to early Cretaceous extensional event. This extension appears to have occurred over crust already slightly thinned beneath the graben, and the seismic data suggest that total extension since the early Permian may have been more than 100km. The data presented here may all be explained using a simple model of uniform extension of the lithosphere.

Shear and compressional velocity models of the mantle from cluster analysis of long-period waveforms

Abstract: SUMMARY We present a new technique for the efficient measurement of the traveltimes of long period body wave phases. The technique is based on the fact that all arrivals of a particular seismic phase are remarkably similar in shape for a single event. This allows the application of cross-correlation techniques that are usually used in a regional context to measure precise global differential times. The analysis is enhanced by the inclusion of a clustering algorithm that automatically clusters waveforms by their degree of similarity. This allows the algorithm to discriminate against unusual or distorted waveforms and makes for an extremely efficient measurement technique. This technique can be applied to any seismic phase that is observed over a reasonably large distance range. Here, we present the results of applying the algorithm to the long-period channels of all data archived at the IRIS DMC from 1976 to 2005 for the seismic phases S and P (from 23 ◦ to 100 ◦ ) and SS and PP (from 50 ◦ to 170 ◦ ). The resulting large data sets are inverted along with existing surface wave and updated differential traveltime measurements for new mantle models of S and P velocity. The resolution of the new model is enhanced, particularly, in the mid-mantle where SS and PP turn. We find that slow anomalies in the central Pacific and Africa extend from the core‐mantle boundary to the upper mantle, but their direct connection to surface hotspots is beyond our resolution. Furthermore, we find that fast anomalies that are likely associated with subducting slabs disappear between 1700 and 2500 km, and thus are not continuous features from the upper to lower mantle despite our extensive coverage and high resolution of the mid-mantle.

Source inversion of W phase: speeding up seismic tsunami warning

Abstract: W phase is a long period phase arriving before S wave. It can be interpreted as superposition of the fundamental, first, second and third overtones of spheroidal modes or Rayleigh waves and has a group velocity from 4.5 to 9 km s^−1 over a period range of 100–1000 s. The amplitude of long period waves better represents the tsunami potential of an earthquake. Because of the fast group velocity of W phase, most of W phase energy is contained within a short time window after the arrival of the P wave. At a distance of 50°, W phase energy is contained within 23 min after the origin time which is the distinct advantage of using W phase for rapid tsunami warning purposes. We use a time domain deconvolution method to extract W phases from the broad-band records of global seismic networks. The bandwidth of W phase is approximately from 0.001 to 0.01 Hz, and we bandpass filter the data from 0.001 to 0.005 Hz in most cases. Having extracted W phase from the vertical component records, we perform a linear inversion using a point source to determine Mw and the source mechanism for several large earthquakes including the 2004 Sumatra–Andaman earthquake, the 2005 Nias earthquake, the 2006 Kuril Is. earthquake and the 2007 Sumatra earthquake. W phase inversion yields reliable solutions and holds promise of the use of W phase for rapid assessment of tsunami potential.

Theoretical background for continental‐ and global‐scale full‐waveform inversion in the time–frequency domain

Abstract: SUMMARY We propose a new approach to full seismic waveform inversion on continental and global scales. This is based on the time–frequency transform of both data and synthetic seismograms with the use of time- and frequency-dependent phase and envelope misfits. These misfits allow us to provide a complete quantification of the differences between data and synthetics while separating phase and amplitude information. The result is an efficient exploitation of waveform information that is robust and quasi-linearly related to Earth's structure. Thus, the phase and envelope misfits are usable for continental- and global-scale tomography, that is, in a scenario where the seismic wavefield is spatially undersampled and where a 3-D reference model is usually unavailable. Body waves, surface waves and interfering phases are naturally included in the analysis. We discuss and illustrate technical details of phase measurements such as the treatment of phase jumps and instability in the case of small amplitudes. The Frechet kernels for phase and envelope misfits can be expressed in terms of their corresponding adjoint wavefields and the forward wavefield. The adjoint wavefields are uniquely determined by their respective adjoint-source time functions. We derive the adjoint-source time functions for phase and envelope misfits. The adjoint sources can be expressed as inverse time–frequency transforms of a weighted phase difference or a weighted envelope difference. In a comparative study, we establish connections between the phase and envelope misfits and the following widely used measures of seismic waveform differences: (1) cross-correlation time-shifts; (2) relative rms amplitude differences; (3) generalized seismological data functionals and (4) the L2 distance between data and synthetics used in time-domain full-waveform inversion. We illustrate the computation of Frechet kernels for phase and envelope misfits with data from an event in the West Irian region of Indonesia, recorded on the Australian continent. The synthetic seismograms are computed for a heterogeneous 3-D velocity model of the Australian upper mantle, with a spectral-element method. The examples include P body waves, Rayleigh waves and S waves, interfering with higher-mode surface waves. All the kernels differ from the more familar kernels for cross-correlation time-shifts or relative rms amplitude differences. The differences arise from interference effects, 3-D Earth's structure and waveform dissimilarities that are due to waveform dispersion in the heterogeneous Earth.

New advances in three‐dimensional controlled‐source electromagnetic inversion

Abstract: SUMMARY New techniques for improving both the computational and imaging performance of the three-dimensional (3-D) electromagnetic inverse problem are presented. A non-linear conjugate gradient algorithm is the framework of the inversion scheme. Full wave equation modelling for controlled sources is utilized for data simulation along with an efficient gradient computation approach for the model update. Improving the modelling efficiency of the 3-D finite difference (FD) method involves the separation of the potentially large modelling mesh, defining the set of model parameters, from the computational FD meshes used for field simulation. Grid spacings and thus overall grid sizes can be reduced and optimized according to source frequencies and source–receiver offsets of a given input data set. Further computational efficiency is obtained by combining different levels of parallelization. While the parallel scheme allows for an arbitrarily large number of parallel tasks, the relative amount of message passing is kept constant. Image enhancement is achieved by model parameter transformation functions, which enforce bounded conductivity parameters and thus prevent parameter overshoots. Further, a remedy for treating distorted data within the inversion process is presented. Data distortions simulated here include positioning errors and a highly conductive overburden, hiding the desired target signal. The methods are demonstrated using both synthetic and field data.

Global upper-mantle tomography with the automated multimode inversion of surface and S-wave forms

Abstract: SUMMARY We apply the Automated Multimode Inversion of surface and S-wave forms to a large global data set, verify the accuracy of the method and assumptions behind it, and compute an S vvelocity model of the upper mantle (crust‐660 km). The model is constrained with ∼51 000 seismograms recorded at 368 permanent and temporary broadband seismic stations. Structure of the mantle and crust is constrained by waveform information both from the fundamentalmode Rayleigh waves (periods from 20 to 400 s) and from S and multiple S waves (higher modes). In order to enhance the validity of the path-average approximation, we implement the automated inversion of surface- and S-wave forms with a three-dimensional (3-D) reference model. Linear equations obtained from the processing of all the seismograms of the data set are inverted for seismic velocity variations also relative to a 3-D reference, in this study composed of a 3-D model of the crust and a one-dimensional (1-D), global-average depth profile in the mantle below. Waveform information is related to shear- and compressional-velocity structure within approximate waveform sensitivity areas. We use two global triangular grids of knots with approximately equal interknot spacing within each: a finely spaced grid for integration over sensitivity areas and a rougher-spaced one for the model parametrization. For the tomographic inversion we use LSQR with horizontal and vertical smoothing and norm damping. We invert for isotropic variations in S- and P-wave velocities but also allow for S-wave azimuthal anisotropy—in order to minimize errors due to possible mapping of anisotropy into isotropic heterogeneity. The lateral resolution of the resulting isotropic upper-mantle images is a few hundred kilometres, varying with data sampling. We validate the imaging technique with a ‘spectral-element’ resolution test: inverting a published global synthetic data set computed with the spectral-element method using a laterally heterogeneous mantle model we are able to reconstruct the synthetic model accurately. This test confirms both the accuracy of the implementation of the method and the validity of the JWKB and path-average approximations as applied in it. Reviewing the tomographic model, we observe that low-S v-velocity anomalies beneath mid-ocean ridges and backarc basins extend down to ∼100 km depth only, shallower than according to some previous tomographic models; this presents a close match to published estimates of primary melt production depth ranges there. In the seismic lithosphere beneath cratons, unambiguous high velocity anomalies extend to ∼200 km. Pronounced low-velocity zones beneath cratonic lithosphere are rare; where present (South America; Tanzania) they are neighboured by volcanic areas near cratonic boundaries. The images of these low-velocity zones may indicate hot material—possibly of mantle-plume origin—trapped or spreading beneath the thick cratonic lithosphere.

Cracks, pulses and macroscopic asymmetry of dynamic rupture on a bimaterial interface with velocity-weakening friction

Abstract: We study in-plane ruptures on a bimaterial fault governed by a velocity-weakening friction with a regularized normal stress response. Numerical simulations and analytical estimates provide characterization of the ranges of velocity-weakening scales, nucleation lengths and background stresses for which ruptures behave as cracks or pulses, decaying or sustained, bilateral or unilateral. With strongly velocity-weakening friction, ruptures occur under a wide range of conditions as large-scale pulses with a preferred propagation direction, that of slip of the more compliant material. Such ruptures have macroscopic asymmetry manifested by significantly larger seismic potency and propagation distance in the preferred direction, and clearly quantified by the directivity ratio derived from the second order moments of the spatio-temporal distribution of slip rate. The macroscopic rupture asymmetry of the large-scale pulses stems from the difference in the criticality conditions for self-sustained propagation in each rupture direction, induced by the asymmetric normal stress changes operating in bimaterial interfaces. In contrast, crack-like ruptures show macroscopic asymmetry under restrictive conditions. The discussed mechanism is robust with respect to regularization parameters, ranges of stress heterogeneities and a proxy for off-fault yielding and should operate similarly for crustal-scale rupture pulses even in the absence of velocity-weakening. Small-scale pulses, driven by the bimaterial normal stress reduction at the scale of the process zone, can detach from the rupture front of the large-scale pulses that propagate in the preferred direction. However, their occurrence depends on the relaxation scale in the regularization of the normal stress response and their development can be hindered by off-fault yielding.

Lithospheric structure and deep earthquakes beneath India, the Himalaya and southern Tibet

Abstract: SUMMARY This paper is concerned with the implications of earthquake depth distributions in the Himalayan‐Tibetan collision zone for the general understanding of lithosphere rheology. In particular, recent studies have argued that microearthquakes in the uppermost mantle beneath Nepal and some earthquakes at 80‐90 km depth, close to the Moho in SE and NW Tibet, reinforce the conventional view of the last 25 yr that the continental lithosphere is well represented by strong seismogenic layers in the upper crust and uppermost mantle, separated by a weak and aseismic lower crust. That view was recently challenged by an alternative one suggesting that the continental lithosphere contained a single strong seismogenic layer, that was either the upper crust or the whole crust, but did not involve the mantle. We re-examine the seismic structure and seismicity of the Himalayan‐Tibetan collision zone, recalculating earthquake depths in velocity structures that are consistent with seismic receiver functions and surface wave dispersion studies, and calculating a geotherm for the Indian Shield consistent with kimberlite nodule geochemistry. Earthquakes occur throughout the crustal thickness of the Indian Shield, where the lower crust is thought to consist of dry granulite, responsible for its seismogenic behaviour and strength as manifested by its relatively large effective elastic thickness. The crust of the Indian Shield is thin (∼35 km) for an Archean shield, and this, in turn, leads to a steady-state Moho temperature that could be as low as ∼500 ◦ C. When this shield is thrust beneath the Himalaya in Nepal, the relatively low mantle temperature, together with the high strain rates associated with it adopting a ‘ramp-and-flat’ geometry, may be responsible for the mantle microearthquakes that accompany other earthquakes in the lower crust. Further north, the upper crust of India south of the Indus Suture Zone has been removed, the uppermost lower crust of India has heated up, and seismicity is restricted to a few earthquakes very close to the Moho at 80‐90 km, where errors in Moho and earthquake depth determinations make it unclear whether these events are in the crust or mantle. A similar situation exists in NW Tibet beneath the Kunlun, where earthquakes at 80‐90 km depth occur very close to the Moho. Both places are about 400 km north of the Himalayan front, and we suspect both represent the minimum distance India has underthrust Tibet, so that India underlies most of the SE and nearly all of the NW Tibetan plateau. The distribution of earthquake depths throughout the region is consistent with a generic global view of seismicity in which earthquakes occur in (1) ‘wet’ upper crustal material to a temperature of ∼350 ◦ C, or (2) higher temperatures in dry granulite-facies lower crust or (3) mantle that is colder than ∼600 ◦ C.

The magnetic structure of convection-driven numerical dynamos

Abstract: The generation of a magnetic field in numerical simulations of the geodynamo is an intrinsically 3-D and time-dependent phenomenon. The concept of magnetic field lines and the frozen-flux approximation can provide insight into such systems, but a suitable visualization method is required. This paper presents results obtained using the Dynamical Magnetic Field line Imaging (DMFI) technique, which is a representation of magnetic field lines accounting for their local magnetic energy, together with an algorithm for the time evolution of their anchor points. The DMFI illustrations are consistent with previously published dynamo mechanisms, and allow further investigation of spatially and temporally complex systems. We highlight three types of magnetic structures: (i) magnetic cyclones and (ii) magnetic anticyclones are expelled by, but corotate with axial flow vortices; (iii) magnetic upwellings are amplified by stretching and advection within flow upwellings, and show structural similarity with helical plumes found in rotating hydrodynamic experiments. While magnetic anticyclones are responsible for the regeneration of a stable axial dipole, herewe showthat excursions and reversals of the dipole axis are caused by the emergence of magnetic upwellings, which amplify and transport a generally multipolar magnetic field from the inner to the outer boundary of the models. Geomagnetic observations suggest the presence of magnetic structures similar to those found in our models; thus,we discuss howour results may pertain to Earth's core dynamo processes. In order to make DMFI a standard tool for numerical dynamo studies, a public software package is available upon request to the authors (supplementary material is available at: http://www.ipgp.jussieu. fr/∼aubert/DMFI.html).

The design of an optimal filter for monthly GRACE gravity models

Abstract: SUMMARY Most applications of the publicly released Gravity Recovery and Climate Experiment monthly gravity field models require the application of a spatial filter to help suppressing noise and other systematic errors present in the data The most common approach makes use of a simple Gaussian averaging process, which is often combined with a ‘destriping’ technique in which coefficient correlations within a given degree are removed As brute force methods, neither of these techniques takes into consideration the statistical information from the gravity solution itself and, while they perform well overall, they can often end up removing more signal than necessary Other optimal filters have been proposed in the literature; however, none have attempted to make full use of all information available from the monthly solutions By examining the underlying principles of filter design, a filter has been developed that incorporates the noise and full signal variance–covariance matrix to tailor the filter to the error characteristics of a particular monthly solution The filter is both anisotropic and nonsymmetric, meaning it can accommodate noise of an arbitrary shape, such as the characteristic stripes The filter minimizes the mean-square error and, in this sense, can be considered as the most optimal filter possible Through both simulated and real data scenarios, this improved filter will be shown to preserve the highest amount of gravity signal when compared to other standard techniques, while simultaneously minimizing leakage effects and producing smooth solutions in areas of low signal

Interseismic Plate coupling and strain partitioning in the Northeastern Caribbean

Abstract: SUMMARY The northeastern Caribbean provides a natural laboratory to investigate strain partitioning, its causes and its consequences on the stress regime and tectonic evolution of a subduction plate boundary. Here, we use GPS and earthquake slip vector data to produce a present-day kinematic model that accounts for secular block rotation and elastic strain accumulation, with variable interplate coupling, on active faults. We confirm that the oblique convergence between Caribbean and North America in Hispaniola is partitioned between plate boundary parallel motion on the Septentrional and Enriquillo faults in the overriding plate and plateboundary normal motion at the plate interface on the Northern Hispaniola Fault. To the east, the Caribbean/North America plate motion is accommodated by oblique slip on the faults bounding the Puerto Rico block to the north (Puerto Rico subduction) and to the south (Muertos thrust), with no evidence for partitioning. The spatial correlation between interplate coupling, strain partitioning and the subduction of buoyant oceanic asperities suggests that the latter enhance the transfer of interplate shear stresses to the overriding plate, facilitating strike-slip faulting in the overriding plate. The model slip rate deficit, together with the dates of large historical earthquakes, indicates the potential for a large (M w7.5 or greater) earthquake on the Septentrional fault in the Dominican Republic. Similarly, the Enriquillo fault in Haiti is currently capable of a M w7.2 earthquake if the entire elastic strain accumulated since the last major earthquake was released in a single event today. The model results show that the Puerto Rico/Lesser Antilles subduction thrust is only partially coupled, meaning that the plate interface is accumulating elastic strain at rates slower than the total plate motion. This does not preclude the existence of isolated locked patches accumulating elastic strain to be released in future earthquakes, but whose location and geometry are not resolvable with the present data distribution. Slip deficit on faults from this study are used in a companion paper to calculate interseismic stress loading and, together with stress changes due to historical earthquakes, derive the recent stress evolution in the NE Caribbean.

Far-field tsunami hazard from mega-thrust earthquakes in the Indian Ocean

Abstract: SUMMARY We evaluate far-field tsunami hazard in the Indian Ocean Basin based on hydrodynamic simulations of ten case studies of possible mega earthquakes at the major seismic zones surrounding the basin. They represent worst-case scenarios of seismic rupture along the full extent of seismogenic faults having supported large earthquakes in the historical record. In a series of numerical experiments in which the source parameters of the 2004 Sumatra tsunami are allowed to vary one by one, while keeping the seismic moment and the fault orientation unchanged, we document that the main patterns of far-field tsunami amplitudes are remarkably robust with respect to nominal variations in such parameters as hypocentral depth, exact centroid location, and slip distribution on the fault plane. These results validate the concept of modelling case scenarios of potential future earthquakes whose source is by definition imprecise. We consider seismic sources located at the extremities of the 2004 Sumatra‐Andaman rupture, namely along the southern coast of Sumatra and in the Andaman‐Myanmar province; along the Makran coast of Pakistan and Iran; and also along the southern coast of Java, where the possibility of a large interplate thrust earthquake cannot be entirely dismissed. The results of our hydrodynamic simulations indicate that the distribution of maximum amplitudes in the Indian Ocean Basin is primarily controlled by the classical effect of source directivity, and additionally by refraction and focusing along bathymetric features. As a result, many provinces in the basin could be threatened by higher tsunami amplitudes than in 2004. This pattern is particularly important along the coast of East Africa, from Somalia to and including South Africa, in Madagascar and the Mascarene Islands, especially under a South Sumatra scenario involving an earthquake comparable to, or even possibly larger than, the 1833 event, whose epicentral area is widely believed to be under enhanced seismic risk as a result of stress transfer from the 2004 and 2005 ruptures to the northwest, possibly even in the wake of the 2007 Bengkulu earthquakes.

The interior penalty discontinuous Galerkin method for elastic wave propagation: grid dispersion

Abstract: SUMMARY Recently, there has been an increased interest in applying the discontinuous Galerkin method (DGM) to wave propagation. In this work, we investigate the applicability of the interior penalty DGM to elastic wave propagation by analysing it’s grid dispersion properties, with particular attention to the effect that different basis functions have on the numerical dispersion. We consider different types of basis functions that naturally yield a diagonal mass matrix. This is relevant to seismology because a diagonal mass matrix is tantamount to an explicit and efficient time marching scheme. We find that the Legendre basis functions that are traditionally used in the DGM introduce numerical dispersion and anisotropy. Furthermore, we find that using Lagrange basis functions along with the Gauss nodes has attractive advantages for numerical wave propagation.

Three‐dimensional modelling of magnetotelluric data from the Rotokawa geothermal field, Taupo Volcanic Zone, New Zealand

Abstract: SUMMARY The resistivity structure of the Rotokawa geothermal system in New Zealand’s Taupo Volcanic Zone has been determined by 3-D modelling of data from a closely spaced (64 measurement sites) magnetotelluric (MT) survey. 3-D conductivity models were constructed using trial and error forward modelling of the phase-tensor data and 3-D inverse modelling of the impedance tensor data. Both the forward and the inverse resistivity models show good consistency. The most interesting feature of these models is a resistive (∼100 � m) zone within the otherwise conductive material of the geothermal system. This zone coincides with the high temperature (300‐335 ◦ C) core of the geothermal system in which seismicity induced by fluid injection occurs and may mark the zone of fracture permeability that is feeding high temperature fluid into the geothermal system from deeper levels.

Seismic interferometry, surface waves and source distribution

Abstract: SUMMARY Seismic interferometry can be used to estimate interreceiver surface wave signals by cross-correlation of signals recorded at each receiver. The quality of the estimated surface waves is controlled by the distribution of sources exciting the cross-correlated wavefields, and it is commonly thought that only sources at or near the surface are required to generate accurate estimates. We study the role of source distribution in surface wave interferometry for both surface and subsurface sources using surface wave Green's functions for laterally homogeneous media. We solve the interferometric integral using a Rayleigh wave orthogonality relationship combined with a stationary phase approach. Contrary to popular opinion we find that sources at depth do indeed play a role in the recovery of surface waves by interferometry. We find that interferometry performs well when surface sources are distributed homogeneously at the surface of the Earth. However, when this homogeneous distribution is not available amplitude errors are introduced, and when multiple modes are present strong spurious events appear and higher mode surface waves may not be correctly estimated. In order to recover higher mode surface waves we propose an additional step in the processing of surface wave data for seismic interferometry: by separating modes and applying interferometry to each mode individually it is possible to recover the interreceiver surface wave modes, without the artefacts introduced by limited source coverage.

A model of plate kinematics in Gondwana breakup

Abstract: SUMMARY An accurate model of relative plate motions in Gondwana breakup is based on visual fitting of seafloor isochrons and fracture zones (FZ) from the Riiser-Larsen Sea and Mozambique Basin. Used predictively, the model precisely locates kinematic markers in the West Somali Basin, which allows the conclusion that the spreading centres in the West Somali and Mozambique basins and the Riiser-Larsen Sea formed parts of the boundary between the same two plates. The locations of FZ and less well-defined isochrons from neighbouring regions are also consistent with their formation on other lengths of this same boundary and with its relocation from the West Somali Basin and northern Natal Valley to the West Enderby Basin and Lazarev Sea during chron M10n. Small independently moving plates thus played no role in the breakup of this core part of Gondwana. In an inversion procedure, the data from these areas yield more precise finite rotations that describe the history of the two plates’ separation. Breakup is most simply interpreted to have occurred in coincidence with Karoo volcanism, and a reconstruction based on the rotations shows the Lebombo and Mateke-Sabi monoclines and the Mozambique and Astrid ridges as two sets of conjugate volcanic margins. Madagascar’s pre-drift position can be used as a constraint to reassess the positions of India and Sri Lanka in the supercontinent.

Quasi-geostrophic flows responsible for the secular variation of the Earth's magnetic field

Abstract: SUMMARY We present core flows constructed from high resolution secular variation (SV) models for the epochs 2001, 2002.5 and 2004, assuming that these flows are quasi-geostrophic in the core interior and that they do not cross the axial cylindrical surface tangent to the inner core. A large jet encircling the inner core and carrying a significant part of the core angular momentum and axial vortices of ∼700 km diameter mainly clustering around the cylinder tangent to the solid inner core, are inferred from geomagnetic SV. New regularizations are suggested from dynamic considerations. It is found that medium and small-scale velocity fields contribute significantly to the large-scale SV. Accordingly, final models of core flows are calculated after an iterative process, whereby the magnetic field variation produced by small-scale stochastic magnetic fields and medium to small-scale computed velocity fields are incorporated into the inversion itself as modelling errors. This study represents a significant step in an effort to join geomagnetic observations and the fluid core dynamics on short timescales.

Testing alarm‐based earthquake predictions

Abstract: SUMMARY Motivated by a recent resurgence in earthquake predictability research, we present a method for testing alarm-based earthquake predictions. The testing method is based on the Molchan diagram—a plot of miss rate and fraction of space–time occupied by alarm—and is applicable to a wide class of predictions, including probabilistic earthquake forecasts varying in space, time, and magnitude. A single alarm can be simply tested using the cumulative binomial distribution. Here we consider the more interesting case of a function from which a continuum of well-ordered alarms can be derived. For such an ‘alarm function’ we construct a cumulative performance measure, the area skill score, based on the normalized area above a trajectory on the Molchan diagram. A score of unity indicates perfect skill, a score of zero indicates perfect non-skill, and the expected score for a random alarm function is 1/2. The area skill score quantifies the performance of an arbitrary alarm function relative to a reference model. To illustrate the testing method, we consider the 10-yr experiment by J. Rundle and others to predict M≥ 5 earthquakes in California. We test forecasts from three models: relative intensity (RI), a simple spatial clustering model constructed using only smoothed historical seismicity; pattern informatics (PI), a model that aims to capture seismicity dynamics by pattern recognition; and the U. S. Geological Survey National Seismic Hazard Map (NSHM), a model that comprises smoothed historical seismicity, zones of ‘background’ seismicity, and explicit fault information. Results show that neither PI nor NSHM provide significant performance gain relative to the RI reference model. We suggest that our testing method can be used to evaluate future experiments in the Collaboratory for the Study of Earthquake Predictability and to iteratively improve reference models for earthquake prediction hypothesis testing.

Capturing magma intrusion and faulting processes during continental rupture: seismicity of the Dabbahu (Afar) rift

Abstract: Continental rupture models emphasize the role of faults in extensional strain accommodation; extension by dyke intrusion is commonly overlooked. A major rifting episode that began in 2005 September in the Afar depression of Ethiopia provides an opportunity to examine strain accommodation in a zone of incipient plate rupture. Earthquakes recorded on a temporary seismic array (2005 October to 2006 April), direct observation of fault patterns and geodetic data document ongoing strain and continued dyke intrusion along the ?60-km long Dabbahu rift segment defined in earlier remote sensing studies. Epicentral locations lie along a ?3 km wide, ?50 km long swath that curves into the SE flank of Dabbahu volcano; a second strand continues to the north toward Gab'ho volcano. Considering the ?8 m of opening in the September crisis, we interpret the depth distribution of microseismicity as the dyke intrusion zone; the dykes rise from ?10 km to the near-surface along the ?60-km long length of the tectono-magmatic segment. Focal mechanisms indicate slip along NNW-striking normal faults, perpendicular to the Arabia–Nubia plate opening vector. The seismicity, InSAR, continuous GPS and structural patterns all suggest that magma injection from lower or subcrustal magma reservoirs continued at least 3 months after the main episode. Persistent earthquake swarms at two sites on Dabbahu volcano coincide with areas of deformation identified in the InSAR data: (1) an elliptical, northwestward-dipping zone of seismicity and subsidence interpreted as a magma conduit, and (2) a more diffuse, 8-km radius zone of shallow seismicity (<2 km) above a shadow zone, interpreted as a magma chamber between 2.5 and 6 km subsurface. InSAR and continuous GPS data show uplift above a shallow source in zone (2) and uplift above the largely aseismic Gab'ho volcano. The patterns of seismicity provide a 3-D perspective of magma feeding systems maintaining the along-axis segmentation of this incipient seafloor spreading segment.

Finite-frequency sensitivity kernels for global seismic wave propagation based upon adjoint methods

Abstract: SUMMARY We determine adjoint equations and Frechet kernels for global seismic wave propagation based upon a Lagrange multiplier method. We start from the equations of motion for a ro- tating, self-gravitating earth model initially in hydrostatic equilibrium, and derive the cor- responding adjoint equations that involve motions on an earth model that rotates in the opposite direction. Variations in the misfit function χ then may be expressed as δχ = � V Km δ ln md 3 r + � � Kd δ ln dd 2 r + � � FS Kd ·∇ � δ ln dd 2 r, where δ ln m = δm/m denotes relative model perturbations in the volume V , δ ln d denotes relative topographic variations on solid-solid or fluid-solid boundaries � , and ∇ � δ ln d denotes surface gradients in relative topographic variations on fluid-solid boundariesFS. The 3-D Frechet kernel K m determines the sensitivity to model perturbations δ ln m, and the 2-D kernels K d and Kd determine the sensitivity to topographic variations δ ln d. We demonstrate also how anelasticity may be in- corporated within the framework of adjoint methods. Finite-frequency sensitivity kernels are calculated by simultaneously computing the adjoint wavefield forward in time and reconstruct- ing the regular wavefield backward in time. Both the forward and adjoint simulations are based upon a spectral-element method. We apply the adjoint technique to generate finite-frequency traveltime kernels for global seismic phases (P, P diff, PKP, S, SKS, depth phases, surface- reflected phases, surface waves, etc.) in both 1-D and 3-D earth models. For 1-D models these adjoint-generated kernels generally agree well with results obtained from ray-based methods. However, adjoint methods do not have the same theoretical limitations as ray-based methods, and can produce sensitivity kernels for any given phase in any 3-D earth model. The Frechet kernels presented in this paper illustrate the sensitivity of seismic observations to structural parameters and topography on internal discontinuities. These kernels form the basis of future 3-D tomographic inversions.

GRIMM: the GFZ Reference Internal Magnetic Model based on vector satellite and observatory data

Abstract: SUMMARY In this paper the new GFZ Reference Internal Magnetic Model (GRIMM) is presented. The model has been derived from nearly 6 yr of CHAMP satellite data and 5 yr of observatory hourly means. At high latitudes, full vector satellite data are used at all local times which allows a separation between, on one hand, the fields generated by ionosphere and field aligned currents, and, on the other hand, the fields generated in the Earth's core and lithosphere. This selection technique leads to a data set without gaps during the polar summers resulting in a core field model that has an unprecedented time resolution. The modelled static core field, secular variation and lithospheric field are all in good agreement with previously published magnetic field models. Order five B-splines are used to model the variation in time of the core field. The energy in the secular acceleration has, therefore, a smooth behaviour in time and increases continuously from 2003.5. Mapping the field acceleration from 2001.5 to 2005.5 reveals its rapid and complex evolution over this time period at the Earth's surface. Due to the applied regularization technique, the acceleration energy in spherical harmonics 6–11 is significantly larger than for other models and we show that such a spectrum is acceptable.

Spectral-element simulations of wave propagation in porous media

Abstract: We present a derivation of the equations describing wave propagation in porous media based upon an averaging technique which accommodates the transition from the microscopic to the macroscopic scale. We demonstrate that the governing macroscopic equations determined by Biot remain valid for media with gradients in porosity. In such media, the well-known expression for the change in porosity, or the change in the fluid content of the pores, acquires two extra terms involving the porosity gradient. One fundamental result of Biot's theory is the prediction of a second compressional wave, often referred to as 'type II' or 'Biot's slow compressional wave', in addition to the classical fast compressional and shear waves. We present a numerical implementation of the Biot equations for 2-D problems based upon the spectral-element method (SEM) that clearly illustrates the existence of these three types of waves as well as their interactions at discontinuities. As in the elastic and acoustic cases, poroelastic wave propagation based upon the SEM involves a diagonal mass matrix, which leads to explicit time integration schemes that are well suited to simulations on parallel computers. Effects associated with physical dispersion and attenuation and frequency-dependent viscous resistance are accommodated based upon a memory variable approach. We perform various benchmarks involving poroelastic wave propagation and acoustic–poroelastic and poroelastic–poroelastic discontinuities, and we discuss the boundary conditions used to deal with these discontinuities based upon domain decomposition. We show potential applications of the method related to wave propagation in compacted sediments, as one encounters in the petroleum industry, and to detect the seismic signature of buried landmines and unexploded ordnance.

Effect of 3-D viscoelastic structure on post-seismic relaxation from the 2004 M= 9.2 Sumatra earthquake

Abstract: SUMMARY The 2004 M = 9.2 Sumatra‐Andaman earthquake profoundly altered the state of stress in a large volume surrounding the ∼1400 km long rupture. Induced mantle flow fields and coupled surface deformation are sensitive to the 3-D rheology structure. To predict the post-seismic motions from this earthquake, relaxation of a 3-D spherical viscoelastic earth model is simulated using the theory of coupled normal modes. The quasi-static deformation basis set and solution on the 3-D model is constructed using: a spherically stratified viscoelastic earth model with a linear stress‐strain relation; an aspherical perturbation in viscoelastic structure; a ‘static’ mode basis set consisting of Earth’s spheroidal and toroidal free oscillations; a “viscoelastic” mode basis set; and interaction kernels that describe the coupling among viscoelastic and static modes. Application to the 2004 Sumatra‐Andaman earthquake illustrates the profound modification of the post-seismic flow field at depth by a slab structure and similarly large effects on the near-field post-seismic deformation field at Earth’s surface. Comparison with postseismic GPS observations illustrates the extent to which viscoelastic relaxation contributes to the regional post-seismic deformation.

Stable inversions for complete moment tensors

Abstract: SUMMARY The seismic moment tensors for certain types of sources, such as volcanic earthquakes and nuclear explosions are expected to contain an isotropic component. Some earlier efforts to calculate the isotropic component of these sources are flawed due to an error in the method of Jost & Herrmann. We corrected the method after Herrmann & Hutchensen and found great improvement in the recovery of non-double-couple moment tensors that include an isotropic component. Tests with synthetic data demonstrate the stability of the corrected linear inversion method, and we recalculate the moment tensor solutions reported in Dreger et al. for Long Valley caldera events and Dreger & Woods for Nevada Test Site nuclear explosions. We confirm the findings of Dreger et al. that the Long Valley volcanic sources contain large statistically significant isotropic components. The nuclear explosions have strikingly anomalous source mechanisms, which contain very large isotropic components, making it evident that these events are not tectonic in origin. This indicates that moment tensor inversions could be an important tool for nuclear monitoring.

The African upper mantle and its relationship to tectonics and surface geology

Abstract: This paper focuses on the upper-mantle velocity structure of the African continent and its relationship to the surface geology. The distribution of seismographs and earthquakes providing seismograms for this study results in good fundamental and higher mode path coverage by a large number of relatively short propagation paths, allowing us to image the SV-wave speed structure, with a horizontal resolution of several hundred kilometres and a vertical resolution of similar to 50 km, to a depth of about 400 km. The difference in mantle structure between the Archean and Pan-African terranes is apparent in our African upper-mantle shear wave model. High-velocity (4-7 per cent) roots exist beneath the cratons. Below the West African, Congo and Tanzanian Cratons, these extend to 225-250 km depth, but beneath the Kalahari Craton, the high wave speed root extends to only similar to 170 km. With the exception of the Damara Belt that separates the Congo and Kalahari Cratons, any high-speed upper-mantle lid below the Pan-African terranes is too thin to be resolved by our long-period surface wave technique. The Damara Belt is underlain by higher wave speeds, similar to those observed beneath the Kalahari Craton. Extremely low SV-wave speeds occur to the bottom of our model beneath the Afar region. The temperature of the African upper mantle is determined from the SV-wave speed model. Large temperature variations occur at 125 km depth with low temperatures beneath west Africa and all of southern Africa and warm mantle beneath the Pan-African terrane of northern Africa. At 175 km depth, cool upper mantle occurs below the West African, Congo, Tanzanian and Kalahari Cratons and anomalously warm mantle occurs below a zone in northcentral Africa and beneath the region surrounding the Red Sea. All of the African volcanic centres are located above regions of warm upper mantle. The temperature profiles were fit to a geotherm to determine the thickness of the African lithosphere. Thick lithosphere exists beneath all of the cratonic areas; independent evidence for this thick lithosphere comes from the locations of diamondiferous kimberlites. Almost all diamond locations occur where the lithosphere is 175-200 km thick, but they are largely absent from the regions of the thickest lithosphere. The lithosphere is thin beneath the Pan-African terranes of northern Africa but appears to be thicker beneath the Pan-African Damara Belt in southern Africa.

Three-dimensional gravity inversion for Moho depth at rifted continental margins incorporating a lithosphere thermal gravity anomaly correction

Abstract: SUMMARY This paper describes a method for determining Moho depth, lithosphere thinning factor (γ = 1 − 1/β) and the location of the ocean‐continent transition at rifted continental margins using 3-D gravity inversion which includes a correction for the large negative lithosphere thermal gravity anomaly within continental margin lithosphere. The lateral density changes caused by the elevated geotherm in thinned continental margin and adjacent ocean basin lithosphere produce a significant lithosphere thermal gravity anomaly which may be in excess of −100 mGal, and for which a correction must be made in order to determine Moho depth accurately from gravity inversion. We describe a method of iteratively calculating the lithosphere thermal gravity anomaly using a lithosphere thermal model to give the present-day temperature field from which we calculate the lithosphere thermal density and gravity anomalies. For continental margin lithosphere, the lithosphere thermal perturbation is calculated from the lithosphere thinning factor (γ = 1 − 1/β) obtained from crustal thinning determined by gravity inversion and breakup age for thermal re-equilibration time. For oceanic lithosphere, the lithosphere thermal model used to predict the lithosphere thermal gravity anomaly may be conditioned using ocean isochrons from plate reconstruction models to provide the age and location of oceanic lithosphere. A correction is made for crustal melt addition due to decompression melting during continental breakup and seafloor spreading. We investigate the sensitivity of the lithosphere thermal gravity anomaly and the predicted Moho depth from gravity inversion at continental rifted margins to the methods used to calculate and condition the lithosphere thermal model using both synthetic models and examples from the North Atlantic.