Geophysical Problems Research Articles

Model-Based Probabilistic Inversion Using Magnetic Data: A Case Study on the Kevitsa Deposit

Structural geological models are widely used to represent relevant geological interfaces and property distributions in the subsurface. Considering the inherent uncertainty of these models, the non-uniqueness of geophysical inverse problems, and the growing availability of data, there is a need for methods that integrate different types of data consistently and consider the uncertainties quantitatively. Probabilistic inference provides a suitable tool for this purpose. Using a Bayesian framework, geological modeling can be considered as an integral part of the inversion and thereby naturally constrain geophysical inversion procedures. This integration prevents geologically unrealistic results and provides the opportunity to include geological and geophysical information in the inversion. This information can be from different sources and is added to the framework through likelihood functions. We applied this methodology to the structurally complex Kevitsa deposit in Finland. We started with an interpretation-based 3D geological model and defined the uncertainties in our geological model through probability density functions. Airborne magnetic data and geological interpretations of borehole data were used to define geophysical and geological likelihoods, respectively. The geophysical data were linked to the uncertain structural parameters through the rock properties. The result of the inverse problem was an ensemble of realized models. These structural models and their uncertainties are visualized using information entropy, which allows for quantitative analysis. Our results show that with our methodology, we can use well-defined likelihood functions to add meaningful information to our initial model without requiring a computationally-heavy full grid inversion, discrepancies between model and data are spotted more easily, and the complementary strength of different types of data can be integrated into one framework.

Bayesian inversion using nested trans-dimensional Gaussian processes

SUMMARYTo understand earth processes, geoscientists infer subsurface earth properties such as electromagnetic resistivity or seismic velocity from surface observations of electromagnetic or seismic data. These properties are used to populate an earth model vector, and the spatial variation of properties across this vector sheds light on the underlying earth structure or physical phenomenon of interest, from groundwater aquifers to plate tectonics. However, to infer these properties the spatial characteristics of these properties need to be known in advance. Typically, assumptions are made about the length scales of earth properties, which are encoded a priori in a Bayesian probabilistic setting. In an optimization setting, appeals are made to promote model simplicity together with constraints which keep models close to a preferred model. All of these approaches are valid, though they can lead to unintended features in the resulting inferred geophysical models owing to inappropriate prior assumptions, constraints or even the nature of the solution basis functions. In this work it will be shown that in order to make accurate inferences about earth properties, inferences can first be made about the underlying length scales of these properties in a very general solution basis. From a mathematical point of view, these spatial characteristics of earth properties can be conveniently thought of as ‘properties’ of the earth properties. Thus, the same machinery used to infer earth properties can be used to infer their length scales. This can be thought of as an ‘infer to infer’ paradigm analogous to the ‘learning to learn’ paradigm which is now commonplace in the machine learning literature. However, it must be noted that (geophysical) inference is not the same as (machine) learning, though there are many common elements which allow for cross-pollination of useful ideas from one field to the other, as is shown here. A non-stationary trans-dimensional Gaussian Process (TDGP) is used to parametrize earth properties, and a multichannel stationary TDGP is used to parametrize the length scales associated with the earth property in question. Using non-stationary kernels, that is kernels with spatially variable length scales, models with sharp discontinuities can be represented within this framework. As GPs are multidimensional interpolators, the same theory and computer code can be used to solve geophysical problems in 1-D, 2-D and 3-D. This is demonstrated through a combination of 1-D and 2-D non-linear regression examples and a controlled source electromagnetic field example. The key difference between this and previous work using TDGP is generalized nested inference and the marginalization of prior length scales for better posterior subsurface property characterization.

Forward modelling of geophysical electromagnetic data on unstructured grids using an adaptive mimetic finite-difference method

While the mimetic finite-difference method shares many similarities with the finite-element and finite-volume methods, it has the advantage of naturally accommodating grids with arbitrary polyhedral elements. In this study, we use this attribute to develop an adaptive scheme for the solution of the geophysical electromagnetic modelling problem on unstructured grids. Starting with an initial tetrahedral grid, our mesh adaptivity implements an iterative h-refinement where a residual- and jump-based goal-oriented error estimator is used to mark a certain portion of the elements. The marked elements are decomposed into new tetrahedra by regular subdivision, creating an octree-like unstructured grid. Since arbitrary polyhedra are naturally permitted in the mimetic finite-difference method, the added nodes are not regarded as hanging nodes and hence any level of non-conformity can be implemented without a modification to the mimetic scheme. In this study, we consider 2-irregularity where two levels of non-conformity between the adjacent elements is permitted. We use a total field approach where the electric field is defined at the edges of the polyhedral elements and the electromagnetic source may have an arbitrary shape and location. The accuracy of the mimetic scheme and the effectiveness of the proposed mesh adaptivity are verified using benchmark and realistic examples that represent various magnetotelluric and controlled-source survey scenarios. The mesh adaptivity generates grids with refinement generally concentrated at the transmitter and receiver locations and the interfaces of materials with contrasting conductivities, and the mimetic finite-difference solutions have good agreement with the reference numerical and real data. We also demonstrate the practicality of our method using an example with an analytical solution and comparison with a standard mesh regeneration technique. The results show that our mesh adaptivity procedure can result in a higher accuracy, with similar numbers of elements, when compared with the mesh regeneration approach. Also, using a generic sparse direct solver, our method is found to be more efficient than the mesh regeneration approach in terms of computation time and memory usage. Moreover, a comparison between 1- and 2-irregularity shows the higher efficiency of the latter in terms of the number of elements required to reach a certain level of accuracy.

Petrov–Galerkin flux upwinding for mixed mimetic spectral elements, and its application to geophysical flow problems

Upwinded mass fluxes are described and analysed for advection operators discretised using mixed mimetic spectral elements. This involves a Petrov–Galerkin formulation by which the mass flux test functions are evaluated at downstream locations along velocity characteristics. As for the original mixed mimetic spectral element advection operator, the upwinded mass flux advection operator is conservative, however unlike the original advection operator, which is purely hyperbolic, the upwinded advection operator adds dissipation which is biased towards high wave numbers. The upwinded advection operator also removes the spectral gaps present in the dispersion relation for the original advection operator. As for the original advection operator, a material form advection operator may be constructed by similarly downwinding the trial functions of the tracer gradients. Both methods allow for the recovery of exact energy conservation for an incompressible flow field via skew-symmetric formulations. However these skew-symmetric formulations are once again purely hyperbolic operators which do not suppress oscillations. The scheme is implemented within a shallow water code on the sphere in order to diagnose and interpolate the potential vorticity. In the absence of other dissipation terms, it is shown to yield more coherent results for a standard test case of barotropic instability.

Nonextensive statistical mechanics for robust physical parameter estimation: the role of entropic index

The problem of inferencing parameters of complex systems from measured data has been extensively studied based on the inverse problem theory. Classically, an inverse problem is formulated as an optimisation task in which the objective function is based on the maximisation of Boltzmann–Gibbs entropy under appropriate constraints. In this classical framework, errors are assumed to follow a Gaussian-behaviour. However, in many situations, the errors are non-Gaussian and therefore the classical approach is predestined to failure. To mitigate the inverse-problem sensitivity to non-Gaussian errors, we have considered in this study a q-generalised objective function within the context of nonextensive statistical mechanics. In this regard, the errors are assumed to follow the q-Gaussian distribution, which arises from the maximisation of the nonadditive Tsallis entropy. We study the robustness properties of the q-objective function, which is a generalisation of the classical objective function, through the so-called influence function. In this work, we analyse and discuss from an analytical and numerical perspective, the role of the entropic index (q) of Tsallis entropy on the effectiveness and robustness in the inference of physical parameters from with strongly noisy-data. In this perspective, we show that exists an optimum value for the entropic index at the limit $$q \rightarrow 3$$ , which implies in $$\delta _q = 3 - q \approx 0$$ in the q-objective function. To validate our proposal, we consider a classic geophysical inverse problem by employing a realistic geological model and a data set contaminated by spikes.

Semi-Supervised Learning for Seismic Impedance Inversion Using Generative Adversarial Networks

Seismic impedance inversion is essential to characterize hydrocarbon reservoir and detect fluids in field of geophysics. However, it is nonlinear and ill-posed due to unknown seismic wavelet, observed data band limitation and noise, but it also requires a forward operator that characterizes physical relation between measured data and model parameters. Deep learning methods have been successfully applied to solve geophysical inversion problems recently. It can obtain results with higher resolution compared to traditional inversion methods, but its performance often not fully explored for the lack of adequate labeled data (i.e., well logs) in training process. To alleviate this problem, we propose a semi-supervised learning workflow based on generative adversarial network (GAN) for acoustic impedance inversion. The workflow contains three networks: a generator, a discriminator and a forward model. The training of the generator and discriminator are guided by well logs and constrained by unlabeled data via the forward model. The benchmark models Marmousi2, SEAM and a field data are used to demonstrate the performance of our method. Results show that impedance predicted by the presented method, due to making use of both labeled and unlabeled data, are better consistent with ground truth than that of conventional deep learning methods.

Linking bulk modulus to an unilateral damage yield criterion: A thermodynamic modeling approach

This work presents a new damage criterion suitable for elastic, elastic-plastic/viscous or elastic-viscous-plastic materials involving rupture effects. Its derivation, made here within a thermodynamic framework, follows previous scalar-valued damage mechanics approaches. Such approaches are appropriate to many geophysical problems involving quasi-brittle materials for which there is no clear physical justification for the level of complexity of a tensorial damage variable. Distinction between the mechanical response to compressive and tensile stresses is therefore not introduced by the damage itself but via a special definition of the Helmholtz free energy. This scheme differs from previous ones in that it combines with an evolution of Poisson’s ratio with the level of damage, which allows expressing the damage criterion in the principal stresses space. Moreover, there is no need to compute the stress eigensystem, which makes it simpler to implement than the Mohr–Coulomb damage criterion. Here we derive this damage criterion and compare it to observations of the variations of the bulk modulus in damaged geomaterials. We also compare it to in-situ stress measurements and find a good agreement in terms of the shape of the criterion in the stress space. We tentatively interpret the results in the context of previous studies of rock and ice mechanics.

Inference of thermodynamic state in the asthenosphere from anelastic properties, with applications to North American upper mantle

Inference of thermodynamic state and full-spectrum mechanical behavior of the lithosphere and asthenosphere is a central problem in geophysics, implicating our understanding of the convection patterns, transient responses and chemical composition of the planet. Anelasticity is responsible for significant relaxation of stress associated with seismic wave propagation in the asthenosphere, while irreversible transient creep may be important in the lithosphere. This paper focuses on the processes that may act at the time scales of seismic wave propagation, and current questions in the effort to determine the dependence of these effects on thermodynamic state. We introduce a free code library, the “Very Broadband Rheology calculator” (VBRc), designed to calculate frequency-dependent mechanical properties and easily compare different constitutive models favored by different laboratories. The methods operate only in the forward sense, starting with arrays of models of thermodynamic state, proceeding to arrays of mechanical properties. These calculations are incorporated into a Bayesian framework to infer variation in mantle thermodynamic state from Vs and Q, applied here to four locations in Western North America. The results demonstrate how well we can constrain the state, given the input models and the measurements. Results for sites in the Basin and Range, Colorado Plateau and interior craton east of the Rio Grande separate into distinct state variable ranges consistent with their tectonic environments.

An alternative amelioration of ME-based DZT-PML for truncating FDTD problems

ABSTRACT In order to further alleviate computational burdens and enhance accuracies, we propose an efficient amelioration on the complex-frequency-shifted perfectly matched layer (CFS- PML) via introducing both direct Z-transform (DZT) and matrix exponential (ME) methods for terminating the finite-difference time-domain (FDTD) region. The method is here referred as the DZT-ME-PML. Due to introducing the DZT method, the proposal can avoid the disadvantages of approximate bilinear and matched Z-transform in modeling fine and course structures, leading to achieving higher accuracies. Moreover, by employing the ME technique, the proposal develops into the compact first-order differential matrix form applied to evading the convolution manipulations, the variable replacements, and the formulation rearrangements during the mathematical derivations. Besides, the proposed DZT-ME-PML can not only flexibly terminate the arbitrary media when adding the DB constitutive relation, but efficiently attenuate evanescent waves, and reduce the late-time reflections due to the CFS scheme. Three-dimensional (3D) numerical examples have been carried out to illustrate the validity and accuracy of the proposed DZT-ME-PML. Eventually, applied to solving geophysical problem with the very low frequency (VLF) case, the DZT-ME-PML shows stronger absorption accuracy than the general ME-PML method.

Lp Norm Smooth Inversion of Magnetic Anomaly Based on Improved Adaptive Differential Evolution

Due to the approved applicability of differential evolution (DE) in geophysical problems, the algorithm has been widely concerned. The DE algorithms are mostly applied to solve the geophysical parametric estimation based on specific models, but they are rarely used in solving the physical property inverse problem of geophysical data. In this paper, an improved adaptive differential evolution is proposed to solve the lp norm magnetic inversion of 2D data, in which the perturbation direction in the mutation strategy is smoothed by using the moving average technique. Besides, a new way of updating the regularization coefficient is introduced to balance the effect of the model constraint adaptively. The inversion results of synthetic models demonstrate that the presented method can obtain a smoother solution and delineate the distributions of abnormal bodies better. In the field example of Zaohuoxi iron ore deposits in China, the reconstructed magnetic source distribution is in good agreement with the one inferred from drilling information. The result shows that the proposed method offers a valuable tool for magnetic anomaly inversion.

Efficient multiscale imaging of subsurface resistivity with uncertainty quantification using ensemble Kalman inversion

SUMMARYElectrical resistivity tomography (ERT) is widely used to image the Earth’s subsurface and has proven to be an extremely useful tool in application to hydrological problems. Conventional smoothness-constrained inversion of ERT data is efficient and robust, and consequently very popular. However, it does not resolve well sharp interfaces of a resistivity field and tends to reduce and smooth resistivity variations. These issues can be problematic in a range of hydrological or near-surface studies, for example mapping regolith-bedrock interfaces. While fully Bayesian approaches, such as those using Markov chain Monte Carlo sampling, can address the above issues, their very high computation cost makes them impractical for many applications. Ensemble Kalman inversion (EKI) offers a computationally efficient alternative by approximating the Bayesian posterior distribution in a derivative-free manner, which means only a relatively small number of ‘black-box’ model runs are required. Although common limitations for ensemble Kalman filter-type methods apply to EKI, it is both efficient and generally captures uncertainty patterns correctly. We propose the use of a new EKI-based framework for ERT which estimates a resistivity model and its uncertainty at a modest computational cost. Our EKI framework uses a level-set parametrization of the unknown resistivity to allow efficient estimation of discontinuous resistivity fields. Instead of estimating level-set parameters directly, we introduce a second step to characterize the spatial variability of the resistivity field and infer length scale hyperparameters directly. We demonstrate these features by applying the method to a series of synthetic and field examples. We also benchmark our results by comparing them to those obtained from standard smoothness-constrained inversion. Resultant resistivity images from EKI successfully capture arbitrarily shaped interfaces between resistivity zones and the inverted resistivities are close to the true values in synthetic cases. We highlight its readiness and applicability to similar problems in geophysics.

Analytical Atmospheric Ekman-Type Solutions with Height-Dependent Eddy Viscosities

We study the geophysical fluid dynamics problem in the steady atmospheric Ekman layer for height-dependent eddy viscosities in the form of some quadratic and rational power functions. We use the appropriate boundary conditions to obtain the wind velocity profiles analytically. For these Ekman-type solutions, the angle between the wind vector at the bottom and the geostrophic wind vector can be smaller or bigger than the traditional value of $$\frac{\pi }{4}$$ .

Solving the Seismic Inverse Problem Using an Ensemble of Models: A Case Study of Inversion of Receiver Functions

This paper proposes a new method for determining the seismic parameters of the subsurface below a seismic station using the waveforms of converted waves. The subsurface below the seismic station is represented by a plane-layered model. In contrast to other methods, the new approach does not involve determining all the elastic parameters of each of the layers. Instead, it is proposed to use a combination of parameters related directly to the object under study. In our case, we estimate the depth dependencies of the seismic velocity and density. The quality of solutions is assessed based on the rms distance between the observed waveforms and the seismograms derived from the averaged depth dependences of the seismic parameters. The approach proposed is applicable to a wide range of inverse problems in geophysics.

Influence of reservoir properties on the dynamics of a migrating current of carbon dioxide

Storage of carbon dioxide (CO2) in saline aquifers is a promising tool to stabilize the anthropogenic CO2 emissions. At the reservoir conditions, injected CO2 is buoyant with respect to the ambient fluid (brine) and spreads as a current laterally and toward the top cap rock of the aquifer, with the potential risk of a leakage into the upper aquifer layers. However, CO2 is partially soluble in brine and the resulting mixture (CO2 + brine) is denser than both starting fluids. This heavy mixture makes the configuration unstable, producing a convective flow that enhances the dissolution of CO2. Motivated by this geophysical problem, we analyze the influence of the porous medium properties on the evolution of a buoyant current that is weakly soluble with the ambient fluid. A time-dependent large-scale model [C. W. MacMinn et al., “Spreading and convective dissolution of carbon dioxide in vertically confined, horizontal aquifers,” Water Resour. Res. 48, W11516 (2012)] is used to analyze the evolution of the flow. In this work, we include additional physical effects to this model, and we investigate the role of horizontal confinement, anisotropy, and dispersion of the porous layer in the dynamics of the fluid injected. The effect of anisotropy and dispersion is accounted by changing the dissolution rate of CO2 in brine, which is obtained from experiments and Darcy simulations and represents a parameter for the model. Our results reveal that while the confinement has a remarkable effect on the long-term dynamics, i.e., on the lifetime of the current, anisotropic permeability and dispersion of the medium influence mainly the short-term behavior of the flow. Finally, we outline possible implications for the CO2 sequestration process.

Building ontologies for solving compute-intensive problems

The aim of the study, results of which are presented in this paper, is to analyse methods and tools for constructing an ontology related to solving compute-intensive problems and to form algorithms of its use. This problem arises from the need to solve problems using modern and future supercomputers, containing millions and, in the long term, billions of simultaneously operating computing cores and having a huge degree of parallelism. In solving such problems, the researcher should be well versed in both computational methods for solving the problem and modern supercomputer technologies, which is not always the case. One of the solutions to this problem is the creation of a knowledge base that includes ontological descriptions of methods for solving compute-intensive problems and architectures of supercomputers that can be used to solve them. The development of ontologies for a given subject area is one of the most important stage in creating an intelligent support system for solving specific compute-intensive problems. The paper discusses the methods and tools that are used to build the ontology. The paper also presents examples of the development of ontologies for astrophysics and geophysics problems.

Mathematical substantiation of pulsed electromagnetic soundings for new problems of petroleum geophysics

Mathematical substantiation of pulsed electromagnetic soundings for new problems of petroleum geophysics

Ближайшие перспективы развития геоинформационного комплекса INTEGRO

The article describes the main directions of development of the INTEGRA geoinformation complex, as well as geological, geophysical and information problems, for the solution of which the corresponding software modules were created. The problems of organizing a multi-user mode, designing software for building a 3D model of the territory, developing a block for storing and processing well information, and some others are discussed. The article also outlines the ways of further development of the software package. Keywords: multi-user mode, 3D model, montage method, parallel computing, data analysis, borehole information.

Dynamic Transitions and Bifurcations for a Class of Axisymmetric Geophysical Fluid Flow

In this article, we aim to classify the dynamic transitions and bifurcations for a family of axisymmetric geophysical fluid problems of a generic fourth-second order structure. A transition theorem is established by reducing the governing partial differential equations to a complex-valued ordinary differential equation, derived by employing approximate invariant manifolds. We develop an algorithm for the numerical determination of the transition/bifurcation types. Finally we apply the transition theorem and algorithm to examine the baroclinic instability in a two-layer quasi-geostrophic system in an annular channel and with different bathymetry profiles. Our numerical results show that with concave bathymetry the transition (bifurcation) is always continuous (supercritical Hopf bifurcation), whereas for convex bathymetry a jump transition (subcritical Hopf bifurcation) may occur in the basic azimuthal currents that rotate in the same direction.

Application of neural networks in problems of determining geometrical properties of objects placed in geophysical elastic media

The use of a neural network approach is associated with high requirements for software implementation. The main criteria is a high degree of accuracy in the recovery of the desired geological-physical model of the environment (GPME) and processing speed. The paper presents a description of the developed neural network architecture for solving the inverse problem of geophysics for determining the geometric properties of GPME objects. The performance of a single neural network and an ensemble of neural networks (NN) has been evaluated. The results are presented comparing the operating time of the NN when restoring models on various computing devices: CPU and GPU. The results of experiments on the restoration of various GPME using the developed NN based on the LSTM layer and the U-Net architecture are presented. This work was supported by the RFBR grant No. 19-07-00170.

The approximation neural-network method for solving nonlinear multi-criteria inverse problems of geophysics

A multi-criteria inverse problem is reduced to a system of operator equations on compact sets. It is shown that some a posteriori error estimates of solutions of the multi-criteria problem decrease with an increasing number of the criteria used. The approximation neural-network method for solving the multi-criteria inverse problem is presented. An example of a numerical solution of the two-criterion problem of geoelectrics is given, and the a posteriori error estimates are calculated.