Three-dimensional massively parallel electromagnetic inversion—I. Theory

An iterative solution to the non-linear 3-D electromagnetic inverse problem is obtained by successive linearized model updates using the method of conjugate gradients. Full wave equation modelling for controlled sources is employed to compute model sensitivities and predicted data in the frequency domain with an efficient 3-D finite-difference algorithm. Necessity dictates that the inverse be underdetermined, since realistic reconstructions require the solution for tens of thousands of parameters. In addition, large-scale 3-D forward modelling is required and this can easily involve the solution of over several million electric field unknowns per solve. A massively parallel computing platform has therefore been utilized to obtain reasonable execution times, and results are given for the 1840-node Intel Paragon. The solution is demonstrated with a synthetic example with added Gaussian noise, where the data were produced from an integral equation forward-modelling code, and is different from the finite difference code embedded in the inversion algorithm     

Two-dimensional non-linear inversion of VLF-R data using simulated annealing

The VLF-R (very low frequency-resistivity) data, i.e. the apparent resistivity ( ρ _a ) and phase ( φ ) data, were inverted individually and jointly using the VFSA (very fast simulated annealing) global inversion approach. Global inversion results for synthetic data without and with various amounts of random and normally distributed Gaussian noise reveal that the inversion of neither the ρ _a nor φ data alone yields the true parameters of the structures. However, the joint inversion of the ρ _a and φ data yields very good estimates of the model parameters. Five models, representing typical subsurface structures in the shield areas, are studied here. Various models achieved after 10 VFSA runs were used to compute the mean model and the corresponding covariance and correlation matrices, which were used to estimate the uncertainties in the mean model parameters and correlations between the model parameters. We observe that these correlations follow the physics associated with the problem. VLF-R field data due to a nearly vertical contact structure and a very thick dyke-like structure were also inverted to demonstrate the efficacy of the approach in the delineation of the parameters of 2-D structures.     

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.     

Comparison of global synthetic seismograms calculated using the spherical 2.5-D finite-difference method with observed long-period waveforms including data from the intra-Antarctic region

We have been developing an accurate and efficient numerical scheme, which uses the finite-difference method (FDM) in spherical coordinates, for the computation of global seismic wave propagation through laterally heterogeneous realistic Earth models. In the field of global seismology, traditional axisymmetric modeling has been used widely as an efficient approach since it can solve the 3-D elastodynamic equation in spherical coordinates on a 2-D cross-section of the Earth, assuming structures to be invariant with respect to the axis through the seismic source. However, it has the severe disadvantages that asymmetric structures about the axis cannot be incorporated and the source mechanisms with arbitrary shear dislocation have not been attempted for a long time. Our scheme is based on the framework of axisymmetric modeling but has been extended to treat asymmetric structures, arbitrary moment-tensor point sources, anelastic attenuation, and the Earth center which is a singularity of wave equations in spherical coordinates. All these types of schemes which solve 3-D wavefields on a 2-D model cross-section are classified as 2.5-D modeling, so we have named our scheme the spherical 2.5-D FDM. In this study, we compare synthetic seismograms calculated using our FDM scheme with three-component observed long-period seismograms including data from stations newly installed in Antarctica in conjunction with the International Polar Year (IPY) 2007–2008. Seismic data from inland Antarctica are expected to reveal images of the Earth's deep interior with enhanced resolution because of the high signal-to-noise ratio and wide extent of this region, in addition to the rarity of sampling paths along the rotation axis of the Earth. We calculate synthetic seismograms through the preliminary reference earth model (PREM) including attenuation using a moment-tensor point source for the November 9, 2009 Fiji earthquake. Our results show quite good agreement between synthetic and observed seismograms, which indicates the accuracy of observations in the Antarctica, as well as the feasibility of the spherical 2.5-D modeling scheme.     

Modeling synthetic aperture radar (SAR) scattering from a seasonally varying snow-covered sea ice volume at 5.3 and 9.25 GHz

A series of sensitivity analyses using dielectric, mixture and microwave scattering models is presented. Data from the Seasonal Sea Ice Monitoring and Modeling Site (SIMMS) in 1990 and 1991 are used to initialize the models. The objective of the research is to investigate the role of various geophysical and electrical properties in specifying the total relative scattering cross section (ρ') of snow covered first-year sea ice during the spring period.
The seasonal transition period from the Winter SAR scattering season to Early Melt was shown to signal a transition in dielectric properties which caused the snow volume to become a factor in the microwave scattering process. The effect of the thermal insulation of a snow cover on sea ice was shown to be significant for both ε' and ε'. Higher atmospheric temperatures caused proportionally greater changes in the dielectric properties of the sea ice at the base of the snow cover. Model ρ⁰ was computed for a range of sensor, sensor-earth geometry, and geophysical properties. In the Winter season the surface roughness terms (o_hand L) were shown to have a significant impact on ρ⁰ when the ice surface was the primary scattering mechanism. Once the snow cover began to warm and water was available in a liquid phase, the ice surface became masked because of the decrease in microwave penetration depths. During this period the water volume variable dominated ρ⁰, both from its impact on ρ_v⁰, and due to its control over the dielectric mismatch created at the air/snow interface.     

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.     

Finite-difference modelling of magnetotelluric fields in two-dimensional anisotropic media

An algorithm for the numerical modelling of magnetotelluric fields in 2-D generally anisotropic block structures is presented. Electrical properties of the individual homogeneous blocks are described by an arbitrary symmetric and positive-definite conductivity tensor. The problem leads to a coupled system of partial differential equations for the strike-parallel components of the electromagnetic field. E _x, and H _x These equations are numerically approximated by the finite-difference (FD) method, making use of the integro-interpolation approach. As the magnetic component H _x, is constant in the non-conductive air, only equations for the electric mode are approximated within the air layer. The system of linear difference equations, resulting from the FD approximation, can be arranged in such a way that its matrix is symmetric and band-limited, and can be solved, for not too large models, by Gaussian elimination. The algorithm is applied to model situations which demonstrate some non-trivial phenomena caused by electrical anisotropy. In particular, the effect of 2-D anisotropy on the relation between magnetotelluric impedances and induction arrows is studied in detail.     

Quadrangle-grid velocity-stress finite-difference method for elastic-wave-propagation simulation

I present a 2-D numerical-modelling algorithm based on a first-order velocity-stress hyperbolic system and a non-rectangular-grid finite-difference operator. In this method the velocity and stress are defined at different nodes for a staggered grid. The scheme uses non-orthogonal grids, thereby surface topography and curved interfaces can be easily modelled in the seismic-wave-propagation stimulation. The free-surface conditions of complex geometry are achieved by using integral equilibrium equations on the surface, and the stability of the free-surface conditions is improved by introducing local filter modification. The method incorporates desirable qualities of the finite-element method and the staggered-grid finite-difference scheme, which is of high accuracy and low computational cost.     

Statistics and geometry of the eigenspectra of three-dimensional second-rank symmetric random tensors

Geoscientists have undertaken mapping of the Earth's crustal strain (or stress) fields using a great variety of field data. The output can be represented by a 3-D second-rank symmetric random strain tensor. The random principal strains-land rotations of the random tensor are frequently computed. The accuracy is calculated using a first-order approximation. The distribution aspects of the random principal strains and rotations have received almost no attention in Earth Sciences. A first-order approximation of accuracy may not be sufficient if the signal-to-noise ratio is small, as is often the case for geodetically derived random strain tensors. Therefore, the purpose of this paper is to investigate the distribution and estimation problems of the general 3-D second-rank tensor equation GΛG ^T= T , where T is a given 3-D second-rank symmetric random tensor, Λ a diagonal (3 × 3) random eigenvalue matrix, and G a (3 × 3) random orientation matrix, which is also orthogonal. Λ and G are to be estimated (or solved) from T . If some eigenvalues coincide, additional conditions are imposed on the eigenvectors so that they can be chosen uniquely. The joint probability density function (pdf) of the random eigenvalues and rotations will be worked out, given a joint pdf of the elements of random tensors T. Because the rotations are of special interest in Earth Sciences, we shall also derive the joint marginal pdf of random rotations. The geometry of eigenspectra will be studied. The biases of random eigenvalues and rotations will be derived, which have been neglected in the past. They can be very crucial in interpreting the pattern of a derived strain field, however, when applied to a real Earth Science problem. The variance-covariance matrices will be computed using a second-order approximation.     

Three-dimensional image of the Moho undulations beneath the Gulf of Bothnia using wide-angle seismic data

The BABEL marine seismic experiment has been carried out to investigate the lithospheric structure and antecedent tectonic signatures of the Baltic Shield, including the Archaean-Proterozoic collisional structure in the northern part of the Gulf of Bothnia.
Lithospheric seismic-reflection streamer data and simultaneously recorded wide-angle reflection and refraction data collected in the Gulf of Bothnia as part of the BABEL project have been used for 3-D modelling. The distribution of land stations around the Gulf provides a good 3-D ray coverage of the PMP reflection data recorded at the eight stations in the area and allows an estimation of strikes and dips of the Moho boundary in the area. The traveltimes of reflected phases are calculated using a method that utilizes the finite-difference solution of the eikonal equation. The Moho wide-angle-reflection (PMP) traveltimes are modelled using an inversion method. A 2-D model from the Gulf of Bothnia extended into the third dimension is used as an initial model. During the inversion the velocity is kept constant and only the Moho boundary is allowed to vary. To estimate the strike of the Moho boundary and the stability of the inversion, two initial models with different strikes are examined.
The results indicate that the Moho depth in the Gulf of Bothnia undulates and has a maximum depth of 55 km in the south, rising to 42 km in the north. The Moho depth variations seem to be step-like. This change in the Moho depth coincides with the location of the presumed fossil subduction zone in the area. The crustal-thickness variations seem to be well approximated by a nearly 2-D structure striking parallel to a postulated subduction zone immediately to the south of the Skellefte area. The presence of the step at the crust/mantle boundary can be interpreted as a result of a plate-collision event at about 2 Ga.     

Frequency response modelling of seismic waves using finite difference time domain with phase sensitive detection (TD–PSD)

This paper describes an efficient approach for computing the frequency response of seismic waves propagating in 2- and 3-D earth models within which the magnitude and phase are required at many locations. The approach consists of running an explicit finite difference time domain (TD) code with a time harmonic source out to steady-state. The magnitudes and phases at locations in the model are computed using phase sensitive detection (PSD). PSD does not require storage of time-series (unlike a fast Fourier transform), reducing its memory requirements. Additionally, the response from multiple sources can be obtained from a single finite difference run by encoding each source with a different frequency. For 2-D models with many sources, this time domain phase sensitive detection (TD–PSD) approach has a higher arithmetic complexity than direct solution of the finite difference frequency domain (FD) equations using nested dissection re-ordering (FD–ND). The storage requirements for 2-D finite difference TD–PSD are lower than FD–ND. For 3-D finite difference models, TD–PSD has significantly lower arithmetic complexity and storage requirements than FD–ND, and therefore, may prove useful for computing the frequency response of large 3-D earth models.     

Effects of galvanic distortion on magnetotelluric data over a three-dimensional regional structure

Although the galvanic distortion due to local, near-surface inhomogeneities is frequency-independent, its effect on the magnetotelluric data becomes, in a 3-D structure, frequency-dependent. Therefore, both the apparent resistivity and the phase responses are disturbed, and a correction should be carried out prior to the 3-D interpretation in order to retrieve the 3-D regional impedance tensor. In many cases, the structure is 2-D for depths corresponding to a first range of periods and 3-D for longer periods (called 2-D/3-D). For these cases, a simple method which allows us to retrieve the 3-D regional impedance tensor (except the static shift) is presented. The method proposed uses the Groom & Bailey decomposition of the distortion matrix for the short periods. Three examples are presented: two using synthetic data and one employing real data. These examples show the effect of the galvanic distortion over a regional 2-D/3-D model and the retrieval of the regional transfer functions from the distorted ones.     

Spectral estimation on a sphere in geophysics and cosmology

We address the problem of estimating the spherical-harmonic power spectrum of a statistically isotropic scalar signal from noise-contaminated data on a region of the unit sphere. Three different methods of spectral estimation are considered: (i) the spherical analogue of the one-dimensional (1-D) periodogram, (ii) the maximum-likelihood method and (iii) a spherical analogue of the 1-D multitaper method. The periodogram exhibits strong spectral leakage, especially for small regions of area   A ≪ 4π  , and is generally unsuitable for spherical spectral analysis applications, just as it is in 1-D. The maximum-likelihood method is particularly useful in the case of nearly-whole-sphere coverage,   A ≈ 4π  , and has been widely used in cosmology to estimate the spectrum of the cosmic microwave background radiation from spacecraft observations. The spherical multitaper method affords easy control over the fundamental trade-off between spectral resolution and variance, and is easily implemented regardless of the region size, requiring neither non-linear iteration nor large-scale matrix inversion. As a result, the method is ideally suited for most applications in geophysics, geodesy or planetary science, where the objective is to obtain a spatially localized estimate of the spectrum of a signal from noisy data within a pre-selected and typically small region.     

Simulating three-dimensional seismograms in 2.5-dimensional structures by combining two-dimensional finite difference modelling and ray tracing

Finite difference (FD) simulation of elastic wave propagation is an important tool in geophysical research. As large-scale 3-D simulations are only feasible on supercomputers or clusters, and even then the simulations are limited to long periods compared to the model size, 2-D FD simulations are widespread. Whereas in generally 3-D heterogeneous structures it is not possible to infer the correct amplitude and waveform from 2-D simulations, in 2.5-D heterogeneous structures some inferences are possible. In particular, Vidale & Helmberger developed an approach that simulates 3-D waveforms using 2-D FD experiments only. However, their method requires a special FD source implementation technique that is based on a source definition which is not any longer used in nowadays FD codes. In this paper, we derive a conversion between 2-D and 3-D Green tensors that allows us to simulate 3-D displacement seismograms using 2-D FD simulations and the actual ray path determined in the geometrical optic limit. We give the conversion for a source of a certain seismic moment that is implemented by incrementing the components of the stress tensor.
Therefore, we present a hybrid modelling procedure involving 2-D FD and kinematic ray-tracing techniques. The applicability is demonstrated by numerical experiments of elastic wave propagation for models of different complexity.     

Gravimetric determination of an intrusive complex under the Island of Faial (Azores): some methodological improvements

We present some improvements of a gravity inversion method to determine the geometry of the anomalous bodies for priori density contrasts. The 3-D method is based on an exploratory process applied, not for the global model, but for the steps of a growth approach. The (positive and/or negative) anomalous structure is described by successive aggregation of cells, while its corresponding gravity field remains nearly proportional to the observed one. Moreover, a simple (e.g. linear) regional trend can be simultaneously adjusted. The corresponding program is applied to new gravity data on the volcanic island of Faial (Azores archipelago). The inversion approach shows a subsurface anomalous structure for the island, the main feature being an elongated high-density body. The body is interpreted as a compact sheeted dyke swarm, emplaced along Faial-Pico Fracture Zone, a leaky transform structure that forms the current boundary between Eurasian and African plates in the Azores area. The new results in this paper are (1) a Bouguer gravity anomaly map, (2) several improvements in the inversion process (robust process, optimal balance fitness/model magnitude), (3) a new gravimetric method for estimating the mean terrain density, (4) a 3-D model for subsurface mass anomalies in Faial and (5) some interpretative conclusions about a main intrusive complex detected under the island as a wall-like structure extending from a depth of 0.5 to 6 km b.s.l., with a N100°E trend and corresponding to an early fissural volcanic episode controlled by the regional tectonics.     

Seismic imaging of the laterally varying D" region beneath the Cocos Plate

We use an axisymmetric, spherical Earth finite difference algorithm to model SH -wave propagation through cross-sections of laterally varying lower mantle models beneath the Cocos Plate derived from recent data analyses. Synthetic seismograms with dominant periods as short as 4 s are computed for several models: (1) a D" reflector 264 km above the core–mantle boundary with laterally varying S -wave velocity increases of 0.9–2.6 per cent, based on localized structures from a 1-D double-array stacking method; (2) an undulating D" reflector with large topography and uniform velocity increase obtained using a 3-D migration method and (3) cross-sections through the 3-D mantle S -wave velocity tomography model TXBW. We apply double-array stacking to assess model predictions of data. Of the models explored, the S -wave tomography model TXBW displays the best overall agreement with data. The undulating reflector produces a double Scd arrival that may be useful in future studies for distinguishing between D" volumetric heterogeneity and D" discontinuity topography. Synthetics for the laterally varying models show waveform variability not observed in 1-D model predictions. It is challenging to predict 3-D structure based on localized 1-D models when lateral structural variations are on the order of a few wavelengths of the energy used, particularly for the grazing geometry of our data. Iterative approaches of computing synthetic seismograms and adjusting model characteristics by considering path integral effects are necessary to accurately model fine-scale D" structure.     

Surface-wave group-delay and attenuation kernels

We derive both 3-D and 2-D Fréchet sensitivity kernels for surface-wave group-delay and anelastic attenuation measurements. A finite-frequency group-delay exhibits 2-D off-ray sensitivity either to the local phase-velocity perturbation  δ c / c   or to its dispersion  ω(∂/∂ω)(δ c / c )  as well as to the local group-velocity perturbation  δ C / C   . This dual dependence makes the ray-theoretical inversion of measured group delays for 2-D maps of  δ C / C   a dubious procedure, unless the lateral variations in group velocity are extremely smooth.     

Earthquake location in strongly heterogeneous media

Traveltime computation methods for strongly heterogeneous 3-D media developed during recent years are well suited for earthquake location. We present here a new method based on the traveltime algorithm of Podvin-Lecomte, related to the inverse problem formulation of Tarantola & Valette. The Podvin-Lecomte method, based on the Huygens principle, is very robust and allows arbitrary surface topography and station placement even for borehole instruments. First arrival traveltimes are computed for each of the recording stations using a fine 3-D velocity mesh (up to 10⁶ cells on a workstation). The traveltime grid allows the use of the Tarantola & Valette formulation, which enables a full non-linear approach. The solution is given as a 3-D probability density function of hypocentre coordinates, which accounts for the arrival time measurements as well as a priori information for the location, the accuracy of both the arrival time readings and the computation of the theoretical traveltimes. This powerful method called 3DGRIDLOC gives the location of the induced seismicity of the gas field of Lacq (France) using 443 520 cells of a 3-D velocity mesh and the observations from nine recording stations, one of which is located at the bottom of a 3880 m deep borehole. Location of synthetic foci as well as more than 500 actual earthquakes shows the real advantages of this new method versus the classical HYPO71. A new insight into the induced seismicity is now possible: induced seismicity may occur as far away as 10 km from the gas reservoir and involve a much greater volume of rock than expected using earlier locations.     

Two-dimensional finite-difference model of seafloor compliance

Seafloor compliance is the measure of seafloor deformation under a pressure signal. Our new 2-D finite-difference compliance modelling algorithm presents several advantages over the existing compliance models, including the ability to handle any gridded subsurface structure with no limitations on the gradients of the material properties, as well as significantly improved performance. Applying this method to some of the problems inaccessible to previously existing methods, demonstrates that lateral variations in subsurface structure must be accounted for to adequately interpret compliance data. In areas with significant lateral variations, the utilization of 1-D modelling and inversion is likely to result in high interpretation errors, even when additional subsurface structure information is available. We find that flattened pure melt bodies have a significantly higher compliance than cylindrical melt bodies with the same cross-sectional area. The compliance created by such bodies often has side peaks over their edges, which are as strong as or stronger than the central peak, requiring a series of measurements to best constrain their size and shear velocity. Finally, we find that the compliance data are far and away most sensitive to the broad, thick, lower-crustal partial melt zone. Our simple data fitting model for the compliance measurements on the East Pacific Rise at 9°48'N required shear velocities as low as 700 m s⁻¹ in the centre of this zone, far below the values previously estimated using 1-D model based inversions, suggesting higher melt percentages than those previously estimated, while small melt bodies in the upper part of the crust were found to have little or no effect on the measured compliance.