International Journal of Geophysics
Volume 2016 (2016), Article ID 1702164, 12 pages
http://dx.doi.org/10.1155/2016/1702164

Research Article

3D Gravity Modeling of Complex Salt Features in the Southern Gulf of Mexico

Mauricio Nava-Flores,¹ Carlos Ortiz-Aleman,² Mauricio G. Orozco-del-Castillo,² Jaime Urrutia-Fucugauchi,³ Alejandro Rodriguez-Castellanos,² Carlos Couder-Castañeda,⁴ and Alfredo Trujillo-Alcantara²

¹Facultad de Ingeniería, Universidad Nacional Autónoma de México (UNAM), Avenida Universidad No. 3000, CU, Coyoacán, 04510 Ciudad de México, DF, Mexico
²Instituto Mexicano del Petróleo, Eje Central Lázaro Cárdenas No. 152, San Bartolo Atepehuacan, Gustavo A. Madero, 07730 Ciudad de México, DF, Mexico
³Programa de Perforaciones en Océanos y Continentes, Instituto de Geofísica, Universidad Nacional Autónoma de México, 04510 Ciudad de México, DF, Mexico
⁴Centro de Desarrollo Aeroespacial, Instituto Politécnico Nacional, Belisario Domínguez No. 22, 06010 Ciudad de México, DF, Mexico

Received 30 November 2015; Accepted 3 February 2016

Academic Editor: Robert Tenzer

Copyright © 2016 Mauricio Nava-Flores et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

Table 1: P-wave velocity, density, and permeability of the rock salt bodies placed in oil and gas prospecting zones (SI units).

Property
Range of values
Reference

Seismic velocity (P-wave)
4270 to 5190 [m/s]
Grant and West, 1967 [8]

Density
2,100 to 2,200 [kg/m³]
Gardner et al., 1974 [9]

Permeability
<10⁻²⁰ [m²]
Carter et al., 1993 [10]

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)

Figure 1: (a) Discretized media formed by rectangular prisms and the vertical component of the gravity response calculated over an observation point at the surface observation plane and (b) prismatic body located in a right rectangular coordinate system. is the observation point, where the gravity response is computed.

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], is

where 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.

Figure 2: Salt bodies interpreted from a 3D PSDM velocity model. Each body is composed of many discrete points with a regular volume grid interval.

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/m³] and is depth [m].

Figure 3: Density versus depth curve, representative of sediments of Gulf of Mexico and rock salt density. Based on Nelson and Fairchild [14].

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/m³].

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.

Figure 4: Gravity anomaly caused by the salt bodies depicted in Figure 2.

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)

Figure 5: Lateral extent of source bodies from a set of depth-to-source estimation methods. (a) Horizontal Gradient, (b) Analytic Signal Amplitude, (c) 1st-order EAS, and (d) 2nd-order Enhanced Analytic Signal.

(a)

(b)

Figure 6: (a) Maxima location estimated from the methods depicted in Figure 5 and (b) sources in a plan view as interpreted from those maxima.

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.

Figure 7: Distribution with depth of gravity anomaly sources estimated by the 3DED algorithm.

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)

Figure 8: 3D initial model built from 3D Euler deconvolution solutions and the maxima of the lateral extent estimation methods previously applied to the gravity anomaly.

Figure 9: Gravity anomaly caused by the 3D initial model shown in Figure 8.

Table 2 summarizes the main characteristics of this 3D initial model.

Table 2: 3D initial model main characteristics.

Model size
Ensemble discretization
Density range

direction: 20,000 [m] 
direction: 20,000 [m]  
direction: 10,000 [m]
direction: 30 prisms 
direction: 30 prisms 
direction: 30 prisms 
Total number of prisms: 27,000
Individual prism size: 666.66 [m] × 666.66 [m] × 333.33 [m]
(, , and directions)
1,900 to 2,590 [kg/m³]

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 by

where 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].

Figure 10: Flux diagram representation of the simulated annealing algorithm applied in this work (modified from Ortiz-Alemán and Martin [15]).

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.

Table 3: Parameter values related with the 3D gravity data inversion procedure.

Parameters
Values

Number of inverted parameters
27,000

Number of observed data points
2,601

Initial temperature
1.0

Final temperature

Energy of the initial model
33.484543

Energy of the final model
0.0178259

Number of temperature reduction steps
1,000

Reduction temperature factor
0.98

Previous cycles to the parameter auto adjustment
10

Cycles of thermal equilibrium
5

Total number of tried models
149,400,000

Number of accepted models
115,281,390

Number of rejected models
34,118,610

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)

Figure 11: 3D Inverted Model generated from the 3D gravity anomaly inversion.

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).

Figure 12: Gravity anomaly caused by the 3DInvM shown in Figure 11.

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.

Figure 13: Spatial distribution of residuals along the grid and their amplitudes, computed as the absolute differences between gravity data before and after inversion.

Figure 14: Relative frequency distribution of residuals shown in Figure 13, including mean and standard deviation values.

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).

Figure 15: Temperature versus energy curve in the 3D gravity inversion process.

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.

References

C. Ortiz-Alemán and J. Urrutia-Fucugauchi, “Aeromagnetic anomaly modeling of central zone structure and magnetic sources in the Chicxulub crater,” Physics of the Earth and Planetary Interiors, vol. 179, no. 3-4, pp. 127–138, 2010.

S. Nagihara and S. A. Hall, “Three-dimensional gravity inversion using simulated annealing: constraints on the diapiric roots of allochthonous salt structures,” Geophysics, vol. 66, no. 5, pp. 1438–1449, 2001.

J. Zhang, C.-Y. Wang, Y. Shi et al., “Three-dimensional crustal structure in central Taiwan from gravity inversion with a parallel genetic algorithm,” Geophysics, vol. 69, no. 4, pp. 917–924, 2004.

L. Roy, M. K. Sen, D. D. Blankenship, P. L. Stoffa, and T. G. Richter, “Inversion and uncertainty estimation of gravity data using simulated annealing: an application over Lake Vostok, East Antarctica,” Geophysics, vol. 70, no. 1, pp. J1–J12, 2005.

R. A. Krahenbuhl, Binary Inversion of Gravity Data for Salt Imaging, Colorado School of Mines, 2005.

R. M. Rene, “Gravity inversion using open, reject, and ‘shape-of-anomaly’ fill criteria,” Geophysics, vol. 51, no. 4, pp. 988–994, 1986.

L. Uieda and V. C. F. Barbosa, “Use of the ‘shape-of-anomaly’ data misfit in 3D inversion by planting anomalous densities,” in SEG Technical Program Expanded Abstracts, vol. 4, pp. 1–6, 2012.

F. S. Grant and G. F. West, Interpretation Theory in Applied Geophysics, vol. 130, McGraw Hill, New York, NY, USA, 1967.

G. H. F. Gardner, L. W. Gardner, and A. R. Gregory, “Formation velocity and density-the diagnostic basics for stratigraphic traps,” Geophysics, vol. 39, no. 6, pp. 770–780, 1974.

N. L. Carter, S. T. Horseman, J. E. Russell, and J. Handin, “Rheology of rocksalt,” Journal of Structural Geology, vol. 15, no. 9-10, pp. 1257–1271, 1993.

D. Plouff, “Gravity and magnetic fields of polygonal prisms and application to magnetic terrain corrections,” Geophysics, vol. 41, no. 4, pp. 727–741, 1976.

R. J. Blakely, Potential Theory in Gravity and Magnetic Applications, Cambridge University Press, Cambridge, UK, 1995.

D. Nagy, G. Papp, and J. Benedek, “The gravitational potential and its derivatives for the prism,” Journal of Geodesy, vol. 74, no. 7-8, pp. 552–560, 2000.

T. H. Nelson and L. Fairchild, “Emplacement and evolution of salt sills in the northern Gulf of Mexico,” Houston Geological Society Bulletin, vol. 32, no. 1, pp. 6–7, 1989.

C. Ortiz-Alemán and R. Martin, “Inversion of electrical capacitance tomography data by simulated annealing: application to real two-phase gas-oil flow imaging,” Flow Measurement and Instrumentation, vol. 16, no. 2-3, pp. 157–162, 2005.

M. R. Hudec and M. P. A. Jackson, “Terra infirma: understanding salt tectonics,” Earth-Science Reviews, vol. 82, no. 1-2, pp. 1–28, 2007.

L. Cordell, “Gravimetric expression of graben faulting in Santa Fe country and the Espanola Basin,” in New Mexico Geological Guidebook, 30th Field Conference, no. 1977, pp. 59–64, 1979.

M. N. Nabighian, “Toward a three-dimensional automatic interpretation of potential field data via generalized Hilbert transforms: fundamental relations,” Geophysics, vol. 49, no. 6, pp. 780–786, 1984.

S. K. Hsu, J. C. Sibuet, and C. T. Shyu, “High-resolution detection of geologic boundaries from potential‐field anomalies: an enhanced analytic signal technique,” Geophysics, vol. 61, no. 2, pp. 373–386, 1996.

A. B. Reid, J. M. Allsop, H. Granser, A. J. Millett, and I. W. Somerton, “Magnetic interpretation in three dimensions using Euler deconvolution,” Geophysics, vol. 55, no. 1, pp. 80–91, 1990.

R. J. Blakely and R. W. Simpson, “Approximating edges of source bodies from magnetic or gravity anomalies,” Geophysics, vol. 51, no. 7, pp. 1494–1498, 1986.

P. Stavrev and A. Reid, “Degrees of homogeneity of potential fields and structural indices of Euler deconvolution,” Geophysics, vol. 72, no. 1, pp. L1–L12, 2007.

J. L. Rodríguez-Zúñiga, C. Ortiz-Alemán, G. Padilla, and R. Gaulon, “Application of genetic algorithms to constrain shallow elastic parameters using in situ ground inclination measurements,” Soil Dynamics and Earthquake Engineering, vol. 16, no. 3, pp. 223–234, 1997.

M. G. Orozco-del-Castillo, C. Ortiz-Alemán, J. Urrutia-Fucugauchi, R. Martin, A. Rodriguez-Castellanos, and P. E. Villaseñor-Rojas, “A genetic algorithm for filter design to enhance features in seismic images,” Geophysical Prospecting, vol. 62, no. 2, pp. 210–222, 2014.

S. Kirkpatrick, C. D. Gelatt Jr., and M. P. Vecchi, “Optimization by simulated annealing,” Science, vol. 220, no. 4598, pp. 671–680, 1983.

M. K. Sen and P. L. Stoffa, Global Optimization Methods in Geophysical Inversion, Cambridge University Press, Cambridge, UK, 2013.

C. Ortiz-Aleman, R. Martin, and J. C. Gamio, “Reconstruction of permittivity images from capacitance tomography data by using very fast simulated annealing,” Measurement Science and Technology, vol. 15, no. 7, pp. 1382–1390, 2004.

J. C. Ortiz-Alemán and R. Martin, “Two-phase oil-gas pipe flow imaging by simulated annealing,” Journal of Geophysics and Engineering, vol. 2, no. 1, pp. 32–37, 2005.

A. Corana, M. Marchesi, C. Martini, and S. Ridella, “Minimizing multimodal functions of continuous variables with the ‘simulated annealing’ algorithm,” ACM Transactions on Mathematical Software, vol. 13, no. 3, pp. 262–280, 1987.

W. L. Goffe, G. D. Ferrier, and J. Rogers, “Global optimization of statistical functions with simulated annealing,” Journal of Econometrics, vol. 60, no. 1-2, pp. 65–99, 1994.

C. Couder-Castañeda, C. Ortiz-Alemán, M. G. Orozco-del-Castillo, and M. Nava-Flores, “TESLA GPUs versus MPI with OpenMP for the forward modeling of gravity and gravity gradient of large prisms ensemble,” Journal of Applied Mathematics, vol. 2013, Article ID 437357, 15 pages, 2013.

C. Couder-Castañeda, J. C. Ortiz-Alemán, M. G. Orozco-del-Castillo, and M. Nava-Flores, “Forward modeling of gravitational fields on hybrid multi-threaded cluster,” Geofísica Internacional, vol. 54, no. 1, pp. 31–48, 2015.

R. Martin, V. Monteiller, D. Komatitsch et al., “Gravity inversion using wavelet-based compression on parallel hybrid CPU/GPU systems: application to southwest Ghana,” Geophysical Journal International, vol. 195, no. 3, pp. 1594–1619, 2013.