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     

A search for the ratio between gravity variation and vertical displacement due to a surface load

A data space approach to magnetotelluric (MT) inversion reduces the size of the system of equations that must be solved from M × M , as required for a model space approach, to only N × N , where M is the number of model parameter and N is the number of data. This reduction makes 3-D MT inversion on a personal computer possible for modest values of M and N . However, the need to store the N × M sensitivity matrix J remains a serious limitation. Here, we consider application of conjugate gradient (CG) methods to solve the system of data space Gauss–Newton equations. With this approach J is not explicitly formed and stored, but instead the product of J with an arbitrary vector is computed by solving one forward problem. As a test of this data space conjugate gradient (DCG) algorithm, we consider the 2-D MT inverse problem. Computational efficiency is assessed and compared to the data space Occam's (DASOCC) inversion by counting the number of forward modelling calls. Experiments with synthetic data show that although DCG requires significantly less memory, it generally requires more forward problem solutions than a scheme such as DASOCC, which is based on a full computation of J .     

A Bayesian approach to inverse modelling of stratigraphy, part 2: Validation tests

A novel inverse modelling method is applied to the problem of constraining the environmental parameters (e.g. relative sea level, sediment supply) that control stratigraphic architecture. This technique links forward modelling of shallow-marine wave/storm-dominated stratigraphy to a combination of inverse methods formulated in a Bayesian framework. We present a number of examples in which relative sea-level and sediment-supply curves were inferred from synthetic vertical successions of grain size (e.g. wells) and synthetic thickness curves (e.g. seismically derived isopachs) extracted from a forward model simulation. These examples represent different scenarios that are designed to test the impact of data distribution, quantity and quality on the uncertainty of the inferred parameters. The inverse modelling approach successfully reproduces the gross stratigraphic architectures and relative sea level and sediment-supply histories of the synthetic forward model simulation, within the constraints of the modelled data quality. The relative importance of the forcing parameters can be evaluated by their sensitivity and impact on the inverted data. Of equal importance, the inverse results allow complete characterisation of the uncertainties inherent to the stratigraphic modelling tool and to the data quality, quantity and distribution. The numerical scheme also successfully deals with the problem of non-uniqueness of the solution of the inverse problem. These preliminary results suggest that the inverse method is a powerful tool in constraining stratigraphic architecture for hydrocarbon reservoir characterisation and modelling, and it may ultimately provide a process-based geological complement to standard geostatistical tools.     

Geophysical inversion with a neighbourhood algorithm—I. Searching a parameter space

This paper presents a new derivative-free search method for finding models of acceptable data fit in a multidimensional parameter space. It falls into the same class of method as simulated annealing and genetic algorithms, which are commonly used for global optimization problems. The objective here is to find an ensemble of models that preferentially sample the good data-fitting regions of parameter space, rather than seeking a single optimal model. (A related paper deals with the quantitative appraisal of the ensemble.)
  The new search algorithm makes use of the geometrical constructs known as Voronoi cells to derive the search in parameter space. These are nearest neighbour regions defined under a suitable distance norm. The algorithm is conceptually simple, requires just two 'tuning parameters', and makes use of only the rank of a data fit criterion rather than the numerical value. In this way all difficulties associated with the scaling of a data misfit function are avoided, and any combination of data fit criteria can be used. It is also shown how Voronoi cells can be used to enhance any existing direct search algorithm, by intermittently replacing the forward modelling calculations with nearest neighbour calculations.
  The new direct search algorithm is illustrated with an application to a synthetic problem involving the inversion of receiver functions for crustal seismic structure. This is known to be a non-linear problem, where linearized inversion techniques suffer from a strong dependence on the starting solution. It is shown that the new algorithm produces a sophisticated type of 'self-adaptive' search behaviour, which to our knowledge has not been demonstrated in any previous technique of this kind.     

New advances in three-dimensional controlled-source electromagnetic inversion

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

Inversion for interface structure using teleseismic traveltime residuals

An inversion method is presented for the reconstruction of interface geometry between two or more crustal layers from teleseismic traveltime residuals. The method is applied to 2-D models consisting of continuous interfaces separating constant-velocity layers. The forward problem of determining ray paths and traveltimes between incident wave fronts below the structure and receivers located on the Earth's surface is solved by an efficient and robust shooting method. A conjugate gradient method is employed to solve the inverse problem of minimizing a least-squares type objective function based on the difference between observed and calculated traveltimes. Teleseismic data do not accurately constrain average vertical structure, so a priori information in the form of layer velocities and average layer thicknesses is required. Synthetic tests show that the method can be used to reconstruct interface geometry accurately, even in the presence of data noise. Tests also show that, if layer velocities and initial interface positions are poorly chosen, lateral structure is still recoverable. The inversion method was applied to previously published teleseismic data recorded by an in-line array of portable seismographs that traversed the northern margin of the Musgrave Block, central Australia. The solution based on interface parametrization is consistent with models given by other studies that used the same data but different methods, most notably the standard tomographic approach that inverts for velocity rather than interface structure.     

Rapid relaxation inversion of CSAMT data

In this paper an inversion algorithm for controlled-source audio frequency magnetotelluric data is presented. This algorithm combines 2.5-D finite element forward modelling with the concepts of rapid relaxation inversion of magnetotelluric data. The inversion uses the same technique to compute sensitivities as the rapid relaxation inversion, and these approximate sensitivities are validated by comparison with exact 2.5-D sensitivities. The comparison shows that the approximate sensitivities are similar in shape to the exact sensitivities when transmitter–receiver offsets are greater than one skin depth in the Earth. The magnitudes of the two sensitivities differ but the variations with depth are similar. Tests of the algorithm on synthetic data and field data provide promising results.     

A Bayesian approach to inverse modelling of stratigraphy, part 1: method

The inference of ancient environmental conditions from their preserved response in the sedimentary record still remains an outstanding issue in stratigraphy. Since the 1970s, conceptual stratigraphic models (e.g. sequence stratigraphy) based on the underlying assumption that accommodation space is the critical control on stratigraphic architecture have been widely used. Although these methods considered more recently other possible parameters such as sediment supply and transport efficiency, they still lack in taking into account the full range of possible parameters, processes, and their complex interactions that control stratigraphic architecture. In this contribution, we present a new quantitative method for the inference of key environmental parameters (specifically sediment supply and relative sea level) that control stratigraphy. The approach combines a fully non‐linear inversion scheme with a ‘process–response’ forward model of stratigraphy. We formulate the inverse problem using a Bayesian framework in order to sample the full range of possible solutions and explicitly build in prior geological knowledge. Our methodology combines Reversible Jump Markov chain Monte Carlo and Simulated Tempering algorithms which are able to deal with variable‐dimensional inverse problems and multi‐modal posterior probability distributions, respectively. The inverse scheme has been linked to a forward stratigraphic model, BARSIM (developed by Joep Storms, University of Delft), which simulates shallow‐marine wave/storm‐dominated systems over geological timescales. This link requires the construction of a likelihood function to quantify the agreement between simulated and observed data of different types (e.g. sediment age and thickness, grain size distributions). The technique has been tested and validated with synthetic data, in which all the parameters are specified to produce a ‘perfect’ simulation, although we add noise to these synthetic data for subsequent testing of the inverse modelling approach. These tests addressed convergence and computational‐overhead issues, and highlight the robustness of the inverse scheme, which is able to assess the full range of uncertainties on the inferred environmental parameters and facies distributions.     

Recovering magnetic susceptibility from electromagnetic data over a one-dimensional earth

While the inversion of electromagnetic data to recover electrical conductivity has received much attention, the inversion of those data to recover magnetic susceptibility has not been fully studied. In this paper we invert frequency-domain electromagnetic (EM) data from a horizontal coplanar system to recover a 1-D distribution of magnetic susceptibility under the assumption that the electrical conductivity is known. The inversion is carried out by dividing the earth into layers of constant susceptibility and minimizing an objective function of the susceptibility subject to fitting the data. An adjoint Green's function solution is used in the calculation of sensitivities, and it is apparent that the sensitivity problem is driven by three sources. One of the sources is the scaled electric field in the layer of interest, and the other two, related to effective magnetic charges, are located at the upper and lower boundaries of the layer. These charges give rise to a frequency-independent term in the sensitivities. Because different frequencies penetrate to different depths in the earth, the EM data contain inherent information about the depth distribution of susceptibility. This contrasts with static field measurements, which can be reproduced by a surface layer of magnetization. We illustrate the effectiveness of the inversion algorithm on synthetic and field data and show also the importance of knowing the background conductivity. In practical circumstances, where there is no a priori information about conductivity distribution, a simultaneous inversion of EM data to recover both electrical conductivity and susceptibility will be required.     

Estimating uncertainties on net erosion from well-log porosity data

Estimating the amount of erosion experienced by a sedimentary basin during its geological history plays a key role in basin modelling. In this paper, we present a novel probabilistic approach to estimate net erosion from porosity–depth data from a single well. Our approach uses a Markov chain Monte Carlo algorithm which readily allows us to deal with imprecise knowledge of the lithology-dependent compaction parameters in a joint inversion scheme using multiple lithologies. The results using synthetic data highlight the advantages of our approach over conventional techniques for net erosion estimation: (a) uncertainties on compaction parameters can be effectively mapped into a probabilistic solution for net erosion; (b) posterior uncertainties are easy to quantify; (c) the joint inversion scheme can automatically reconcile porosity data from different lithologies. Our results also underscore the critical role of prior assumptions on controlling the retrieved estimates for net erosion. Using real data from a well in the Barents Sea, we simulate three possible scenarios of variable prior assumptions on compaction parameters to demonstrate the general applicability of our approach. Strong prior assumptions on the compaction parameters led to unrealistic estimates of net erosion for the target well, indicating the assumptions are probably inappropriate. Our preferred strategy for this dataset is to include additional data to constrain the normal compaction trend of the sediments. This provides a net erosion estimate for the target well of about 2300 m with a standard deviation of 140 m which is in line with previous studies. Finally, we discuss potential guidelines to deal with real applications in which data from normally compacted sediments are not available. One is to use our algorithm as a hypothesis-testing tool to evaluate the results under a large set of assumed compaction parameters. A second is to infer compaction parameters and net erosion simultaneously from the target well porosity data. Although appealing and successful with synthetic data, this strategy provides results which are strongly dependent on the calibration data and the geological history of the sediments sampled by the target well.     

Inversion of time domain three-dimensional electromagnetic data

We present a general formulation for inverting time domain electromagnetic data to recover a 3-D distribution of electrical conductivity. The forward problem is solved using finite volume methods in the spatial domain and an implicit method (Backward Euler) in the time domain. A modified Gauss–Newton strategy is employed to solve the inverse problem. The modifications include the use of a quasi-Newton method to generate a pre-conditioner for the perturbed system, and implementing an iterative Tikhonov approach in the solution to the inverse problem. In addition, we show how the size of the inverse problem can be reduced through a corrective source procedure. The same procedure can correct for discretization errors that inevidably arise. We also show how the inverse problem can be efficiently carried out even when the decay time for the conductor is significantly larger than the repetition time of the transmitter wave form. This requires a second processor to carry an additional forward modelling. Our inversion algorithm is general and is applicable for any electromagnetic field  ( E , H , d B / dt )  measured in the air, on the ground, or in boreholes, and from an arbitrary grounded or ungrounded source. Three synthetic examples illustrate the basic functionality of the algorithm, and a result from a field example shows applicability in a larger-scale field example.     

Effects of topography and crustal heterogeneities on the source estimation of LP event at Kilauea volcano

The main goal of this study is to improve the modelling of the source mechanism associated with the generation of long period (LP) signals in volcanic areas. Our intent is to evaluate the effects that detailed structural features of the volcanic models play in the generation of LP signal and the consequent retrieval of LP source characteristics. In particular, effects associated with the presence of topography and crustal heterogeneities are here studied in detail. We focus our study on a LP event observed at Kilauea volcano, Hawaii, in 2001 May. A detailed analysis of this event and its source modelling is accompanied by a set of synthetic tests, which aim to evaluate the effects of topography and the presence of low velocity shallow layers in the source region. The forward problem of Green's function generation is solved numerically following a pseudo-spectral approach, assuming different 3-D models. The inversion is done in the frequency domain and the resulting source mechanism is represented by the sum of two time-dependent terms: a full moment tensor and a single force. Synthetic tests show how characteristic velocity structures, associated with shallow sources, may be partially responsible for the generation of the observed long-lasting ringing waveforms. When applying the inversion technique to Kilauea LP data set, inversions carried out for different crustal models led to very similar source geometries, indicating a subhorizontal cracks. On the other hand, the source time function and its duration are significantly different for different models. These results support the indication of a strong influence of crustal layering on the generation of the LP signal, while the assumption of homogeneous velocity model may bring to misleading results.     

Simultaneous velocity and interface tomography of normal-incidence and wide-aperture seismic traveltime data

A tomographic inversion technique that inverts traveltimes to obtain a model of the subsurface in terms of velocities and interfaces is presented. It uses a combination of refraction, wide-angle reflection and normal-incidence data, it simultaneously inverts for velocities and interface depths, and it is able to quantify the errors and trade-offs in the final model. The technique uses an iterative linearized approach to the non-linear traveltime inversion problem. The subsurface is represented as a set of layers separated by interfaces, across which the velocity may be discontinuous. Within each layer the velocity varies in two dimensions and has a continuous first derivative. Rays are traced in this medium using a technique based on ray perturbation theory, and two-point ray tracing is avoided by interpolating the traveltimes to the receivers from a roughly equidistant fan of rays. The calculated traveltimes are inverted by simultaneously minimizing the misfit between the data and calculated traveltimes, and the roughness of the model. This 'smoothing regularization' stabilizes the solution of the inverse problem. In practice, the first iterations are performed with a high level of smoothing. As the inversion proceeds, the level of smoothing is gradually reduced until the traveltime residual is at the estimated level of noise in the data. At this point, a minimum-feature solution is obtained, which should contain only those features discernible over the noise.
The technique is tested on a synthetic data set, demonstrating its accuracy and stability and also illustrating the desirability of including a large number of different ray types in an inversion.     

Wavelet modelling of broad-band receiver functions

We present a wavelet modelling approach to invert for S -wave velocities from broad-band receiver functions. Taking spline function as the basic wavelet, the broad-band receiver function is decomposed into five resolution scales by Mallat's pyramid algorithm. The linearized least-squares inversion procedure is applied to every resolution scale. The fifth-scale approximation of receiver function is first inverted to recover the slowly varying background velocity variations with respect to a reference model. This solution is then taken as the initial model for fitting the fourth-scale wavelet coefficients of receiver function to further tune the solution to resolve sharper variations. This procedure is iteratively carried out up to the first-scale wavelet coefficients of receiver function. In this manner, the model neighbourhood containing the global minimum is first searched from the coarsest-scale receiver function, and the search gradually focuses on the global minimum by introducing finer-scale information of receiver function. Noise-free synthetic receiver function tests show that wavelet modelling of receiver functions can guide a certain range of initial models to converge to the true velocity distribution. Tests on actual data indicate that wavelet modelling can provide results very similar to those inferred by joint inversion of receiver function and surface wave dispersion.     

Interpretation of long-offset transient electromagnetic data from the Odenwald area, Germany, using two-dimensional modelling

The standard 1-D inversion approach for the interpretation of transient electromagnetic (TEM) data usually fails in the presence of near-surface conductivity anomalies. Since multidimensional inversion codes are not routinely available, the only alternative to discarding the data may be trial-and-error forward modelling. We interpret data from a long-offset transient electromagnetic (LOTEM) survey which was carried out in 1995 in the Odenwald area, using 2-D finite-difference modelling. We focus on a subsegment of the LOTEM profile, which was shot with two different electric dipole transmitters. A model is found which consistently explains the electric and magnetic field data at eight locations for both transmitters. First, we introduce a conductive dyke under the receiver spread to explain sign reversals in the magnetic field transients. A conductive slab under one of the transmitters is required to obtain a reasonable quantitative fit for that transmitter. Consideration of the electric field data then requires a modification of the layered earth background. Finally, we study the response of a crustal conductor, which was the original target of the survey. The data are sensitive to the conductor, and for the investigated subset of the data the fits are slightly better without the conductive layer.     

Spectral-finite-element approach to two-dimensional electromagnetic induction in a spherical earth

We present a spectral-finite-element approach to the 2-D forward problem for electromagnetic induction in a spherical earth. It represents an alternative to a variety of numerical methods for 2-D global electromagnetic modelling introduced recently (e.g. the perturbation expansion approach, the finite difference scheme). It may be used to estimate the effect of a possible axisymmetric structure of electrical conductivity of the mantle on surface observations, or it may serve as a tool for testing methods and codes for 3-D global electromagnetic modelling. The ultimate goal of these electromagnetic studies is to learn about the Earth's 3-D electrical structure.
Since the spectral-finite-element approach comes from the variational formulation, we formulate the 2-D electromagnetic induction problem in a variational sense. The boundary data used in this formulation consist of the horizontal components of the total magnetic intensity measured on the Earth's surface. In this the variational approach differs from other methods, which usually use spherical harmonic coefficients of external magnetic sources as input data. We verify the assumptions of the Lax-Milgram theorem and show that the variational solution exists and is unique. The spectral-finite-element approach then means that the problem is parametrized by spherical harmonics in the angular direction, whereas finite elements span the radial direction. The solution is searched for by the Galerkin method, which leads to the solving of a system of linear algebraic equations. The method and code have been tested for Everett & Schultz's (1995) model of two eccentrically nested spheres, and good agreement has been obtained.     

Forward modelling receiver functions for crustal structure beneath station TBZ (Trabzon, Turkey)

We use teleseismic three-component digital data from the Trabzon, Turkey broadband seismic station TBZ to model the crustal structure by the receiver function method. The station is located at a structural transition from continental northeastern Anatolia to the oceanic Black Sea basin. Rocks in the region are of volcanic origin covered by young sediments. By forward modelling the radial receiver functions, we construct 1-D crustal shear velocity models that include a lower crustal low-velocity zone, indicating a partial melt mechanism which may be the source of surfacing magmatic rocks and regional volcanism. Within the top 5 km, velocities increase sharply from about 1.5 to 3.5 km s⁻¹. Such near-surface low velocities are caused by sedimentation, extending from the Black Sea basin. Velocities at around 20 km depth have mantle-like values (about 4.25 km s⁻¹ ), which easily correlate to magmatic rocks cropping out on the surface. At 25 km depth there is a thin low-velocity layer of about 4.0 km s⁻¹. The average Moho velocity is about 4.6 km s⁻¹, and its depth changes from 32 to 40 km. Arrivals on the tangential components indicate that the Moho discontinuity dips approximately southwards, in agreement with the crustal thickening to the south. We searched for the solution of receiver functions around the regional surface wave group velocity inversion results, which helped alleviate the multiple solution problem frequently encountered in receiver function modelling.
Station TBZ is a recently deployed broadband seismic station, and the aim of this study is to report on the analysis of new receiver function data. The analysis of new data in such a structurally complex region provides constraining starting models for future structural studies in the region.     

Joint inversion of temperature,vitrinite reflectance and fission tracks in apatite with examples from the eastern North Sea area

As sediment accumulation indicates basin subsidence, erosion often is understood as tectonic uplift, but the amplitude and timing may be difficult to determine because the sedimentary record is missing. Quantification of erosion therefore requires indirect evidence, for example thermal indicators such as temperature, vitrinite reflectance and fission tracks in apatite. However, as always, the types and quality of data and the choice of models are important to the results. For example, considering only the thermal evolution of the sedimentary section discards the thermal time constant of the lithosphere and essentially ignores the temporal continuity of the thermal structure. Furthermore, the types and density of thermal indicators determine the solution space of deposition and erosion, the quantification of which calls for the use of inverse methods, which can only be successful when all models are mutually consistent. Here, we use integrated basin modelling and Markov Chain Monte Carlo inversion of four deep boreholes to show that the erosional pattern along the Sorgenfrei–Tornquist Zone (STZ) in the eastern North Sea is consistent with a tectonic model of tectonic inversion based on compression and relaxation of an elastic plate. Three wells in close proximity SW of the STZ have different data and exhibit characteristic differences in erosion estimates but are consistent with the formation of a thick chalk sequence, followed by minor Cenozoic erosion during relaxation inversion. The well on the inversion ridge requires ca. 1.7 km Jurassic-Early Cretaceous sedimentation followed by Late Cretaceous–Palaeocene erosion during inversion. No well demands thick Cenozoic sedimentation followed by equivalent significant Neogene exhumation. When data are of high quality and models are consistent, the thermal indicator method yields significant results with important tectonic and geodynamic implications.     

Integration of linear programming and GIS for land-use modelling

Abstract

Geographical Information Systems (GIS) are becoming basic tools for a wide variety of earth science and land-use applications. This article presents linear programming (LP) as a promising tool for spatial modelling within a GIS. Although LP is not properly a spatial technique, it may be used to optimize spatial distributions or to guide the integration of variables. An example of the use of LP in land-use planning is described, with minimizing rural unemployment as the main goal. Technical, financial and ecological constraints are established to show the influence of several limitations on achieving the optimal solution. LP makes it possible to achieve optimal land-use, where the objective is maximized and the constraints respected. LP can also be used to simulate different planning scenarios, by modifying both the objective function coefficients and the constraints. The integration of LP and GIS is presented in two phases: (i) acquisition of attribute data for the LP model, and (ii) modelling and mapping the results.     

Numerical scheme for the inversion of acoustical impedance profile based on the Gelfand-Levitan method

Summary. An exact method for the solution of the inverse problem in plane wave propagation modelled after the Gelfand-Levitan technique is reviewed and refined. A numerical scheme for the solution of the integral equation that arises in the method is proposed. A discussion on the stability and an error analysis of the numerical approximation are presented. The applicability of the inversion algorithm is demonstrated in a numerical experiment.