Multiple scattering and energy transfer of seismic waves – separation of scattering effect from intrinsic attenuation – I. Theoretical modelling

Summary. In order to separate the scattering effect from the intrinsic attenuation, we need a multiple scattering model for seismic wave propagation in random heterogeneous media. In this paper, we apply radiative transfer theory to seismic wave propagation and formulate in the frequency domain the energy density distribution in space for a point source. We consider the cases of isotropic scattering and strong forward scattering. Some numerical examples are shown. It is seen that the energy density–distance curves have quite different shapes depending on the values of medium seismic albedo B ₀=η_s/(η_s+η_a) where η_s is the scattering coefficient and η_a is the absorption coefficient of the medium. For a high albedo ( B > 0.5) medium, the energy–distance curve is of arch shape and the position of the peak is a function of the extinction coefficient of the medium η_e=η_s+η_a. Therefore it is possible to separate the scattering effect and the absorption based on the measured energy density distribution curves.     

Focal mechanisms of small earthquakes in the southeastern Brazilian shield: a test of stress models of the South American plate

Focal mechanisms of small earthquakes with magnitudes of about 3 in the SE Brazilian shield are calculated using S / P amplitude ratios. Low attenuation ( Q _p from 400 to 800) in the shield upper-crustal layers allowed sharp S arrivals to be recorded up to distances of 100 km. Besides P -wave polarities, SH -wave first motions were also used to constrain the nodal-plane orientations. Normal and reverse faulting mechanisms with strike-slip components were found. The inversion of four mechanisms to estimate the stress tensor indicated a strike-slip stress regime with roughly E–W-orientated σ ₁ and N–S σ ₃. Both the orientations and the shape factor ( φ =0.7) of the inverted stress are in excellent agreement with theoretical predictions for that part of Brazil from the driving-force model of Coblentz & Richardson (1996) . Good agreement with the nature of the stress, as well as its orientation, was also found for the model of Meijer (1995) . Both of these theoretical models include spreading stresses along the continent/ocean lithospheric transition. Because the earthquakes are more than 300 km from the continental shelf they should not be affected by the local flexural forces caused by sediment load in the marginal basins. The agreement between observed and theoretical stresses then confirms the importance of continental spreading forces in modelling intraplate stresses.     

Magnetotelluric response of a uniformly stratified earth containing a magnetized layer

Summary The magnetotelluric (MT) response is studied of a uniformly stratified earth which contains a magnetized layer. The impedance as a function of the layer parameters (resistivity, ρ permeability, μ and thickness h ) is discussed. The MT response from a layer (μ, ρ, h ) is equivalent to that from a layer (μ/μ_r=μ₀, μ_rρ, μ_r h ) where μ_r is the relative permeability of the layer. Thus the effect of a magnetized layer is to make it apppear μ_r times more resistive and μ_r times thicker than an unmagnetized equivalent layer. Master curves of apparent resistivity and phase are computed for three-layer models with varying permeability associated with varying resistivity in each layer. An example of MT field data is presented in which the most reasonable interpretation is that a magnetized layer exists beneath the observatory site.     

Calculating the two-dimensional magnetotelluric Jacobian in finite elements using reciprocity

To speed up the calculation of the field Jacobian for 2-D magnetoteliuric inversion using finite elements, the principle of electromagnetic reciprocity is applied. The governing relationship for the Jacobian of the field along strike is obtained by differentiating the Helmholtz equation with respect to the resistivity of each region in the finite-element mesh. The result is a similar Helmholtz equation for the Jacobian, with new sources distributed over all nodes within the parameter medium. However, according to the principle of electromagnetic reciprocity, the roles of sources and receivers are interchangeable. Utilizing reciprocity, the field values obtained from the original forward problem and for new unit sources imposed at the receivers are then utilized in the calculation of the Jacobian by simple multiplication and summation with finite-element terms at each rectangle in the mesh. For the auxiliary (across-strike) fields, the Jacobian terms are obtained by solving source vectors loaded with parabola coefficients used in the approximation to Maxwell's equations. Jacobian terms for the apparent resistivity ( p _a), the impedance phase (φ) and the vertical magnetic field ( K _zy) are then calculated utilizing the parallel- and auxiliary-field Jacobians. Comparison of Jacobian values obtained from reciprocity calculations and by differencing two forward solutions show that the reciprocity method is accurate and can be used to decrease the number of calculations required to obtain sensitivities by one to two orders of magnitude.     

The solution of the inverse gravity potential problem for cylindrical bodies

Summary. The inverse gravity potential problem consists in the determination of the form and the density of the body by its exterior gravity potential. We describe two similar classes of bodies for which this problem has a unique constructive solution.
(1) The first class contains the cylindrical bodies with finite length, arbitrary form of section and ρ( R , ø, z) =ρ₁( z )ρ₂( R , ø) density distribution, where z is the cylindrical coordinate; R , ø are the polar coordinates in a section plane. This class is important for prospecting geophysics in that it allows us to determine in a unique and constructive way, the function ρ₁( R , ø), the length, form and orientation of the cylinder if we know the function ρ₁( z ) and the exterior potential. The classical moment problem of functions is the basis for the solution of this problem.
(2) The analogous problem for the class of the spherical cylinders, or bodies bounded by arbitrary similar sections of two different concentric spheres and the radial lateral surface, appears when bodies of planetary size are studied. (An example of these bodies would be the Moon mascons.) The density distribution of these cylinders is ρ(τ, θ, ø) =ρ₁(τ)ρ₂(θ, ø) where τ, θ, ø are the spherical coordinates. The function ρ₁(θ, ø), length and form of spherical sections can be uniquely determined by exterior potential if we know the function ρ₁(τ). We propose a new constructive method for harmonic continuation of the gravity potential into the region containing the perturbing masses for the solution of the problem.     

A pseudo-Thellier relative palaeointensity record, and rock magnetic and geochemical parameters in relation to climate during the last 276 kyr in the Azores region

In the pseudo-Thellier method for relative palaeointensity determinations (Tauxe et al. 1995) the slope of the NRM intensity left after AF demagnetization versus ARM intensity gained at the same peak field is used as a palaeointensity measure. We tested this method on a marine core from the Azores, spanning the last 276  kyr. We compared the pseudo-Thellier palaeointensity record with the conventional record obtained earlier by Lehman et al . (1996 ), who normalized NRM by SIRM. The two records show similar features: intensity lows with deviating palaeomagnetic directions at 40–45  ka and at 180–190  ka. The first interval is associated with the Laschamps excursion, while the 180–190  ka low represents the Iceland Basin excursion (Channell et al. 1997). The pseudo-Thellier method, in combination with a jackknife resampling scheme, provides error estimates on the palaeointensity.
  Spectral analysis of the rock magnetic parameters and the palaeointensity estimates shows orbitally forced periods, particularly 23  kyr for climatic precession. This suggests that palaeointensity is still slightly contaminated by climate. Fuzzy c -means cluster analysis of rock magnetic and geochemical parameters yields a seven-cluster model of predominantly calcareous clusters and detrital clusters. The clusters show a strong correlation with climate, for example samples from detrital clusters predominantly appear during rapid warming. Although both the pseudo-Thellier palaeointensity m _a and fuzzy clusters show climatic influences, we have not been able to find an unambiguous connection between the clusters and m _a .     

On the connection between the critical value of the "vertical" electrical anisotropy and earthquake prediction

Summary. It is shown that at any orientation of a dipole of direct current at the surface of a uniform anisotropic halfspace there exists a unique critical value η_cr of the "vertical" electric anisotropy η. For η > η_cr not all field lines are closed, and three straight field lines pass through the centre of the dipole. Also there exists a unique value η_eq > η_cr, where equatorial η is equal to zero. If η > η_eq then equatorial ρ_app is negative since the electric vector is turned more than through 90° at the equator of dipole, and it is shown that η_cr= 1.79, η_eq= 1.93 if the axis of dipole is directed symmetrically with respect to the axes of anisotropy. Routine observations of variations of ρ_app with the aim of earthquake prediction should be carried out taking into account the possibility of great fluctuations on η far from the epicentre of an earthquake.     

A normal approximation to the Fisher distribution

Summary. A simple normal approximation is given for the joint probability density function of the polar coordinates (θ, φ) of a random vector following the Fisher distribution with arbitrary mean direction (θ₀, φ₀). Two important applications of this approximation are discussed, and conditions for the validity of the approximation are investigated and summarized.     

Transient and impulse responses of a one-dimensional linearly attenuating medium—II. A parametric study

Summary. We investigate one-dimensional waves in a standard linear solid for geophysically relevant ranges of the parameters. The critical parameters are shown to be T*= t_u/Q_m where t _u is the travel time and Q_m the quality factor in the absorption band, and τ⁻¹ _m , the high-frequency cut-off of the relaxation spectrum. The visual onset time, rise time, peak time, and peak amplitude are studied as functions of T* and τ _m. For very small τ _m , this model is shown to be very similar to previously proposed attenuation models. As τ _m grows past a critical value which depends on T* , the character of the attenuated pulse changes. Seismological implications of this model may be inferred by comparing body wave travel times with a'one second'earth model derived from long-period observations and corrected for attenuation effects assuming a frequency independent Q over the seismic band. From such a comparison we speculate that there may be a gap in the relaxation spectrum of the Earth's mantle for relaxation times shorter than about one second. However, observational constraints from the attenuation of body waves suggest that such a gap might in fact occur at higher frequencies. Such a hypothesis would imply a frequency dependence of Q in the Earth's mantle for short-period body waves.     

Non-linear altimetric geoid inversion for lithospheric elastic thickness and crustal density

The inversion of high-resolution geoid anomaly maps derived from satellite altimetry should allow one to retrieve the lithospheric elastic thickness, T _e , and crustal density, _c . Indeed, the bending of a lithospheric plate under the load of a seamount depends on both parameters, and the associated geoid anomaly is correspondingly dependent on the two parameters. The difference between the observed and modelled geoid signatures is estimated by a cost function, J , of the two variables, T _e and _c . We show that this cost function forms a valley structure along which many local minima appear, the global minimum of J corresponding to the true values of the lithospheric parameters. Classical gradient methods fail to find this global minimum because they converge to the first local minimum of J encountered, so that the final parameter estimate strongly depends on the starting pair of values ( T _e ,   _c ). We here implement a non-linear optimization algorithm to recover these two parameters from altimetry data. We demonstrate from the inversion of synthetic data that this approach ensures robust estimates of T _e and _c by activating two search phases alternately: a gradient phase to find a local minimum of J , and a tunnelling phase through high values of the cost function. The accuracy of the solution can be improved by a search in an iteratively restricted parameter subspace. Applying our non-linear inversion to the Great Meteor Seamount geoid data, we further show that the inverse problem is intrinsically ill-posed. As a consequence, minute geoid (or gravity) data errors can induce large changes in any recovery of lithospheric elastic thickness and crustal density.     

Focal mechanisms of recent earthquakes in the Southern Korean Peninsula

We evaluate the stress field in and around the southern Korean Peninsula with focal mechanism solutions, using the data collected from 71 earthquakes ( M_L = 1.9–5.2) between 1999 and 2004. For this, the hypocentres were relocated and well-constrained fault plane solutions were obtained from the data set of 1270 clear P -wave polarities and 46 SH / P amplitude ratios. The focal mechanism solutions indicate that the prevailing faulting types in South Korea are strike-slip-dominant-oblique-slip faultings with minor reverse-slip component. The maximum principal stresses (σ₁) estimated from fault-slip inversion analysis of the focal mechanism solutions show a similar orientation with E–W trend (269°–275°) and low-angle plunge (10°–25°) for all tectonic provinces in South Korea, consistent with the E–W trending maximum horizontal stress (σ_Hmax) of the Amurian microplate reported from in situ stress measurements and earthquake focal mechanisms. The directions of the intermediate (σ₂) and minimum (σ₃) principal stresses of the Gyeongsang Basin are, however, about 90 deg off from those of the other tectonic provinces on a common σ₂–σ₃ plane, suggesting a permutation of σ₂ and σ₃. Our results incorporated with those from the kinematic studies of the Quaternary faults imply that NNW- to NE-striking faults (dextral strike-slip or oblique-slip with a reverse-slip component) are highly likely to generate earthquakes in South Korea.     

Source tomography by simulated annealing using broad-band surface waves and geodetic data: application to the Mw= 8.1 Chile 1995 event

We invert surface-wave and geodetic data for the spatio-temporal complexity of slip during the M _w =8.1 Chile 1995 event by simulated annealing. This quasi-global inversion method allows for a wide exploration of model space, and retains the non-linearity of the source tomography problem. Complex source spectra are obtained from 5 to 45 mHz from first- and second-orbit fundamental-mode Rayleigh waves using an empirical Green's function cross-correlation technique. Coseismic displacement vectors were measured at 10 GPS sites near Antofagasta. They are part of a French-Chilean experiment which monitors the Northern Chile seismic gap. The spectra, together with the geodetic data, are inverted for the moment distribution on a 2-D dipping fault, under the physical constraints of slip positivity and causality. Marginal a posteriori distributions of the model parameters are obtained from several independently inverted solutions. In general, features of the slip model are well resolved. Data are well fitted by a purely unilateral southward rupture with a nearly uniform velocity around 2.5–3.0 km s⁻¹, and a total duration of 65 s. Several regions of moment release were imaged, one near the hypocentre, a major one 80 km south of it and a minor one 160 km south of it. The major patch of moment release seemed to have propagated to relatively shallow depths near the trench, 100 km SSW of the epicentre. The region of major slip is located updip of the 1987, M _w =7.5 earthquake, suggesting a causal relationship. Most of the slip occurred updip of the hypocentre (36 km), but the entire coupled plate interface (20–40 km) ruptured during the Chile 1995 event.     

Palaeomagnetism of Lower Cambrian sediments from the Olenek River section (northern Siberia): palaeopoles and the problem of magnetic polarity in the Early Cambrian

New palaeomagnetic data from the Lower and Middle Cambrian sedimentary rocks of northern Siberia are presented. During stepwise thermal demagnetization the stable characteristic remanence (ChRM) directions have been isolated for three Cambrian formations. Both polarities have been observed, and mean ChRM directions (for normal polarity) are: Kessyusa Formation (Lower Cambrian) D = 145°, I = -40°, N = 12, α₉₅= 12.8°; pole position: φ= 38°S, A = 165°E; Erkeket Formation (Lower Cambrian, stratigraphically highly) D = 152°, I = - 47°, N = 23, α₉₅= 6.8°; pole position: φ= 45°S, A = 159°E; Yunkyulyabit-Yuryakh Formation (Middle Cambrian) D = 166°, I = - 33°, N = 38, α₉₅= 4.6°; pole position: φ= 36°S, L = 140°E. These poles are in good agreement with the apparent polar wander path based on the bulk of existing Cambrian palaeomagnetic data from the Siberian platform. In Cambrian times, the Siberian platform probably occupied southerly latitudes stretching from about 35° to 0°, and was oriented 'reversely' with respect to its present position. Siberia moved northwards during the Cambrian by about 10° of latitude. This movement was accompanied by anticlockwise rotation of about 30°. The magnetostratigraphic results show the predominance of reversed polarity in the Early Cambrian and an approximately equal occurrence of both polarities in the part of the Middle Cambrian studied. These results are in good agreement with the palaeomagnetic polarity timescale for the Cambrian of the Siberian platform constructed previously by Khramov et al. (1987).     

Long-period geomagnetic variations and mantle conductivity: an inversion using Bailey's method

Summary. We have implemented an algorithm which is based on Bailey's solution of the inverse problem of electromagnetic induction in the Earth. The study was motivated by recent determinations of very long period data and also benefited from recent redeterminations of high frequency data. The algorithm has been successfully tested to provide reliable estimates of conductivity down to a depth of 2000 km, using synthetic data in the period range from 4 days to 11 years. Smooth data sets, which are required for the inversion, were constructed from various sources. At a given depth, the range of inverted models is less than one order of magnitude. Due to the lack of high frequency data, the conductivity of the upper 600 km of the mantle, which is found to be of the order of 10⁻¹Ω⁻¹ m⁻¹, may be overestimated. The algorithm performs well in the middle mantle, where conductivity rises steadily from 1 to 50 Ω⁻¹ m⁻¹. The lack of very low frequency data and limitations of the algorithm prevent one from obtaining meaningful estimates in the lower mantle. However, the study of the propagation of the late 1960s secular variation acceleration provides an estimate of the mean conductivity of the whole mantle. Thus, a complete mantle profile can be constructed. It is found that deep mantle conductivity probability does not exceed a few hundred Ω⁻¹ m⁻¹.     

Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature

We present a regional surface waveform tomography of the Pacific upper mantle, obtained using an automated multimode surface waveform inversion technique on fundamental and higher mode Rayleigh waves, to constrain the   V_SV   structure down to ∼400 km depth. We have improved on previous implementations of this technique by robustly accounting for the effects of uncertainties in earthquake source parameters in the tomographic inversion. We have furthermore improved path coverage in the South Pacific region by including Rayleigh wave observations from the French Polynesian Pacific Lithosphere and Upper Mantle Experiment deployment. This improvement has led to imaging of vertical low-velocity structures associated with hotspots within the South Pacific Super-Swell region. We have produced an age-dependent average cross-section for the Pacific Ocean lithosphere and found that the increase in   V_SV   with age is broadly compatible with a half-space cooling model of oceanic lithosphere formation. We cannot confirm evidence for a Pacific-wide reheating event. Our synthetic tests show that detailed interpretation of average   V_SV   trends across the Pacific Ocean may be misleading unless lateral resolution and amplitude recovery are uniform across the region, a condition that is difficult to achieve in such a large oceanic basin with current seismic stations.     

Point sources of elastic deformation: elementary sources, dynamic displacements

We construct a catalogue of all the possible elementary point sources of seismic waves. There are three general classes of sources, two spheroidal and one toroidal. We consider excitation functions for these point-like sources as well as for sources of finite size in far-, intermediate- and near-field for an infinite homogeneous isotropic medium. The sources corresponding to seismic-moment tensors for the second-, third- and fourth-ranks are considered in more detail; we identify 10 different seismic sources in this range: one monopole, two or three dipoles, three quadrupoles, etc. For the step-function of the scalar seismic-moment release, the amplitude spectrum for the third-rank sources is proportional to the angular frequency ω in the region below the corner frequency ω _cr. The fourth-rank sources have an ω ² spectrum in the same range. The possibility of separate and simultaneous inversion of seismic body-wave data and static deformation data for sources of different order is discussed. Some equivalent-force moment higher-rank sources are 'shielded' by lower-rank sources of the same order; the former sources cannot be inverted from seismic data without additional assumptions. Because of their simple radiation pattern, the lower order multipoles, i.e. the monopole and dipoles, are the first sources other than the double-couple which should be considered for inversion.     

A three-dimensional radially anisotropic model of shear velocity in the whole mantle

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   V_SH > V_SV   starting at ∼80 km under oceanic regions and ∼200 km under stable continental lithosphere, suggesting horizontal flow beneath the lithosphere. We also observe a   V_SV > V_SH   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   V_SV > V_SH   , while the ridge signal decreases. While the mid-mantle has lower amplitude anisotropy (<1 per cent), we also confirm the observation of radially symmetric   V_SH > V_SV   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 under the central Pacific and Africa, suggesting that   V_SH > V_SV   is generated in the predominant horizontal flow of a mechanical boundary layer, with a change in signature related to transition to upwelling at the superplumes.     

The deformation and compaction of partial molten zones

Summary. The segregation of melt from a partially molten source region requires a corresponding deformation of the unmelted residue ('matrix'). The role of matrix deformation during melt segregation is examined using simple one-dimensional models, for which the deformation consists only of bulk compression or 'compaction'. In model I, a volume fraction φ₀ of ascending mantle material undergoes pressure-release melting at a depth z = 0 (localized melting). Compaction of the matrix occurs in a boundary layer whose thickness (reduced compaction length δ_R) is proportional to the square root of the matrix viscosity. In the Earth's mantle, δ_R∼ 10–100 m, indicating that compaction cannot be important over large distances. Model II examines the case in which melting occurs over a depth range of order h (distributed melting). In the limit h ≪δ_R, the solution is the same as for the case of localized melting, except in a 'melting layer' of thickness ∼ h near z = 0. In the more realistic limit h ≫δ_R, compaction makes a negligible contribution to the balance of forces associated with melt segregation. This result is also valid for the more general case of two-dimensional flow. Compaction is therefore likely to be of negligible importance in the Earth's mantle, with the consequence that melt segregation can be accurately described by Darcy's law.     

Inversion of Love-wave data for velocity and anelasticity using exact kernels

Summary. The general problem of inverting Love-wave dispersion and amplitude data to obtain a velocity and Q_s structure is considered. A formulation is used which incorporates attenuation into the Haskell-Thompson matrix method in an exact manner and thus retains the inherent non-linearity in the anelasticity. The resulting exact inversion kernels allow simultaneous inversion for velocity and intrinsic attenuation parameters. The method is applied to synthetic data which allows a comparison to be made with inexact kernels. The results indicate that the use of inexact kernels may introduce spurious oscillations into the Q_s structure and that a simultaneous inversion can be more stable than inverting for velocity alone.     

Inversion of sonic waveforms with unknown source and receiver functions

A seismogram that is several times the length of the source-receiver wavelet is windowed into two parts—these may overlap—to obtain two seismograms with approximately the same source function but different Green's functions. A similarly windowed synthetic seismogram gives two corresponding synthetic seismograms. The spectral product of the window 1 data with the window 2 synthetic is equal to the spectral product of the window 1 synthetic with the window 2 data only if the correct earth model is used to compute the synthetic. This partition principle is applied to well-log sonic waveform data from Ocean Drilling Project hole 806B, a slow formation, and used there to estimate Poisson's ratio from a single seismogram whose transmitter and receiver functions are unknown. A multichannel extension of the algorithm gives even better results. The effective borehole radius R _b, was included in the inversion procedure, because of waveform sensitivity to R _b. Inversion results for R _b agreed with the sonic caliper, but not the mechanical caliper; thus if R _b is not included in the inversion its value should be taken from the sonic caliper.