Earth's Subsurface Research Articles

Deep Learning Assisted Extremely Low-Dimensional Representation of Subsurface Earth

Deep Learning Assisted Extremely Low-Dimensional Representation of Subsurface Earth

Forward modeling of magnetotellurics using Comsol Multiphysics

Magnetotellurics is an electromagnetic geophysical method that has been widely used to study structures in Earth's subsurface. Numerical modeling of magnetotellurics is important for survey design, inversion, geological interpretation and many other aspects of geophysical studies. For example, modeling a subsurface conductive body in terms of its conductivity, geometry and dipping angle would yield substantial information on the phase response and sensitivity in an MT survey. While there are many different modeling techniques, the finite element method is most commonly used. In this effort, we present magnetotelluric models of layered Earth, uplift structures, auroral electrojets, and geomagnetically induced currents in power-line skywires using the commercial finite-element package Comsol Multiphysics. The AC/DC module in Comsol can be used to solve Maxwell's equations in the quasi-static limit for modeling the magnetotelluric response. One of the advantages of Comsol modeling is its Graphical User Interface (GUI), which allows users to solve complex single or multi-physics problems in a meshed domain. The use of Comsol also reduces the need for sophisticated computer coding when solving partial differential equations such as Maxwell's equations. In the effort presented here, we first discuss model validation for layered Earth geometries. We then present two examples of magnetotellurics modeling in impact crater and geomagnetically induced current studies. Numerical results were compared with analytical solutions or benchmark results whenever possible.

Temperature and Pressure Dependence of Salt-Brine Dihedral Angles in the Subsurface.

Elevated temperature and pressure in the earth's subsurface alters the permeability of salt formations, due to changing properties of the salt-brine interface. Molecular dynamics (MD) simulations are used to investigate the mechanisms of temperature and pressure dependence of liquid-solid interfacial tensions of NaCl, KCl, and NaCl-KCl brines in contact with (100) salt surfaces. Salt-brine dihedral angles vary between 55 and 76° across the temperature (300-450 K) and pressure range (0-150 MPa) evaluated. Temperature-dependent brine composition results in elevated dihedral angles of 65-80°, which falls above the reported salt percolation threshold of 60°. Mixed NaCl-KCl brine compositions increased this effect. Elevated temperatures excluded dissolved Na+ ions from the interface, causing the strong temperature dependence of the liquid-solid interfacial tension and the resulting dihedral angle. Therefore, at higher temperature, pressure, and brine concentrations Na-Cl systems may underpredict the dihedral angle. Higher dihedral angles in more realistic mixed brine systems maintain low permeability of salt formations due to changes in the structure and energetics of the salt-brine interface.

A Geochemical Comparison of Three Terrestrial Sites of Serpentinization: The Tablelands, the Cedars, and Aqua de Ney

AbstractAlthough the Earth's subsurface hosts an abundance of microbial life, the influence of geochemistry on these communities remains poorly constrained. Ophiolites, sites where oceanic ultramafic minerals can be hydrated to serpentine minerals and metal oxides, create unique conditions capable of sustaining life. The fluid geochemistry of the Tablelands (NL, CAN), The Cedars (CA, USA), and Aqua de Ney (CA, USA), were studied to better characterize the range of fluid compositions observed at terrestrial sites of serpentinization. Fluids from these sites shared many commonalities including being ultra‐basic and reducing as well as having elevated levels of Cl−, Na+, K+, and Br− and depleted concentrations of Mg2+. They also exhibited a wide range of geochemistry. Isotopic and compositional data suggested the CH4 from The Cedars was a mixture of microbial and non‐microbial sources while the CH4 from the Tablelands was non‐microbial in origin. Aqua de Ney was the only site where the CH4 plotted in the abiogenic field. Despite being a known product of serpentinization, no H2 was detected at Aqua de Ney, likely due to the formation of abiogenic CH4 as well as the reaction of H2 and SO42− in the system to produce H2S. These unique sites of terrestrial serpentinization help to better understand the range of geochemistry at sites of serpentinization and its influence on the microbial communities in the subsurface.

Ions in the Deep Subsurface of Earth, Mars, and Icy Moons: Their Effects in Combination with Temperature and Pressure on tRNA-Ligand Binding.

The interactions of ligands with nucleic acids are central to numerous reactions in the biological cell. How such reactions are affected by harsh environmental conditions such as low temperatures, high pressures, and high concentrations of destructive ions is still largely unknown. To elucidate the ions’ role in shaping habitability in extraterrestrial environments and the deep subsurface of Earth with respect to fundamental biochemical processes, we investigated the effect of selected salts (MgCl2, MgSO4, and Mg(ClO4)2) and high hydrostatic pressure (relevant for the subsurface of that planet) on the complex formation between tRNA and the ligand ThT. The results show that Mg2+ salts reduce the binding tendency of ThT to tRNA. This effect is largely due to the interaction of ThT with the salt anions, which leads to a strong decrease in the activity of the ligand. However, at mM concentrations, binding is still favored. The ions alter the thermodynamics of binding, rendering complex formation that is more entropy driven. Remarkably, the pressure favors ligand binding regardless of the type of salt. Although the binding constant is reduced, the harsh conditions in the subsurface of Earth, Mars, and icy moons do not necessarily preclude nucleic acid–ligand interactions of the type studied here.

Early Results of Time Domain Receiver Function Data Processing in Mt Merapi and Mt Merbabu

The receiver function method is a method to image the earth subsurface by utilizing Ps conversion waves that are formed due to the velocity boundary. In general, the receiver function method estimates depth of structures such as the mantle-crust boundary by deconvoluting the vertical component from the horizontal component. Typical receiver function data processing is done in the frequency domain where the deconvolution process can be seen as a division between two components. In this study, we tried to reprocess the data using a deconvolution technique in time domain, popularly known as iterative time-domain deconvolution. The principle of iterative time domain deconvolution consists of iterative cross-correlation between the horizontal and vertical component. We use data from the DOMERAPI seismic station network located in the vicinity of Mt Merapi and Mt Merbabu. Mt Merapi is one of the most active volcanoes in the world with frequent eruptions and located at the ring of fire chain volcano in Indonesia. Note that the previous receiver function study in this region showed complex signals at some stations that may be related to sediment at shallow sediment and possible layers of low velocity zone that interfering main signal for a crust-mantle boundary. Our current results show iterative time domain RFs have clearer and smoother signal than the frequency domain that help interpreting the waveform signals. We estimate a range of crust thickness between 26-31 km near Mt Merapi. However, we noticed that iterative time domain calculation requires longer computation time and input signal.

Machine Learning Tools for Fossil and Geothermal Energy Production and Carbon Geo-sequestration—a Step Towards Energy Digitization and Geoscientific Digitalization

Eighty percent of current world energy consumption is satisfied by subsurface resources. In future, billions of watts of electrical power will be generated from geothermal energy sources. Subsurface earth can store the energy produced from renewable sources, such as wind and solar, and could provide safe storage of contaminants and hazardous nuclear waste. Engineering of subsurface earth is critical for fossil energy production, geothermal energy production, and carbon geo-sequestration. However, the subsurface is opaque, inaccessible, and heterogenous with nano-scale to kilo-scale processes that limits our understanding (characterization) and constrains the efficient engineering of the complex subsurface earth resources. There has been rapid increase in sensor deployment, data acquisition, data storage, and data processing for purposes of better engineering and characterizing the subsurface earth. This has promoted large-scale deployment of data-driven methods, machine learning, and data analytics workflows. Subsurface data ranges from nano-scale to kilometer-scale passive as well as active measurements. Subsurface data is acquired in the form of physical fluid/solid samples, images, 3D scans, time-series data, waveforms, and depth-based multi-modal signals, to name a few. Subsurface data sense various physical phenomena, e.g., transport, chemical, mechanical, electrical, and thermal properties. Integration of such varied data sources being acquired at varying scales, rates, resolutions, and volumes mandates robust machine learning methods to better characterize and engineer the subsurface earth. Machine learning has shown to improve the efficiency, efficacy, and productivity of subsurface engineering and characterization efforts required to produce fossil/geothermal energy and to sequester carbon dioxide.

Opportunistic magnetotelluric transects from CSEM surveys in the Barents Sea

Summary Magnetotelluric (MT) data allow for electrical resistivity probing of the Earth's subsurface. Integration of resistivity models in passive margin studies could help disambiguate non-unique interpretations of crustal composition derived from seismic and potential field data, a recurrent issue in the distal domain. In this contribution, we present the first marine MT data in the Barents Sea, derived from industrial controlled-source electromagnetic (CSEM) surveys. We characterize data quality, dimensionality, depth penetration and elaborate an analysis strategy. The extensive MT database consists of 337 receivers located along 7 regional transects, emanating from ∼70,000 km2 of 3D CSEM surveys acquired for hydrocarbon exploration from 2007 to 2019. High-quality MT data are extracted for periods ranging from 0.5 s to 5000 s. The data show no apparent contamination by the active source nor effects related to large time-gaps in data collection and variable solar activity. Along receiver profiles, abrupt lateral variations of apparent resistivity and phase trends coincide with major structural boundaries and underline the geological information contained in the data. Dimensionality analysis reveals a dichotomy between the western domain of the SW Barents Sea, dominated by a single N-S electromagnetic strike, and the eastern domain, with a two-fold, period-dependent strike. 35 receivers show 3D distortion caused by nearby bathymetric slopes, evidenced by elevated skew values. We delineate geographical areas where the 2D assumption is tenable and lay the foundation for future MT modelling strategies in the SW Barents Sea. We performed 2D MT inversion along one of the regional transects, a ∼220 km-long, E-W profile encompassing a major structural high and sedimentary basin approaching the continent-ocean transition. The resistivity model reveals low crustal resistivity values (1–10 Ω.m) beneath the deep sedimentary basins, in marked contrast with high resistivity values (1000–5000 Ω.m) of the thick crystalline crust on the structural high. We interpret this abrupt lateral resistivity variation as a rapid transition from a thick, dry continental crust to a hyperextended and hydrated crustal domain. Integration of resistivity with seismic velocity, density and magnetic susceptibility models may further refine these structural models and the underlying tectonic processes in the SW Barents Sea margin. Our methodology is applicable globally where 3D CSEM surveys are acquired and has a large potential for harvesting new knowledge on the electrical resistivity properties of the lithosphere.

An illumination‐compensated Gaussian beam migration for enhancing subsalt imaging

ABSTRACTSeismic images of Earth's subsurface are essential for oil and gas exploration. Gaussian beam migration is popular for seismic imaging owing to its flexibility and efficiency of implementation. However, the practical use of classic Gaussian beam migration is limited for complicated structures such as subsalt because of poor seismic illumination in those areas. We propose an illumination‐compensated Gaussian beam migration under the framework of least‐squares migration, enhancing the subsalt imaging effectively. A novel scheme based on the Born modelling and migration of Gaussian beam is developed to estimate the point spreading function, with which the illumination compensation can be efficiently implemented in the local wavenumber domain. As the aim of least‐squares migration, the proposed illumination‐compensated Gaussian beam migration produces true‐amplitude images in subsalt areas, which facilitates the seismic interpretation of subsalt structures. The total computational cost of the proposed method includes one Born modelling process and two conventional Gaussian beam migrations, and thus it is much more efficient than the classic least‐squares migration which requires multiple iterations. Numerical examples have demonstrated the effectiveness of the adaptive strategies, manifesting applicable potential in oil and gas exploration of subsalt targets.

RBF-FD analysis of 2D time-domain acoustic wave propagation in heterogeneous media

Radial Basis Function-generated Finite Differences (RBF-FD) is a popular variant of local strong-form meshless methods that do not require a predefined connection between the nodes, making it easier to adapt node-distribution to the problem under consideration. This paper investigates an RBF-FD solution of time-domain acoustic wave propagation in the context of seismic modeling in the Earth's subsurface. Through a number of numerical tests, ranging from homogeneous to highly-heterogeneous velocity models including non-smooth irregular topography, we demonstrate that the present approach can be further generalized to solve large-scale seismic modeling and full waveform inversion problems in arbitrarily complex models enabling more robust interpretations of geophysical observations.

Evaluating the Earth Subsurface for Civil Engineering Site Characterization in Agege, Southwest Nigeria Using Integrated Geoelectrical and Multichannel Analysis of Surface Wave (MASW)

A geophysical investigation involving 1D Vertical Electrical Sounding (VES), 2D Electrical Resistivity Imaging (2D ERI) and Multichannel Analysis of Surface Wave (MASW) has been carried out at Agege, Lagos, Nigeria with a view to delineating the subsurface stratigraphy and locate some competent strata/stratum for founding civil engineering structures. Six (6) 200 m long traverses were established within the study area. Along these traverses, 2D ERI were carried out adopting Wenner electrode configuration. Vertical Electrical Sounding (VES) adopting Schlumberger electrode array were carried out at selected points along profiles 1, 2 and 3 to determine the lithological sequence at depth. MASW data also were acquired along traverses 1, 2 and 3. The data were processed and the result yielded interpretable 2D resistivity structure and geoelectrical parameters (layer resistivity, thicknesses and depth) from the VES. The interpreted VES results were used to generate geoelectric section while the MASW resulted in 2D velocity sections. Three subsoils including topsoil, clay and clayey sand/sand were delineated beneath the study area. The resistivity and thickness range of the layers are; topsoil (34.0-54.6 ohm-m, 0.9 – 1.7 m), clay (10.3 – 17.7 ohm-m, 8.9 – 12.3 m) and clayey sand/sand (48.9 – 323 ohm-m) while the S-wave velocity range for the subsoil falls between 40 – 500 m/sec.

Distribution of metallogenic zones of the Caucasus region originated as a result of the subduction of the lithosphere of the Tethys Paleo-Oceanic plate under the East-European Paleo-Continental plate

On the assumption that the 2D heat flux anomaly in the territory of Caucasus is provided by the convective heat transfer from the mantle wedge, the angle and velocity of subduction of the paleo-lithospheric plate Tethys (Mediterranean and East Black Sea paleo-lithospheric plates) under the East-European paleo-lithospheric plate (Caucasus region) are estimated numerically. Within the framework of the geodynamical model constructed the horizontal extent of the 2D heat flux anomaly observed in the rear of the Caucasus mountain belt corresponds to subduction velocity of ~ 40 mm per year which is close to that observed with the help of geodetic means. The upwelling convective mantle flows can transport the mantle calc-alkali magmas (containing metals) to the subsurface Earth’s crust layers, thus making the ore deposits to be attached to the upper lithospheric locations above the ascending flows of the Karig-Richter vortices.

Subsurface Investigation of Basement Structures Using Electrical Resistivity Methods at Zainawa Village Kano, Nigeria

Electrical Resistivity Methods involving Schlumberger Vertical Electrical Sounding (VES) and Wenner Electrical Profiling (EP) were carried out to map the Geological features of the earth subsurface in Zainawa Area of Kano State, Nigeria. Five profiles were established; consist of six (6) VES points at each profile. GEOPULSE resistivity meter (SAS 300) was used for the data acquisition. The field data obtained have been analyzed using computer software (IPI2win) which gives an automatic interpretation of the apparent resistivity. A maximum of three geoelectric subsurface layers were delineated from the VES master curves. The geoelectric section beneath the study area was composed of top soil (clayey-sandy and sandy-lateritic), weathered layer, partly weathered (fractured basement) and fresh basement. The resistivity value for the topsoil layer varies from 20 Ωm to 600 Ωm with thickness ranging from 0.5 to 7.2 m. The weathered basement has resistivity values ranging from 15 Ωm to 593 Ωm and thickness of between 2.75 to 33.04 m. The fractured basement has resistivity values ranging from 201 Ωm to 835 Ωm and thickness of between 11 to 20.4 m. The fresh basement (bedrock) has resistivity values ranging from 1161 Ωm to 3115 Ωm with infinite depth. The depth to basement map was produced to give a good picture of the basement topography within the study area. The depth to basement ranges from 11 m around VES 01 to 85 m around VES 25 m. The map also reveals linear structures (VES 05, 21, 22 and VES 23) which trends in the NE-SW direction. These structures suggest a basement depression at these points. However, the depth from the topsoil to the bedrock surface varies between 2.5 to 37.75 m.

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.

Deep Learning for Seismic Inverse Problems: Toward the Acceleration of Geophysical Analysis Workflows

Seismic inversion is a fundamental tool in geophysical analysis, providing a window into Earth. In particular, it enables the reconstruction of large-scale subsurface Earth models for hydrocarbon exploration, mining, earthquake analysis, shallow hazard assessment, and other geophysical tasks.

Reciprocity-gap misfit functional for distributed acoustic sensing, combining data from passive and active sources

Quantitative imaging of subsurface earth properties in elastic media is performed from distributed acoustic sensing data. A new misfit functional based upon the reciprocity gap is designed, taking crosscorrelations of displacement and strain, and these products further associate an observation with a simulation. In comparison with other misfit functionals, this functional has the advantage of only requiring little a priori information on the exciting sources. In particular, the misfit criterion enables the use of data from regional earthquakes (teleseismic events can be included as well), followed by exploration data to perform a multiresolution reconstruction. The data from regional earthquakes contain the low-frequency content that is missing in the exploration data, allowing for the recovery of the long spatial wavelength, even with very few sources. These data are used to build prior models for the subsequent reconstruction from the higher frequency exploration data. This results in the elastic full reciprocity-gap waveform inversion method, and we illustrate its performance with a pilot experiment for elastic isotropic reconstruction.

A data mining approach for improved interpretation of ERT inverted sections using the DBSCAN clustering algorithm

SUMMARY Geophysical imaging using the inversion procedure is a powerful tool for the exploration of the Earth's subsurface. However, the interpretation of inverted images can sometimes be difficult, due to the inherent limitations of existing inversion algorithms, which produce smoothed sections. In order to improve and automate the processing and interpretation of inverted geophysical models, we propose an approach inspired from data mining. We selected an algorithm known as DBSCAN (Density-Based Spatial Clustering of Applications with Noise) to perform clustering of inverted geophysical sections. The methodology relies on the automatic sorting and clustering of data. DBSCAN detects clusters in the inverted electrical resistivity values, with no prior knowledge of the number of clusters. This algorithm has the advantage of being defined by only two parameters: the neighbourhood of a point in the data space, and the minimum number of data points in this neighbourhood. We propose an objective procedure for the determination of these two parameters. The proof of concept described here is applied to simulated ERT (electrical resistivity tomography) sections, for the following three cases: two layers with a step, two layers with a rebound, and two layers with an anomaly embedded in the upper layer. To validate this approach, sensitivity studies were carried out on both of the above parameters, as well as to assess the influence of noise on the algorithm's performance. Finally, this methodology was tested on real field data. DBSCAN detects clusters in the inverted electrical resistivity models, and the former are then associated with various types of earth materials, thus allowing the structure of the prospected area to be determined. The proposed data-mining algorithm is shown to be effective, and to improve the interpretation of the inverted ERT sections. This new approach has considerable potential, as it can be applied to any geophysical data represented in the form of sections or maps.

The Weight of Cities: Urbanization Effects on Earth's Subsurface

AbstractAcross the world, people increasingly choose to live in cities. By 2050, 70% of Earth's population will live in large urban areas. Upon considering a large city, questions arise such as, how much does that weigh? What are its effects on the landscape? Does it cause measurable subsidence? Here I calculate the weight of San Francisco Bay region urbanization, where 7.75 million people live at, or near the coast. It is difficult to account for everything that is in a city. I assume that most of the weight is buildings and their contents, which allows the use of base outline and height data to approximate their mass, which is cumulatively 1.6·1012 kg. I build a series of finite element models to study effects of pressure exerted by the weight distribution. Within the elastic realm, I look at compression, flexure, isostatic compensation, stress change, dilatation, and fluid flow changes. Within the nonlinear realm I show example calculations of primary and secondary settlement of soils under load. The combined modeled subsidence from building loads is at least 5–80 mm, with the largest contributions coming from nonlinear settlement and creep in soils. A general result is closing of pore space and redirection of pore fluids. While the calculated subsidence of the Bay Area is relatively small compared with other sources of elevation change such as pumping and recharge of aquifers, all sources of subsidence are concerning given an expected 200–300 mm sea level rise at San Francisco by the year 2050.

On hydrogen wettability of basaltic rock

One of the major challenges pertaining to a full hydrogen economy is hydrogen storage in the Earth's subsurface formations, such as deep saline aquifers and depleted oil and gas reservoirs. Generally, sedimentary basins are considered for this purpose due to their potential for large-scale storage. However, basaltic formations are also in comparable proportions having voluminous volcanic rocks, vast expansions, and high thickness. We, thus, in this work, evaluate brine/hydrogen/rock wettability of the real Basalt rock samples from CarbFix site in Iceland. Wettability is one of the significant parameters that highly influences the spreading and migration of the gas (as seen for CO2 capillary trapping investigations), interfacial area, mass transfer, and mineralization. We carried out brine/gas/Basalt three-phase contact angle measurements for some gases (CO2, N2, and He) under storage conditions of high temperature (323 K) and pressures (5, 10, 15, and 20 MPa). The analogy drawn from developed correlation for H2 from the contact angle trend of He (density close to H2 density) and empirical correlations reveals that a strong water-wetting state of Basalt can be expected for H2. The high water-wetting state of Basalt in the presence of H2 at storage conditions may be suitable for capillary trapping of H2 for long term storage in Basaltic formation.

Pseudo 3D seismic using kriging interpolation

Seismic method is one of geophysical method that is used to determine subsurface condition of the earth. There are two types of seismic data, 2D seismic and 3D seismic. The advantages of using 3D seismic data is it has less uncertainty for imaging the earth subsurface. However, the cost to acquire 3D data is higher than producing 2D data, thus some company needs to invest more money to be able to get a better data quality. One of the ways to have 3D seismic data without having to spend a lot of money is by creating a 3D pseudo seismic data from 2D seismic data. The purpose of this research is to create a 3D pseudo seismic data by using a geostatistical method. The geostatistical method that is used in this research is Kriging interpolation. By the help of Petrel software, kriging interpolation is done to create the 3D seismic. The pseudo 3D seismic results show a good image such as difference in amplitude, also structural geology features of subsurface. In conclusion, the use of kriging interpolation in pseudo 3D seismic still needs more work and development so the results will be good for the next interpretation steps.