OpenMP Implementation of a Novel Potential-Field-Data Source-Growth-Based Inversion Approach for 3D Salt Imaging in Deepwater Gulf of Mexico

Abstract

Potential-field-data imaging of complex geological features in deepwater salt-tectonic regions in the Gulf of Mexico remains an open active research field. There is still a lack of resolution in seismic imaging methods below and in the surroundings of allochthonous salt bodies. In this work, we present a novel three-dimensional potential-field-data simultaneous inversion method for imaging of salt features. This new approach incorporates a growth algorithm for source estimation, which progressively recovers geological structures by exploring a constrained parameter space; restrictions are posed from a priori geological knowledge of the study area. The algorithm is tested with synthetic data corresponding to a real complex salt-tectonic geological setting commonly found in exploration areas of deepwater Gulf of Mexico. Due to the huge amount of data involved in three-dimensional inversion of potential field data, the use of parallel computing techniques becomes mandatory. In this sense, to alleviate computational burden, an easy to implement parallelization strategy for the inversion scheme through OpenMP directives is presented. The methodology was applied to invert and integrate gravity, magnetic and full tensor gradient data of the study area.

1. Introduction

Imaging of salt-tectonic-related geological features is still a challenging research matter in geophysical exploration, particularly in the oil exploration industry. Subsalt exploration is fundamental as salt domes are geological structures that often act as hydrocarbon traps, due to their impermeability to gas and oil [1,2,3,4]. Although three-dimensional seismic imaging has experienced a quite significant development in recent years [5,6,7,8,9,10], it is well known that it still tends to fail in subsalt regions due to the high-density contrasts between salt structures and their surrounding-encapsulating-rocks, yielding varying degrees of uncertainty in estimation of base and highly dipping flanks of salt bodies [11]. In this context, potential field methods are often excellent supporting tools for geophysical subsalt exploration and can provide accurate estimations of depth distribution and geometry of allochthonous salt features [12,13]. In addition, integration of potential field methods as an aid to seismic-structural interpretation is useful for the construction of initial velocity and density models for seismic RTM and FWI migration processes [13,14].

Among the class of potential field methods for geophysical exploration are gravity, magnetics and, more recently, Full Tensor Gradient method (FTG) [15]. Although gravity provides information of anomalies in the vertical component of deeper structures due to its a large bandwidth, FTG allows a better estimation of shallow geological structures due to the short wavelength signals which are not present in gravity data. Since the gradiometric tensorGis symmetric and traceless, its independent components allow determination of the locations of the mass center (_Gxz,_Gyz), edges (_Gyy,_Gxx) and corners (_Gxy) of the anomaly;_Gzz provides the correct spatial position of the anomaly, above the geological mass center, as conventional gravity [16,17,18,19,20,21,22,23,24,25,26]. On the other side, magnetometry is sensitive to anomalies of parts-per-million, making of it a very sensitive tool of geophysical surveying [27].

Inversion of depth and geometry of geological structures through potential field methods is an ill-posed problem, yielding to non-unique non-linear solutions [28]. Under some practical considerations, the problem can be linearized, but this approach can yield to numerical instabilities and non-unique solutions [29,30]. Thus, optimization techniques for the non-linear problem must be considered to be an alternative. Optimization-based data inversion consists in the extensive search over the parameter space [31], in a trial and error basis. This process continues until the difference between observed and synthetic data is minimized. In the case of three-dimensional inversions, allocation of huge amounts of memory for the sensitivity matrices, and a very high computational burden are required to solve the large unstructured linear systems. Among the techniques to alleviate computational burden in linear optimization, there are: (1) the use of Fourier transforms for the matrix-vector multiplication [20,32], (2) the compression of the sensitivity matrix [19,33], (3) the computation of only the elements of the matrix that are relevant to a determined search space [20,34,35], (4) solution methods that do not solve the linear system, but search solutions only in a particular space [36,37], and (5) systematic search methods [38,39,40,41]. Many systematic search methods are based on the idea of placing a seed somewhere within the search space from which the algorithm begins to optimize a solution, like those presented in [40,42,43]. The inversion algorithm in [40] and further refined in [44], is a systematic search method to invert gravity data. In this algorithm, the distribution of the mass density is estimated by growing iteratively the solution with respect to a prismatic element per iteration from a null prisms distribution. In each iteration, a new prism is chosen from the total set of prisms which remains with a value of zero for density. Its incorporation to the estimation minimizes an objective function which is the misfit of real and simulated data plus the_L2 norm of the weighted distribution of contrasts of density [31,45].

In this work, we generalize the algorithm reported in [40,44] in two ways: first, we modify it to also invert magnetic and gradiometric data; and second, a modification is introduced to be able to invert several variables simultaneously. These modifications were introduced in a way in which the algorithm is easily parallelizable. In this sense, in order to optimize the data processing burden, a parallelization strategy based on OpenMP directives is cautiously implemented in order to take advantage of current powerful multi-threading (MT) computers with shared memory architectures. First, parallelization and inversion results were analyzed for synthetic models. Finally, the resulting algorithm was applied to a real geological setting of intense salt tectonics in the northern Gulf of Mexico. For such purpose, an inversion of gravity, magnetic and FTG data is performed over a geological frame including salt tectonics, to determine the top and base of salt features.

The structure of this paper is as follows: Section 2, devoted to the presentation of materials and methods, includes some theoretical basis of potential-field-data forward and inverse modelling, a detailed description of the aforementioned source-growth algorithm and its parallel implementation, as well as the setup of both two reference synthetic models (T and S models) and a specific deepwater complex geologic scenario involving intense salt tectonics. In Section 3, results of the integrated potential-field-data inversion, along with the obtained performance are shown for both the T and S models, as well as for the real geological setup. Finally, Section 4 presents a thorough discussion about the herein proposed algorithm as well as the obtained parallel performance, and it poses some final remarks as well as significant recommendations and guidelines for future work.

2. Materials and Methods 2.1. Theoretical Background

2.1.1. Potential Field Theory

Theα-th component of the gravitational field vectorGmeasured in an observation point P at a positionrdue to a source body with mass densityρis given by

_Gα(r)=γ_{∫∫∫source}ρ(^r′)α−^α′|r−^{r′ |3}d^v′,

whereα=x,y,z,^r′is a vector pointing to the differential volume elementd^v′of the source body andγ=6,67428×^{10−11 m3s−2 kg−1}is the gravitational constant. In geophysical surveys, the commonly used units for gravitational field are mGals (1 Gal =1×¹⁰⁻⁵m/s²).

As an evolution of the gravity method, FTG method (full tensor gradient method) directly measures the derivatives of the gravitational vector field, which lead to a second-rank (covariant) symmetric tensorG, whoseα, β-th component is given by

_Gαβ(r)=∂_Gβ(r)∂α=γ_{∫∫∫source}3α−^α′β−^β′−^r−r′2 _δαβ^r−r′5ρ^r′d^v′.

The gravity field gradient tensor also has the property that the sum of its trace is zero, satisfying Laplace equation_Gxx+_Gyy+_Gzz=0. In practice, some combinations of the components of the gradiometric tensor are also used, particularly_Guν , which aids to determine the limits of the source body in N-S and E-W directions [20,46],

_Guν=_Gxx−_Gyy2.

FTG provides a greater sensitivity to low wavelengths with respect to gravity data. In geophysics, gradiometric anomalies are commonly measured in Eötvös (Eo), where 1Eo=1×^{10−9 s−2}.

On the other side, geophysical magnetic anomalies are commonly measured in nanoTeslas (1 nT =¹⁰⁻⁹T). The magnetic fieldBat a point of spacerdue to a material with a magnetizationMis given by

B(r)=−_Cm∇_{∫∫∫source}M·^∇′1r−^r′d^v′,

where_Cm=¹⁰⁻⁷is a scale factor, the prime at each coordinate denotes the position of the differential volume elementd^v′of magnetized material, and the components ofM are referred with respect to the International Geomagnetic Reference Field (IGRF) [27].M=χH=χ(_H0+_Ha)≈χ_H0, whereχis the magnetic susceptibility of a linear, isotropic and homogeneous material,_H0is the auxiliary Earth’s magnetic field (_H0=_B0/_μ0, where_μ0is vacuum’s magnetic permeability), and_Ha is the auxiliary magnetic field generated by rocks, which is usually negligible for most magnetic rocks in Earth’s crust. Equation (4) also neglects remnant magnetization and self-demagnetization effects, which results in a weaker magnetization.

As the magnitude of the magnetic field, the Total Magnetic Intensity (TMI),_BTMI=|B|is acquired during geophysical surveys, it ought to be corrected by subtracting the regional magnetic fieldFin order to obtain the total field anomalyΔT, which yields the following magnetic response

ΔT(r)=−_Cmk∇_{∫∫∫source}M·∇1r−^r′d^v′,

where k is the unit vector which involves the direction cosines of the geomagnetic field.

2.1.2. Numerical Approximation for Forward Modelling

The superposition principle of magnetic and gravitational fields always allows the decomposition of any domain in a collection of rectangular prisms which, if small enough, might be regarded as possessing a constant mass density/magnetic susceptibility. Such principle guarantees that the total magnetic/gravitational response of the domain is just the sum of the responses of each prism in the set. Considering rectangular differential prisms of volumed^V′=d^x′d^y′d^z′ , with constant mass densities/magnetic susceptibilities, and if we set a reference Cartesian frame such that the sides of such prisms are always parallel to the coordinate planes, then each of such differential prisms has initial and final limits shown in Figure 1.

In the case of gravimetry, from Equation (1), each prism contribution to the observed value of the vertical component ofGin the observation pointr can be numerically approximated by [47]

_Gz=γρ∑i,j,k=12_μijkzkarctan_{xi yjzk rijk}−_xiln(_rijk+_yj)−_yjln(_rijk+_xi),

where_xi=_x0−_ξi, _yi=_y0−_ηi, _zi=_z0−_ζi,_rijk=_xi2+_yi2+_zi2and_μijk=^{(−1)i (−1)j (−1)k} [48].

On the other side, the independent components of the gradiometric tensor (Equation (2)) can be numerically computed by [48,49]

_Gxx=γρ∑i,j,k=12_μijkarctan_{yj zkxi rijk},_Gxy=γρ∑i,j,k=12_μijkln(_zk−_rijk),_Gyy=γρ∑i,j,k=12_μijkarctan_{xi zkyj rijk},_Gxz=γρ∑i,j,k=12_μijkln(_yj−_rijk),_Gzz=γρ∑i,j,k=12_μijkarctan_{xi yjzk rijk},_Gyz=γρ∑i,j,k=12_μijkln(_xi−_rijk),

where_xi=_x0−_ξi, _yi=_y0−_ηi, _zi=_z0−_ζi,_rijk=_xi2+_yi2+_zi2and_μijk=^{(−1)i (−1)j (−1)k}.

Finally, in order to numerically compute the total magnetic field anomaly produced in an observation pointP=(x,y,0)at the surfacez=0 by a prism with arbitrary constant magnetization, we follow the algorithm presented in [50,51] which consists in computing the magnetic anomaly due to an infinitely extended in-depth prism (Figure 2).

Equation (5) in the schema proposed by [50] takes the form

ΔT(x,y,0)=_{Cm Mp}[_α232log_r0−^x′_r0+^x′+_α132log_r0−^y′_r0+^y′−_α12log(_r0+_z1)−_{F^x M^x} ^{tan−1x′ y′x′2}+_{r0 z1}+_z12−_{F^y M^y} ^{tan−1x′ y′}_r02+_{r0 z1}+^x′2+_{F^z M^z} ^{tan−1x′ y′}_{r0 z1}]_{|^x′=ξ1^x′=ξ2 |^y′=η1^y′=η2},

where_Mp=M(_Mx^i^+_My^j^+_Mz^k^)is the magnetization vector of the prism,F^=(_Fx^i^+_Fy^j^+_Fz^k^)is a unit vector in the environment field direction,_r0=^(x′−x)2+^(y′−y)2+_z12is the distance from the observation point to the differential prism of volumedv, and_α12=_Mx^_Fy^+_My^_Fx^,_α13=_Mx^_Fz^+_Mz^_Fx^,_α23=_My^_Fz^+_Mz^_Fy^.

2.1.3. Forward and Inverse Modelling

Forward modelling refers to the computing of the field anomalies with respect to the physical parameters. Measurements_miand parameters_pjare related through functions of the form_mi=f(_pj)∀i=1,⋯,N;j=1,⋯,M , where N and M are the total number of observations and parameters, respectively. Considering the linearized form of f, which has shown to provide acceptable results [31,45], the problem takes the form

Δm=KΔp,

whereKis the Jacobian matrix and it is also known as sensitivity kernel, andΔm=^m(C)−^m(O)gives the distance between the observed and the computed values of the physical propertymin terms of the distance between the measured and the computational model parametersΔp=^p(C)−^p(O).

On the other hand, the inverse problem consists in finding the values of the M parameterspof the model that provides the observed values of the field measurementsmin the N nodes of the observations mesh. In practice, a good approximation is obtained by minimizing an error, cost, objective or misfit function e. In this work, we used_L2 as e because it provides a similar weight to all data as well as stability [31,45],

e=Δm=^{∑absm(O)−m(C)12}.

A suitable inversion procedure ought to be stable (convergent) and it is generally constructed following the steps in [28].

2.1.4. Bodies Growth Inversion Algorithm

The bodies growth algorithm is a non-linear inversion method which estimates the geometry of the anomalous source by the addition of previously known density/magnetic susceptibility contrast prisms to the model. The main advantages of this algorithm are that: (1) it admits irregular prisms ensembles, (2) it allows determination of the regional tendency of the observed data, (3) it does not limit the growth to determined region of the domain, (4) it allows fitting of the positive and negative contrasts, (5) it allows variations on the contrast’s values during the process, and (6) it takes into account the measurement’s errors. The initial model_p0 can be a set of empty (null contrasts) prisms, or some prisms can be initialized with non-null contrasts, so the anomalous body tends to grow in the periphery of the initial model [40,44].

The inversion process consists in an iterative procedure to find the adequate prisms so to assign to it a value of the physical property. In the k-th iteration of the algorithm, k prisms have been adjusted andM−kprisms remain unchanged. If the forward response of the model in the k-th iteration is given by^m[k], and if in the next iteration thej=M−k+1-th prism is updated, the forward response of the model in thek+1-th can be easily updated as

^m[k]=K_pλ,

where_pλis a zero vector except in itsλ-th position, where it has the fitted contrast value of the added prism. For an isolated geological structure, the difference between the computed forward response of the model in thek+1-th iteration and the observed (measured) values^m(O)is given by

Δ^m[k+1]=^m(O)−^{f[k+1] m[k+1]}−_dreg[k+1],

where^f[k+1]≥1is a scale factor used to fit^m(O)and^m[k+1], and_dreg[k+1]is the regional field at thek+1-th iteration, linearly approximated by

_dreg[k+1]=_c0[k+1]+_cx[k+1](x−_xM)+_cy[k+1](y−_yM),

where(_xM,_yM)are the coordinates of the mid (central) point of the domain, and_c0,_cx,_cys are coefficients to fit the regional tendency. Thus, the coefficients to be optimized are_fk,_c0k,_cxk,_cyk. For gravimetric and magnetic data, the cost or objective function to fit such coefficients is given by

ϕ^p[K+1]=e^p[k+1],λ2=Δ^m[k+1]T _Qm−1Δ^m[k+1]+λ^{f[k+1]2 p[k+1]T} _Qp−1 ^p[k+1],

where T stands for transpose,_Qmis the measurements (response) covariance matrix (N×N) and_Qp=^KT _Qd−1Kis the covariance matrix (M×M ) parameter. For gradiometric data, the total cost function is the sum of the cost functions of each component of the tensor being inverted. The first term in Equation (14) represents the misfit of the response normalized by_Qm, the second term is the sum of the parameter vector weighted by theλparameter (positive), which value should be chosen carefully. Low values ofλprovide good fitting of data, but the model may grow excessively, creating fictitious structures, while high values ofλyield simple and shallow models that poorly adjust measurements. The coefficients to be optimized can be computed as

_cy[k+1]=_{Gxx Gyd}−_{Gxy Gxd}/_{Gyy Gxx}−_{Gxy Gxy},

_cx[k+1]=_Gxd−_{cy[k+1] Gxy}/_Gxx,

_c0[k+1]=_Fud−_{cx[k+1] Fux}−_{cy[k+1] Fuy}/_Fuu,

^f[k+1]=_srd−_{c0[k+1] sru}−_{cx[k+1] srx}−_{cy[k+1] sry}/_srr+_spp,

where_spp=^pT _Qp−1p, the subindexes r, d, u, x and y refer to the vectorsr=^m[k],d=T,u=1,x=x−_xMandy=y−_yM, and where the coefficients_srg,_sru,_srx,_sry,_srr,_smm,_Fug,_Fux,_Fuy,_Fuu,_Gxg,_Gxy,_Gxx,_Gygand_Gyyare computed as

_Gab=_{Fuu Fab}−_{Fua Fub},_Fab=_sab(_srr+_smm)−_{sar sbr},_sab=^aT _QD−1b,

where a and b can take the values of the subindexes r, g, u, x and y. In each step,^f[k]is diminished, and the inversion process stops when^f[k]≈1and the inverted model is accepted.

If the parameters (contrasts) are changed during the inversion process, the election of excessively high contrasts shall fit the observed anomaly, yielding to very simple models, while very small contrasts would fit not only the real anomalous bodies but also the observed noise, growing fictitious bodies in the model. In this way, a criterion for the variation of contrasts along the iterative process is necessary. Taking the maximum possible (positive or negative) contrast_pmax, it will evolve according to

Δp=_pmax1−1^f+0.1τ,

wheref≥1is a scale factor used to describe the growth instant, andτis the smoothing parameter that indicates the desired contrast’s variations. High values ofτproduce homogeneous models with angulate geometries and contrasts close to_pmax, while low values yield models with very varied contrasts which decrease to the exterior of the grown body, with smoothed geometries.

2.2. Parallel Implementation of Forward Modelling It is a well-known fact that geophysical inversion problems imply a very high computational burden, which exponentially increases with the size of the geological setup under study. Several strategies to accelerate geophysical inversion algorithms in very diverse computing architectures have been developed. The selection of a particular approach to parallelize an inversion scheme heavily relies in both the nature of the particular inversion algorithm to be implemented as well as the available computing architecture in which the algorithm should be executed, and the programming paradigm.

With respect to the computing architecture, it is well known that a deep acquaintance of the architecture leads to a very efficient parallel solution to the problem, which in turn becomes complicated because it requires programming skills at low level and a large development (coding) time. This is the case of using CUDA C for Graphics Processing Units (GPU) or Message Passing Interface (MPI) for memory distributed architectures like clusters or Xeon Phi Co-processors [48,52,53]. As one of the goals of this paper is to make this research accessible to the greatest possible part of the geophysical and related communities, the acceleration of the herein presented algorithms are based on shared memory CPU-based architectures.

In this sense, shared processing architectures, as multi-socket-multi-core machines, are more common each day, even in home personal computers. To increase the computing power, the trend during the last 20 years has been to increase the number of integrated cores in the same chip. This is the case, for instance, of Intel Xeon W-3175X^®, with 28 physical cores able to manage up to 56 processing threads, as well as the AMD Ryzen Threadripper 3990X^®, with 64 physical cores, able to manage 128 execution threads. These vertiginous technological advancements in shared memory architectures should be exploited to solve problems of high computational burden as those of geophysics. In this sense, it is necessary to use programming paradigms which allow the community to center the attention in the geophysical problem rather than in the low-level programming implementation of the algorithm. In this sense, OpenMP has practically become a standard for the parallel programming of Symmetric Multi-Processing (SMP) systems in order to execute a code without the necessity to write low-level instructions, as POSIX Threads (pthreads).

OpenMP is a set of functions and directives for C, FORTRAN and C++, which delegates the generation of the low-level parallel instructions to the compiler. Although OpenMP alleviates the code development burden for the programmer, he is still in charge of the design of the parallel algorithm. Thus, a deficient design of the parallel strategy is prone to a deficient performance. OpenMP has benefited several geophysical problems [48,52]. OpenMP also provides an implicit parallelism model which produces medium granularity tasks for MT applications, providing a higher computational abstraction than MPI. It is highly portable between platforms and its ease of usage allows significant parallelism with some directives inserted in the serial code, diminishing (sometimes dramatically) the development costs of parallelization [54]. In exchange of these features, OpenMP does not provide code execution speed-ups as those provided by a carefully parallel implementation of the algorithm in MPI.

The inversion of potential data by the growth of bodies is an algorithm which explores different individual solutions (prisms) that gradually diminish the misfit between observed and computed data, finally yielding the best configuration of models in the search space. The main time drawback of this algorithm is that it fits just one prism per iteration. The algorithm’s parallelization steps are:

1 OpenMP performs the division of the total search space by the available processing threads. This operation is carried out in an automatic way according to an OpenMP schedule; in general, the schedule selected by the compiler is optimal.

2 Every thread independently finds the best prism in its own search subspace.

3 The algorithm finds the prism which reduces the most the misfit between observed and computed data, which is selected as the best prism in the parallel iteration.

4 This process is repeated until a maximum number of iterations is reached.

This process is illustrated in Figure 3, and the design of the algorithm is performed using Foster’s methodology [55]: Partitioning, Communication, Agglomeration, Mapping. Once this cycle is finished, this optimal prism is taken as part of the grown body, restarting the cycle with updated search subspaces for each thread, considering the last updated prism as not modifiable.

Performance Tests

Performance tests of the parallelization scheme were performed in both the synthetic and real data models. We present performance results for gravimetric, magnetic and gradiometric inversions, considering the prisms and observations mesh configurations taken in the previous sections. The technical specifications of the workstation where the tests were performed are:

• Two X5650 Intel Xeon (R) processors with six 2.66 GHz physical cores per socket (12 physical cores), 24 logical cores with HyperThreading (HT) activated,
• 48 GB of RAM memory,
• Red Hat 7.7 Operating System,
• GNU FORTRAN Devtoolset-6.

In the field of SMP parallel computing, thread affinity restricts execution of certain threads to a subset of the physical processing units in a multiprocessor computer. Depending upon the topology of the machine, thread affinity is always important in terms of performance, particularly when the HT technology is enabled, which basically means that one CPU physical core can attend two threads at the same time, sharing the same Arithmetic-Logical Unit (ALU). Two common affinity setups for the workstation with two sockets and the HT enabled used in this research, are depicted in the Figure 4.

For the performance tests, threads from 1 to 24 were spawned to cover the total number of logical processors; creating more than 24 threads would have consequences as a drop in performance.

The results of the tests show that the distribution of the threads using the scatter affinity configuration is slightly better; in this configuration the threads are assigned first to a physical core and then to the corresponding logical cores in a consecutive way up to reaching (in this case) 24 threads (see Figure 4b). In such way, even though we tested both configurations, from here on we only show performance results for the inversion processes obtained with a scatter affinity setup.

With HT technology enabled, it is necessary to create two threads and assign them to the same core to obtain a performance gain. According to the scatter affinity, the threads 0 and 12 are assigned to the same physical core, sharing the resources of the core, but switching between them very quickly. 2.3. Synthetic Test Models Configuration

To test the proposed parallel algorithm for inversion of geophysical data, in this section we provide the characteristics of two regular well-known synthetic models: T and S models, which can be observed in Figure 5.

The parameters used to compute field anomalies for each model can be observed in Table 1.

Observations mesh for the two models consists in a grid of equispaced nodes of 180 [m] in X and Y directions, with a total of 2601 nodes (51 in each direction ), considering a Cartesian reference system with X, Y and Z in the east, north and in-depth directions, respectively. The mesh is in theZ=0plane. The magnetic anomalies were computed at the pole (I = 90°).

For the inversion of synthetic models, we considered an initial model of null (density or magnetic susceptibility) contrasts. Thus, the algorithm was set to test only positive ones, not allowing prisms with negative contrasts. The search space was constrained to the limits shown in Figure 5, and consists in a set of 31,104 prisms.

For the inversion of gradiometric data for synthetic models, in order to speed up and guarantee convergence of the results, we simultaneously invert the following three components of the gradiometric tensor (_Guν,_Gxy,_Gzz); knowing_Guνand_Gzz,_Gxxand_Gyycan be computed considering thatGis traceless. This choice also diminishes the required RAM memory for the three used sensitivity matrices. For the results of the inversion for both models, the initial response is shown in the color map and the inverted anomalies for_Gxx,_Gyy,_Gzzand_Gxy, are shown in black contour lines.

2.4. Geological Salt Dome Setup

To further evaluate the performance of the herein presented algorithm, it was tested on a real-case scenario. The region under study corresponds to a real geological formation at the southern Gulf of Mexico, which features a salt body of a previously determined thickness, so the depth of its top and base is known a priori. The initial geologically feasible model was built with data from stratigraphic columns, interpreted images from geological horizons (up to the top of basement) and the sediment compaction curve of the Gulf of Mexico which is shown in Figure 6 [36]. The exact location and details of the zone under study cannot be disclosed in this paper due to business confidentiality.

2.4.1. Initial Model Configuration

An interpolation algorithm was applied to the XYZ dataset representing the top of each horizon where necessary to build equispaced meshes, in order to guarantee that the obtained meshes per horizon are mutually coherent with the identified tops, and that they were not morphologically modified. The interpolated geological horizons can be observed in Figure 7. Furthermore, to avoid an overlapping of points, the lithological sequence of the stratigraphic column is followed. The size of the domain is 40.5 × 36.5 × 14.1 [km] in the east-west, north-south and depth directions, respectively, which also considers a 5 [km] region in each extreme to keep border effects out of the region of interest.

The geometry of the model to compute the field anomalies consists in a prisms ensemble of 500 × 500 × 50 [m] in the X, Y and Z directions respectively, forming 1,667,466 prisms. The observations mesh consists in 35 × 31 nodes in the X and Y direction (1085 total nodes), at a height of 100 [m] above mean sea level.

The density range for the model is between 1300 [kg/m³] and 2700 [kg/m³ ] from the top to the basement. Constant density layers were adjusted from geological horizons shown in Figure 7. A graphical representation can be observed in Figure 8. In this work, we only consider the large salt body emplaced in the center of the domain, with respect to other bodies in the zone [52].

For the magnetic properties of the model, we considered that there is no variation of the susceptibility with depth. For the values of magnetic susceptibility, we considered a mean value for surrounding sediments of350×¹⁰⁻⁶ [1], while for the salt a value of−10×¹⁰⁻⁶ [1,56]. Geomagnetic vector and media magnetization were taken vertical and parallel with respect to the pole. Gravity, FTG and magnetic data were especially acquired through airplane for this research work.

2.4.2. Forward Gravimetric Modelling

With the synthetic model built from the real data (Figure 8), we proceeded to compute its forward (air-free) gravimetric response, which can be observed in Figure 9a. To remove the regional field effect so to isolate the salt body of interest, and considering that the geology around the salt body is known, the prisms pertaining to the sediments surrounding the salt body are considered with a null density contrast. The gravimetric response of the isolated salt body can be observed in Figure 9b.

2.4.3. Forward Gradiometric Modelling

The forward (air-free) gradiometric response of the real synthetic model (Figure 8) can be observed in Figure 10. Following the same strategy to remove the regional field, the gravimetric response of the isolated salt body can be observed in Figure 11.

2.4.4. Forward Magnetic Modelling

Finally, the forward magnetic response of the isolated salt body can be observed in Figure 12.

3. Results 3.1. Inversion of Synthetic Models

3.1.1. T Model

Table 2 shows the configuration of the inversion performed for the T model. In Figure 13, the initial and inverted gravimetric vertical anomalies corresponding to the T model can be observed, while Figure 14 shows the obtained (inverted) configuration of prisms.

In the case of magnetic data, Table 3 shows the configuration for the inversion performed to the initial T model. In Figure 15, the initial and inverted magnetic anomalies can be observed, while Figure 16 shows different perspectives of the inverted (with magnetic data) prisms configuration for the T model.

In Table 4, the configuration of the inversion performed to the T model is shown. Figure 17 shows the initial and inverted gradiometric anomalies, while Figure 18 shows the inverted (gradiometric) prisms configuration for the T model.

3.1.2. S Model

Table 5 shows the configuration of the gravimetric inversion performed for the S model. In Figure 19, the initial and inverted gravimetric vertical anomalies corresponding to the S model can be observed, while Figure 20 shows the obtained (inverted) configuration of prisms for the S model.

In the case of magnetic data, Table 6 shows the configuration for the inversion performed to the initial S model. In Figure 21, the initial and inverted magnetic anomalies can be observed, while Figure 22 shows different perspectives of the inverted magnetic prisms configuration for the S model.

In Table 7, the configuration of the gradiometric inversion performed to the S model is shown. Figure 23 shows the initial and inverted gradiometric anomalies, while Figure 24 shows the inverted (gradiometric) prisms configuration for the S model.

With respect to the performance, we measured the time obtained for such threads to assess the performance of the parallel code for the three inversion cases: gravimetric, gradiometric and magnetic. Figure 25 and Figure 26 show the computing time and speed-up results for the T and S models, respectively.

3.2. Inversion of Real Salt Structures in the Gulf of Mexico

For the real data-based model, the strategy of beginning with a null contrasts model that was applied to synthetic S and T models is changed. In this case, as the salt body dimensions are previously known, the initial search space is configured with a rectangular block of 24.5 × 10.5 × 6.2 [km] in the X, Y and Z directions respectively, representing the lateral and vertical limits of the salt body. The block is divided in prisms of constant (density or susceptibility) contrast of 250 [m] × 250 [m] × 100 [m] in the X, Y and Z directions respectively. The initial model for the inversion can be observed in Figure 27. This yields a total number of prisms of 72,269. The observations mesh remains unchanged from the one used in the previous section.

The growth algorithm is limited to test parameters pertaining to prisms below the interpolated horizon corresponding to the salt top, and without restrictions in depth.

3.2.1. Gravimetric Inversion of Real Data

Table 8 shows the configuration of the inversion performed to the salt body. In Figure 28, its initial and inverted gravimetric vertical anomalies can be observed, while Figure 29 shows the final configuration of prisms after the inversion.

3.2.2. Magnetometric Inversion of Real Data

In the case of magnetic data, Table 9 shows the configuration for the inversion performed to the salt body. In Figure 30, the initial and inverted magnetic anomalies can be observed, while Figure 31 shows different perspectives of the final prisms configuration after magnetic data inversion for the salt body.

3.2.3. Gradiometric Inversion of Real Data

Finally, in Table 10, the configuration of the inversion performed to the salt body is shown. Figure 32 shows the initial and inverted gradiometric anomalies, while Figure 33 shows the prisms configuration for the salt body after the gradiometric data inversion.

With respect to the performance of the parallel algorithm for the salt structure, execution times and their corresponding speed-up values for 1 to 24 threads are shown in Figure 34, for the three inversion cases: gravimetric, gradiometric and magnetic.

4. Discussion 4.1. Inversion of Synthetic Models From the exercise of inverting synthetic models, it can be observed that from the point of view of the observed and inverted field responses, there is a good agreement; the most notorious differences appear in the gradiometric dataset. From the point of view of the prisms configuration, the agreement between initial and inverted models is particularly notable in the shallow parts of the models. The sum of two terms in the objective function, along with the bodies growth algorithm, yields compact models which acceptably define the top of the models. The obtained prisms model with the joint inversion of gradiometric datasets tends to be more compact and better laterally determined than the ones inverted from magnetic or gravity data. This is likely due to the dependence of gradiometric data with the cube of the distance from the source to the observation point, which allows better definition of shallow structures. As the growth algorithm was unrestricted in the synthetic cases, these models show the necessity of favoring the growth in the vertical direction in order for the algorithm to explore deeper regions of the search space. For real datasets, these constraints must come from geologic/geophysical information of the region under study.

As previously mentioned, the smoothing factor(τ≥0) is a parameter which allows a priori the fitting of the variation of the density contrast during the inversion process. Figure 35 shows the behavior of the inversion process as a function of different values of smoothing factorτ for the T model. As it can be observed in Figure 35a, a small value ofτproduces an ensemble with more prisms with a lower density, as well as a greater scattering of the elements. On the other side, large values of the coefficientτ produce greater density contrasts, yielding more compact ensembles with few prisms out of the main ensemble, with grater values of density, as it can be observed in Figure 35c. There is not an a priori method to precisely determine the smoothing factor, which is generally done on a trial and error basis. In this work we chose the valueτ=8 (see Figure 35b), because after several experiments, this value fitted the best the observed values of the fields.

4.2. Performance Results for Synthetic Models

According to the obtained results shown in Figure 25a, execution times were reduced 10.73×, 11.72× and 10.21× for the gravimetric, magnetometric and gradiometric inversions, respectively, in the best case, i.e., for 12 threads. This is a satisfactory result considering that the workstation consists in 12 physical cores. Once the number of physical cores is reached, variations in performance are introduced in the logical cores. This result is common because the threads are sharing the physical cores [48,57]. Furthermore, with HT enabled, a small gain of 0.46×, 2.38× and 0.33× in the computing times is observed for each inversion process for 24 spawned threads, with respect to the results of 12 threads. It is important to point out that search algorithms by their nature have time variations in their execution, because searches are not deterministic. In this sense, the results presented in Figure 25a are the result of a statistical average of 10 different executions under the same circumstances. As a general result, the execution time is reduced around 12× when executing it with 24 spawned threads.

The maximum speed-up values for the T model were 10.73, 11.72 and 10.21 with 24 threads, for gravimetric, magnetic and gradiometric inversions, respectively (see Figure 25b). Considering that the obtained speed-up values with 12 threads are 10.32, 11.39 and 9.91, the difference in performance between 24 and 12 threads is not significant. It can even be observed from Figure 25b that after 12 threads, a significant drop in the performance occurs up to 14 threads. After the thread 15, the performance is recovered slowly. This behavior is common with HT enabled, and it shows to bring an oscillating behavior in this case, although in certain applications and circumstances it can achieve up to 25% in performance [48,58,59]. For this case, only with 24 threads, the performance slightly improves in about 4% with respect to the results with 12 threads.

On the other hand, it is important to note that for 1 to 12 threads, the general behavior of the OpenMP parallelization for the bodies growth algorithm is close to the perfect speed-up (notably, for the gravimetric inversion). This fact demonstrates that the herein presented OpenMP parallelization strategy of the algorithm is well developed and coded, so the remaining fraction of serial code is minimum.

For the S model case, 1 to 24 threads were also spawned for the three inversion cases: gravimetric, gradiometric and magnetic data inversion. Such performance results for the S model are presented in Figure 26, where it can be observed that the results are analogous to those of the T model (Figure 25). We can observe the same behavior in computing times, providing time reductions of 11.26×, 11.10× and 11.05×, respectively, for the 24 threads case (see Figure 26a). 24 threads case yields an improvement of only 4.8% with respect to the results of 12 threads in the best case (magnetometric inversion). On the other side, the speed-up results for the S model are remarkable, in the sense that the three inversion cases are very close to the ideal theoretic speedup. This is likely to the fact that a configuration of prisms that optimized the misfit was found in more consistent execution times for the S model than for the T model.

4.3. Performance Results for the Real Salt Structure Inversion

For the real data-based inversion of the considered salt body, it is confirmed that the final configuration obtained from the joint inversion of the tensor gradiometric data is more compact with respect to the configurations obtained by gravimetric/magnetic data, defining its limits in a better way. This is likely due to the fact that the variation of densities corresponding to the smoothing factorτis less considerable in the gradiometric inverted model, concentrating these variations in the borders of the salt body to a greater extent.

As can be observed, the results are analogous to those obtained for the synthetics models. Again, as a search-based parallel algorithm is not deterministic, execution times can vary slightly among multiple executions of the codes, so results in Figure 34 are statistical between 10 different executions of the code under the same conditions. Nevertheless, the general behavior is always the same.

With respect to the computing time, it is notable from Figure 34a that the execution time for the magnetometric inversion improves 12.34% for 24 threads (11.75×) with respect to 12 threads (10.30×), while the improvement for the same threads of gradiometric and gravimetric inversions is only of 2.94×. With respect to the speed-up results (see Figure 34b), it is noticeable that the gravimetric case follows an almost perfect speed-up in the 1 to 12 threads range, while the results of magnetometric and gradiometric inversions are still close to the perfect speed-up in the same range of threads. Once the logical cores are reached (after thread 13), performance drops and HT offers no more than a small increase in performance which is only important for 24 threads. Finally, in the workstation in which the tests were performed, it must be emphasized that in the most resources-demanding case (gradiometric inversion) for the real model involving a salt structure, the parallelization scheme reduced the execution time from 16 h:56 m:15 s (in serial mode) to 01 h:30 m:41 s (for 24 threads), which means that the design and implementation of the parallel algorithm based on OpenMP is satisfactory.

4.4. Concluding Remarks The source-growth algorithm-based inversion scheme herein presented has shown to be an efficient method to invert the geometry of geobodies in highly complex zones such as the deepwater Gulf of Mexico. Part of the success is the incorporation of constrictions to mass density or magnetic susceptibility (pre-fixed) contrasts, which optimizes the inversion of large-scale bodies. On the other hand, the incorporation of FTG data provided an excellent approximation to the geometry of salt structures and their surroundings within the study region. While the velocities models of the zone under study obtained by seismics presented severe limitations in their interpretation due to the existence of a highly complex geology with intense allochthonous salt and clay tectonics—in particular because the lateral and vertical anisotropic changes are very abrupt—the geometric elements provided by potential data inversion allowed to reconstruct the interest horizons with a good approximation.

The selection of an adequateλparameter and contrasts for the inversion of data coming from the geological acquaintance of the region under study are fundamental to estimate accurate results. In these cases, this information aids in establishing adequate and more concise search spaces which highly reduces the non-unicity problems or the appearance of fictitious sources. The significant reduction of the space search allows the generation of non-prohibitive solutions from the computational burden point of view. This, because the initial model was built with appropriate distributions of mass densities/magnetic susceptibilities according to the current geology of the zone under study, which provided an initial model which is close to the global optimum of the objective function. In this sense, the pre-fixed mass density and magnetic susceptibility for the inversion algorithm were of 2200 kg/m³ and −0.00036 [1], respectively. In all cases, the total number of parameters for the inversion algorithm is 72,269 prisms; the observations mesh was located 100 m above the origin, and the misfit between computed and observed values was evaluated through the_L2norm.

In the case of FTG data, the joined inversion of the_Guν,_Gxyand_Gzzcomponents of the gradiometric tensor allowed to diminish the memory allocation for the Jacobian used in the simultaneous inversion, which in turn also speeds up the convergence of the process. As gradiometric data represent a better sensitivity to characterize the limits of bodies than gravimetric or magnetic data, their inversion provides results of great aid for the interpreter of reservoir engineering, due to the importance of possessing clear information to define geometries of bodies which volumetries are later to be estimated if they are of economic interest.

Although simultaneous inversion of gravity, magnetic and FTG were not performed due to the capacity of the available workstation, it is important to point out that the proposed method for the simultaneous inversion of FTG components can be easily extended to provide a simultaneous inversion of all potential field data in order to yield to a univocal and better determined imaging of 3D geological structures. With respect to the parallelization strategy, as OpenMP could be considered a standard for MT programming on SMP systems, sustained by a combination of functions and compiler directives, it reduces the lines of code and complexity with respect to low-level coding (like with the usage of pthreads), easing the migration of serial codes to a parallel shared memory environment. Additionally, the performance is really acceptable; in this regard, for the performance results in this work, we can observe an almost perfect speed-up until 12 threads; beyond this limit, the performance no longer scales on logical cores, this means that the parallel OpenMP-based implementation was well developed. In this sense, the source-growth inversion approach parallelized in this work is intended to be executed in multi-socket-multi-core shared memory computers, with or without MT technology activated. These results are encouraging with respect to the possibility of characterizing massive geological structures through the bodies growth algorithm in more affordable execution times. It is important to note that with the methodology herein developed, the work with sensitivity kernels is almost null, because the algorithm does not require any conditioning or solution to the system of equations. We can finally remark that this inversion approach can be considered to be a support in the definition of geobodies of interest in areas with high allochthonous and clay tectonics in deepwaters, where the definition of the geological structures is extremely complex. In this sense, this algorithm poses as an excellent alternative to invert potential field data in less complex geological scenarios than the defiant deepwater Gulf of Mexico. 5. Recommendations and Future Work Future work, as well as some important recommendations are:

• The performance of the proposed parallelization strategy in other CPU-based architectures than the one used in this work, remains to be tested.
• For extremely large structures, compression techniques for sensitivity matrices can be considered to alleviate the memory allocation.

• It is also possible to follow [60] to obtain a more accurate initial density model, considering also information from a velocities model of the area under study.

• In this work we did not consider other parallel computing architectures like GPUs, mainly for their limited memory size. Nevertheless, in a future when GPUs have more memory, this algorithm can be migrated to such architecture to compare performance.
• As previously mentioned, we did not invert simultaneously magnetic, gravimetric and gradiometric data due to memory limitations in our workstation. To further constrain the possible solutions of the model to be estimated, it remains to assess the behavior of the growth inversion algorithm when inverting such datasets jointly.

(a) For gravity-based anomalies.
Model	Δρ=[Kgm3]			Depth [km]
T	300			1
S	300			1
(b) For magnetic anomalies.
Model	Δm=Am	FI, FD	MI, MD	Depth [km]
T	0.01	90°, 0°	90°, 0°	1
S	0.01	90°, 0°	90°, 0°	1

λ parameter [1]	2.04
Initial maximum densityρ[kg/m3]	300
Smoothing coefficientτ [1]	8
Initial model misfit [mGal]	116.251
Minimum reached misfit [mGal]	0.373
Total iterations	1314

Initial maximum magnetization [A/m]	0.01
λ parameter [1]	4.04
Smoothing coefficientτ [1]	8
Initial model misfit [nT]	29,304
Minimum reached misfit [nT]	0.366
Total iterations	1323

λparameter	7.04		Initial maximum densityρ[kg/m3]	300
Smoothing coefficientτ [1]	8		Total iterations	1354
Guνinitial model misfit [Eo]	271.326		Guνminimum reached misfit [Eo]	2.102
Gxyinitial model misfit [Eo]	183.120		Gxyminimum reached misfit [Eo]	1.032
Gzzinitial model misfit [Eo]	581.733		Gzzminimum reached misfit [Eo]	4.223

λ parameter [1]	2.04
Initial maximum densityρ[kg/m3]	300
Smoothing coefficientτ [1]	8
Initial model misfit [mGal]	84.281
Minimum reached misfit [mGal]	0.273
Total iterations	1299

Initial maximum magnetization [A/m]	0.01
λ parameter [1]	5.04
Smoothing coefficientτ [1]	8
Initial model misfit [nT]	18.645
Minimum reached misfit [nT]	0.276
Total iterations	1212

λparameter	7.04		Initial maximum densityρ[kg/m3]	300
Smoothing coefficientτ [1]	8		Total iterations	1146
Guνinitial model misfit [Eo]	200.341		Guνminimum reached misfit [Eo]	2.583
Gxyinitial model misfit [Eo]	113.830		Gxyminimum reached misfit [Eo]	1.771
Gzzinitial model misfit [Eo]	370.240		Gzzminimum reached misfit [Eo]	5.043

λ parameter [1]	2.04
Initial maximum densityρ[kg/m3]	2200
Smoothing coefficientτ [1]	8
Initial model misfit [mGal]	28.411
Minimum reached misfit [mGal]	0.061
Total iterations	31,677

Initial maximum magnetization [A/m]	−0.00036
λ parameter [1]	8.04
Smoothing coefficientτ [1]	8
Initial model misfit [nT]	10.547
Minimum reached misfit [nT]	0.111
Total iterations	33,878

λparameter	8.04		Initial maximum densityρ[kg/m3]	2200
Smoothing coefficientτ [1]	8		Total iterations	30,611
Guνinitial model misfit [Eo]	20.607		Guνminimum reached misfit [Eo]	0.451
Gxyinitial model misfit [Eo]	16.446		Gxyminimum reached misfit [Eo]	0.180
Gzzinitial model misfit [Eo]	50.209		Gzzminimum reached misfit [Eo]	0.624

Author Contributions

Conceptualization, C.C. and C.O.; data curation, A.T. and C.O.; formal analysis, J.H. and N.V.; funding acquisition, C.C. and A.T. investigation, J.H., A.T. and M.O.; methodology, C.C. and N.V.; project administration, C.O.; resources, M.O. and C.C.; software, N.V., C.C. and J.H.; supervision, C.O.; validation, C.C.; visualization, N.V. and M.O.; writing-original draft preparation, J.H., A.T. and M.O.; writing-review and editing, J.H., C.C., N.V., M.O., C.O. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Secretaría de Investigación y Posgrado, Instituto Politécnico Nacional grant numbers 20200378, 20201040 and 20200259, by Centro de Desarrollo Aeroespacial, Instituto Politécnico Nacional computational time grant number FS0008, by project 292637 "Litoteca Nacional Sede Yucatán", from fund Secretaría de Energía-Consejo Nacional de Ciencia y Tecnología-Hidrocarburos, as well as by Tecnológico Nacional de México/IT de Mérida project # 8285.20-P. The APC was funded by Comisión de Operación y Fomento de Actividades Académicas del Instituto Politécnico Nacional and Secretaría de Investigación y Posgrado, Instituto Politécnico Nacional.

Acknowledgments

Authors acknowledge Secretaría de Energía-Consejo Nacional de Ciencia y Tecnología fund for the data as well as the scholarship provided to perform this study, as well as EDI grant by Secretaría de Investigación y Posgrado, Instituto Politécnico Nacional.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.