Modeling in Geophysics is a permanent exercise, trying to reproduce, as close as possible, the Earth systems and processes that scientists try to understand. Here we are concerned with one of the Potential Methods used for this purpose: the gravity method. As technology has progressed, determination of the gravity field of the Earth at given locations has also changed. For most of last century gravimetric determinations were mainly performed on the surface of the Earth, with oceanic measurements occurring along ships’ navigation trajectories, including a few aerial acquisitions. The process was usually slow, expensive, and limited in extent. By the end of the century a radical change occurred, and satellites were incorporated to gravimetric data acquisition, introducing global data coverage     . Presently there is a wide variety of satellites from which gravity data can be obtained, however, there is not much information on how much scope they can have as regional study tools, and it is not clear among them, which may be the most appropriate for certain type of study. We will focus on two of the most recent data sets, having greater spatial definition; they are the EIGEN-6C4 model (ICGEM,  ) the EGM2008 model  and the GGMplus model .
The EIGEN-6C4 gravimetric satellite data model is the highest resolution model with the highest coverage worldwide. Many studies have proven the value of their data for numerous regional tectonic studies, e.g., . Studies that otherwise would have been impossible or immensely costly for the area of acquisition. Recently the GGMplus model has been available, with an unprecedented resolution (200 m). However, there is no systematic comparison of its performance versus the EIGEN-6C4 model, or whether it manages to determine the regional structures as well as the EIGEN-6C4 model. Its high-resolution is generated mainly using the interpolations of a finer topography, combining it with a database of a previous gravimetric model  achieving an effective use in the determination of geological structures of shorter wavelength. To resolve this doubt, we will compare both models in a real situation and observe their behavior.
2. Materials and Methods
The gravimetric satellite models used for evaluation in this study are: the EIGEN-6C4 model, made available through International Centre for Global Earth Models (ICGEM,  ) it was sampled to 0.009˚ (approximately 1 km; the maximum resolution for the EIGEN-6C4;); the EGM2008 model , made available from https://topex.ucsd.edu/cgi-bin/get_data.cgi (Topex data; Scripps Institution of Oceanography, University of California San Diego), with 1 arc-minute resolution (1' = 1.85 km); and finally, the GGMplus, available fromhttp://ddfe.curtin.edu.au/gravitymodels/GGMplus/ with 7.5 arc-sec (220 m) resolution .
For the calculation of the Bouguer Anomaly (AB) at a given location, it is necessary to use the elevation at the given point. The elevation is provided by a digital elevation model (DEM). In this work we used 3 DEMs with different resolutions. The ETOPO1 topography model (National Center for Environmental Information; NOAA; https://ngdc.noaa.gov/mgg/global/global.html), with 1 arcminute resolution, sampled to 0.009˚ resolution (approximately 1 km), the maximum resolution for the EIGEN-6C4. The SRTM15 topography model https://topex.ucsd.edu/cgi-bin/get_srtm15.cgi; , has a resolution 15 arc-second (approximately 450 m); finally, the SRTM90 (https://srtm.csi.cgiar.org/srtmdata/;  ) with a 7.5 arc-sec resolution. This was done to make the comparison more realistic, and to diminish the DEM’s influence to a minimum. The DEM used was SRTM15, which has an intermediate resolution between the EIGEN-6C4 and the GGMplus.
2.1. Free Air Anomalies
For the direct comparison of the GGMplus model it is necessary to obtain the Free Air (FA) anomaly of the EIGEN-6C4; to this end, we will use the ETOPO DEM, and the EIGEN-6C4 data. We calculate the FA anomaly of Model EIGEN-6C4 with:
2.2. Obtaining Gobs
For a direct comparison of the Gobs of the EIGEN-6C4 model with the GGMplus model, it is necessary to convert the GGMplus model (FA anomaly) to an observed gravity value at ground surface (Gobs); this is done using the elevation of the station or, in general, the elevation provided by a DEM. The precision of the Gobs will depend, thus, on the precision of the DEM used to obtain it. In the case of the GGMplus model, the DEM has a resolution of approximately 7.5 arc-sec (or 200 m approximately  ), for this reason we will use the DEM of the SRTM90. In this way, we get Gobs
Gravity on Earth is distributed in a similar way as it would be in a sphere, but with a slight increase towards the geographical poles, due to the flattening of the terrestrial spheroid  . In this way, the gravity of the Earth varies according to the latitude, it is what is known as theoretical gravity (Gtheo). To calculate the theoretical gravity at any point on the planet we need Equation (2) (Geodetic Reference System of 1980, GRS80).
The altitude correction (Calt) is applied to correct the distance of the measurement to the reference level; therefore, it has a + or - sign depending on whether the altitude is less or greater than the reference level. The equation for calculating the simplified height correction is :
where h is the station elevation of the DEM. From (3) and (4) in (2) we obtain the value of Gobs:
2.3. Calculation of Bouguer Anomaly
Since satellites acquire data at several hundred kilometers above the ground surface, when computing the gravity field, a new type of correction must be included, corresponding to the weight of the column of air between ground and the satellite position, and it is also necessary to consider the sphericity of the Earth in the calculation. In  the procedure to obtain the Bouguer anomaly from EIGEN-6C4 data was explained, where such a correction was included according to the new gravimetric standard of the USGS; e.g., .
We perform 3D gravity inversions using the method described by , based in turn on the theoretical considerations of . The inversion results are densities in g/cm3. The code is implemented in the Oasis Montaj program of Seequent. To represent geologic volumes, the program uses a Cartesian Cut Cell algorithm (CCC); to match the observed result with the calculated one, within established error limits, the inversion program uses an Iterative Reweighting Inversion algorithm (IRI) . The program can also perform magnetic data inversions, although this is not relevant to the present problem. The model requires a DEM; we use the same geographic area and the same DEM (SRTM15) to perform the GGMplus and EIGEN-6C4 Bouguer anomaly inversions to better evaluate the inherent differences between the gravity data of each set. The term voxel, derived from pixel, but representing a volume, contains the results of the 3D inversion. We selected a volcanic area in central Mexico to compare 3D inversions with the two data sets. Additional examples of the inversion process can be found in  - .
3.1. Free Air Anomalies
Since GGMplus is a Free Air anomaly model (FA), we shall start by comparing it with other FA models, for a preliminary evaluation. Figure 1(a) shows the one corresponding to GGMplus, Figure 1(b) shows the EGM2008 FA anomaly, and Figure 1(c) the EIGEN-6C4 FA anomaly. The higher resolution of the former
Figure 1. Free Air Anomaly maps obtained over the same geographic region, illustrating (a) the GGMplus, (b) the EGM2008, and (c) the EIGEN-6C4, FA corresponding anomalies. The better definition of the FA anomaly of the GGMplus data is neatly observed.
can be readily appreciated comparing the low-gravity regions on the west side of the maps.
The general aspect of the anomaly (Figure 1(c)) is quite similar to that of model EGM2008 (Figure 1(b)), although the anomaly range (125.8 to −55.3 mGal) is somewhat smaller than that of EMG2008 model (148.1 to -55.3 mGal). In Fig 1 the resolution of the anomalies is clearly observed, reflecting smoother contours in the anomalies of the EIGEN-6C4 and EGM2008 models. Between the EIGEN-6C4 model and the EGM20087 we observed differences in the values of the anomalies, the negative values of FA being the most affected. In the EGM2008 model, more marked transitions between anomalies are observed, while in the EIGEN-6C4 model registers a less marked transition between the anomalies, better resembling the transition obtained in the GGMplus model. The range of GGMplus varies from 125.5 to −66.6 mGal, enhancing details in the negative portion (Figure 1(a)), and being more similar to the behavior of the data distribution of the model EIGEN-6C4 (Figure 2). The distribution histograms (Figure 2), show a greater correlation between the GGMplus models and the EIGEN-6C4 model, in both we observed a bimodal distribution with modes close to −50 and 50 mGal, and with a minimum between modes less than 50% of the accumulated data. In addition, the range of the histogram data varies from −90 mGal to 215 mGal, while the range for the EGM2008 model is −87 to 315 mGal, extending the upper range by about 70 mGal.
The Observed Gravity (Gobs) is the gravity that is measured in the surface of the Earth, which is the value of the GGMplus model we compared with the EIGEEN-6C4 model. We calculate the Gobs of the GGMplus model according to Equation (5) (Figure 3(a)) and compare it to the corresponding EIGEN-6C4 (Figure 3(b)). A general similitude of results is observed although differences appear at closer inspection. The statistical data is very similar and the behavior of the data distribution in the histograms is quite similar.
Figure 2. Frequency histograms for FA models (a) GGMplus, (b) EGM2008, and (c) EIGEN-6C4. There are subtle but relevant differences among them, revealing a greater similitude between EIGEN and GGMplus distributions and their cumulative frequencies (black line).
Figure 3. Comparison of Gobs between (a) GGMplus and, (b) EIGEN-6C4. Notice the great similarity between the gravity ranges in both maps.
3.3. Radial Power Spectrum
A comparison of the radial power spectrum of GGMplus and EIGEN-6C4 on a land region, shows additional differences between the depth sources of the anomalies as shown in Figure 4. And although it is not the purpose of this study to clarify the implications of these variations, this result shows that there is an important difference between them that must be considered when defining the objective of the study. This change occurs between the arrows of Figure 4, in the range of the spectrum 5 - 3.5 CYL/K_unit, which is the range in these areas for structures of less than 10 km.
Figure 5 Comparison between simple Bouguer anomalies of GGMplus and EIGEN-6C4 of the same land region, where no filters have been applied. Whilst the latter shows high-frequency variations, the former shows smooth regions. These differences affect 2D and 3D gravity models. This high-frequency variation has a low impact in the distribution of the histograms, with similar behavior, even though there is a greater content of high-frequencies in the GGMplus model.
Inversions corresponding to EIGEN-6C4 and GGMplus data sets are in the black
Figure 4. Power spectrum of GGMplus and EIGEN-6C4. The arrows indicate where the spectrum behavior changes for structures located shallower than 10 km.
Figure 5. Comparison between simple Bouguer anomalies of GGMplus and EIGEN-6C4. The full anomaly range is comparable, the difference in standard deviation is 3 mGal, and the difference between their means is 5 mGal.
polygon of Figure 6 and Figure 7 are displayed in Figure 8. All inversion parameters were the same, with exception of the BA data sets; thus, differences observed are to be attributed to differences between those data sets. The inversion resolution is 1000 m. We use the region containing Nevado de Toluca volcano to exemplify the differences between the two data sets; it was chosen since it contains a variety of geologic structures  that help enhance the differences between them. The recent activity of this Quaternary volcano indicates that additional activity cannot be discarded . Although the present study provides valuable insights about the structures in this volcanic area, it is not intended as a volcanic study.
Figure 6. Topography of the geographic area where 3D gravity inversions are performed. Nevado de Toluca volcano is located at 99˚45'W, 19˚07'N and Sierra Las Cruces runs NW-SE along the NE side of the map. The black polygon represents the region where inversions are performed. The magenta lines correspond to the cross-sections of the inversions, to be presented below. Figure made with GeoMapApp (www.geomapapp.org) /CC BY/CC BY (Ryan et al., 2009).
Figure 7. GGMplus Bouguer Anomaly map showing the region in which 3D inversions are performed (black polygon). Cross-sections are obtained from the respective inversions along lines N-S and E-W (magenta).
A preliminary assessment can be made considering the front sections along the X-axis (Figure 8). Both sections show two regions of higher densities (red); the higher definition of the section corresponding to GGMplus is quite clear. Whilst a diffuse red portion occupying about half of the EIGEN section is observed, the corresponding area in the GGMplus section shows the region divided in two, by a gap of material of lower density; that is, a diffuse region is substituted by a sharp boundary, made possible along the complete vertical range by the higher resolution of the GGMplus data set.
Another difference is appreciated in the surface distribution of density not shadowed by the displayed DEM. The EIGEN model shows rather uniform portions of green, whilst the GGMplus model shows a speckled surface with a distribution of cells with lower density values (blue).
Sections along the Y-axis of EIGEN-6C4 and GGMplus appear in Figure 9, where the greater resolution of the GGMplus inversion again shows important structural differences with respect to that of EIGEN-6C4. In the former, the low-density anomaly associated with Nevado de Toluca volcano, in the north portion of the section (red arrows), is continuous and of larger dimensions, whilst
Figure 8. 3D inversions of the EIGEN-6C4 and GGMplus BA maps. The DEM used in the calculation is the same in both cases . The voxel contains cells in the X and Y directions of 1 km, respectively; the Z-dimensions are variable, increasing in size downwards. The respective voxels have the same vertical extent: −1185 to 3340 m, as well as the same cell distribution.
Figure 9. Same inversions as those in Figure 8, showing a section through Nevado de Toluca volcano (yellow arrows). To the N of the line, EIGEN-6C4 model displays a faint low (blue) at the bottom of the section, separated from a wider, shallower one to the north (red arrows). GGMplus shows instead a wide, continuous region of low densities not reaching the surface, associated with the same location. A much larger anomalous region of low density is observed to the south (yellow arrows) corresponding to NT volcano. Regarding the low-density anomaly at the volcano surface EIGEN-6C4 shows an elliptic anomaly, whereas GGMplus shows a divided upper portion and more detail in the anomaly reaching the surface. The inversion is performed at 1000 m resolution.
in the latter it appears of smaller size and fragmented. In the south portion of the sections appears an unexpected, low-density anomaly of large dimensions (yellow arrows); we shall discuss the volcanic implications of this anomaly elsewhere. In this study, we will only compare its manifestations in the two cases under comparison. The EIGEN-6C4 section shows the anomaly as an ellipse, whilst the GGMplus section shows a wider anomaly, divided in the surface by a higher density region.
To inspect in more detail the results of the inversions, we compare cross-sections along lines N-S and E-W of Figure 7; the N-S line intersects the edifice of Nevado de Toluca volcano ; the EIGEN-6C4 cross-section (Figure 10) shows a single negative anomaly in the southern portion of the section; the center of this negative anomaly is located at −2000 m approximately, corresponding to NT volcano. On the north side, a deeper anomaly of higher density than that observed under the volcano appears with its center at approximately −3000 m; it may correspond to intruded material.
The cross-section corresponding to GGMplus exhibits important differences. The volcano anomaly shows its center closer to the surface by almost 2 km (i.e., at sea level) and a high-density gradient underneath the center; this may have important implications for the location of a magma chamber. The verticality of the anomaly is modified to a south-dipping anomaly. In addition, a bifurcation is observed, divided by a small, higher-density region in the surface, in agreement with the observation made in Figure 9. The elongated anomaly to the north has now been displaced southwardly and its center is shallower. The northernmost density-low diminished its value, although the region maintains its general shape.
Figure 10. Cross-sections extracted from the voxels corresponding to EIGEN-6C4 and GGMplus along line N-S in Figure 7. Contours enhance the differences between the two sections. The vertical axis ranges from 5000 to −5000 m. Notice that the averaged topographic profiles are the same in both cross-sections, and so are the voxel characteristics in the vertical direction. The inversion resolution is 1000 m. The coordinates are UTM northings.
Cross-sections corresponding to line E-W appear in Figure 11, they are perpendicular to line N-S and go across the structure of Nevado de Toluca volcano. The EIGEN-6C4 section presents a low-density anomaly at the center, corresponding to the position of the volcanic structure. The result is similar to that obtained for the N-S section above. Together they suggest that the volcanic anomaly has the shape of a conus, with a rounded apex on top, and its center at −2 km. The section corresponding to GGMplus shows instead a more compact structure, with a division at the top of the low-density anomaly, like the one observed in the N-S section. The differences observed on the flanking portions of the low-density anomaly are more drastic than those observed in the N-S section, tending to emphasize the presence of finer geologic structures.
5.1 Importance of DEM
In the use of the gravity satellite model, the correct choice of the DEM allows to
Figure 11. Cross-sections extracted from the voxels corresponding to EIGEN-6C4 and GGMplus along line E-W (Figure 7). The central, low-density anomaly corresponds to the location of Nevado de Toluca. Contours enhance the differences between the two sections. Notice that the averaged topographic profiles are the same in both cross-sections, and so are the voxel characteristics in the vertical direction.
get the most out of the data. For example, the ideal case is to use a DEM that has the same spatial resolution as the gravimetric model. Figure 12 shows a comparison between simple Bouguer anomalies of GGMplus with DEM 7.5 arc-sec (the same resolution of model GGMplus), and simple Bouguer anomalies of GGMplus with DEM 15 arc-sec (half the resolution of model GGMplus). The full anomaly range is comparable, and the statistical data is very close, nonetheless the anomaly is slightly distorted, a small distortion difficult to appreciate which blurs the focus of the anomalies without changing its behavior.
In the use of the gravity satellite model, the correct choice of the DEM allows to get the most out of the data. For example, the ideal case is to use a DEM that has the same spatial resolution as the gravimetric model. In Figure 13 show how comparison between simple Bouguer anomalies of EIGEN-6C4 model with DEM with the same resolution of model, and simple Bouguer anomalies of EIGEN-6C4 model with DEM with 15 arc-sec (the double resolution of model). The full anomaly range is comparable, and the statistical data is very close, nonetheless the anomaly presented a high frequency noise, because when using a higher resolution DEM, the C_alt and CB corrections are sampled at twice the resolution of the gravimetric model. Although these discrepancies can be corrected with filters that remove high frequencies.
Figure 12. Comparison between simple Bouguer anomalies of the model GGMplus processed with two DEM: in case (a) with 7.5 arc-sec (220 m) resolution; in case (b) with a 15 arc-sec (450 m) resolution.
Figure 13. Comparison between simple Bouguer anomalies of the model EIGEN-6C4 processed with two DEMs: (a) with a 60 arc-sec (1 km; ETOPO1) resolution; (b) with a 15 arc-sec (450 m; SRTM15) resolution.
A brief discussion of the BA in Figure 7 helps explain some of the observations made regarding the cross-sections. We first note that the shape of the BA in the vicinity of Nevado de Toluca volcano has a horseshoe shape; Nevado de Toluca and an adjacent anomaly to the N, are in its western branch. Although the main traits are maintained between the two cross-sections under comparison (Figure 10 and Figure 11), the anomalies are more sharply exposed by GGMplus. The BA map in Figure 7 shows a low-density anomaly on the eastern portion of the E-W line; the cross-section of EIGEN-6C4 in Figure 11 shows only a small decrement of the prevailing positive anomaly, the GGMplus more precisely reflects the low-density region, showing a clearly defined negative region. Additional differences are evident on the west end of the line; a strong difference is observed between the single anomaly on the EIGEN-6C4 with its center at −2000 m, and an anomaly of a quite different shape in the GGMplus, where the surface is divided into two low-density anomalies, with density increasing with depth.
We selected a region in which geologic differences are present, expecting that they would be reflected in the rock-density variations associated with the intervening geologic sources. We selected a voxel for the inversion with cells of 1 km on the side. The results obtained belong to this spatial resolution; they can be increased or diminished varying those dimensions. We would expect that the larger the cell dimensions the smaller the differences between the 3D inversions with the BA of both data sets. If cell dimensions are diminished, the better resolution of GGMplus is expected to yield more accurate results.
After evaluation of the frequency histograms, the Free Air anomalies, the radial power spectrum, and the simple Bouguer anomalies of the EIGEN-6C4 and GGMplus data sets, we concluded that the latter has a better spatial resolution. We infer that for wavelengths of 5 - 3.5 CYL/K_unit, the former can produce better results with respect to the location of the deep sources of the gravitational field, while the GGMplus model could represent better results for shallower sources. The effects on models built from 3D inversions were evaluated under conditions of complete equality, except for the BA for each data set. The GGMplus model indicated that its resolution advantages are maintained in the modeling process.
The GGMplus model indisputably demonstrated that it achieves a higher spectral resolution in shallow cortical elements, which is reflected in a better identification of local elements. In addition, it is observed that it maintains the regional tectonic trends presented by the EIGEN-6C4 model, which makes it ideal for regional and local gravimetric studies. The biggest problem with the GGMplus model is still its limited coverage, since it only presents data up to 10 km from the coastline, which is why it is suggested to be combined with the EIGEN-6C4 model for study areas that include marine regions.
The use of DEMs with a higher sampling resolution than the resolution of the gravimetric model, allows the overestimation of AB of the model to be reduced by a fraction. However, they generate a random and chaotic high-frequency noise pattern, for this reason their use is discouraged, although it can be corrected using band-pass filters. The use of DEMs of lower resolution than the gravimetric model is discouraged since it generates blurring of the anomalies.
A comparison of the radial power spectrum of GGMplus and EIGEN-6C4, on a land region, shows additional differences between the depth sources of the anomalies; the difference is particularly relevant for the depth sources between approximately 15 and 1 km. This change between spectra is of interest in studies that try to model the depths of sources according to the power spectrum radial average method ; for this reason, we suggest the GGMplus model for shallow structures and the EIGEN-6C4 for deep structures.
During development of this work, MC received support from Consejo Nacional de Ciencia y Tecnología (CONACYT, México. This study has been supported by IIMAS, UNAM; we acknowledge material support from both institutions. Comments of one anonymous reviewer helped improve the original manuscript. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.
 Becker, J.J., Sandwell, D.T., Smith, W.H.F. and Braud, J. (2009) Global Bathymetry and Elevation Data at 30 Arc Seconds Resolution?: Global Bathymetry and Elevation Data at 30 Arc Seconds Resolution: SRTM30_PLUS. Marine Geodesy, 32, 355-371.
 Sandwell, D.T. and Smith, W.H.F. (2009) Global Marine Gravity from Retracked Geosat and ERS-1 Altimetry: Ridge Segmentation versus Spreading Rate. Journal of Geophysical Research: Solid Earth, 114, 1-18.
 Sandwell, D., Garcia, E., Soofi, K., Wessel, P., Chandler, M. and Smith, W.H.F. (2013) Toward 1-mGal Accuracy in Global Marine Gravity from CryoSat-2, Envisat, and Jason-1. Leading Edge, 32, 892-899. https://doi.org/10.1190/tle32080892.1
 Sandwell, D.T., Müller, R.D., Smith, W.H.F., Garcia, E. and Francis, R. (2014) New Global Marine Gravity Model from CryoSat-2 and Jason-1 Reveals Buried Tectonic Structure. Science, 346, 65-67. https://doi.org/10.1126/science.1258213
 Ince, E.S., Barthelmes, F., Reißland, S., Elger, K., Förste, C., Flechtner, F. and Schuh, H. (2019) ICGEM—15 Years of Successful Collection and Distribution of Global Gravitational Models, Associated Services, and Future Plans. Earth System Science Data, 11, 647-674. https://doi.org/10.5194/essd-11-647-2019
 Pavlis, N.K., Holmes, S.A., Kenyon, S.C. and Factor, J.K. (2012) The Development and Evaluation of the Earth Gravitational Model 2008 (EGM2008). Journal of Geophysical Research: Solid Earth, 117, 1-38. https://doi.org/10.1029/2011JB008916
 Hirt, C., Claessens, S., Fecher, T., Kuhn, M., Pail, R. and Rexer, M. (2013) New Ultrahigh-Resolution Picture of Earth’s Gravity Field. Geophysical Research Letters, 40, 4279-4283. https://doi.org/10.1002/grl.50838
 Tozer, B., Sandwell, D.T., Smith, W.H.F., Olson, C., Beale, J.R. and Wessel, P. (2019) Global Bathymetry and Topography at 15 Arc Sec: SRTM15+. Earth and Space Science, 6, 1847-1864. https://doi.org/10.1029/2019EA000658
 Götze, H.J. and Li, X. (1996) Topography and Geoid Effects on Gravity Anomalies in Mountainous Areas as Inferred from the Gravity Field of the Central Andes. Physics and Chemistry of the Earth, 21, 295-297.
 Camacho, M. and Alvarez, R. (2020) Gravimetric Analysis of the Rifts and Volcanic Fields of the Jalisco Block, Mexico. Tectonophysics, 791, Article ID: 228577.
 Hildenbrand, T.G., Briesacher, A., Flanagan, G. and Hinze, W.J. (2002) Rationale and Operational Plan to Upgrade the U.S. Gravity Database. USGS Open File Report, 12 p. https://doi.org/10.3133/ofr02463
 Macleod, I.N. and Ellis, R.G. (2013) Magnetic Vector Inversion, a Simple Approach to the Challenge of Varying Direction of Rock Magnetization. 2013 Australian Society of Exploration Geophysicists, Petroleum Exploration Society of Australia (ASEG-PESA) 23rd International Geophysical Conference and Exhibition, Melbourne, 11-14 August 2013, 6 p. https://doi.org/10.1071/PVv2012n159news
 Ellis, R.G., de Wet, B. and Macleod, I.N. (2012) Inversion of Magnetic Data from Remnant and Induced Sources. 2012 Australian Society of Exploration Geophysicists (ASEG) 22nd International Geophysical Conference and Exhibition, Brisbane, 26-29 February 2012, 4 p.
 Alvarez, R. and Yutsis, V. (2015) Southward Migration of Magmatic Activity in the Colima Volcanic Complex, Mexico: An Ongoing Process. International Journal of Geosciences, 6, 1077-1099. https://doi.org/10.4236/ijg.2015.69085
 Alvarez, R. (2017) Mapping Geologic Interfaces that May Alter Seismic Wave Propagation in the Mexico City Basin. Geofísica Internacional, 56, 37-56.
 Guevara, R., Yutsis, V, Varley, N., Almaguer, J., Calderón-Moctezuma, A. and Guevara-Mansilla, O. (2021) Geophysical Determination of the Jalisco and Michoacán Blocks Boundary along the Colima Graben. Journal of South American Earth Sciences, 109, Article ID: 103208. https://doi.org/10.1016/j.jsames.2021.103208
 García-Palomo, A., Macías, J.L., Arce, J.L., Capra, L., Garduño, V.H. and Espindola, J.M. (2002) Geology of Nevado de Toluca Volcano and Surrounding Areas, Central Mexico. Geological Society of America: Map and Chart Series MCH089.
 Martínez-Serrano, R.G., Schaaf, P., Solís-Pichardo, G., Hernández-Bernal, M.S., Hernández-Treviño, T., Morales-Contreras, J.J. and Macías, J.L. (2004) Sr, Nd and Pb Isotope and Geochemical Data from the Quaternary Nevado de Toluca Volcano, a Source of Recent Adakitic Magmatism, and the Tenango Volcanic Field, Mexico. Journal of Volcanology and Geothermal Research, 138, 77-110.
 Ryan, W.B.F., Carbotte, S.M., Coplan, J., O’Hara, S., Melkonian, A., Arko, R., Weissel, R.A., Ferrini, V., Goodwillie, A., Nitsche, F., Bonczkowski, J. and Zemsky, R. (2009) Global Multi-Resolution Topography (GMRT) Synthesis Data Set. Geochemistry Geophysics Geosystems, 10, Article ID: Q03014.