Abstract
We present a three-dimensional (3D) gravity modeling and inversion approach and its application to complex geological settings characterized by several allochthonous salt bodies embedded in terrigenous sediments. Synthetic gravity data were computed for 3D forward modeling of salt bodies interpreted from Prestack Depth Migration (PSDM) seismic images. Density contrasts for the salt bodies surrounded by sedimentary units are derived from density-compaction curves for the northern Gulf of Mexico’s oil exploration surveys. By integrating results from different shape- and depth-source estimation algorithms, we built an initial model for the gravity anomaly inversion. We then applied a numerically optimized 3D simulated annealing gravity inversion method. The inverted 3D density model successfully retrieves the synthetic salt body ensemble. Results highlight the significance of integrating high-resolution potential field data for salt and subsalt imaging in oil exploration.
1. Introduction
Hydrocarbon exploration is largely based on geophysical methods among which seismic reflection is the most intensely employed. Increased interest in subsalt related plays in the Gulf of Mexico and in other sedimentary basins around the world has turned oil and gas prospecting within these regions into a major challenge. Physical property contrasts of salt features such as highly contrasting seismic velocities relative to the surrounding media lead to complex wave diffraction patterns and lack of illumination near and below them.
In this context, gravity methods are well suited to support seismic prospection and improve subsalt imaging by taking full advantage of the density contrasts between salt bodies and surrounding sedimentary targets. Salt bodies retain low densities, whereas upon burial sediments compact increasing the density contrast. Table 1 shows typical ranges for seismic wave velocity, density, and permeability of salt. The seismic wave high velocity contrasts at salt-sediment interface result in strong refractions and reflections, making it difficult to image the bottom and structures beneath salt bodies. Major advances with improved images of subsalt plays have resulted from prestack imaging, with velocity-depth modeling and Prestack Depth Migration (PSDM). Nevertheless, subsalt structural complexities present major barriers, requiring new approaches and integrative analyses of seismic and potential field data.
Complex geological imaging using modeling and inversion of potential field anomalies has been examined in recent studies. Ortiz-Alemán and Urrutia-Fucugauchi [1] applied 3D magnetic field modeling to the study of the central zone of the Chicxulub impact structure, which had proved difficult to image from seismic reflection data. Nagihara and Hall [2] applied the simulated annealing (SA) global optimization method to invert synthetic gravity data due to a simplified salt diapir. They constrained the shape of the diapir at depth through inverse modeling in 3D. Zhang et al. [3] determined the crustal structure of central Taiwan through gravity inversion with a parallel genetic algorithm, using an initial model derived from 3D P-wave velocity tomography. Roy et al. [4] inverted gravity data using SA over Lake Vostok in East Antarctica in order to estimate the water-sediment and sediment-basement interfaces. Krahenbuhl [5] introduced an approach called binary inversion, which uses an assemblage of equal-volume prisms as model space, and density contrasts as model parameters, which only could take one of two possible values: zero for prisms located in salt-free zones and one for prisms placed in salt areas. Rene [6] developed a method of gravity inversion by iteratively applying open, reject, and fill criteria within a modeling procedure based on the use of prisms ensembles with density contrasts previously assigned and the “shape-of-anomaly” fill criterion. Uieda and Barbosa [7] performed a 3D gravity inversion method by planting anomalous densities. They used the “shape-of-anomaly” data misfit in conjunction with the -norm data-misfit function achieving better delineation of elongated sources and the recovery of compact geologic bodies.
In this work we built a 3D gravity model including several allochthonous salt bodies as interpreted from a Prestack Depth Migration seismic volume, integrating results from different potential field techniques such as edge-source detection, depth-to-source estimation, and 3D gravity inversion.
2. Materials and Methods
2.1. Forward Gravity Modeling of Salt Structures
The whole computational domain, including several salt features, was discretized into an ensemble of rectangular prismatic elements. Its gravity response, that is, , was calculated by summing the individual responses of every single prismatic element, on all points belonging to the observation grid. Figure 1 shows a sketch of the gravity response generated by an ensemble formed by prisms, computed over an observation plane.
(a)
(b)
The total gravity response calculated at some observation point was the sum of the gravity contributions generated by the prisms of the ensemble: where represents the gravity response on the observation point and is the gravity response due to the particular prism , observed on the position .
Now, the vertical component of the gravity vector, due to each individual rectangular prism, with constant density , according to Plouff [11], Blakely [12], and Nagy et al. [13], iswhere is the universal gravitational constant, , , and . Figure 1(b) shows the elements involved in (2).
As illustrated in Figure 1(b), and according to (2), to calculate the gravity response of a rectangular prismatic body it is necessary to determine the density and coordinates of its extreme vertices. In this approach it is just necessary to know the density, the coordinates of only one point, and the volume grid interval of the ensemble to which it belongs.
Taking into account the above and the fact that salt bodies embedded in terrigenous sediments are geometrically irregular, they can be modeled as ensembles of regular rectangular prisms, formed by discrete points. Analyzed salt bodies (Figure 2) were imaged from digital 3D Prestack Depth Migration seismic velocity modeling, using real data from areas located in the southeastern sector of the Gulf of Mexico.
On computing the gravity response of the salt bodies illustrated in Figure 2, a background density, representative of the surrounding sedimentary rocks, is needed. Curves of sedimentary rocks densities, representative of the Gulf of Mexico, were published by the Applied Geodynamics Laboratory (AGL) by Hudec and Jackson [16], based on the work of Nelson and Fairchild [14], in which they propose that the density of those sedimentary units could be modeled as a function of depth (Figure 3) by the exponential curve:where is the sediments density [kg/m3] and is depth [m].
The gravity response is therefore calculated by first obtaining a density contrast between the salt bodies and the surrounding sediments, in the position of each point source (prism), subtracting (3) from each single salt prism, assuming that its density is constant, with the value 2,180 [kg/m3].
Figure 4 shows the gravity anomaly calculated from the ensemble formed by all the prismatic sources simulating the salt bodies depicted in Figure 3. The number of rectangular prisms considered in the calculation of the anomaly was 201,540; each prism size is 50 [m] × 50 [m] × 25 [m] (, and directions). The grid interval for the observation grid in both directions, and , was 0.4 km, and the number of observation points in both directions was 51, that is, 2,601 grid observation points.
2.2. Shape and Depth Estimation of Salt Structures
Several methods especially suited to enhance anomalies and estimate depth to source are commonly applied to potential field data. While there are methods that use systematic search algorithms to find a solution of the distribution of the densities of the model [7], and others which use a great amount of rectangular prisms to obtain a good approximation of the gravimetric anomalies [6], in this work we applied a series of those approaches to gravity gridded data (Figure 4), in order to infer an initial 3D density distribution for inverse modeling. These methods include the Horizontal Gradient (HG), as proposed by Cordell [17], the 3D Analytic Signal Amplitude (AS), developed by Nabighian [18], the Enhanced Analytic Signal (EAS) introduced by Hsu et al. [19], and the 3D Euler deconvolution (3DED) developed by Reid et al. [20].
Figure 5 shows the results obtained after applying the HG, AS, and EAS methods to the gravity anomaly of Figure 4. Location of the maxima in those grids was estimated by the method of Blakely and Simpson [21] and roughly corresponds to the lateral extent of the gravity field sources, that is, the salt bodies in a plan view. The edges of salt bodies as interpreted from these maxima are shown in Figure 6.
(a)
(b)
(c)
(d)
(a)
(b)
After estimation of projected source location on a horizontal plane (plant-view of sources), the corresponding distribution of the sources with depth was inferred in order to build an initial 3D structural model. For this purpose, we applied the 3DED algorithm to the gravity anomaly grid, considering a structural index , assuming a geologic contact-type of source [22], and a 4 km size square window. The computed solutions are shown in Figure 7.
We built a 3D structural model including two huge salt diapirs surrounded by sedimentary rocks with relative density contrast assigned as a function of depth (3), by considering source location in plant (Figure 6) and the 3D Euler deconvolution solutions (Figure 7). We then computed the gravity anomaly for this 3D model, in order to evaluate how well it correlates with the gravity anomaly response from the originally postulated salt bodies, shown in Figure 4.
Figure 8 shows this 3D model, labeled as 3D initial model (3DIniM), in three orthogonal projections, and a 3D perspective oblique view. The salt masses were presented gray colored for display purposes, and their gravity anomaly is shown in Figure 9.
(a)
(b)
(c)
(d)
Table 2 summarizes the main characteristics of this 3D initial model.
The computed gravity anomaly for the 3D initial model qualitatively resembles the shape of the anomaly produced by the salt bodies interpreted from a PSDM volume (Figure 4) but in terms of the amplitude of relative error remains still quite large. To minimize such error, we inverted the gravity anomaly by using a numerically optimized simulated annealing algorithm, as discussed in the next section.
2.3. 3D Inversion of Gravity Data by Simulated Annealing
According to (2), the gravity anomaly of the entire prism ensemble is the sum of the gravity anomalies generated by each of the individual prisms, so we can rewrite (2) as
This is a linear system of equations, where denotes the amount of observation points and is the number of prisms in the ensemble. The linear system can be also represented as
Here, are the elements in the sensitivity kernel or sensitivity matrix. Each one of its elements accounts for the contribution to the complete gravity anomaly in the observation point, due to the prism located on the position inside the ensemble.
To solve the inverse problem, we chose the simulated annealing global optimization method. A main drawback of global optimization lies in the excessive amount of forward problem computations required to solve the inverse problem. In the past decades, global optimization has been successfully applied to several geophysical exploration issues, where dimensionality of the inverse problem did not represent a bottleneck [23, 24].
The simulated annealing method was conceived as a mathematical analogy with the natural optimization process of crystal formation from a mineral fluid at high temperature. Its basic concepts were taken from the statistical mechanics.
The simulated annealing optimization process emulates the evolution of a physical system as it slowly cools down and crystallizes at a state of minimum energy. If temperature, , is gradually reduced after the thermal equilibrium has been reached, then in the limit, as , the minimum energy state becomes predominantly likely, as well as the crystal formation, and therefore the parameter configuration could be considered as optimum model.
Following Kirkpatrick et al. [25], we used the Metropolis algorithm as central part of our simulated annealing method, taking advantage of its ability to escape from local minima, which increases the chances to reach the global optimum, and at the same time approaches to the Boltzmann probability density function asymptotically [26]. It consists basically in disturbing some initial model, , which already has the energy content , and getting a new model, , with energy , and then calculating the energy level change due to the disturbance applied, , and accepting or rejecting on the basis of the value of calculated: if , then will be unconditionally accepted, but in the case that , will be accepted with the probability .
This acceptance-rejection procedure is repeated several times for a fixed temperature, , until the thermal equilibrium is reached, which is characterized by not exhibiting substantial changes in the energy level before the temperature reduction.
To compute the energy level in each stage, we used a normalized norm [27, 28], given bywhere is the observed gravity anomaly and is the gravity anomaly calculated for the model.
The cooling schedule we choose reduces the temperature in an exponential fashion by multiplying the actual temperature by some parameter in each temperature cycle pass and, according to Nagihara and Hall [2], is characterized by ensuring convergence to the global minimum:where is the initial temperature of the system, is the temperature at th stage, and is the reduction temperature parameter .
Finally, this process is repeated until reaching the limit or controlled by some stop criterion given as a tolerance error with respect to and, at the same time, by a maximum number of temperature reductions.
One first improvement made to the basic simulated annealing method in this work was to accelerate the product , as proposed by Ortiz-Alemán and Martin [15] by using a previously computed forward problem, and using it to update the actual one, summing it to the product ( is the disturbance applied to the model parameter ). This improvement justifies the employment of a global optimization method in the inversion of a quite large linear system, as after the first iteration it is no longer required to compute the forward problem for a complete ensemble, and hence the numerical burden will be dramatically reduced.
The final improvement made to the SA method consisted in applying an auto adjustment to the amplitude control parameter, from iteration to , proposed by Corana et al. [29]. Let be the relation between the numbers of accepted () and rejected models () by the Metropolis criterion. If , then , and if , then , where is a constant value fixed at 2.0.
The final SA inversion algorithm is summarized in the diagram shown in Figure 10. Its three-nested-loop structure is based on the algorithm presented by Goffe et al. [30].
This modified SA algorithm was applied to the gravity anomaly data generated by the postulated set of salt bodies, with the following restrictions:(1)The lateral extent of all models generated by the inversion procedure was restricted to the interpreted source borders (Figure 7).(2)The model space was bounded according to the salt and sediments density contrast (Figure 4).Table 3 shows the 3D gravity data inversion parameters of the SA algorithm method as applied in this work.
3. Results and Discussion
The 3D density model resulting from the inversion procedure (Figure 11) is shown in the same projections and perspective angles as displayed in the 3DIniM. In order to differentiate this inverted model from the initial one, it is labeled as 3D Inverted Model (3DInvM).
(a)
(b)
(c)
(d)
The gravity anomaly grid generated by the 3DInvM shows that, despite the apparent differences in the central part of the grid corresponding to the gravity minimum, the amplitudes are similar to the observed gravity (Figure 12).
In order to quantify the quality of the 3DInvM, we calculated the differences between the gridded gravity anomaly data of the salt bodies (Figure 4) and the 3DInvM (Figure 12). The absolute values of these differences are shown in Figure 13, illustrating the spatial distribution of residuals along the grid and their amplitudes, whose maximum difference (0.005624884 [mGal]) is 4.26% of Figure 12 total range (0.132137625 [mGal]). The histogram and the mean and standard deviation values of the residuals are shown in Figure 14, where the mean value (−0.000186 [mGal]) and low standard deviation (0.001043 [mGal]) indicate that the inverted gravity anomaly successfully resembles the observed gravity anomaly. Based on this last fact we find our inversion results as very encouraging.
The misfit curve, representing the relationship between temperature and energy parameters along the inversion process, exhibits three different kinds of convergence rates: a gradual decay in the beginning of the inversion, an intermediate region of sharp decreasing misfit, and a zone of progressively slower convergence rates until final entrapment (Figure 15).
4. Conclusions
In this study we applied 3D gravity modeling and inversion in a complex geological setting involving several allochthonous salt features embedded in terrigenous sediments, representing a challenging and quite realistic scenario commonly found in the southern Gulf of Mexico.
Several methods especially suited to enhance anomalies and estimate depth to source are used in this work to determine an initial 3D density distribution for inverse modeling. These methods include the Horizontal Gradient (HG), the 3D Analytic Signal Amplitude (AS), the Enhanced Analytic Signal (EAS), and the 3D Euler deconvolution (3DED). We built an initial density model by integrating results from this set of shape- and depth-source estimation algorithms. We finally applied a numerically optimized three-dimensional simulated annealing gravity inversion approach. As the total amount of evaluated forward models in this study case was quite large (~150 million), application to other realistic gravity modeling efforts should consider the use of high performance computing for the forward and inverse problems, as recently introduced by Couder-Castañeda et al. [31, 32] and Martin et al. [33]. Results highlight the significance of integrating high-resolution potential field data to imaging complex salt tectonics media for oil exploration.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors of this work acknowledge the financial support provided by SENER-CONACyT Project 128376 and Mexican Institute of Petroleum Projects D.00475 and H.61006.