Showing papers in "Geophysical Journal International in 2006"

A review of the adjoint-state method for computing the gradient of a functional with geophysical applications

Abstract: SUMMARY Estimating the model parameters from measured data generally consists of minimizing an error functional. A classic technique to solve a minimization problem is to successively determine the minimum of a series of linearized problems. This formulation requires the Frechet derivatives (the Jacobian matrix), which can be expensive to compute. If the minimization is viewed as a non-linear optimization problem, only the gradient of the error functional is needed. This gradient can be computed without the Frechet derivatives. In the 1970s, the adjoint-state method was developed to efficiently compute the gradient. It is now a well-known method in the numerical community for computing the gradient of a functional with respect to the model parameters when this functional depends on those model parameters through state variables, which are solutions of the forward problem. However, this method is less well understood in the geophysical community. The goal of this paper is to review the adjoint-state method. The idea is to define some adjoint-state variables that are solutions of a linear system. The adjoint-state variables are independent of the model parameter perturbations and in a way gather the perturbations with respect to the state variables. The adjoint-state method is efficient because only one extra linear system needs to be solved. Several applications are presented. When applied to the computation of the derivatives of the ray trajectories, the link with the propagator of the perturbed ray equation is established.

Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis: I - Phase velocity maps

Abstract: SUMMARY Empirical Green’s functions (EGFs) between pairs of seismographs can be estimated from the time derivative of the long-time cross-correlation of ambient seismic noise. These EGFs reveal velocity dispersion at relatively short periods, which can be used to resolve structures in the crust and uppermost mantle better than with traditional surface-wave tomography. We combine Rayleigh-wave dispersion estimates from EGFs and from traditional two-station (TS) analysis into a new approach to surface-wave array tomography with data from dense receiver arrays. We illustrate the methodology with continuous broad-band recordings from a temporary seismographic network on the southeastern part of the Tibetan plateau, in Sichuan and Yunnan provinces, SW China. The EGFs are robust under temporal changes in regional seismicity and the use of either ambient noise (approximated by records without signal from events with magnitude mb ≥ 5 or 4) or surface wave coda produces similar results. The EGFs do not strongly depend on the presence of large earthquakes, but they are not reciprocal for stations aligned in the N‐S direction. This directionality reflects the paucity of seismicity to the north of the array. Using a far-field representation of the surface-wave Green’s function and an image transformation technique, we infer from the EGFs the Rayleigh-wave phase velocity dispersion in the period band from 10‐30 s. A classical TS approach is used to determine Rayleigh-wave phase velocity dispersion between 20‐120 s. Together, they constrain phase velocity variations for T = 10‐120 s, which can be used to study the structure from the crust to the upper mantle. Beneath SE Tibet, short and intermediate period (10‐80 s) phase velocities are prominently low, suggesting that the crust and upper mantle beneath SE Tibet is characterized by slow shear wave propagation.

Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields

Abstract: SUMMARY We study numerically an extensive set of dynamo models in rotating spherical shells, varying all relevant control parameters by at least two orders of magnitude. Convection is driven by a fixed temperature contrast between rigid boundaries. There are two distinct classes of solutions with strong and weak dipole contributions to the magnetic field, respectively. Non-dipolar dynamos are found when inertia plays a significant role in the force balance. In the dipolar regime the critical magnetic Reynolds number for self-sustained dynamos is of order 50, independent of the magnetic Prandtl number Pm. However, dynamos at low Pm exist only at sufficiently low Ekman number E. For dynamos in the dipolar regime we attempt to establish scaling laws that fit our numerical results. Assuming that diffusive effects do not play a primary role, we introduce non-dimensional parameters that are independent of any diffusivity. These are a modified Rayleigh number based on heat (or buoyancy) flux Ra ∗ , the Rossby number Ro measuring the flow velocity, the Lorentz number Lo measuring magnetic field strength, and a modified Nusselt number Nu ∗ for the advected heat flow. To first approximation, all our dynamo results can be collapsed into simple power-law dependencies on the modified Rayleigh number, with approximate exponents of 2/5, 1/2 and 1/3 for the Rossby number, modified Nusselt number and Lorentz number, respectively. Residual dependencies on the parameters related to diffusion (E, Pm, Prandtl number Pr) are weak. Our scaling laws are in agreement with the assumption that the magnetic field strength is controlled by the available power and not necessarily by a force balance. The Elsasser number � , which is the conventional measure for the ratio of Lorentz force to Coriolis force, is found to vary widely. We try to assess the relative importance of the various forces by studying sources and sinks of enstrophy (squared vorticity). In general Coriolis and buoyancy forces are of the same order, inertia and viscous forces make smaller and variable contributions, and the Lorentz force is highly variable. Ignoring a possible weak dependence on the Prandtl numbers or the Ekman number, a surprising prediction is that the magnetic field strength is independent both of conductivity and of rotation rate and is basically controlled by the buoyancy flux. Estimating the buoyancy flux in the Earth’s core using our Rossby number scaling and a typical velocity inferred from geomagnetic secular variations, we predict a small growth rate and old age of the inner core and obtain a reasonable magnetic field strength of order 1 mT inside the core. From the observed heat flow in Jupiter, we predict an internal field of 8 mT, in agreement with Jupiter’s external field being 10 times stronger than that of the Earth.

The detection of low magnitude seismic events using array-based waveform correlation

Abstract: It has long been accepted that occurrences of a known signal are most effectively detected by cross-correlating the incoming data stream with a waveform template. Such matched signal detectors have received very little attention in the field of detection seismology because there are relatively few instances in which the form of an anticipated seismic signal is known a priori. Repeating events in highly confined geographical regions have been observed to produce very similar waveforms and good signals from events at a given site can be exploited to detect subsequent co-located events at lower magnitudes than would be possible using traditional power detectors. Even greater improvement in signal detectability can be achieved using seismic arrays; running correlation coefficients from single sensors can be stacked over an array or network to result in a network correlation coefficient displaying a significant array gain. If two events are co-located, the time separating the corresponding patterns in the wave train as indicated by the cross-correlation function is identical for all seismic stations and this property means that the correlation coefficient traces are coherent even when the waveforms are not. We illustrate the power of array-based waveform correlation using the 1997 August 16 Kara Sea event. The weak event that occurred 4 hr after the main event was barely detected using an STA/LTA detector on the SPITS array but is readily detected by signal matching on a single channel. The main event was also recorded by the far more distant NORSAR array but no conventional detection can be made for the second event. A clear detection is, however, made when the correlation coefficient traces are beamformed over all sensors of the array. We estimate the reduction in detection threshold of a test signal on a regional seismic array using waveform correlation by scaling down a master signal and immersing it into seismic noise. We show that, for this case, waveform correlation using a single channel detects signals of approximately 0.7 orders of magnitude lower than is possible using an STA/LTA detector on the array beam. Waveform matching on the full array provides an additional improvement of approximately 0.4 magnitude units. We describe a case study in which small seismic events at the Barentsburg coal mine on Spitsbergen were detected using the signals from a major rockburst as master waveforms. Many spurious triggers occurred in this study whereby short sections of signal exhibited coincidental similarity with unrelated incoming wave fronts. We demonstrate how such false alarms can almost always be identified and screened out automatically by performing frequency–wavenumber analysis upon the set of individual correlation coefficient traces.

Three‐dimensional modelling and inversion of dc resistivity data incorporating topography – II. Inversion

Abstract: SUMMARY We present a novel technique for the determination of resistivity structures associated with arbitrary surface topography. The approach represents a triple-grid inversion technique that is based on unstructured tetrahedral meshes and finite-element forward calculation. The three grids are characterized as follows: A relatively coarse parameter grid defines the elements whose resistivities are to be determined. On the secondary field grid the forward calculations in each inversion step are carried out using a secondary potential (SP) approach. The primary fields are provided by a one-time simulation on the highly refined primary field grid at the beginning of the inversion process. We use a Gauss‐Newton method with inexact line search to fit the data within error bounds. A global regularization scheme using special smoothness constraints is applied. The regularization parameter compromising data misfit and model roughness is determined by an L-curve method and finally evaluated by the discrepancy principle. To solve the inverse subproblem efficiently, a least-squares solver is presented. We apply our technique to synthetic data from a burial mound to demonstrate its effectiveness. A resolution-dependent parametrization helps to keep the inverse problem small to cope with memory limitations of today’s standard PCs. Furthermore, the SP calculation reduces the computation time significantly. This is a crucial issue since the forward calculation is generally very time consuming. Thus, the approach can be applied to large-scale 3-D problems as encountered in practice, which is finally proved on field data. As a by-product of the primary potential calculation we obtain a quantification of the topography effect and the corresponding geometric factors. The latter are used for calculation of apparent resistivities to prevent the reconstruction process from topography induced artefacts.

H/V ratio: a tool for site effects evaluation. Results from 1-D noise simulations

Abstract: SUMMARY Ambient vibration techniques such as the H/V method may have the potential to significantly contribute to site effect evaluation, particularly in urban areas. Previous studies interpret the so-called Nakamura’s technique in relation to the ellipticity ratio of Rayleigh waves, which, for a high enough impedance contrast, exhibits a pronounced peak close to the fundamental S-wave resonance frequency. Within the European SESAME project (Site EffectS assessment using AMbient Excitations) this interpretation has been tested through noise numerical simulation under well-controlled conditions in terms of source type and distribution and propagation structure. We will present simulations for a simple realistic site (one sedimentary layer over bedrock) characterized by a rather high impedance contrast and low quality factor. Careful H/V and array analysis on these noise synthetics allow an in-depth investigation of the link between H/V ratio peaks and the noise wavefield composition for the soil model considered here: (1) when sources are near (4 to 50 times the layer thickness) and surficial, H/V curves exhibit one single peak, while the array analysis shows that the wavefield is dominated by Rayleigh waves; (2) when sources are distant (more than 50 times the layer thickness) and located inside the sedimentary layer, two peaks show up on the H/V curve, while the array analysis indicates both Rayleigh waves and strong S head waves; the first peak is due to both fundamental Rayleigh waves and resonance of head S waves, the second is only due to the resonance of head S waves; (3) when sources are deep (located inside the bedrock), whatever their distance, H/V ratio exhibit peaks at the fundamental and harmonic resonance frequencies, while array analyses indicate only non-dispersive body waves; the H/V is thus simply due to multiple reflections of S waves within the layer. Therefore, considering that experimental H/V ratio (i.e. derived from actual noise measured in the field) exhibit in most cases only one peak, we conclude that H/V ratio is (1) mainly controlled by local surface sources, (2) mainly due to the ellipticity of the fundamental Rayleigh waves. Then the amplitude of H/V peak is not able to give a good estimate of site amplification factor.

An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes – II. The three-dimensional isotropic case

Abstract: SUMMARY We present a new numerical method to solve the heterogeneous elastic wave equations formulated as a linear hyperbolic system using first-order derivatives with arbitrary high-order accuracy in space and time on 3-D unstructured tetrahedral meshes. The method combines the Discontinuous Galerkin (DG) Finite Element (FE) method with the ADER approach using Arbitrary high-order DERivatives for flux calculation. In the DG framework, in contrast to classical FE methods, the numerical solution is approximated by piecewise polynomials which allow for discontinuities at element interfaces. Therefore, the well-established theory of numerical fluxes across element interfaces obtained by the solution of Riemann-Problems can be applied as in the finite volume framework. To define a suitable flux over the element surfaces, we solve so-called Generalized Riemann-Problems (GRP) at the element interfaces. The GRP solution provides simultaneously a numerical flux function as well as a time-integration method. The main idea is a Taylor expansion in time in which all time-derivatives are replaced by space derivatives using the so-called Cauchy–Kovalewski or Lax–Wendroff procedure which makes extensive use of the governing PDE. The numerical solution can thus be advanced for one time step without intermediate stages as typical, for example, for classical Runge–Kutta time stepping schemes. Due to the ADER time-integration technique, the same approximation order in space and time is achieved automatically. Furthermore, the projection of the tetrahedral elements in physical space on to a canonical reference tetrahedron allows for an efficient implementation, as many computations of 3-D integrals can be carried out analytically beforehand. Based on a numerical convergence analysis, we demonstrate that the new schemes provide very high order accuracy even on unstructured tetrahedral meshes and computational cost and storage space for a desired accuracy can be reduced by higher-order schemes. Moreover, due to the choice of the basis functions for the piecewise polynomial approximation, the new ADER–DG method shows spectral convergence on tetrahedral meshes. An application of the new method to a well-acknowledged test case and comparisons with analytical and reference solutions, obtained by different well-established methods, confirm the performance of the proposed method. Therefore, the development of the highly accurate ADER–DG approach for tetrahedral meshes provides a numerical technique to approach 3-D wave propagation problems in complex geometry with unforeseen accuracy.

An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes - I. The two-dimensional isotropic case with external source terms

Abstract: SUMMARY We present a new numerical approach to solve the elastic wave equation in heterogeneous media in the presence of externally given source terms with arbitrary high-order accuracy in space and time on unstructured triangular meshes. We combine a discontinuous Galerkin (DG) method with the ideas of the ADER time integration approach using Arbitrary high-order DERivatives. The time integration is performed via the so-called Cauchy-Kovalewski procedure using repeatedly the governing partial differential equation itself. In contrast to classical finite element methods we allow for discontinuities of the piecewise polynomial approximation of the solution at element interfaces. This way, we can use the well-established theory of fluxes across element interfaces based on the solution of Riemann problems as developed in the finite volume framework. In particular, we replace time derivatives in the Taylor expansion of the time integration procedure by space derivatives to obtain a numerical scheme of the same high order in space and time using only one single explicit step to evolve the solution from one time level to another. The method is specially suited for linear hyperbolic systems such as the heterogeneous elastic wave equations and allows an efficient implementation. We consider continuous sources in space and time and point sources characterized by a Delta distribution in space and some continuous source time function. Hereby, the method is able to deal with point sources at any position in the computational domain that does not necessarily need to coincide with a mesh point. Interpolation is automatically performed by evaluation of test functions at the source locations. The convergence analysis demonstrates that very high accuracy is retained even on strongly irregular meshes and by increasing the order of the ADER‐DG schemes computational time and storage space can be reduced remarkably. Applications of the proposed method to Lamb’s Problem, a problem of strong material heterogeneities and to an example of global seismic wave propagation finally confirm its accuracy, robustness and high flexibility.

A Three-Dimensional Radially-Anisotropic Model of Shear Velocity in the Whole Mantle

Abstract: SUMMARY We present a 3-D radially anisotropic S velocity model of the whole mantle (SAW642AN), obtained using a large three component surface and body waveform data set and an iterative inversion for structure and source parameters based on Non-linear Asymptotic Coupling Theory (NACT). The model is parametrized in level 4 spherical splines, which have a spacing of ! 8 " . The model shows a link between mantle flow and anisotropy in a variety of depth ranges. In the uppermost mantle, we confirm observations of regions with VSH > VSV starting at ! 80 km under oceanic regions and ! 200 km under stable continental lithosphere, suggesting horizontal flow beneath the lithosphere. We also observe a VSV > VSH signature at ! 150‐300 km depth beneath major ridge systems with amplitude correlated with spreading rate for fast-spreading segments. In the transition zone (400‐700 km depth), regions of subducted slab material are associated with VSV > VSH, while the ridge signal decreases. Whilethemid-mantlehasloweramplitudeanisotropy( VSV in the lowermost 300 km, which appears to be a robust conclusion, despite an error in our previous paper which has been corrected here. The 3-D deviations from this signature are associated with the large-scale low-velocity superplumes underthecentralPacificandAfrica,suggestingthatVSH >VSV isgeneratedinthepredominant horizontal flow of a mechanical boundary layer, with a change in signature related to transition to upwelling at the superplumes.

Relocation and assessment of seismicity in the Iran region

Abstract: SUMMARY More than 2000 instrumentally recorded earthquakes occurring in the Iran region during the period 1918‐2004 have been relocated and reassessed, with special attention to focal depth, using an advanced technique for 1-D earthquake location. A careful review of starting depths, association of teleseismic depth phases, and the effects of reading errors on these phases are made and, when necessary, waveforms have been examined to better constrain EHB focal depths. Uncertainties in EHB epicentres are on the order of 10‐15 km in the Iran region, owing to the Earth’s lateral heterogeneity and uneven station distribution. Uncertainties of reviewed EHB focal depth estimates are on the order of 10 km, as compared to about 4 km for long-period P and SH body-waveform inversions. Nevertheless, these EHB depth estimates are sufficiently accurate to resolve robust differences in focal depth distribution throughout the Iran region and, within their errors, show patterns that are in agreement with the smaller number of earthquakes whose depths have been confirmed by body-wave modelling or local seismic networks. The importance of this result is that future earthquakes with apparently anomalous depths can easily be identified, and checked, if necessary. Most earthquakes in the Iranian continental lithosphere occur in the upper crust, with the crustal shortening produced by continental collision accommodated entirely by thickening and distributed deformation. In the Zagros Mountains nearly all earthquakes are confined to the upper crust (depths <20 km), and there is no evidence for a seismically active subducted slab dipping NE beneath central Iran. By contrast, in southeastern Iran, where the Arabian seafloor is being subducted beneath the Makran coast, low-level earthquake activity occurs in the upper crust as well as to depths of at least 150 km within a northward-dipping subducting slab. Near the Oman Line, a region transitional between the Zagros and the Makran, seismicity extends to depths of up to 30‐45 km in the crust, consistent with low-angle thrusting of Arabian basement beneath central Iran. In north-central Iran, along the Alborz mountain belt, seismic activity occurs primarily in the upper crust but with some infrequent events in the lower crust, particularly in the western part of the belt (the Talesh), where the South Caspian basin underthrusts NW Iran. Earthquakes that occur in a band across the central Caspian, following the Apscheron‐Balkhan sill between Azerbaijan and Turkmenistan, have depths in the range 30‐100 km, deepening northwards. These are thought to be connected with either incipient or remnant northeast subduction of the South Caspian basin basement beneath the east-west trending Apscheron‐Balkhan sill. Curiously, in this region of genuine mantle seismicity, there is no evidence for earthquakes shallower than 30 km.

Large igneous provinces generated from the margins of the large low-velocity provinces in the deep mantle

Abstract: SUMMARY There is a clear correlation between downward projected large igneous province (LIP) eruption sites of the past 200 Myr and the margins of the large low-velocity provinces (LLVPs) at the base of the mantle. We established this correlation by using palaeomagnetic as well as fixed and moving hotspot reference frames. Our finding indicates that the majority of the LIPs have been generated by plumes that rose from the D �� zone at the edges of the LLVPs. Most LIP eruption sites project radially downwards to the core‐mantle boundary (CMB) within ±10 ◦ of the 1 per cent slow shear wave velocity contour in the SMEAN tomographic model. Steep shear wave velocity gradients have been mapped near the CMB along much of the lengths of the LLVP margins close to that contour which marks a faster/slower boundary (FSB) within the D �� zone. The observation that eruption sites of LIPs as old as 200 Myr can be linked to this prominent present day seismic structure shows that the FSBs of the two LLVPs have occupied their current positions for at least as long and that the process that leads to the generation of deep-seated plumes has been localized on the FSBs at the margins of the African and Pacific LLVPs for the same interval. The persistence of the LLVPs over 200 Myr is consistent with independent evidence that they are compositionally distinct and are not just simply hotter than the material making up the rest of the D �� zone.

Trans-dimensional inverse problems, model comparison and the evidence

Abstract: SUMMARY In most geophysical inverse problems the properties of interest are parametrized using a fixed number of unknowns. In some cases arguments can be used to bound the maximum number of parameters that need to be considered. In others the number of unknowns is set at some arbitrary value and regularization is used to encourage simple, non-extravagant models. In recent times variable or self-adaptive parametrizations have gained in popularity. Rarely, however, is the number of unknowns itself directly treated as an unknown. This situation leads to a transdimensional inverse problem, that is, one where the dimension of the parameter space is a variable to be solved for. This paper discusses trans-dimensional inverse problems from the Bayesian viewpoint. A particular type of Markov chain Monte Carlo (MCMC) sampling algorithm is highlighted which allows probabilistic sampling in variable dimension spaces. A quantity termed the evidence or marginal likelihood plays a key role in this type of problem. It is shown that once evidence calculations are performed, the results of complex variable dimension sampling algorithms can be replicated with simple and more familiar fixed dimensional MCMC sampling techniques. Numerical examples are used to illustrate the main points. The evidence can be difficult to calculate, especially in high-dimensional non-linear inverse problems. Nevertheless some general strategies are discussed and analytical expressions given for certain linear problems.

Models of large‐scale viscous flow in the Earth's mantle with constraints from mineral physics and surface observations

Abstract: SUMMARY Modelling the geoid has been a widely used and successful approach in constraining flow and viscosity in the Earth’s mantle. However, details of the viscosity structure cannot be tightly constrained with this approach. Here, radial viscosity variations in four to five mantle layers (lithosphere, upper mantle, one to two transition zone layers, lower mantle) are computed with the aid of independent mineral physics results. A density model is obtained by converting s-wave anomalies from seismic tomography to density anomalies. Assuming both are of thermal origin, conversion factors are computed based on mineral physics results. From the density and viscosity model, a model of mantle flow, and the resulting geoid and radial heat flux profile are computed. Absolute viscosity values in the mantle layers are treated as free parameters and determined by minimizing a misfit function, which considers fit to geoid, ‘Haskell average’ determined from post-glacial rebound and the radial heat flux profile and penalizes if at some depth computed heat flux exceeds the estimated mantle heat flux 33 TW. Typically, optimized models do not exceed this value by more than about 20 per cent and fit the Haskell average well. Viscosity profiles obtained show a characteristic hump in the lower mantle, with maximum viscosities of about 10 23 Pa s just above the D �� layer— several hundred to about 1000 times the lowest viscosities in the upper mantle. This viscosity contrast is several times higher than what is inferred when a constant lower mantle viscosity is assumed. The geoid variance reduction obtained is up to about 80 per cent—similar to previous results. However, because of the use of mineral physics constraints, a rather small number of free model parameters is required, and at the same time, a reasonable heat flux profile is obtained. Results are best when the lowest viscosities occur in the transition zone. When viscosity is lowest in the asthenosphere, variance reduction is about 70‐75 per cent. Best results were obtained with a viscous lithosphere with a few times 10 22 Pa s. The optimized models yield a core-mantle boundary excess ellipticity several times higher than observed, possibly indicating that radial stresses are partly compensated due to non-thermal lateral variations within the lowermost mantle.

The influence of fluid-sensitive dispersion and attenuation on AVO analysis

Abstract: SUMMARY Analysis of seismic data suggests that hydrocarbon deposits are often associated with higher than usual values of attenuation, but this is generally ignored during amplitude-versus-offset (AVO) analysis. The effect can be modelled with equivalent medium theory based on the squirt flow concept, but the excess attenuation is associated with strong velocity dispersion. Consequently, when we study reflections from the interface between such an equivalent medium and an elastic overburden we find that the reflection coefficient varies with frequency. The impact of this variation depends on the AVO behaviour at the interface; class I reflections tend to be shifted to higher frequency while class III reflections have their lower frequencies amplified. We calculate synthetic seismograms for typical models using the reflectivity method for materials with frequency dependent velocities and attenuations, and find that these effects are predicted to be detectable on stacked data. Two field data sets show frequency anomalies similar to those predicted by the analysis, and we suggest that our modelling provides a plausible explanation of the observations.

CHAOS—a model of the Earth's magnetic field derived from CHAMP, Ørsted, and SAC‐C magnetic satellite data

Abstract: SUMMARY We have derived a model of the near-Earth magnetic field (up to spherical harmonic degree n = 50 for the static field, and up to n = 18 for the first time derivative) using more than 6.5 yr of high-precision geomagnetic measurements from the three satellites Orsted, CHAMP and SAC-C taken between 1999 March and 2005 December. Our modelling approach goes in several aspects beyond that used for recent models: (i) we use different data selection criteria and allow for higher geomagnetic activity (index Kp ≤ 2o), thus we include more data than previous models; (ii) we describe the temporal variation of the core field by splines (for n ≤ 14); (iii) we take magnetometer vector data in the instrument frame and co-estimate the Euler angles that describe the transformation from the magnetometer frame to the star imager frame, avoiding the inconsistency of using vector data that have been aligned using a different (pre-existing) field model; (iv) we account for the bending of the CHAMP optical bench connecting magnetometer and star imager by estimating Euler angles in 10 day segments and (v) we co-estimate degree-1 external fields separately for every 12 hr interval. The model provides a reliable representation of the static (core and crustal) field up to spherical harmonic degree n = 40, and of the first time derivative up to n = 15.

Impact of fabric, microcracks and stress field on shale anisotropy

Abstract: SUMMARY Few data are available on shales in terms of seismic to ultrasonic properties and anisotropy, although all are important with regards to imaging problems often encountered in such lithologies. Additionally, mechanisms causing changes in these properties are not well documented due to the fine grain size of such materials and time required for testing under controlled pore pressure conditions. The results presented here are derived from a set of experiments run on Muderong Shale with pore pressure control in order to evaluate the effect of stress magnitude and stress anisotropy on ultrasonic response. This shale was noted to have a linear velocitymean effective stress response and extremely high anisotropy, both likely the result of the presence of fluid-filled cracks in a low- permeability porous medium. Changes in velocity and Vp/Vs ratios are dependent on both stress and smectite content. S-wave velocity is significantly affected by the presence of smectite in this and other shales and at low stress (<20 MPa) is less sensitive to stress change than P-wave velocity. Vp/Vs ratios are noted to increase in this shale up to 20 MPa effective stress, then decrease slightly due to stress-induced loss of interlayer water in smectite. Intrinsic anisotropy comes from composition, a strong compaction fabric and the presence of microfractures; changes to ultrasonic anisotropy are the result of the magnitude of the stresses, their orientation with respect to the fractures and the degree of stress anisotropy.

Imaging subduction from the trench to 300 km depth beneath the central North Island, New Zealand, with Vp and Vp/Vs

Abstract: SUMMARY Recent dense deployments of portable digital seismographs have provided excellent control on earthquakes beneath the central North Island of New Zealand. Here we use a subset of the best-recorded earthquakes in an inversion for the 3-D Vp and Vp/Vs structure. The data set includes 39 123 P observations and 18 331 S observations from 1239 earthquakes and nine explosions. The subducted plate is imaged as a high Vp, low Vp/Vs feature. Vp within the mantle of the subducted slab is almost always >8.5 km s−1, which requires the ca. 120 Myr slab to be unusually cold. The low Vp/Vs within the subducted plate closely parallels the lower plane of the dipping seismic zone. It most likely indicates fluid resulting from dehydration of serpentine in the slab mantle, and the earthquakes themselves are likely to be promoted by dehydration embrittlement. We identify a region with Vp 8.0 km s−1 directly above the dipping seismic zone can be interpreted as sinking, entrained with the motion of the subducted slab and forming a viscous blanket that insulates the slab from the high-temperature mantle wedge. Material in the overlying low Vp region can be interpreted as rising within a return flow within the wedge. The volcanic front appears to be controlled by where this dipping low Vp region meets the base of the crust. The thickness of the backarc crust also shows significant variation along strike. In the central Taupo Volcanic Zone (TVZ) the crust is ca. 35 km thick, while southwest of Mt Ruapehu the crust thickens by ca. 10 km. There is no significant low Vp zone in the mantle wedge in this southwestern region, suggesting that this thicker crust has choked off mantle return flow. The seismic tomography results, when combined with constraints on mantle flow from previous shear-wave splitting results, provide a plausible model for both the distribution of volcanism in the central North Island, and the exceptional magmatic productivity of the central TVZ.

Traction image method for irregular free surface boundaries in finite difference seismic wave simulation

Abstract: SUMMARY In this study, we propose a new numerical method, named as Traction Image method, to accurately and efficiently implement the traction-free boundary conditions in finite difference simulation in the presence of surface topography. In this algorithm, the computational domain is discretized by boundary-conforming grids, in which the irregular surface is transformed into a ‘flat’ surface in computational space. Thus, the artefact of staircase approximation to arbitrarily irregular surface can be avoided. Such boundary-conforming gridding is equivalent to a curvilinear coordinate system, in which the first-order partial differential velocity-stress equations are numerically updated by an optimized high-order non-staggered finite difference scheme, that is, DRP/opt MacCormack scheme. To satisfy the free surface boundary conditions, we extend the Stress Image method for planar surface to Traction Image method for arbitrarily irregular surface by antisymmetrically setting the values of normal traction on the grid points above the free surface. This Traction Image method can be efficiently implemented. To validate this new method, we perform numerical tests to several complex models by comparing our results with those computed by other independent accurate methods. Although some of the testing examples have extremely sloped topography, all tested results show an excellent agreement between our results and those from the reference solutions, confirming the validity of our method for modelling seismic waves in the heterogeneous media with arbitrary shape topography. Numerical tests also demonstrate the efficiency of this method. We find about 10 grid points per shortest wavelength is enough to maintain the global accuracy of the simulation. Although the current study is for 2-D P-SV problem, it can be easily extended to 3-D problem.

Lithospheric structure of the Arabian and Eurasian collision zone in eastern Turkey from S-wave receiver functions

Abstract: SUMMARY Crustal and upper-mantle seismic discontinuities beneath eastern Turkey are imaged using teleseismic S-to-P converted phases. Three crustal phases are observed: the Moho with depth ranging between 30 and 55 km, indicating variable tectonic regimes within this continental collision zone; an upper-crustal discontinuity at approximately 10 km depth; and various crustal low-velocity zones, possibly associated with recent Quaternary volcanism. Imaging of the upper mantle is complicated by the 3-D geometry of the region, in particular due to the Bitlis‐Zagros suture zone. However, several upper-mantle S-to-P converted phase are identified as being the signature of the lithosphere‐asthenosphere boundary (LAB). The inferred LAB for the Eastern Anatolian Accretionary Complex indicates that eastern Turkey has an anomalously thin (between ∼60 and 80 km) lithosphere which is consistent with an oceanic slab detachment model. The observed LAB phases for the Arabian shield and Iranian plateau indicate that lithospheric thickness for these stable regions is on the order of 100 to 125 km thick, which is typical of continental margins.

Regional tomographic inversion of the amplitude and phase of Rayleigh waves with 2-D sensitivity kernels

Abstract: SUMMARY In this study, we test the adequacy of 2-D sensitivity kernels for fundamental-mode Rayleigh waves based on the single-scattering (Born) approximation to account for the effects of heterogeneous structure on the wavefield in a regional surface wave study. The calculated phase and amplitude data using the 2-D sensitivity kernels are compared to phase and amplitude data obtained from seismic waveforms synthesized by the pseudo-spectral method for plane Rayleigh waves propagating through heterogeneous structure. We find that the kernels can accurately predict the perturbation of the wavefield even when the size of anomaly is larger than one wavelength. The only exception is a systematic bias in the amplitude within the anomaly itself due to a site response. An inversion method of surface wave tomography based on the sensitivity kernels is developed and applied to synthesized data obtained from a numerical simulation modelling Rayleigh wave propagation over checkerboard structure. By comparing recovered images to input structure, we illustrate that the method can almost completely recover anomalies within an array of stations when the size of the anomalies is larger than or close to one wavelength of the surface waves. Surface wave amplitude contains important information about Earth structure and should be inverted together with phase data in surface wave tomography.

Difference in the GPS deformation pattern of North and Central Zagros (Iran)

Abstract: SUMMARY Measurements on either side of the Kazerun fault system in the Zagros Mountain Belt, Iran, show that the accommodation of the convergence of the Arabian and Eurasian Plates differs across the region. In northwest Zagros, the deformation is partitioned as 3‐6 mm yr −1 of shortening perpendicular to the axis of the mountain belt, and 4‐6 mm yr −1 of dextral strikeslip motion on northwest‐southeast trending faults. No individual strike-slip fault seems to slip at a rate higher than ∼ 2m m yr −1 . In southeast Zagros, the deformation is pure shortening of 8 ± 2m m yr −1 occurring perpendicular to the simple folded belt and restricted to the Persian Gulf shore. The fact that most of the deformation is located in front of the simple folded belt, close to the Persian Gulf, while seismicity is more widely spread across the mountain belt, confirms the decoupling of the surface sedimentary layers from the seismogenic basement. A comparison with the folding and topography corroborates a southwestward propagation of the surface deformation. The difference in deformation between the two regions suggests that right-lateral shear cumulates on the north‐south trending Kazerun strike-slip fault system to 6 ± 2m m yr −1 .

Seismological evidence for crustal‐scale thrusting in the Zagros mountain belt (Iran)

Abstract: SUMMARY Crustal receiver functions (RFs) computed from the records of 45 temporary seismological stations installed on a 620-km-long profile across central Zagros provide the first direct evidence for crustal thickening in this mountain belt. Due to a rather short 14 km average station spacing, the migrated section computed from radial RFs displays the Moho depth variations across the belt with good spatial resolution. From the coast of the Persian Gulf to 25 km southwest of the Main Zagros Thrust (MZT), the Moho is almost horizontal with slight depth variations around 45 km. Crustal thickness then increases abruptly to a maximum of ∼70 km beneath the Sanandaj‐Sirjan metamorphic zone, between 50 and 90 km northeast of the surface exposure of the MZT. Further northeast, the Moho depth decreases to ∼42 km beneath the UrumiehDokhtar magmatic assemblage and the southern part of the Central Iranian microcontinent. The region of thickest crust is located ∼75 km to the northeast of the Bouguer anomaly low at −220 mGals. Gravity modelling shows that the measured Moho depth variations can be reconciled with gravity observations by assuming that the crust of Zagros underthrusts the crust of central Iran along the MZT considered as a crustal-scale structure. This hypothesis is compatible with shortening estimates by balanced cross-sections of the Zagros folded belt, as well as with structural and petrological studies of the metamorphic Sanandaj‐Sirjan zone.

The S receiver functions: synthetics and data example

Abstract: SUMMARY Recently, the S receiver function method has been successfully developed to identify upper mantle interfaces. S receiver functions have the advantage of being free of S-wave multiple reflections and can be more suitable than P receiver functions for studying mantle lithosphere. However, because of specific ray geometry and interference of diverse phases, the S receiver function method has some technical difficulties and limitations. We use synthetic seismograms to demonstrate the feasibility and limitations of S receiver functions for studying mantle structures. Full-wavefield seismograms were calculated using the reflectivity method and processed to generate synthetic S receiver functions for S, SKS and ScS waves. Results show that S receiver functions can be obtained from waveforms of S, SKS and ScS waves. The synthetic S receiver functions for these incident waves show S-to-P converted phases at all discontinuities in the crust and upper mantle. Useful ranges of epicentral distances for calculation of S receiver functions are: 55 ◦ ‐85 ◦ for S, >85 ◦ for SKS and 50 ◦ ‐75 ◦ for ScS waves. We apply both the S and P receiver function methods to data recorded at broadband station YKW3 in Northwest Canada. The study shows that there is significant agreement among different receiver function methods, and demonstrates the usefulness of S receiver functions for imaging the mantle lithosphere.

Post-seismic relaxation following the great 2004 Sumatra-Andaman earthquake on a compressible self-gravitating Earth

Abstract: SUMMARY The Mwγ 9.0 2004 December 26 Sumatra-Andaman and Mw= 8.7 2005 March 28 Nias earthquakes, which collectively ruptured approximately 1800 km of the Andaman and Sunda subduction zones, are expected to be followed by vigorous viscoelastic relaxation involving both the upper and lower mantle. Because of these large spatial dimensions it is desirable to fully account for gravitational coupling effects in the relaxation process. We present a stable method of computing relaxation of a spherically-stratified, compressible and self-gravitating viscoelastic Earth following an impulsive moment release event. The solution is cast in terms of a spherical harmonic expansion of viscoelastic normal modes. For simple layered viscoelastic models, which include a low-viscosity oceanic asthenosphere, we predict substantial post-seismic effects over a region several 100s of km wide surrounding the eastern Indian Ocean. We compare observed GPS time-series from ten regional sites (mostly in Thailand and Indonesia), beginning in 2004 December, with synthetic time-series that include the coseismic and post-seismic effects of the 2004 December 26 and 2005 March 28 earthquakes. A viscosity structure involving a biviscous (Burgers body) rheology in the asthenosphere explains the pattern and amplitude of post-seismic offsets remarkably well.

The relationship between depth, age and gravity in the oceans

Abstract: We reassess the applicability of the thermal plate cooling model to the subsidence of the North Pacific, Atlantic and North Indian Ocean Basins. We use a new numerical plate model in which the thermophysical parameters of the lithosphere vary with temperature according to the results of laboratory experiments, and the ridge temperature structure is consistent with the thickness of the oceanic crust. We first attempt to exclude thickened crust from our data set, and then to exclude swells and downwellings by masking regions of the data that remains that have significant gravity anomalies when there exists a clear regional correlation between intermediate-wavelength gravity and topography. We find that the average variation of depth with age is consistent with conventional half-space models until about 90 Myr. Thereafter, the departure from the half-space cooling curve is more rapid than predicted using simple conductive plate cooling models. The depth–age curves in the Pacific and Atlantic show ∼250 m of temporary shallowing between the ages of 90–130 Myr, a result consistent with the outcome of experiments on the initiation of small-scale boundary layer convection. The results do not change significantly if the estimated component of the gravity arising from plate cooling is subtracted prior to calculation of the correlation between gravity and topography. A 90-km-thick conductive plate is nevertheless a reasonable model for the average temperature structure of the oldest part of the Pacific ocean lithosphere. In the Pacific, the broad topographic undulations associated with the Line Island Swell, the Hawaiian Swell and surrounding basins have correlated gravity anomalies and an admittance of approximately 30 mGal km−1 and are likely to result from convective circulation in the upper mantle. In the Northeast Atlantic, the intermediate-wavelength admittance over the Cape Verde swell is similar; in the Northwest Atlantic over the Bermuda Swell it is slightly larger but not as well constrained.

Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology

Abstract: SUMMARY We perform analytical and numerical studies of aftershock sequences following abrupt steps of strain in a rheologically layered model of the lithosphere. The model consists of a weak sedimentary layer, over a seismogenic zone governed by a viscoelastic damage rheology, underlain by a viscoelastic upper mantle. The damage rheology accounts for fundamental irreversible aspects of brittle rock deformation and is constrained by laboratory data of fracture and friction experiments. A 1-D version of the viscoelastic damage rheology leads to an exponential analytical solution for aftershock rates. The corresponding solution for a 3-D volume is expected to be sum of exponentials. The exponential solution depends primarily on a material parameter R given by the ratio of timescale for damage increase to timescale for accumulation of gradual inelastic deformation, and to a lesser extent on the initial damage and a threshold strain state for material degradation. The parameter R is also inversely proportional to the degree of seismic coupling across the fault. Simplifying the governing equations leads to a solution following the modified Omori power-law decay with an analytical exponent p = 1. In addition, the results associated with the general exponential expression can be fitted for various values of R with the modified Omori law. The same holds for the decay rates of aftershocks simulated numerically using the 3-D layered lithospheric model. The results indicate that low R values (e.g. R ≤ 1) corresponding to cold brittle material produce long Omori-type aftershock sequences with high event productivity, while high R values (e.g. R ≥ 5) corresponding to hot viscous material produce short diffuse response with low event productivity. The frequency-size statistics of aftershocks simulated in 3-D cases with low R values follow the Gutenberg–Richter power law relation, while events simulated for high R values are concentrated in a narrow magnitude range. Increasing thickness of the weak sedimentary cover produces results that are similar to those associated with higher R values. Increasing the assumed geothermal gradient reduces the depth extent of the simulated earthquakes. The magnitude of the largest simulated after

Lithospheric and upper mantle structure of central Chile and Argentina

Abstract: SUMMARY The tectonics of central Chile and Argentina have been greatly affected by the shallow dips of the subducting Nazca plate, which controlled patterns of magmatism and deformation nearly 1000 km away from the plate boundary. We calculate receiver functions from data recorded by the CHARGE array, which transected the Andes and Sierras Pampeanas in central Chile and Argentina, to better constrain the crustal structure of this region. Beneath the northern transect of the CHARGE array, where the Nazca slab flattens near 100 km, we find the crust is over 60 km thick beneath the Andes and thin to the east. The thick crust, however, extends ∼200 km to the east of the high elevations. Estimates of VP/VS obtained from receiver functions vary along ancient terrane boundaries exhibiting higher values to the west. Interestingly, we observe that the amplitude of the phase corresponding to the Moho on receiver functions diminishes to the west, complicating our images of crustal structure. We proposed that the observations presented here of thickened crust within a region of low elevations, diminished receiver function arrivals, and reports of high shear-wave speeds atop of the mantle wedge overlying the shallowly subducted Nazca slab, can be explained by partial eclogitization of the lower crust. The Moho appears simpler across the southern transect where it can be identified near 50 km depth at its deepest point beneath the Andes and shallows eastwards. Volcanism remains active near the latitudes of our southern transect and we observe multiple crustal lowvelocity zones indicative of regions of partial melt near the centres of volcanism. Signals related to the Nazca slab remain more elusive, suggestive of a small impedance contrast between the slab and overlying mantle.

Improved strategies for the automatic selection of optimized sets of electrical resistivity tomography measurement configurations

Abstract: Two strategies are presented for obtaining the maximum spatial resolution in electrical resistivity tomography surveys using a limited number of four-electrode measurement configurations. Both methods use a linearized estimate of the model resolution matrix to assess the effects of including a given electrode configuration in the measurement set. The algorithms are described in detail, and their execution times are analysed in terms of the number of cells in the inverse model. One strategy directly compares the model resolution matrices to optimize the spatial resolution. The other uses approximations based on the distribution and linear independence of the Jacobian matrix elements. The first strategy produces results that are nearer to optimal, however the second is several orders of magnitude faster. Significantly however, both offer better optimization performance than a similar, previously published, method. Realistic examples are used to compare the results of each algorithm. Synthetic data are generated for each optimized set of electrodes using simple forward models containing resistive and/or conductive prisms. By inverting the data, it is demonstrated that the linearized model resolution matrix yields a good estimate of the actual resolution obtained in the inverted image. Furthermore, comparison of the inversion results confirms that the spatial distribution of the estimated model resolution is a reliable indicator of tomographic image quality.

Regression problems for magnitudes

Abstract: SUMMARY Least-squares linear regression is so popular that it is sometimes applied without checking whether its basic requirements are satisfied. In particular, in studying earthquake phenomena, the conditions (a) that the uncertainty on the independent variable is at least one order of magnitude smaller than the one on the dependent variable, (b) that both data and uncertainties are normally distributed and (c) that residuals are constant are at times disregarded. This may easily lead to wrong results. As an alternative to least squares, when the ratio between errors on the independent and the dependent variable can be estimated, orthogonal regression can be applied. We test the performance of orthogonal regression in its general form against Gaussian and non-Gaussian data and error distributions and compare it with standard leastsquare regression. General orthogonal regression is found to be superior or equal to the standard least squares in all the cases investigated and its use is recommended. We also compare the performance of orthogonal regression versus standard regression when, as often happens in the literature, the ratio between errors on the independent and the dependent variables cannot be estimated and is arbitrarily set to 1. We apply these results to magnitude scale conversion, which is a common problem in seismology, with important implications in seismic hazard evaluation, and analyse it through specific tests. Our analysis concludes that the commonly used standard regression may induce systematic errors in magnitude conversion as high as 0.3‐0.4, and, even more importantly, this can introduce apparent catalogue incompleteness, as well as a heavy bias in estimates of the slope of the frequency‐magnitude distributions. All this can be avoided by using the general orthogonal regression in magnitude conversions.

VS30 mapping and soil classification for seismic site effect evaluation in Dinar region, SW Turkey

Abstract: SUMMARY The Dinar earthquake (M S = 6.1) of 1995 October 1 killed 90 people and destroyed more than 4000 buildings. Despite the moderate size of the earthquake, the level of damage was extremely high, which led to many studies that were carried out in the region. The majority of these studies concluded that the main reasons for the damage were the construction errors and the poor soil conditions. However, at that time no appropriate soil condition map based on extended, high density measurements was available. Shear wave velocity is an important parameter for evaluating the dynamic behaviour of soil in the shallow subsurface. Thus site characterization in calculating seismic hazards is usually based on the near surface shear wave velocity values. The average shear wave velocity for the top 30 m of soil is referred to as V S 30 . For earthquake engineering design purposes, both the Uniform Building Code (UBC) and Eurocode 8 (EC8) codes use V S 30 to classify sites according to the soil type. The V s 30 values calculated by using multichannel analysis of surface waves (MASW) were used to create a new soil classification map of the Dinar region. Surface seismic measurements were carried out at 50 locations mostly in Dinar city and its surroundings. The dispersion data of the recorded Rayleigh waves were inverted using a Genetic Algorithm (GA) method to obtain shear wave velocity profiles of the investigated sites. Thus the derived V s 30 map of the Dinar region was transformed to the UBC and EC8 standards. Soil classification results show that most parts of the region, located in alluvial basin, have low shear wave velocity values. These values are within the range of 160–240 m s −1 and thus fall into the S D and S E categories according to the UBC and the C and D categories according to EC8. Within the region, some parts located on the hill zone and the transition zone have better soil conditions [corresponding to S C (UBC) and B (EC8) categories] and have comparatively high shear wave velocities in the range of 500–740 m s −1 and 350–450 m s −1 , respectively. V S 30 and soil classification maps were compared with the damage distribution associated with the earthquake. In possession of a detailed shear wave velocity map of Dinar City, in general, the results show that there is a correlation between the V S 30 values and the damage distribution of the region. In addition to the low V S 30 values, the likely causes of the damage were investigated, and it is observed that one of the major factors for high levels of damage is 3-D variations of geological structures.