· · · Institutional Repository

The performance of a multigrid solver for time-harmonic electromagnetic problems in geophysical settings was investigated. With the low frequencies used in geophysical surveys for deeper targets, the light-speed waves in the earth can be neglected. Diffusion of induced currents is the dominant physical effect. The governing equations were discretised by the Finite-Integration Technique. The resulting set of discrete equation was solved by a multigrid method. The multigrid method provided excellent convergence with constant grid spacings, but not on stretched grids. The slower convergence rate of the multigrid method could be compensated by using bicgstab2, in which case multigrid acted as a preconditioner. Still, the overall performance was less than satisfactory with substantial grid stretching.

Migration velocity analysis with the wave equation can be accomplished by focusing of extended migration images, obtained by introducing a subsurface offset or shift. A reflector in the wrong velocity model will show up as a curve in the extended image. In the correct model, it should collapse to a point. The usual approach to obtain a focused image involves a cost functional that penalizes energy in the extended image at non-zero shift. Its minimization by a gradient-based method should then produce the correct velocity model. Here, asymptotic analysis and numerical examples show that this method may be too sensitive to amplitude peaks at large shifts at the wrong depth and to artifacts. A more robust alternative is proposed that can be interpreted as a generalization of stack power and maximizes the energy at zero subsurface shift. A real-data example is included.

Seismic data enable imaging of the Earth, not only of velocity and density but also of attenuation contrasts. Unfortunately, the Born approximation of the constant-density visco-acoustic wave equation, which can serve as a forward modelling operator related to seismic migration, exhibits an ambiguity when attenuation is included. Different scattering models involving velocity and attenuation perturbations may provide nearly identical data. This result was obtained earlier for scatterers that did not contain a correction term for causality. Such a term leads to dispersion when considering a range of frequencies. We demonstrate that with this term, linearized inversion or iterative migration will almost, but not fully, remove the ambiguity. We also investigate if attenuation imaging suffers from the same ambiguity when using non-linear or full waveform inversion. A numerical experiment shows that non-linear inversion with causality convergences to the true model, whereas without causality, a substantial difference with the true model remains even after a very large number of iterations. For both linearized and non-linear inversion, the initial update in a gradient-based optimization scheme that minimizes the difference between modelled and observed data is still affected by the ambiguity and does not provide a good result. This first update corresponds to a classic migration operation. In our numerical experiments, the reconstructed model started to approximate the true model only after a large number of iterations.

Wave-equation traveltime tomography tries to obtain a subsurface velocity model from seismic data, either passive or active, that explains their traveltimes. A key step is the extraction of traveltime differences, or relative phase shifts, between observed and modelled finite-frequency waveforms. A standard approach involves a correlation of the observed and measured waveforms. When the amplitude spectra of the waveforms are identical, the maximum of the correlation is indicative of the relative phase shift. When the amplitude spectra are not identical, however, this argument is no longer valid. We propose an alternative criterion to measure the relative phase shift. This misfit criterion is a weighted norm of the correlation and is less sensitive to differences in the amplitude spectra. For practical application it is important to use a sensitivity kernel that is consistent with the way the misfit is measured. We derive this sensitivity kernel and show how it differs from the standard banana–doughnut sensitivity kernel. We illustrate the approach on a cross-well data set.

We discuss some computational aspects of resistivity imaging by inversion of offshore controlled-source electromagnetic data. We adopt the classic approach to imaging by formulating it as an inverse problem. A weighted least-squares functional measures the misfit between synthetic and observed data. Its minimization by a quasi-Newton algorithm requires the gradient of the functional with respect to the model parameters. We compute the gradient with the adjoint-state technique. Preconditioners can improve the convergence of the inversion. Diagonal preconditioner based on a Born approximation are commonly used. In the context of CSEM inversion, the Born approximation is not really accurate, this limits the possibility of estimating a correct approximation of the Hessian in a smooth medium or, in fact, in any reference background that does not roughly account for the resistors. We hence rely on the limited memory BFGS approximation of the inverse of the Hessian and we improve the inversion convergence with the help of a heuristic data and depth weighting. Based on a numerical example, we show that a simple exponential depth weighting combined with an offset or frequency data weighting significantly improves the convergence rate of a deep-water controlled-source electromagnetic data inversion.

In seismic imaging, one tries to infer the medium properties of the subsurface from seismic reflection data. These data are the result of an active source experiment, where an explosive source and an array of receivers are placed at the surface. Due to the absence of low frequencies in the data, the corresponding inverse problem is strongly non-linear in the slowly varying component of the velocity. The least-squares misfit functional typically exhibits local minima and has a small basin of attraction. The usual approach of fitting the data in a least-squares sense by employing a gradient-based optimisation method will therefore most likely result in a wrong velocity model. In the geophysical community, this problem has long been recognised and alternative formulations of the inverse problem have been developed. We review several of these formulations and analyse the sensitivity to the error in the smooth velocity component. This analysis is carried out for laterally homogeneous velocities using an asymptotic solution of the wave equation. The analysis suggests that formulations which are geared towards fitting the phases of the data, rather than the amplitudes, have smooth corresponding misfit functionals with a large basin of attraction.

The classic least-squares cost functional for full waveform inversion suffers from local minima due to loop skipping in the absence of low frequencies in the seismic data. Velocity model building based on subsurface spatial or temporal shifts may break down in the presence of multiples in the data. Cost functionals that translate this idea to the data domain, with offset- or time-shifts, can handle multiples. An earlier data-domain formulation suffered from cross-talk between events. Here, we present a multishot extension that should be less sensitive to cross-talk. It has the property of an annihilator, similar to the functional used for velocity analysis with extended images based on subsurface shifs. However, since it operates in the data domain, it should be able to handle multiples. For 2-D models with line acquistion, the proposed functional applies a singular-value decomposition on the observed data and uses the eigenvectors to build data panel that should be diagonal in the correct velocity model, but has significant off-diagonal entries in the wrong model. By minimizing these offdiagonal entries or maximizing the main diagonal, the correct model should be found. We therefore named it the diagonalator. We present initial tests on a simple, horizontally layered velocity model.

We measured the electric parameters for four different configurations of unconsolidated homogeneous and layered sands as a function of frequency, water saturation, and salinity under fluid flow conditions. Our objective is to determine if the effect of heterogeneities at scales much smaller than the skin depth can be captured by introducing effective frequency‐dependent electrical values whose behavior can be described by simple functions. We employed the parallel plate capacitor technique to measure the complex impedance over a broad frequency range, from 100 kHz up to 3 MHz. We conducted main drainage and secondary imbibition cycles at atmospheric pressure and temperatures between 21°C and 22°C. The hysteretic effect in the real part of the effective complex permittivity at higher concentrations of NaCl is more pronounced for the homogeneous configurations than for the heterogeneous samples. Effective medium theory works well for dry and saturated layered sand, when the NaCl solution concentration is 1 mmol/l. It fails for fully saturated layered sands at salinities of 10 mmol/l or more. It also does not work for partially saturated sands, independent of salinity. A description of the electric properties of a layered sand at all saturation levels by means of an effective homogeneous medium will therefore require a dependence on frequency, saturation level, and salinity of the pore fluid. An extended version of the Cole‐Cole model fits the nonmonotonic behavior of the real part of permittivity versus saturation.

The spreading adoption of computationally intensive techniques such as Reverse Time Migration and Full Waveform Inversion increases the need of efficiently solving the three-dimensional wave equation. Common finite-difference discretization schemes lose their accuracy and efficiency in complex geological settings with discontinuities in the material properties and topography. Finite elements on tetrahedral meshes follow the interfaces while maintaining their accuracy and can have smaller meshes if the elements are scaled with the velocity. Here, we consider two higher-order finite element methods that allow for explicit time stepping: the continuous mass-lumped finite-element method (CMLFE) and the symmetric interior penalty discontinuous Galerkin method (SIPDG). The price paid for the ability to perform explicit time stepping is an increase in computational cost: CMLFE requires a larger number of discretization nodes to preserve accuracy, whereas SIPDG needs additional fluxes to impose the continuity of the solution. Therefore, it is not obvious which one is more efficient. We compare the two methods in terms of accuracy, stability and computational cost. Experiments on a three-dimensional problem with a dipping interface show that CMLFE and SIPDG have similar stability conditions, accuracy and efficiency, the last being measured as the computational time required to reach a given accuracy of the result.

We modeled time-domain EM measurements of induction currents for marine and land applications with a frequency-domain code. An analysis of the computational complexity of a number of numerical methods shows that frequency-domain modeling followed by a Fourier transform is an attractive choice if a sufficiently powerful solver is available. A recently developed, robust multigrid solver meets this requirement. An interpolation criterion determined the automatic selection of frequencies. The skin depth controlled the construction of the computational grid at each frequency. Tests of the method against exact solutions for some simple problems and a realistic marine example demonstrate that a limited number of frequencies suffice to provide time-domain solutions after piecewise-cubic Hermite interpolation and a fast Fourier transform.

In the application of controlled source electromagnetics for reservoir monitoring on land, repeatability errors in the source will mask the time-lapse signal due to hydrocarbon production when recording surface data close to the source. We demonstrate that at larger distances, the airwave will still provide sufficient illumination of the target. The primary airwave diffuses downward into the earth and then is scattered back to the surface. The time-lapse difference of its recorded signal reveals the outline on the surface of the resistivity changes in a hydrocarbon reservoir under production. However, repeatability errors in the primary airwave can destroy the signal-to-noise ratio of the time-lapse data. We present a simple and effective method to remove the primary airwave from the data, which we call partial airwave removal. For a homogeneous half space and a delta-function type of source, the surface expression of the airwave does not depend on frequency. For this reason, the primary airwave can be subtracted from the data using recordings at two frequencies, one low enough with a skin depth of the order of the reservoir depth that is sensitive to the reservoir, the other high enough to only sense the near surface. The method does not affect secondary airwave components created by signals that have propagated through the earth and returned to the surface. We show that the method provides a direct indicator of production-related time-lapse changes in the reservoir. We illustrate this for several models, including a general 3D heterogeneous model and one with strong surface topography, for situations where survey repeatability errors are large.

Macroscopic measurements of electrical resistivity require frequency-dependent effective models that honor the microscopic effects observable in macroscopic measurements. Effective models based on microscopic physics exist alongside with empirical models. We adopted an empirical model approach to modify an existing physical model. This provided a description of electrical resistivity as a function of not only frequency, but also water saturation. We performed two-electrode laboratory measurements of the complex resistivity on a number of fine and medium-grained unconsolidated sand packs saturated with water of three different salinities. For frequencies between 0.1 and 1 MHz, the data were fitted with the new model and compared to fits with Archie’s law. Our model described the relaxation times and DC resistivity values as negative exponential functions with increasing water saturation. All data could be accurately described as a function of frequency and water saturation with nine parameters.

In the application of controlled source electromagnetics for reservoir monitoring on land, the timelapse signal measured with a surface-to-surface acquisition can reveal the lateral extent on the surface of resistivity changes at depth in a hydrocarbon reservoir under production. However, a direct interpretation of the time-lapse signal may generally be difficult and biased. We investigated if non-linear inversion can use time-lapse responses to characterize the subsurface resistivity changes. We examined two different strategies, using a full non-linear inversion algorithm as the interpretation tool: inverting the reference and monitor data independently or in sequence. In the second case, the inversion result of the reference data set serves as an initial guess for the inversion of the monitor data set. Numerical examples show that independent inversion of the data sets can provide an estimate of the depth and lateral extent of the resistivity changes. The second strategy of sequential inversion produces less satisfactory results. We illustrate the independent inversion approach for an example with large survey repeatability errors are large and another one with a complex overburden.

Modeling and imaging techniques for geophysics are extremely demanding in terms of computational resources. Seismic data attempt to resolve smaller scales and deeper targets in increasingly more complex geologic settings. Finite elements enable accurate simulation of time-dependent wave propagation in heterogeneous media. They are more costly than finite-difference methods, but this is compensated by their superior accuracy if the finite-element mesh follows the sharp impedance contrasts and by their improved efficiency if the element size scales with wavelength, hence with the local wave velocity. However, 3D complex geologic settings often contain details on a very small scale compared to the dominant wavelength, requiring the mesh to contain elements that are smaller than dictated by the wavelength. Also, limitations of the mesh generation software may produce regions where the elements are much smaller than desired. In both cases, this leads to a reduction of the time step required to solve the wave propagation and significantly increases the computational cost. Local time stepping (LTS) can improve the computational efficiency and speed up the simulation. We evaluated a local formulation of an LTS scheme with second-order accuracy for the discontinuous Galerkin finite-element discretization of the wave equation.We tested the benefits of the scheme by considering a geologic model for a North-Sea-type example.

The time-harmonic wave equation, also known as the Helmholtz equation, is obtained if the constant-density acoustic wave equation is transformed from the time domain to the frequency domain. Its discretization results in a large, sparse, linear system of equations. In two dimensions, this system can be solved efficiently by a direct method. In three dimensions, direct methods cannot be used for problems of practical sizes because the computational time and the amount of memory required become too large. Iterative methods are an alternative. These methods are often based on a conjugate gradient iterative scheme with a preconditioner that accelerates its convergence. The iterative solution of the time-harmonic wave equation has long been a notoriously difficult problem in numerical analysis. Recently, a new preconditioner based on a strongly damped wave equation has heralded a breakthrough. The solution of the linear system associated with the preconditioner is approximated by another iterative method, the multigrid method. The multigrid method fails for the original wave equation but performs well on the damped version. The performance of the new iterative solver is investigated on a number of 2D test problems. The results suggest that the number of required iterations increases linearly with frequency, even for a strongly heterogeneous model where earlier iterative schemes fail to converge. Complexity analysis shows that the new iterative solver is still slower than a time-domain solver to generate a full time series. We compare the time-domain numeric results obtained using the new iterative solver with those using the direct solver and conclude that they agree very well quantitatively. The new iterative solver can be applied straightforwardly to 3D problems.

Full waveform inversion suffers from local minima, due to a lack of low frequencies in the data. A reflector below the zone of interest may, however, help in recovering the long-wavelength components of a velocity perturbation, as demonstrated in a paper by Mora. With the Born approximation for the perturbation in a reference model consisting of two homogeneous isotropic acoustic halfspaces, analytic expressions can be found that describe the spatial spectrum of the recorded seismic signal as a function of the spatial spectrum of the inhomogeneity. We study this spectrum in more detail by separately considering direct, reflected and head waves. Taking the reflection coefficient of the deeper reflector into account, we obtain sensitivity estimates for each of these types of waves. Although the head waves have a relatively small contribution to the reconstruction of the velocity perturbation, compared to the other waves, they contain reliable long-wavelength information that can be beneficial for full waveform inversion.