Building uncertainty-aware subsurface models with 3D magnetotelluric inversion

Magnetotelluric (MT) inversion algorithms are an essential tool in exploration geophysics because they provide us with resistivity models of the subsurface. Inverse problems are ill-posed and non-unique. Resistivity models produced by geologically unconstrained MT inversions are typically smooth, which means that the geological interpretations obtained based on the resistivity models will be ambiguous and uncertain. To solve the above problems, geophysicists have made many efforts over the years to introduce any prior information (such as geological structure information and physical property data) into models, which could help to better constrain the inverted model with the inversion algorithms. In the last few years, there has been renewed interest in machine learning techniques. Various machine learning methods have been applied to inversions. Using the fuzzy c-means (FCM) clustering method, geophysicists proposed new inversion algorithms that are capable of building petrophysical information into the inversion. Although the FCM has been used and advanced in previous work, geophysicists concluded that a priori information of the correct number and value of cluster centers is very important for a successful FCM inversion. In fact, it is difficult to obtain the appropriate clustering information in some geophysical survey areas. In this study, we present an effective way to build the petrophysical information for the MT inversion based on FCM clustering algorithm. When the actual petrophysical information is insufficient, considering the characteristics of MT data, we perform a one-dimensional blocky inversion that can clearly identify the interface with distinct electrical property contrast; then, we obtain the number and value of cluster centers from 1D blocky inversion for the further MT inversion with FCM clustering. The algorithm uses guided FCM clustering to improve the model within the iterative minimization during MT inversion. In our MT inversion method, we integrate the geophysical inversion and geological differentiation into a unified scheme, which interact and enhance each other. The resistivity models obtained from the inversion respect the geophysical and petrophysical data and are easy to geologically interpret. We tested the algorithm using two synthetic examples and a field data example.

Summary The Simulated Annealing (SA) algorithm has been widely used in geophysical data processing. However, the seri al SA algorithm is limited by low computing efficiency. In this paper, the annealing schedule has been improved in the SA algorithm, and a constrained inversion of magnetotelluric (MT) data using the parallel SA algorithm has been applied. The test of synthetic data shows that the improved Very Fast Simulated Annealing (VFSA) method is significantly better than the traditional method. Additionally, the computing efficiency has been improved using the parallel SA inversion algorithm. This method has also been applied in MT data processing for geothermal exploration and good geological results have been achieved.

An extensive magnetotelluric (MT) survey comprised of 85 sites has been acquired over the Darajat geothermal field in Indonesia to map the geothermal reservoir and the overlying clay cap. The rouged topography and the geometry of the margin of the clay cap makes the resistivity structure 3D at reservoir depth. Although 3D MT inversion is now commonly used in geothermal interpretations 1D and 2D MT inversions are still effective tools for a variety of tasks such as quality assurance. Lower dimensional inversion can also play two critical roles in determining and assessing the resistivity model derived by 3D inversion: 1) by providing a good starting model to reduce the computational cost of the 3D inversion, and 2) by providing a computationally feasible path to stochastic inversion of the data that provides realistic parameter’s standard deviations for use in assessing reliability of the resistivity model. Using a spatially constrained 1D stochastic inversion of the MT data, we investigate the common claim that 1D inversion can provide a pseudo 3D model which closely matches the 3D inversion for the overburden and clay cap layers. The discrepancy between the pseudo 3D and true 3D inverse models increases with depth, however the presence of the core resistive feature of the field is still indicated at approximately the same depth as found in the true 3D model. Analysis of the 1D model parameter probability density functions shows that layer thicknesses are better determined than layer resistivities.

The geoelectrical features of the Travale geothermal field (Italy), one of the most productive geothermal fields in the world, have been investigated by means of three-dimensional (3D) magnetotelluric (MT) data inversion. This study presents the first resistivity model of the Travale geothermal field derived from derivative-based 3D MT inversion. We analyzed MT data that have been acquired in Travale over the past decades in order to determine its geoelectrical dimensionality, directionality, and phase tensor properties. We selected data from 51 MT sites for 3D inversion. We carried out a number of 3D MT inversion tests by changing the type of data to be inverted, the inclusion of static-shift correction at some sites where new time-domain electromagnetic soundings (TDEM) were acquired, the grid rotation, as well as the starting model in order to assess the connection between the inversion model and the geology. The final 3D model herein presents deep elongated resistive bodies between the depths of 1.5 and 8 km. They are transverse to the Apennine structures and suggest a correlation with the strike-slip tectonics. Comparison with a seismic velocity model and well log data suggests a highly-fractured volume of rocks with vapor-dominated circulation. The outcome of this study provides new insights into the complex geothermal system of Travale.

Summary The large Kevitsa Ni-Cu-PGE (platinum group elements) deposit in Finland is part of a complex geoelectrical setting. The ore zone is electrically conductive compared to the high resistivity olivine pyroxenite host rock. The ore body and host are surrounded by thin sheet like high electrical conductive carbonaceous phyllites. Unconstrained 2D and 3D inversion of magnetotelluric (MT) data was performed using seven densely sampled MT transverses ranging in length from 2500m to 7000m. Two dimensional inversions included a total of 312 MT stations along the seven transverses, while 3D inversion was completed on a subset of 270 MT stations focused around the ore zone. The inversions were conducted on the Pawsey Centre, Magnus Cray XC30 supercomputer. The challenge for this site was to resolve the ore body in the presence of the carbonaceous phyllites using sparse transverses. We compare the outcome of 2D and 3D MT inversion and find that 3D inversion appears to provide superior definition of the ore zone. This research and the resulting electrical conductivity distributions will inform our future work on full 3D cooperative inversion of seismic and MT data over the Kevitsa site.

The Sumatran Fault runs from the southeast (SE) to the northwest (NW) of Sumatra Island, with the highest slip rates reaching about 3.0 cm per year in the northwestern part. There is a seismic gap along this fault, including the northern Aceh domain, which consists of the Aceh and Seulimeum fault segments. Previous studies have used various methods to investigate the Sumatran Fault system, including seismic, geoelectric, gravity anomaly, and magnetotellurics (MT). The MT method has proven advantageous as it can non-destructively image a wide range of depths. However, previous studies using the two-dimensional (2D) MT inversion did not represent realistic information of the subsurface conditions. Therefore, a three-dimensional (3D) MT data inversion was used in this study to obtain more realistic information about the resistivity structure of the Aceh and Seulimeum segments. The results confirmed that the Sumatran Fault is a strike-slip fault, with a relatively northwest (NW)–southeast (SE) direction of conductivity strike with an angle of S 71.61° E from Groom–Bailey decomposition of MT data. The 3D resistivity distribution model from 33 stations showed that the Aceh Fault Segment is 20–30 km away, while the Seulimeum Fault Segment is 55–60 km away based on the MT data. The results also indicated a creeping zone at a depth of 2 km beneath the Aceh Fault Segment. Different rock formations were identified beneath the fault system, with the western part of the Aceh Segment dominated by high-resistivity metamorphic rocks (150–1000 Ωm) from the Triassic–Cretaceous age. The zone between the Aceh and Seulimeum fault segments exhibited low resistivity, characterized by volcanic rocks (1–15 Ωm) from the Lam Teuba Volcanic Formation and the Indrapuri Formation. Beneath the eastern part of the Seulimeum Fault Segment was found to consist of low-resistivity quaternary volcanic rocks (1–15 Ωm) and high-resistivity andesite rocks (4.5 × 104–1.7 × 105 Ωm). These findings correlated well with the geological map.

Inversion is a fundamental step in magnetotelluric (MT) data routine analysis to retrieve a subsurface geoelectrical model that can be used to inform geological interpretations. To reduce the effect of non-uniqueness and local minimum trapping problems and improve calculation speeds, a data-driven mathematical method with a deep neural network was developed to estimate the subsurface resistivity. In this study, a deep learning (DL) inversion technique using a revised multi-head convolutional neural network (CNN) architecture was investigated for MT data analysis. We created synthetic datasets consisting of 100,000 random samples of resistivity layers to train the network’s parameters. The trained model was validated with independent noised datasets, and the predicted results displayed reasonable accuracy and reliability, which demonstrates the potential application of DL inversion for real-world MT data. The trained model was used to analyze MT data collected in the southwestern Athabasca Basin, Canada. The calculated results from the DL method displayed a detailed subsurface resistivity distribution compared to traditional iterative inversion. Since this approach can predict a resistivity model without multiple forward modeling operations after the CNN model is created, this framework is suitable to speed up the computation of multidimensional MT inversion for subsurface resistivity.

This research activity aims a developing a better understanding of the geophysical response of hydrothermally-altered rocks associated with volcanogenic massive sulphide (VMS) deposits. Integrated modelling of physical rock property, whole rock geochemistry and 3D seismic data from the Lalor VMS deposit, Snow Lake, Manitoba enhanced insight into the seismic responses of hydrothermally-altered rocks in the footwall of the massive sulphide ore zone. Multivariate regression of P-wave velocity and density on lithogeochemical proxies for protolith composition, mineralization and hydrothermal alteration suggests that hydrothermal alteration in the footwall of the Lalor deposit enhanced acoustic impedance contrast between mafic to felsic volcanic units. Compositional variations between mafic, intermediate and felsic volcanic rocks, as modelled by the log(Zr/TiO2) immobile element ratio, remains the predominant factor that controls seismic response, independent from the intensity of hydrothermal alteration. The results obtained from the multivariate regression model are further corroborated by seismic forward modelling experiments, yielding enhanced seismic reflectivity in the intensely-altered footwall when alteration proxies are included in the 3D seismic rock property models. In parallel, geophysical response modelling of magnetotelluric (MT) data will be developed based on forward and inverse modelling routines on tetrahedral meshes, which yield results that are compatible with the general geological structure of the deposit. Initial unconstrained 3D MT inversions show a conductive zone between 600 to 800 m depth corresponding to the Lalor deposit, as well as, a highly conductive zone between depths of 500 to 1500 m, about 0.5 to 1.2 km to the southwest of the deposit.

SUMMARYSimulation-based probabilistic inversions of 3-D magnetotelluric (MT) data are arguably the best option to deal with the nonlinearity and non-uniqueness of the MT problem. However, the computational cost associated with the modelling of 3-D MT data has so far precluded the community from adopting and/or pursuing full probabilistic inversions of large MT data sets. In this contribution, we present a novel and general inversion framework, driven by Markov Chain Monte Carlo (MCMC) algorithms, which combines (i) an efficient parallel-in-parallel structure to solve the 3-D forward problem, (ii) a reduced order technique to create fast and accurate surrogate models of the forward problem and (iii) adaptive strategies for both the MCMC algorithm and the surrogate model. In particular, and contrary to traditional implementations, the adaptation of the surrogate is integrated into the MCMC inversion. This circumvents the need of costly offline stages to build the surrogate and further increases the overall efficiency of the method. We demonstrate the feasibility and performance of our approach to invert for large-scale conductivity structures with two numerical examples using different parametrizations and dimensionalities. In both cases, we report staggering gains in computational efficiency compared to traditional MCMC implementations. Our method finally removes the main bottleneck of probabilistic inversions of 3-D MT data and opens up new opportunities for both stand-alone MT inversions and multi-observable joint inversions for the physical state of the Earth’s interior.

A prominent TE and TM mode split is observed in magnetotelluric (MT) data below 0.5 Hz collected in the Perth Basin over the Harvey Ridge, Western Australia. We investigate the causes of mode splitting and consider implications on inversion of the MT data to subsurface electrical conductivity distribution. Twenty-five broad-band MT stations were acquired and remote reference processing was completed to arrive at a data set located midway between the Darling Fault and the Indian Ocean. We used forward modelling to test our strong suspicion that the Indian Ocean, Darling fault and architecture of the Granitic Basement were indeed the major contributors to mode splitting that we observed. Forward modelling of synthetic data was completed for comparison with the Harvey MT data. We were surprised at the match between synthetic and field data given the simplicity of the forward model and the considerable lateral distance between the MT soundings and the Indian Ocean or Darling fault. We were then able to make significant improvements to the MT inversion outcome by introducing a large scale geo-electrical architecture as the seed model for inversion. Our work demonstrates that large scale geo-electrical contrasts at considerable lateral distance from an MT transect, or the target zone need to be systematically introduced to the inversion if a quality outcome is to be achieved.

ABSTRACTWe present results from the resolution and sensitivity analysis of 1D DC resistivity and IP sounding data using a non‐linear inversion. The inversion scheme uses a theoretically correct Metropolis–Gibbs' sampling technique and an approximate method using numerous models sampled by a global optimization algorithm called very fast simulated annealing (VFSA). VFSA has recently been found to be computationally efficient in several geophysical parameter estimation problems. Unlike conventional simulated annealing (SA), in VFSA the perturbations are generated from the model parameters according to a Cauchy‐like distribution whose shape changes with each iteration. This results in an algorithm that converges much faster than a standard SA. In the course of finding the optimal solution, VFSA samples several models from the search space. All these models can be used to obtain estimates of uncertainty in the derived solution. This method makes no assumptions about the shape of an a posteriori probability density function in the model space. Here, we carry out a VFSA‐based sensitivity analysis with several synthetic and field sounding data sets for resistivity and IP. The resolution capability of the VFSA algorithm as seen from the sensitivity analysis is satisfactory. The interpretation of VES and IP sounding data by VFSA, incorporating resolution, sensitivity and uncertainty of layer parameters, would generally be more useful than the conventional best‐fit techniques.

Summary Estimation of reservoir parameters is im portant for prediction and optimization of oil production. Time-lapse seismic data, commonly acquired in mature oil fields to identify profitable bypassed zones, can be used to estimate reservoir parameters by posing the history matching process as a joint inversion problem. Due to the size of the problem and existence of multiple minima (of the error function), only global optimization techniques, such as Very Fast Simulated Annealing (VFSA) have an opportunity to find the global minimum. However, even for the VFSA technique, the computational cost of the forward modeling (particularly, the reservoir simulation) can make the inversion process impractical in real applications if we consider that many inversions are required to estimate the solution uncertainty. To improve the convergence characteristics of the VFSA technique, we have developed a modification of this technique where a regulation of the local temperature associated with the model variables is driven by the local changes of the error function. We call this modification Local Thermal Regulation (LTR) and show that the time (or number of iterations) spent for convergence can be reduced by 25 to 50%, while maintaining the same degree of ergodicity of the conventional VFSA technique.

SUMMARYBecause the magnetotelluric (MT) method uses natural sources, the electric and magnetic fields recorded in the field acquisition are not directly used but usually converted into other MT response functions for interpretation such as inversion. Considering that inversion results are dependent on types of input data, it can be helpful to analyse different characteristics of MT response functions for inversion. In this study, we examine sensitivity patterns of MT response functions used commonly in MT inversion, which are the impedance tensor, apparent resistivity, phase, tipper, effective impedance and phase tensor; and investigate how their sensitivity patterns affect inversion results. We first describe overall tendencies of 3-D sensitivity patterns of the MT response functions, and then classify the MT response functions into six groups based on 2-D sensitivity patterns computed at the surface, which are briefly called ‘surface-sensitivity patterns’ in this study. The ’diagonal components of the impedance’ and ‘off-diagonal components of the phase tensor’, which have four petals-shaped surface-sensitivity patterns along the diagonal directions, belong to Group 1, and contribute to imaging 3-D subsurface structures from receivers installed evenly at the surface. Group 2 contains the ‘xy-components of the impedance, apparent resistivity and phase’ and ‘yy-component of the phase tensor’ whose surface-sensitivity patterns are linear in the y-axis. The ‘yx-components of the impedance, apparent resistivity and phase’ and ‘xx-component of the phase tensor’ that have strong linear surface-sensitivity patterns along the x-axis are classified into Group 3. The MT response functions of Groups 2 and 3 are useful for inversion of structures close to 2-D, whose strike extends along the y- and x-axes, respectively. Groups 4 and 5 include the ‘x- and y-components of tipper’ that possess linearly aligned two petals-shaped surface-sensitivity patterns in the x- and y-axes, respectively. The tipper can be helpful in imaging both 2-D and 3-D structures. The ‘effective impedance’ belongs to Group 6, whose surface-sensitivity patterns appear as a small circle. The surface-sensitivity patterns allow the effective impedance to have an advantage in interpretation of 1-D structures. By using several MT response functions for specific cases of 1-D, 2-D and 3-D interpretation of MT data, we investigate whether characteristics of the sensitivity patterns are reflected in modelling (simulating field data) and inversion results, and then suggest optimal MT response functions for those cases. In doing so, we show how to utilize the characteristics of the sensitivity patterns in inversion, and recommend the input MT response functions for inversion according to MT exploration situations. Our study provides basic information on similarities and differences of major MT response functions for inversion and insights on which MT response functions are suitable to increase the feasibility of MT inversion for different field situations based on the sensitivity patterns.

WSJointInv2D-MT-DCR: An efficient joint two-dimensional magnetotelluric and direct current resistivity inversion

Magnetotelluric data modeling was used for oil and gas exploration. The parameters are apparent resistivity and phase. The magnetotelluric data of this research were line JBS4 with 71 sounding sites and line JBS8 with 66 sounding sites. The magnetotelluric data was completed by remote reference data which used for reducing the noise so that good data could be gotten. The magnetotelluric data was proceed in stage as follow: time series selection, Fourier transform, robust processing, cross power selection, and 2 dimension inversion. The result of this research showed that magnetotelluric method can figured out the sub-surface model well. From the result of interpretation 2D inversion magnetotelluric and other geophysical data (gravity anomaly map), we can conclude that magnetotelluric modeling figured out 2 sub-basins. Those sub basins had thick sediment layer (4000 meter for Bantarkalong Basin and 3000 meter for Citanduy Basin) and the dominant resistivity value was 22-119 Ohm.meter.