Showing papers in "Geophysical Journal International in 2003"

GPS constraints on Africa (Nubia) and Arabia plate motions

Abstract: SUMMARY We use continuously recording GPS (CGPS) and survey-mode GPS (SGPS) observations to determine Euler vectors for relative motion of the African (Nubian), Arabian and Eurasian plates. We present a well-constrained Eurasia‐Nubia Euler vector derived from 23 IGS sites in Europe and four CGPS and three SGPS sites on the Nubian Plate (−0.95 ± 4.8 ◦ N, −21.8 ± 4.3 ◦ E, 0.06 ± 0.005 ◦ Myr −1 ). We see no significant (> 1m m yr −1 ) internal deformation of the Nubian Plate. The GPS Nubian‐Eurasian Euler vector differs significantly from NUVEL-1A (21.0 ± 4.2 ◦ N, −20.6 ± 0.6 ◦ E, 0.12 ± 0.015 ◦ Myr −1 ), implying more westward motion of Africa relative to Eurasia and slower convergence in the eastern Mediterranean. The Arabia‐ Eurasia and Arabia‐Nubia GPS Euler vectors are less well determined, based on only one CGPS and three SGPS sites on the Arabian Plate. The preliminary Arabia‐Eurasia and Arabia‐ Nubia Euler vectors are 27.4 ± 1.0 ◦ N, 18.4 ± 2.5 ◦ E, 0.40 ± 0.04 ◦ Myr −1 , and 30.5 ± 1.0 ◦ N, 25.7 ± 2.3 ◦ E, 0.37 ± 0.04 ◦ Myr −1 , respectively. The GPS Arabia‐Nubia Euler vector differs significantly from NUVEL-1A (24.1 ± 1.7 ◦ N, 24.0 ± 3.5 ◦ E, 0.40 ± 0.05 ◦ Myr −1 ), but is statistically consistent at the 95 per cent confidence level with the revised Euler vector reported by Chu & Gordon based on a re-evaluation of magnetic anomalies in the Red Sea (31.5 ± 1.2 ◦ N, 23.0 ± 2.7 ◦ E, 0.40 ± 0.05 ◦ Myr −1 ). The motion implied in the Gulf of Aqaba and on the Dead Sea fault (DSF) by the new GPS Nubia‐Arabia Euler vector (i.e. ignoring possible Sinai block motion and possible internal plate deformation) grades from pure left lateral strike-slip in the Gulf and on the southern DSF with increasing compression on the central and northern DSF with relative motion increasing from 5.6 to 7.5 mm yr −1 (± 1m m yr −1 ) from south to north. Along the northern DSF (i.e. north of the Lebanon restraining bend) motion is partitioned between 6 ± 1m m yr −1 left-lateral motion parallel to the fault trace and 4 ± 1m m yr −1 faultnormal compression. Relative motions on other plate boundaries (including the Anatolian and Aegean microplates) derived from the GPS Euler vectors agree qualitatively with the sense of motion indicated by focal mechanisms for large crustal earthquakes (M > 6). Where data are available on fault-slip rates on plate bounding faults (North Anatolian fault, East Anatolian fault, Dead Sea fault, Red Sea rift), they are generally lower than, but not significantly different from, the full plate motion estimates suggesting that the majority of relative plate motion is accommodated on these structures.

An integrated global model of present-day plate motions and plate boundary deformation

Abstract: SUMMARY In this paper we present a global model (GSRM-1) of both horizontal velocities on the Earth’s surface and horizontal strain rates for almost all deforming plate boundary zones. A model strain rate field is obtained jointly with a global velocity field in the process of solving for a global velocity gradient tensor field. In our model we perform a least-squares fit between model velocities and observed geodetic velocities, as well as between model strain rates and observed geological strain rates. Model velocities and strain rates are interpolated over a spherical Earth using bi-cubic Bessel splines. We include 3000 geodetic velocities from 50 different, mostly published, studies. Geological strain rates are obtained for central Asia only and they are inferred from Quaternary fault slip rates. For all areas where no geological information is included a priori constraints are placed on the style and direction (but not magnitude) of the model strain rate field. These constraints are taken from a seismic strain rate field inferred from (normalized) focal mechanisms of shallow earthquakes. We present a global solution of the second invariant of the model strain rate field as well as strain rate solutions for a few selected plate boundary zones. Generally, the strain rate tensor field is consistent with geological and seismological data. With the assumption of plate rigidity for all areas other than the plate boundary zones we also present relative angular velocities for most plate pairs. We find that in general there is a good agreement between the present-day plate motions we obtain and longterm plate motions, but a few significant differences exist. The rotation rates for the Indian, Arabian and Nubian plates relative to Eurasia are 30, 13 and 50 per cent slower than the NUVEL1A estimate, respectively, and the rotation rate for the Nazca Plate relative to South America is 17 per cent slower. On the other hand, Caribbean‐North America motion is 76 per cent faster than the NUVEL-1A estimate. While crustal blocks in the India‐Eurasia collision zone move significantly and self-consistently with respect to bounding plates, only a very small motion is predicted between the Nubian and Somalian plates. By integrating plate boundary zone deformation with the traditional modelling of angular velocities of rigid plates we have obtained a model that has already been proven valuable in, for instance, redefining a no-net-rotation model of surface motions and by confirming a global correlation between seismicity rates and tectonic moment rates along subduction zones and within zones of continental deformation.

A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation

Abstract: SUMMARY The perfectly matched layer absorbing boundary condition has proven to be very efficient for the elastic wave equation written as a first-order system in velocity and stress. We demonstrate how to use this condition for the same equation written as a second-order system in displacement. This facilitates use in the context of numerical schemes based upon such a system, e.g. the finite-element method, the spectral-element method and some finite-difference methods. We illustrate the efficiency of this second-order perfectly matched layer based upon 2-D benchmarks with body and surface waves.

Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method

Abstract: SUMMARY In this paper wavelet transforms and a logarithmic barrier method are applied to the inversion of large-scale magnetic data to recover a 3-D distribution of magnetic susceptibility. The fast wavelet transform is used, along with thresholding the small wavelet coefficients, to form a sparse representation of the sensitivity matrix. The reduced size of the resultant matrix allows the solution of large problems that are otherwise intractable. The compressed matrix is used to carry out fast forward modelling by performing matrix-vector multiplications in the wavelet domain. The reduction in CPU time is directly proportional to the compression ratio of the matrix. A second important feature of the algorithm used here is the use of an interior-point method of optimization to enforce positivity constraints. In this approach, the positivity is incorporated into the inversion by a sequence of non-linear optimizations approximated by truncated Newton steps. At the heart of the algorithm, a linear system of equations is solved. The conjugate gradient technique has been used as the basic solver to take the advantage of the efficient forward modelling offered by the sparse matrix representation. Overall, the combination of wavelet transforms, interior point optimization and conjugate gradient solutions readily allows us to solve magnetic inverse problems that have a few hundred thousand parameters and tens of thousands of data.

Determination of shallow shear wave velocity profiles in the Cologne, Germany area using ambient vibrations

Abstract: SUMMARY We have used both single-station and array methods to determine shallow shear velocity site profiles in the vicinity of the city of Cologne, Germany from ambient vibration records. Based on fk-analysis we assume that fundamental-mode Rayleigh waves dominate the analysed wavefield in the frequency range of 0.7‐2.2 Hz. According to this view a close relation exists between H/V spectral ratios and the ellipticity of the contributing Rayleigh waves. The inversion of the shape of H/V spectral ratios then provides quantitative information concerning the local shear wave velocity structure. However, based on tests with synthetic data believed to represent a typical situation in the Lower Rhine Embayment, dispersion curves were found to provide stronger constraints on the absolute values of the velocity‐depth model than the ellipticities. The shape of the ellipticities was found to be subject to a strong trade-off between the layer thickness and the average layer velocity. We have made use of this observation by combining the inversion schemes for dispersion curves and ellipticities such that the velocity‐depth dependence is essentially constrained by the dispersion curves while the layer thickness is constrained by the ellipticities. In order to test this method in practice, we have used array recordings of ambient vibrations from three sites where the subsurface geology is fairly well known and geotechnical information is at least partially available. In order to keep the parameter space as simple as possible we attempted to fit only a single layer over a half-space model. However, owing to earlier studies from the region, we assume a power-law depth dependence for sediment velocities. For all three sites investigated, the inversion resulted in models for which the shear wave velocity within the sediment layer both in absolute value at the surface and in depth dependence are found to be remarkably similar to the results obtained by Budny from downhole measurements. This is strong support for the interpretation of H/V spectral ratios as Rayleigh wave ellipticities. For all three sites the predicted SH-wave site amplification factors at the fundamental frequency are of the order of 5‐6 with a slightly smaller value south of Cologne.

On post-glacial sea level: I. General theory

Abstract: SUMMARY Modern analyses of sea level changes due to glacial isostatic adjustment (GIA) are based on the classic sea level equation derived by Farrell & Clark (1976, Geophys. J.R. astr. Soc., 46, 647–667). The connection between global sea level variations and changes to ocean height that is assumed within this equation breaks down in the presence of a time-varying shoreline geometry. We present a generalized sea level equation that overcomes this difficulty. We also derive analytic expressions for, and present schematic illustrations of, the error in the ocean height change over finite time intervals introduced in published efforts to incorporate shoreline evolution into the theory of GIA-induced sea level change. This comparison includes studies of shoreline migration due to either local sea level changes or the growth and ablation of marine-based ice. We conclude that the theories applied by Johnston (1993, Geophys. J. Int., 114, 615–634), Milne (1998, PhD thesis, University of Toronto, Toronto) and co-workers are more accurate than the procedure advocated by Peltier (1994, Science, 265, 195–201; 1998a, Geophys. Res. Lett., 25, 3955–3958; 1998b, Rev. Geophys., 114, 615–634), although an improvement in the latter has recently been reported (Peltier & Drummond (2002, Geophys. Res. Lett., 29, 10.1029/2001GL014273). Our generalized theory is valid for any Earth model. In a companion paper we derive the equations necessary to treat the special case of a spherically symmetric, linear viscoelastic and rotating Earth, and we quantify errors associated with previous work.

Elastic properties of dry clay mineral aggregates, suspensions and sandstones

Abstract: The presence of clay minerals can alter the elastic behaviour of rocks significantly. Although clay minerals are common in sedimentary formations and seismic measurements are our main tools for studying subsurface lithologies, measurements of elastic properties of clay minerals have proven difficult. Theoretical values for the bulk modulus of clay are reported between 20 and 50 GPa. The only published experimental measurement of Young's modulus in a clay mineral using atomic force acoustic microscopy (AFAM) gave a much lower value of 6.2 GPa. This study has concentrated on using independent experimental methods to measure the elastic moduli of clay minerals as functions of pressure and saturation. First, ultrasonic P- and S-wave velocities were measured as functions of hydrostatic pressure in cold-pressed clay aggregates with porosity and grain density ranging from 4 to 43 per cent and 2.13 to 2.83 g cm - 3 , respectively. In the second experiment, P- and S-wave velocities in clay powders were measured under uniaxial stresses compaction. In the third experiment, P-wave velocity and attenuation in a kaolinite-water suspension with clay concentrations between 0 and 60 per cent were measured at ambient conditions. Our elastic moduli measurements of kaolinite, montmorillonite and smectite are consistent for all experiments and with reported AFAM measurements on a nanometre scale. The bulk modulus values of the solid clay phase (K s ) lie between 6 and 12 GPa and shear (μ s ) modulus values vary between 4 and 6 GPa. A comparison is made between the accuracy of velocity prediction in shaley sandstones and clay-water and clay-sand mixtures using the values measured in this study and those from theoretical models. Using K s = 12 GPa and μ s = 6 GPa from this study, the models give a much better prediction both of experimental velocity reduction due to increase in clay content in sandstones and velocity measurements in a kaolinite-water suspension.

A shallow fault-zone structure illuminated by trapped waves in the Karadere–Duzce branch of the North Anatolian Fault, western Turkey

Abstract: SUMMARY We discuss the subsurface structure of the Karadere‐Duzce branch of the North Anatolian Fault based on analysis of a large seismic data set recorded by a local PASSCAL network in the 6 months following the Mw = 7.4 1999 Izmit earthquake. Seismograms observed at stations located in the immediate vicinity of the rupture zone show motion amplification and long-period oscillations in both P- and S-wave trains that do not exist in nearby off-fault stations. Examination of thousands of waveforms reveals that these characteristics are commonly generated by events that are well outside the fault zone. The anomalous features in fault-zone seismograms produced by events not necessarily in the fault may be referred to generally as fault-zone-related site effects. The oscillatory shear wave trains after the direct S arrival in these seismograms are analysed as trapped waves propagating in a low-velocity fault-zone layer. The time difference between the S arrival and trapped waves group does not grow systematically with increasing source‐receiver separation along the fault. These observations imply that the trapping of seismic energy in the Karadere‐Duzce rupture zone is generated by a shallow fault-zone layer. Traveltime analysis and synthetic waveform modelling indicate that the depth of the trapping structure is approximately 3‐4 km. The synthetic waveform modelling indicates further that the shallow trapping structure has effective waveguide properties consisting of thickness of the order of 100 m, a velocity decrease relative to the surrounding rock of approximately 50 per cent and an S-wave quality factor of 10‐15. The results are supported by large 2-D and 3-D parameter space studies and are compatible with recent analyses of trapped waves in a number of other faults and rupture zones. The inferred shallow trapping structure is likely to be a common structural element of fault zones and may correspond to the top part of a flower-type structure. The motion amplification associated with fault-zone-related site effects increases the seismic shaking hazard near fault-zone structures. The effect may be significant since the volume of sources capable of generating motion amplification in shallow trapping structures is large.

Inversion of shallow-seismic wavefields: I. Wavefield transformation

Abstract: SUMMARY I calculate Fourier–Bessel expansion coefficients for recorded shallow-seismic wavefields using a discrete approximation to the Bessel transformation. This is the first stage of a full-wavefield inversion. The transform is a complete representation of the data, recorded waveforms can be reconstructed from the expansion coefficients obtained. In a second stage (described in a companion paper) I infer a 1-D model of the subsurface from these transforms and P-wave arrival times by fitting them with their synthetic counterpart. The whole procedure avoids dealing with dispersion in terms of normal modes, but exploits the full signal-content, including the dispersion of higher modes, leaky modes and their true amplitudes. It is robust even in the absence of a priori information. I successfully apply it to the near field of the source. And it is more efficient than direct inversion of seismograms. I have developed this new approach because the inversion of shallow-seismic Rayleigh waves suffers from the interference of multiple modes that are present in the majority of our field data sets. Since even the fundamental-mode signal cannot be isolated in the time domain, conventional phase-difference techniques are not applicable. The potential to reconstruct the full waveform from the transform is confirmed by two field-data examples, which are recorded with 10 Hz geophones at effective intervals of about 1 m and spreads of less than 70 m length and are excited by a hammer source. Their transforms are discussed in detail, regarding aliasing and resolution. They reveal typical properties of shallow surface waves that are at variance with assumptions inherent to conventional inversion techniques: multiple modes contribute to the wavefield and overtones may dominate over the fundamental mode. The total wavefield may bear the signature of inverse or anomalous dispersion, although the excited modes have regular and normal dispersion. The resolution at long wavelengths (and thus the penetration depth of the survey) is limited by the length of the profile rather than by the signal-to-noise ratio at low frequencies. Finally, this approach is compared with conventional techniques of dispersion analysis. This illustrates the advantage of conserving the full wavefield in contrast to the reduction to one dispersion curve.

Partial melt or aqueous fluid in the mid-crust of Southern Tibet? Constraints from INDEPTH magnetotelluric data

Abstract: SUMMARY The INDEPTH project has applied modern geophysical techniques to the study of the crustal structure and tectonic evolution of the Tibetan Plateau. In the Lhasa Block, seismic reflection surveys in 1994 detected a number of bright-spots at 15‐20 km depths that indicate zones of crustal fluids (aqueous fluids or partial melt). Coincident magnetotelluric (MT) data collected in 1995 detected a major zone of high electrical conductivity at the same depth as the brightspots. Using constrained inversion, the MT data require a minimum crustal conductance of 6000 S. This abnormally high electrical conductance can be best explained by a layered model with fluids: partial melt, aqueous fluids or a combination of partial melt and aqueous fluids. The non-uniqueness of the MT method means that a wide range of melt fraction‐thickness combinations for the above models could all explain the 6000 S conductance. To distinguish between these three models, other geophysical and geological data are required. Reflection seismic data suggest that a high fluid content (>15 per cent) is present at the top of the layer. The amplitude-versus-offset data suggest that the top of this layer may be aqueous fluids rather than partial melt. Passive seismic data imaged a 20 km thick layer of lower fluid content that is probably partial melt. Petrological studies suggest that concentrations of aqueous fluids above 0.1 per cent at mid-crustal depth cannot be sustained. Taken together, these data show that the high conductivity in Southern Tibet is most probably the result of a relatively thin layer of aqueous fluids (100‐200 m) overlying a thicker zone of partial melt (>10 km).

Micro and macroscopic models of rock fracture

Abstract: SUMMARY The anelastic deformation of solids is often treated using continuum damage mechanics. An alternative approach to the brittle failure of a solid is provided by the discrete fiber-bundle model. Here we show that the continuum damage model can give exactly the same solution for material failure as the fiber-bundle model. We compare both models with laboratory experiments on the time dependent failure of chipboard and fiberglass. The power-law scaling obtained in both models and in the experiments is consistent with the power-law seismic activation observed prior to some earthquakes.

Crustal velocity field of western Europe from permanent GPS array solutions, 1996–2001

Abstract: We derive a new geodetic velocity field for western Europe and the Western Mediterranean by rigorously combining (1) a selection of 36 ITRF2000 sites, (2) a solution from a subset of sites of the European Permanent GPS Network (EUREF-EPN), (3) a solution of the French national geodetic permanent GPS network (RGP) and (4) a solution of a permanent GPS network in the western Alps (REGAL). The resulting velocity field describes horizontal crustal motion at 64 sites in Western Europe with an accuracy of the order of I mm yr - 1 or better. Its analysis shows that Central Europe (defined as east of the Rhine Graben and north of the Alps and the Carpathians) behaves rigidly at a 0.4 mm yr - 1 level and defines a stable Europe reference frame. In that reference frame, we find no significant motion at sites located west of the Rhine graben and on the Iberian peninsula, which sets an upper bound of 0.6 mm yr - 1 on horizontal motion across the Rhine graben and the Pyrenees. We find that the current strain pattern in the western Alps combines E-W extension and right-lateral shear. We confirm a counterclockwise rotation of the Adriatic microplate, which appears to control the strain pattern along its boundaries in the Friuli area, the Alps and the Apennines. Our results also suggest that the Africa-Eurasia plate motion in the Western Mediterranean may be 40-50 per cent slower that the NUVEL1A plate motion model and rotated 20°-30° counterclockwise.

Seismological constraints on the crustal structure beneath the Zagros Mountain belt (Iran)

Abstract: SUMMARY The Zagros Mountain belt of western Iran results from the collision of the Arabian and Central Iran continental blocks. The stage of the collision is unclear and the crustal structure of the Zagros is rather poorly known. In this study we investigate the velocity structure of the crust and upper mantle beneath the Ghir region located in the Central Zagros using data collected by a temporary local seismological network including a broad-band instrument. The structures of the sedimentary cover and the upper crystalline crust are estimated from the inversion of P and S traveltimes of local earthquakes recorded on a dense seismological network. The upper crust consists of an ∼11 km thick sedimentary layer (V p ∼ 4.70 km s −1 ) above a ∼8 km thick upper crystalline crust (V p ∼ 5.85 km s −1 ). The velocity of the lower crust and the depth of the Moho are found using receiver function analysis of teleseismic earthquakes. The lower crystalline crust is unusually slow (V p ∼ 6.5 km s −1 ) and ∼27 km thick. The upper bound for the total crustal thickness beneath the Ghir region is 46 ± 2 km. A comparison of the thickness of the crystalline crust of the Zagros with available information for the thickness of the crystalline crust of the Arabian Platform shows that, at present, the Zagros has a thinner crust. The current crustal thickness beneath the Zagros is comparable to the pre-collision crustal thickness of the Arabian Platform, suggesting that the Zagros is now in a very early stage of continental collision.

A deterministic seismic hazard map of India and adjacent areas.

Abstract: A seismic hazard map of the territory of India and adjacent areas has been prepared using a deterministic approach based on the computation of synthetic seismograms complete with all main phases. The input data set consists of structural models, seismogenic zones, focal mechanisms and earthquake catalogues. There are few probabilistic hazard maps available for the Indian subcontinent, however, this is the first study aimed at producing a deterministic seismic hazard map for the Indian region using realistic strong ground motion modelling with the knowledge of the physical process of earthquake generation, the level of seismicity and wave propagation in anelastic media. Synthetic seismograms at a frequency of 1 Hz have been generated at a regular grid of 0.2° x 0.2° by the modal summation technique. The seismic hazard, expressed in terms of maximum displacement (D m a x ), maximum velocity (V m a x ), and design ground acceleration (DGA), has been extracted from the synthetic signals and mapped on a regular grid over the studied territory. The estimated values of the peak ground acceleration are compared with the observed data available for the Himalayan region and are found to be in agreement. Many parts of the Himalayan region have DGA values exceeding 0.6 g. The epicentral areas of the great Assam earthquakes of 1897 and 1950 in northeast India represent the maximum hazard with DGA values reaching 1.2-1.3 g. The peak velocity and displacement in the same region is estimated as 120-177 cm s - 1 and 60-90 cm, respectively.

Resolution, stability and efficiency of resistivity tomography estimated from a generalized inverse approach

Abstract: SUMMARY Inconsistencies between an object and its image delivered by tomographical methods are inevitable. Loss of information occurs during the survey through incomplete and inaccurate data sampling and may also be introduced during the inverse procedure by smoothness constraints inadequate to the resolving power of the experimental setup. A quantitative appraisal of image quality (spatial resolution and image noise) is therefore not only required for successful interpretation of images but can be used together with measures of efficiency of the experimental design to optimize survey and inverse procedures. This paper introduces a low-contrast inversion scheme for electrical resistivity tomography that supports the reconstructed image with estimates of model resolution, model covariance and data importance. The algorithm uses a truncated pseudo-inverse and a line search approach to determine the maximum number of degrees of freedom necessary to fit the data to a prescribed target misfit. Though computationally expensive, the virtue of the method is that it reduces subjectivity by avoiding any empirically motivated model smoothness constraints. The method can be incorporated into a full non-linear inversion scheme for which a posteriori quality estimates can be calculated. In a numerical 2-D example the algorithm yielded reasonable agreement between object and image even for moderate resistivity contrasts of 10:100:1000. On the other hand, the resolving power of an exemplary four-electrode data set containing classical dipole–dipole and non-conventional configurations was shown to be severely affected by data inaccuracy. Insight into the resolving power as a function of space and data accuracy can be used as a guideline to designing optimized data sets, smoothness constraints and model parametrization.

Seismotectonics of the Sinai subplate – the eastern Mediterranean region

Abstract: SUMMARY We define the Sinai subplate, from a seismotectonic perspective, as a distinct component in the plate tectonics of the eastern Mediterranean region. This is based on the tectonic characteristics of a comprehensive listing of all ML≥ 4 recorded seismicity in the region during the 20th century, on newly calculated and recalculated fault plane mechanisms of first P-wave arrivals and on published solutions based on waveform inversion of broad-band data. The low seismicity level and scarcity of strong events in the region required a thorough search for useful data and a careful examination of the reliability of the focal solutions. We gathered all available records of first P-wave onsets from the ISS and ISC Bulletins and the local seismic networks. Altogether, we were able to calculate 48 new focal mechanisms and 33 recalculated ones of events that occurred during the years 1940–1992. With the increasing number of teleseismic and regional broad-band stations in the later years, we added 37 solutions based on teleseismic and regional waveform inversions of events that occurred during 1977–2001. These mechanisms enabled us to examine the seismotectonic character of the Sinai subplate. The strike and rake directions of the calculated mechanisms usually reflect the geometry and the large-scale type of deformation observed along the boundaries of the Sinai subplate—the Dead Sea Transform, the Cypriot Arc convergent zone and the Suez Rift. Nevertheless, along each of these boundaries we found anomalous solutions that attest to the complexity of the deformation processes along plate margins. Earthquakes along the Dead Sea Transform exhibit mainly sinistral transtension and transpression, reflecting its leaky manner and local change in the transform geometry. The presence of other unexpected mechanisms near the transform, however, reflects the heterogeneous deformation it induces around. As expected, thrust mechanisms along the Cypriot Arc mirror its convergent nature and typical curved geometry. Transtension and transpressional solutions in the eastern segment of the arc reflect the sinistral shear motion between Anatolia and Sinai there. However, shear mechanisms found between Cyprus and the Eratosthenes Seamount pose a problem regarding its collision process. Most intriguing of all are ML≥ 4 thrust and shear solutions found in the Gulf of Suez. They are associated with predominantly normal mechanisms within a rift zone and therefore constitute a unique phenomenon, yet to be deciphered.

Assessment of crustal velocity models using seismic refraction and reflection tomography

Abstract: Two tomographic methods for assessing velocity models obtained from wide-angle seismic traveltime data are presented through four case studies. The modelling/inversion of wide-angle traveltimes usually involves some aspects that are quite subjective. For example: (1) identifying and including later phases that are often difficult to pick within the seismic coda, (2) assigning specific layers to arrivals, (3) incorporating pre-conceived structure not specifically required by the data and (4) selecting a model parametrization. These steps are applied to maximize model constraint and minimize model non-uniqueness. However, these steps may cause the overall approach to appear ad hoc, and thereby diminish the credibility of the final model. The effect of these subjective choices can largely be addressed by estimating the minimum model structure required by the least subjective portion of the wide-angle data set: the first-arrival times. For data sets with Moho reflections, the tomographic velocity model can be used to invert the PmP times for a minimum-structure Moho. In this way, crustal velocity and Moho models can be obtained that require the least amount of subjective input, and the model structure that is required by the wide-angle data with a high degree of certainty can be differentiated from structure that is merely consistent with the data. The tomographic models are not intended to supersede the preferred models, since the latter model is typically better resolved and more interpretable. This form of tomographic assessment is intended to lend credibility to model features common to the tomographic and preferred models. Four case studies are presented in which a preferred model was derived using one or more of the subjective steps described above. This was followed by conventional first-arrival and reflection traveltime tomography using a finely gridded model parametrization to derive smooth, minimum-structure models. The case studies are from the SE Canadian Cordillera across the Rocky Mountain Trench, central India across the Narmada-Son lineament, the Iberia margin across the Galicia Bank, and the central Chilean margin across the Valparaiso Basin and a subducting seamount. These case studies span the range of modern wide-angle experiments and data sets in terms of shot‐receiver spacing, marine and land acquisition, lateral heterogeneity of the study area, and availability of wide-angle reflections and coincident near-vertical reflection data. The results are surprising given the amount of structure in the smooth, tomographically derived models that is consistent with the more subjectively derived models. The results show that exploiting the complementary nature of the subjective and tomographic approaches is an effective strategy for the analysis of wide-angle traveltime data.

A review of shear wave splitting in the crack-critical crust

Abstract: SUMMARY Over the last 15 years, it has become established that crack-induced stress-aligned shear wave splitting, with azimuthal anisotropy, is an inherent characteristic of almost all rocks in the crust. This means that most in situ rocks are pervaded by fluid-saturated microcracks and consequently are highly compliant. The evolution of such stress-aligned fluid-saturated grain-boundary cracks and pore throats in response to changing conditions can be calculated, in some cases with great accuracy, using anisotropic poro-elasticity (APE). APE is tightly constrained with no free parameters, yet dynamic modelling with APE currently matches a wide range of phenomena concerning anisotropy, stress, shear waves and cracks. In particular, APE has allowed the anisotropic response of a reservoir to injection to be calculated (predicted with hindsight), and the time and magnitude of an earthquake to be correctly stress-forecast. The reason for this calculability and predictability is that the microcracks in the crust are so closely spaced that they form critical systems. This crack-critical crust leads to a new style of geophysics that has profound implications for almost all aspects of pre-fracturing deformation of the crust and for solid-earth geophysics and geology. We review past, present and speculate about the future of shear wave splitting in the crack-critical crust. Shear wave splitting is seen to be a dynamic measure of the deformation of the rock mass. There is some good news and some bad news for conventional geophysics. Many accepted phenomena are no longer valid at high spatial and temporal resolution. A major effect is that the detailed crack geometry changes with time and varies from place to place in response to very small previously negligible changes. However, at least in some circumstances, the behaviour of the rock in the highly complex inhomogeneous Earth may be calculated and the response predicted, opening the way to possible control by feedback. The need is to devise ways to exploit these new opportunities in the crack-critical crust. Recent observations from the SMSITES Project at Husavik in Northern Iceland, gathered while this review was being written, display the extraordinarily sensitivity of in situ rock to small changes at great distances. The effects are far too large to occur in a conventional elastic brittle crust, and their presence confirms the highly compliant nature of the crack-critical crust.

The S-velocity structure of the East Asian mantle from inversion of shear and surface waveforms

Abstract: SUMMARY The structure of the mantle beneath East Asia down to 800 km depth is investigated using full waveforms of seismic shear and surface waves. Epicentral distances are limited to less than 40°. In contrast with previous waveform inversions, we avoid ray-theoretical or path-integral approaches. Instead, we use (1) exact 3-D waveform sensitivity kernels that correctly reflect off-path sensitivity and the existence of Fresnel zones. We apply (2) an accurate 3-D forward modelling technique based on a coupled-mode, multiple-forward-scattering approach, allowing us (3) to iterate the inversion procedure through several 3-D models and (4) to evaluate the true misfit between the data and the synthetics for the 3-D model. Average lateral resolution of the model in regions with good path coverage is 400 km throughout the upper mantle. In the depth range from 100 to 250 km the lateral resolution even approaches 200 km. Since wave front healing is taken into account, the amplitudes of velocity perturbations are larger than in other tomographic models. Moreover, the waveform sensitivity kernels provide an intrinsic physical smoothing of the model, stabilizing the inversion. Finally, owing to the use of exact sensitivities, better resolution can be achieved with a given set of seismograms. Major features of the model are high-velocity subducting slabs along the West Pacific subduction zones stagnating at the 660 km discontinuity, a strong low-velocity zone between 80 and 250 km depth in the West Pacific backarc regions, a plume-like low-velocity feature beneath the southern tip of the Baikal rift zone extending into the transition zone, a low-velocity region under the Tien Shan with connection into the transition zone and thick crust under Tibet reaching its maximum thickness of approximately 80 km close to the 35th parallel. The lithospheric mantle underneath southern Tibet is very fast indicating underthrusting of the Indian lithosphere. In contrast, the upper mantle beneath northern Tibet exhibits average-to-slow values from the Moho down to the 660 km discontinuity. Very high S velocities are again observed beneath the Tarim Basin. The high-velocity mantle lithospheres under southern Tibet and the Tarim join beneath the Karakorum range at the western tip of the Tibet Plateau.

Ritter Island Volcano—lateral collapse and the tsunami of 1888

Abstract: SUMMARY In the early morning of 1888 March 13, roughly 5 km 3 of Ritter Island Volcano fell violently into the sea northeast of New Guinea. This event, the largest lateral collapse of an island volcano to be recorded in historical time, flung devastating tsunami tens of metres high on to adjacent shores. Several hundred kilometres away, observers on New Guinea chronicled 3 min period waves up to 8 m high, that lasted for as long as 3 h. These accounts represent the best available first-hand information on tsunami generated by a major volcano lateral collapse. In this article, we simulate the Ritter Island landslide as constrained by a 1985 sonar survey of its debris field and compare predicted tsunami with historical observations. The best agreement occurs for landslides travelling at 40 m s −1 , but velocities up to 80 m s −1 cannot be excluded. The Ritter Island debris dropped little more than 800 m vertically and moved slowly compared with landslides that descend into deeper water. Basal friction block models predict that slides with shorter falls should attain lower peak velocities and that 40+ ms −1 is perfectly compatible with the geometry and runout extent of the Ritter Island landslide. The consensus between theory and observation for the Ritter Island waves increases our confidence in the existence of mega-tsunami produced by oceanic volcano collapses two to three orders of magnitude larger in scale.

T-matrix approach to shale acoustics

Abstract: SUMMARY The T-matrix approach of quantum scattering theory is used here to place many long-wavelength equivalent-medium approximations for porous composites, polycrystals and cracked media on a common footing and to indicate their limitations, but also to derive some new results based on two-point statistics. In this way, we have obtained an insight into the difficult problem of elastic inclusions at finite concentrations, which is of foremost relevance when estimating the effective material parameters of porous/cracked shales, involving stacks of more or less horizontally aligned clay platelets, mixed together with more rounded silt minerals, and with fluid filling the spaces. A rather involved perturbative analysis of the effects of interactions between (or structural correlations among) the various inclusions (minerals and cavities) making up a real shale of hexagonal symmetry was performed in an attempt to obtain a better match between theoretical predictions (based on a combination of coherent and optical potential approximations) and experimental results (recovered from ultrasonic wave speeds) for the effective elastic stiffness tensor. For the particular data set considered in this study, the T-matrix approach was able to match the data better than the approach of Hornby et al., but the match was not completely satisfactory. Further progress in theoretical shale modelling may come from a better knowledge of the elastic properties of pure clay minerals, a more detailed knowledge of the microstructure of shales, the incorporation of constraints obtained from comparisons between theoretical predictions and experimental results, as well as a continuing development of the T-matrix approach. Numerical results (also for the effect of bedding parallel microcracks on the elasticity of such a real shale) have value in illustrating the importance of taking into account the effects of spatial distribution when trying to deal with non-dilute mixtures of highly contrasting material properties.

Multiscale dynamics of the Tonga—Kermadec subduction zone

Abstract: Our understanding of mantle convection and the motion of plates depends intimately on our understanding of the viscosity structure of the mantle. While geoid and gravity observations have provided fundamental constraints on the radial viscosity structure of the mantle, the influence of short-wavelength variations in viscosity is still poorly understood. We present 2-D and 3-D finite-element models of mantle flow, including strong lateral viscosity variations and local sources of buoyancy, owing to both thermal and compositional effects. We first use generic 2-D models of a subduction zone to investigate how different observations depend on various aspects of the viscosity structure, in particular, the slab and lower-mantle viscosity and the presence of a low-viscosity region in the mantle wedge above the slab. We find that: (1) the strain rate provides a strong constraint on the absolute viscosity of the slab (10^(23) Pa s); (2) stress orientation within the slab is sensitive to the relative viscosity of the slab, lower mantle and the wedge; and (3) the stress state and topography of the overriding plate depend on the wedge viscosity and local sources of buoyancy. In particular, the state of stress in the overriding plate changes from compression to extension with the addition of a low-viscosity wedge. We then use observations of strain rate, stress orientation, dynamic topography and the geoid for the Tonga–Kermadec subduction zone as simultaneous constraints on the viscosity and buoyancy in a 3-D regional dynamic model. Together these observations are used to develop a self-consistent model of the viscosity and buoyancy by taking advantage of the sensitivity of each observation to different aspects of the dynamics, over a broad range of length-scales. The presence of a low-viscosity wedge makes it possible to match observations of shallow dynamic topography and horizontal extension within the backarc, and down-dip compression in the shallow portion of the slab. These results suggest that a low-viscosity wedge plays an important role in controlling the presence of backarc spreading. However, for a model with a low-viscosity and low-density region that provides a good fit to the observed topography, we find that a reduction of the slab density by a factor of 1.3 relative to the reference density model, is required to match the observed geoid. These results suggest that compensation of the slab by dynamic topography may be a much smaller effect at short to intermediate wavelengths than predicted by long-wavelength modelling of the geoid.

Inversion of shallow-seismic wavefields: II. Inferring subsurface properties from wavefield transforms

Abstract: SUMMARY The problem of inferring subsurface properties from shallow-seismic data is solved by a two-stage scheme that fits the full wavefield by its synthetic counterpart. In a first stage (described in a companion paper) I derive Fourier–Bessel expansion coefficients for the recorded data through a wavefield transformation. The present paper describes the joint inversion of these coefficients together with P-wave arrival times to infer subsurface properties. In this way we exploit the full signal-content including the dispersion of higher-modes, leaky-modes, their true amplitudes, and, at least partly, body waves. Owing to the multi-mode character of shallow-seismic field-data conventional techniques of dispersion analysis are not applicable. Since an initial model appropriate for inversion of full seismograms is rarely available in shallow seismics, the direct inversion of waveforms is not feasible. The wavefield transformation removes a remarkable amount of non-linearity from the data. In consequence the proposed method is robust even in the absence of a priori information. In distinction to the inversion of dispersion curves, it does not require the identification of normal-modes prior to inversion. The method performs well when applied to the multi-mode wavefields present in most shallow-seismic data sets. Compared to waveform fitting, it can be more efficient by about a factor of ten, because we need not evaluate the Bessel-expansion and therefore need less calculations of the forward problem. Subsurface properties are derived for the two sets of field-data that were already presented in the first paper. One of them includes a pronounced low-velocity channel. For both we observe a remarkably good resolution of S-velocity down to the bedrock, which is found in 6 and 16 m depth, respectively. In both cases it would not be possible to infer the depth of bedrock from P-wave data alone. Synthetic seismograms calculated from the final model match the recorded waveforms surprisingly well, although no waveform fitting was applied. A subsequent waveform inversion becomes feasible with initial models taken from the results of this method. Finally it is shown by example, that conventional techniques of dispersion-curve fitting are likely to give misleading results when applied to our field-data.

Comparison of azimuthal seismic anisotropy from surface waves and finite strain from global mantle-circulation models

Abstract: SUMMARY We present global models of strain accumulation in mantle flow to compare the predicted finite-strain ellipsoid (FSE) orientations with observed seismic anisotropy. The geographic focus is on oceanic and young continental regions where we expect our models to agree best with azimuthal anisotropy from surface waves. Finite-strain-derived models and alignment with the largest FSE axes lead to better model fits than the hypothesis of alignment of fast propagation orientation with absolute plate motions. Our modelling approach is simplified in that we are using a linear viscosity for flow and assume a simple relationship between strain and anisotropy. However, results are encouraging and suggest that similar models can be used to assess the validity of assumptions inherent in the modelling of mantle convection and lithospheric deformation. Our results substantiate the hypothesis that seismic anisotropy can be used as an indicator for mantle flow; circulation-derived models can contribute to the establishment of a quantitative functional relationship between the two.

Energetics of a two-phase model of lithospheric damage, shear localization and plate-boundary formation

Abstract: SUMMARY The two-phase theory for compaction and damage proposed by Bercovici et al. (2001a, J. Geophys. Res.,106, 8887‐8906) employs a nonequilibrium relation between interfacial surface energy, pressure and viscous deformation, thereby providing a model for damage (void generation and microcracking) and a continuum description of weakening, failure and shear localization. Here we examine further variations of the model which consider (1) how interfacial surface energy, when averaged over the mixture, appears to be partitioned between phases; (2) how variability in deformational-work partitioning greatly facilitates localization; and (3) how damage and localization are manifested in heat output and bulk energy exchange. Microphysical considerations of molecular bonding and activation energy suggest that the apparent partitioning of surface energy between phases goes as the viscosity of the phases. When such partitioning is used in the two-phase theory, it captures the melt-compaction theory of McKenzie (1984, J. Petrol., 25, 713‐765) exactly, as well as the void-damage theory proposed in a companion paper (Ricard & Bercovici, submitted). Calculations of 1-D shear localization with this variation of the theory still show at least three possible regimes of damage and localization: at low stress is weak localization with diffuse slowly evolving shear bands; at higher stress strong localization with narrow rapidly growing bands exists; and at yet higher shear stress it is possible for the system to undergo broadly distributed damage and no localization. However, the intensity of localization is strongly controlled by the variability of the deformational-work partitioning with dilation rate, represented by the parameter γ .F orγ � 1, extreme localization is allowed, with sharp profiles in porosity (weak zones), nearly discontinuous separation velocities and effectively singular dilation rates. Finally, the bulk heat output is examined for the 1-D system to discern how much deformational work is effectively stored as surface energy. In the high-stress, distributed-damage cases, heat output is reduced as more interfacial surface energy is created. Yet, in either the weak or strong localizing cases, the system always releases surface energy, regardless of the presence of damage or not, and thus slightly more heat is in fact released than energy is input through external work. Moreover, increased levels of damage (represented by the maximum work-partitioning f ∗ ) make the localizing system release surface energy faster as damage enhances phase separation and focusing of the porosity field, thus yielding more rapid loss of net interfacial surface area. However, when cases with different levels of damage are compared at similar stages of development (say, the peak porosity of the localization) it is apparent that increased damage causes smaller relative heat release and retards loss of net interfacial surface energy. The energetics and energy partitioning of this damage and shear-localization model are applied to estimating the energy costs of forming plate boundaries and generating plates from mantle convection.

Can the Earth's dynamo run on heat alone?

Abstract: SUMMARY The power required to drive the geodynamo places significant constraints on the heat passing across the core‐mantle boundary and the Earth’s thermal history. Calculations to date have been limited by inaccuracies in the properties of liquid iron mixtures at core pressures and temperatures. Here we re-examine the problem of core energetics in the light of new firstprinciples calculations for the properties of liquid iron. There is disagreement on the fate of gravitational energy released by contraction on cooling. We show that only a small fraction of this energy, that associated with heating resulting from changes in pressure, is available to drive convection and the dynamo. This leaves two very simple equations in the cooling rate and radioactive heating, one yielding the heat flux out of the core and the other the entropy gain of electrical and thermal dissipation, the two main dissipative processes. This paper is restricted to thermal convection in a pure iron core; compositional convection in a liquid iron mixture is considered in a companion paper. We show that heat sources alone are unlikely to be adequate to power the geodynamo because they require a rapid secular cooling rate, which implies a very young inner core, or a combination of cooling and substantial radioactive heating, which requires a very large heat flux across the core‐mantle boundary. A simple calculation with no inner core shows even higher heat fluxes are required in the absence of latent heat before the inner core formed.

Stress, fluid pressure and structural permeability in seismogenic crust, North Island, New Zealand

Abstract: SUMMARY Stress and fluid-pressure conditions within seismogenic crust are compared for two subparallel belts of active deformation and fluid redistribution associated with the obliquely convergent Pacific‐Australia plate boundary in the North Island of New Zealand. Whereas seismic activity on extensional normal faults in the arc-backarc Taupo volcanic zone is restricted to <8 km depth in a high heat-flow, near-hydrostatic fluid-pressure regime undergoing vigorous hydrothermal convection, rupturing along the thrust interface of the contractional Hikurangi subduction margin and in its hangingwall extends to ∼25 km depth in crust with fluids overpressured towards lithostatic values. The contrast in fluid-pressure levels stems partly from the abundance of low-permeability mudrocks in the forearc and partly from superior containment of overpressures by a compressional thrust-fault regime. Maximum supportable levels of differential stress and fluid pressure are critically interdependent in the overpressured regime of the Hikurangi subduction margin. Frictional instability leading to fault rupture in such settings may be triggered by increasing fluid pressure as well as by accumulating shear stress, so that nucleation and recurrence of earthquake ruptures are likely to be affected by cycling of fluid pressure through fault-valve action as well as by stress accumulation. Coupling across the subduction interface is also likely to be highly sensitive to the degree of overpressuring. Different factors are responsible for the localization of active deformation within the two crustal seismic belts. Within the magmatically active Taupo volcanic zone, thermal weakening is clearly responsible for concentrating seismicity and deformation with respect to the surrounding crust. However, in the hangingwall of the Hikurangi subduction margin, where heat flow has been reduced by subduction refrigeration and frictional interaction extends to ∼25 km depth, relative weakening arises principally from reduction of frictional resistance by fluid overpressuring.

Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography

Abstract: SUMMARY Mantle convection models require an initial condition some time in the past. Because this initial condition is unknown for Earth, we cannot infer the geological evolution of mantle flow from forward mantle convection calculations even for the most recent Mesozoic and Cenozoic geological history of our planet. Here we introduce a fluid dynamic inverse problem to constrain unknown mantle flow back in time from seismic tomographic observations of the mantle and reconstructions of past plate motions using variational data assimilation. We derive the generalized inverse of mantle convection and explore the initial condition problem in high-resolution, 3-D spherical mantle circulation models for a time period of 100 Myr, roughly comparable to half a mantle overturn. We present a synthetic modelling experiment to demonstrate that mid-Cretaceous mantle structure can be inferred accurately from fluid dynamic inverse modelling, assuming present-day mantle structure is well-known, even if an initial first guess assumption about the mid-Cretaceous mantle involved only a simple 1-D radial temperature profile. We also demonstrate that convecting present-day mantle structure back in time by reversing the time-stepping of the energy equation is insufficient to model the mantle structure of the past. The difficulty arises, because such backward convection calculations ignore thermal diffusion effects, and therefore cannot account for the generation of thermal buoyancy in boundary layers as we go back in time. Inverse mantle convection modelling should make it possible to infer a number of flow parameters from observational constraints of the mantle.

Analytical solutions for deformable elliptical inclusions in general shear

Abstract: SUMMARY Using Muskhelishvili’s method, we present closed-form analytical solutions for an isolated elliptical inclusion in general shear far-field flows. The inclusion is either perfectly bonded to the matrix or, as in the case of a circular inclusion, to a possible intermediate layer. The solutions are valid for incompressible all-elastic or all-viscous systems. The actual values of the shear modulus or viscosity in the inclusion, mantle and matrix can be different and no limits are imposed on the possible material property contrasts. The solutions presented are complete 2-D solutions and the parameters that can be analysed include all kinematic (e.g. stream functions, velocities, strain rates, strains) and dynamic parameters (e.g. pressure, maximum shear stress, etc.). We refrain from giving the tedious derivation of the presented solutions, and instead we focus on how to use the solutions, how to extract the parameters of interest and how to apply and verify them. In order to demonstrate the usefulness of Muskhelishvili’s method for slow viscous flow problems, we apply our results to a clast in a shear zone and obtain important new insights, even for simple, circular inclusions. Another important application is the benchmarking of numerical codes for which the presented solutions are most suitable due to the infinite range of viscosity contrasts and the strong local gradients of properties and results. To stimulate a broader use of Muskhelishvili’s method, all solutions are implemented in MATLAB and downloadable from the web.

Source parameters of large historical (1917-1961) earthquakes, North Island, New Zealand

Abstract: We have studied 18 earthquakes of M w > 5.5 within the North Island, New Zealand region to determine fault parameters for earthquake hazards and regional tectonic studies. Since most large (M w > 6.5) earthquakes in the North Island occurred prior to 1961, these events provide important information to supplement studies of more recent seismicity. Our results indicate that no large plate interface earthquakes have occurred in the central and southern North Island for the past 80 years, although six strike-slip earthquakes of M w ≥ 6.8 occurred within the Australian (upper) Plate during an extremely active 25 year period that began in 1917. Only the M w > 7.2 strike-slip events were associated with surface faulting, indicating the difficulties that may arise in attempting to identify active faults within this region. Two earthquakes (M w = 6.9-7.1) off the northeastern North Island in 1947 appear to have occurred along the plate interface and were associated with local tsunamis having runup heights of up to 10 m. Two M w = 6.8 events also occurred within the Pacific (lower) Plate, highlighting the hazards related to intraslab events. Slip vectors for the earthquakes studied suggest that the majority of transcurrent motion along the plate margin is accommodated within the Australian Plate, similar to the results obtained from studies of more recent, smaller earthquakes. Pure thrusting occurs along the plate interface and T axes of intraslab events indicate downdip tension in the Pacific Plate.