New Bouguer Gravity Maps of Venezuela: Representation and Analysis of Free-Air and Bouguer Anomalies with Emphasis on Spectral Analyses and Elastic Thickness

A new gravity data compilation for Venezuela was processed and homogenized. Gravity was measured in reference to the International Gravity Standardization Net 1971, and the complete Bouguer anomaly was calculated by using the Geodetic Reference System 1980 and 2.67 Mg/m3. A regional gravity map was computed by removing wavelengths higher than 200 km from the Bouguer anomaly. After the anomaly separation, regional and residual Bouguer gravity fields were then critically discussed in term of the regional tectonic features. Results were compared with the previous geological and tectonic information obtained from former studies. Gravity and topography data in the spectral domain were used to examine the elastic thickness and depths of the structures of the causative measured anomaly. According to the power spectrum analysis results of the gravity data, the averaged Moho depths for the massif, plains, and mountainous areas in Venezuela are 42, 35, and 40 km, respectively. The averaged admittance function computed from the topography and Free-Air anomaly profiles across Mérida Andes showed a good fit for a regional compensation model with an effective elastic thickness of 15 km.


Introduction (Earlier Gravity Mapping in Venezuela)
Gravity surveys have been carried out in Venezuela since 1945 with the intensification of oil exploration.Thus, gravity surveying was first confined to the oil-producing sedimentary basins (i.e., Maracaibo and Eastern Venezuela basins).The first precise gravimetric survey covering the whole country was carried out as late as 1970 under the framework of the Latin American Gravity Standardization Network [1], which established a national network.The Venezuelan National Geographic Institute published the first gravity map of the entire country after 1988 [2].Since then, several Venezuelan universities and official institutes have cooperated with international institutes to improve the coverage of the national network.In addition, the Venezuelan National Oil Company (PDVSA) has released a large amount of data for educational purposes.In this context, Izarra et al. [3] presented the last gravity data compilation from Simon Bolívar University in 2005.However, due to the extent of the country and the presence of inaccessible areas such as the Amazon forest and Mérida Andes, general improvement in the data coverage has been slow.The purpose of this work is to present a new Bouguer anomaly map of Venezuela using the data available thus far.This study has four main aims: (1) to study the correlation between the Bouguer anomalies and the known regional geology by means of regional and residual gravity anomaly separation; (2) to estimate the continental crust thickness (Moho) using a spectral technique and to correlate it with Moho estimations derived from independent geophysical techniques (i.e., P-wave velocity models, which stem from wide-angle refraction seismic and the results from receiver function methods); (3) to use the coherence-admittance International Journal of Geophysics techniques to determine the spatial variations in elastic thickness (Te) in the mountainous areas; (4) to provide a new updated gravity database for future investigations.

Geological Setting
Different theories have been proposed about the origin of the Caribbean Plate [4,5].Nowadays, it is commonly accepted that extensive volcanism in the mid to late Cretaceous (∼90-110 Ma) resulted in the formation of the Caribbean Large Igneous Province and a thickening of large areas of the Caribbean crust.Around this time, the Caribbean plate (CA) initiated an eastward migration relative to South America (SA) [6][7][8][9][10].The resulting variety in geologic features between the CA and SA margins includes accretionary wedges, fold and thrust belts, and extensional and foreland basins [11].Current relative plate motions make this margin predominantly strike-slip; it is limited on the east by subduction of oceanic (Atlantic) SA beneath the Lesser Antilles island arc [12] and in the west by partial subduction of the Caribbean plate under South America (Figure 1).On the continental deformation front observed along the Venezuela coastline, the important (large oilproducing) Barinas-Apure and Eastern basins have been formed.Barinas-Apure is mainly a foreland basin generated by the flexural response to the Mérida-Andes Mountain load [13].On the other hand, the Eastern basin can be considered to be the result of (a) flexural loading of the Cordillera de la Costa range, (b) large and continuous deposition of continental material from the Guyana shield, and (c) subduction dynamics in the east [14].The Barinas-Apure and Eastern Venezuelan sedimentary basins are separated by a geomorphologic-structural high called El Ba úl High, formed by an igneous-metamorphic complex.This massif was described by Bellizzia [15], who differentiated several granitic, volcanic, and metasedimentary units and subunits.For instance, the Maracaibo basin is a foreland basin separated from Barinas-Apure by the emplacement of the Mérida Andes.This basin is bounded on the north by the Oca-Ancon fault system, Sierra de Perijá to the west, and the Mérida Andes mountains to the south and east.Mérida Andes (MA) is a NE2-SW trending mountain chain about 420 km long and with a maximum height of 5 km.This uplifted block was formed as a consequence of the convergence of the Panama arc and western South America [16][17][18][19].According to Schubert [20] and Kohn et al., [21] this convergence has the main period of shortening during the Oligocene-Miocene, and evidence suggest that there is still significant present-day deformation.
Due to the large amount of interactions and its complexity between the CA-SA plates, the location of the plate boundary is controversial.Thus, this boundary has been interpreted as a 300 km wide plate boundary zone [22] that forms an orogenic float.This orogenic float is represented by the Mérida Andes; it is linked to the north with the Northern Cordilleras of Venezuela (i.e., Cordillera de la Costa range) and goes as far as the Trinidad-Tobago islands to the east.The boundary zone is limited by the Southern Caribbean thrust system on the north and the Mérida Southern Foothills, Guárico Frontal thrust belt, and Serranía thrust belt on the south [19].The Cordillera de la Costa range is formed by different geological provinces containing (from north to south) late Jurassic to early Cretaceous basic and ultrabasic rocks, Precambrian and Paleozoic basement rocks, Jurassic to Cretaceous lower crust-upper mantle fragments, volcanosedimentary sequences and basaltic to rhyolitic rocks, and late Cretaceous-Paleocene molassic sediments and flysh sequences [23][24][25][26][27][28].
The stable South American crust is represented by the Guyana Shield.This massif outcrops as sialic Precambrian continental crust composed mainly of metasedimentary and metaigneous rocks at amphibolite to granulite facies that have been intruded by granites [29].The reported ages of these crystalline rocks range from 3600 to 800 Ma [30].

Gravity Database
A new gravity data compilation is presented here that includes data from early compilations by Simon Bolivar University [3] and data from the National Geophysical Data Center (NGDC).In order to improve the quality and coverage of onshore data, new high-resolution gravity surveys from the Venezuelan Foundation for Seismological Research (FUNVISIS) were included.Offshore data coverage was improved by new measurements from several marine surveys provided by the Marine Geoscience Data System (MGDS) [31].
(1) The Simon Bolivar University database consists of datasets from different Venezuelan institutions, oil companies, universities, and international surveys.Most of the data are evenly spaced; however, the compilation includes data with high accuracy (1 × 10 −5 m/s 2 or less for gravity and 0.5 m or less for height) and an average station spacing of less than 100 m.This dataset has coordinates measured through the use of precise leveling methods, which were mainly collected by oil companies.In contrast, for data collected by universities, in some cases, the stations coordinates and heights were derived by reading from topographic maps.These data have elevation errors of several meters and errors of higher than 1 × 10  (http://www.marine-geo.org/).The collected data include research cruises from 1977 to 2004.
(4) The FUNVISIS dataset includes more than 4000 observations collected by a Scintrex CG-5 gravimeter combined with GPS system.The resolution and accuracy of this gravity data are 0.5 × 10 −5 m/s 2 for gravity and ±1 m for station heights with average station spacing of 500 m.
Presently, the gravity compilation contains about 80,000 onshore observations and more than 40,000 offshore stations (Figure 2).The average station interval is less than 1 km, which results in an average station density of 1 station/km 2 or higher.

Data Processing.
A comparison of the different datasets showed that the gravity datasets and surveys in the compilation refer to different datum levels and exhibit variable quality and accuracy.Therefore, data homogenization focusing on the gravity datum and calibration and the coordinate determination and anomaly equations was required for anomaly reduction.After processing, the data were manually edited to remove erroneous measurements.Stations with outlier gravity values were removed after interpolating a high resolution gravity map.Additionally, stations with erroneous coordinates and/or heights have been removed after comparing the station heights with those obtained by a high resolution digital terrain model (90 m spacing).Stations with heights differences of 50 m or higher were also removed.In total, more than 2000 stations were eliminated from the original database [3].
The Bouguer anomaly was calculated using the following assumptions.
(i) The horizontal coordinates and elevations of the gravity stations based on the Geodetic Reference System 1980.
(ii) Absolute gravity datum is referred to International Gravity Standardization Net 1971 (IGSN71).
International Journal of Geophysics (iv) Topography correction calculated for a spherical cap of up to 167 km radius (Hayford zone O2) [33] assuming a constant density of 2670 kg/m 3 .This value is close to the mean density of the surface rocks in the investigated area and the standard value used in the Bouguer anomaly correction.The digital terrain model used was based on the Shuttle Radar Topography Mission (SRTM) with a grid spacing of 90 m (onshore), and Gtopo30 (offshore).
(vi) Height correction estimated by a Taylor series expansion of normal gravity up to 2nd order [35].
Taking into account all the different sources of errors in the databases, the accuracy of the computed Bouguer anomaly values was estimated to be about ±1−5 × 10 −5 m/s 2 (1 × 10 −5 m/s 2 = 1 mGal).Additional details about the Venezuelan gravity network and the dataset can be found in Drewes et al. [36] and Izarra et al. [3].

Gravity Maps
The calculated anomaly map consists of Bouguer anomalies (BA) onshore (correction density of 2.67 Mg/m 3 ) and Free-Air anomalies (FAA) offshore (Figure 3).Anomaly

Regional and Residual Gravity
Maps.In order to analyze the anomalies, a set of wavelength filters were applied to progressively separate local effects from regional effects within the gravity field.Here, the concept of the Butterworth band-pass filter was applied in the frequency domains to separate regional and residual fields.The Butterworth filter is a spectral domain filter with a roll-off and requires an order (n) to implement the transition between the passing and rejected portions of the data spectrum; the higher the order, the steeper the transition.Based on the analysis of the radial power spectrum (Figure 4), the wavelength cutoff adopted here retained short-wavelength components for the observed data of 200 km, and the order was set at 3. Figure 5(a) shows a map of the regional gravity field obtained by applying this filter to the Fourier transform of the demeaned and detrended observed gravity data after conversion to the spatial domain taking the inverse Fourier transform.The residual gravity field (Figure 5(b)) was the result of subtracting the regional gravity field from the observed Bouguer anomaly.The gravity map of the regional component of the gravity field (Figure 5(a)) showed anomalies between +125 and −125 ×10 −5 m/s 2 .These anomalies are characterized by a WNW-ESE relative gravity low.The residual gravity map (Figure 5(b)) ranges from +125 to −120 × 10 −5 m/s 2 .The distribution of anomalies is similar to the observed Bouguer anomaly map, and their main features are the strong gradients associated with the Caribbean-South America boundary fault system.The maximum gravity lows are related to basins in the study area (e.g., Granada and Eastern basins), and the maximum gravity highs related to several local geological features such as volcanic and metamorphic outcrops in the ABC islands.These gravity highs and lows show predominant orientation following the Antilles arc region.

Power Spectral Analysis
The application of Fourier analysis to the interpretation of potential field data is common and is frequently used to obtain the regional/residual field components of the gravity field.Here, the complete Bouguer gravity anomaly was used to get estimated depths for the structures that cause the measured anomaly.The methodology assumed that for large samples, the logarithm of the power spectrum (E) of the gravity field of a monopole source versus the wave number (radial/distance) may show a linear relationship.The slope of the straight line is proportional to the depth to the top of the corresponding body causing the gravity anomaly.Thus, if k denotes the wavenumber and S(k) denotes the power spectrum of the gravity field, the depth (d) to the source can be estimated from the relation The power spectrum analysis was carried out through the 2D fast Fourier transformation of the gravity field.Due to the two-dimensional character of the dataset, radial averaging of the power spectrum was performed to obtain a onedimensional representation [39][40][41][42][43]. Confidence limits for the depth estimations were calculated from the standard errors of the slopes of the best fitting lines for the linear segments.
The gravity data of the study area, covering a surface area of 900 × 900 km 2 , was interpolated to produce a grid with a node spacing of 4 km.Results of the power spectrum analysis for this dataset showed four tendencies for the correlation between the energy (E) and wavenumber (radial/distance) (Figure 6(a)).The most regional part of the spectrum resulted in greater depths of about 79 km.The local part of the spectrum resulted in depths of ∼16 km.The main aim of the spectral analysis was focused on the intermediate depths In order to perform the same analysis over all anomalies in the study area, a data block (window), measuring 300 × 300 km each, was used for the calculation of Moho depths across the N-S and E-W directions.The window size of 300 km 2 corresponds to six-times the expected source depth (Moho) to assure a depth-estimation error of <10% according to Regan and Hinze [44].Figures 6(a

Effective Elastic Thickness Determination
The correlation between gravity observations and sea floor topography has been well documented.This relationship can be analyzed by using a linear transfer function called admittance.The basic techniques for determining the admittance function between the gravity and bathymetry data have been discussed in detail by McKenzie and Bowin [45] and Watts [46].
The method assumes that Free-Air gravity anomalies are caused by topography, and its compensation and attempts to determine a function when convolved with topography produce the gravity anomaly.The advantage of this method is that the admittance function may be derived from the observed data independent of a particular isostatic model and can be interpreted in terms of simple models of isostasy.
The isostatic response method simply involves deriving a filter that, convolved with the bathymetry in space domain,   produces a series that resembles the observed gravity, again in the space domain.This process can be represented in the space domain by using the convolution operator * : where h(x), f (x), and g(x) are the series representing the topography, filter, and Free-Air gravity, respectively.The above convolution in the space domain is equivalent to multiplication in the spatial frequency domain, where G(k), Z(k), and H(k) are the discrete Fourier transforms of g(x), f (x), and h(x), respectively.Z(k) is known as the admittance or transfer function and the wavenumber k = 2π/λ = n2π/L, where n = 0, 1, 2, . . ., L/2Δζ, L is the length of the Fourier series, Δζ is the distance between two consecutive sampling points, and λ is the wavelength.Equation (3) can be rewritten as: However, the function Z(k) obtained in this way is influenced by noise in the gravity field, particularly at short wavelengths.The noise present in any data may be due to measurement errors, the data reduction procedure, or variations in the structure of the linear system under consideration.In the presence of noise, a better estimate of Z(k) [45] can be obtained from where N is the number of profiles, and C(k) is the complex cross-spectrum of bathymetry and gravity.Eb(k) is the power spectrum of bathymetry.* denotes the complex conjugate.
The quality and reliability of the admittance amplitudes were controlled by means of four additional parameters: the coherence, phase of admittance, and the coherent and incoherent energies: where ϕ(k) is phase of admittance.Im[Z(k)] and Re[Z(k)] are the real and imaginary parts of Z(k).γ 2 (k) is the coherence.Eg(k) is the power spectrum of gravity.The observed admittance curve was compared with a set of theoretical admittance curves for the Airy and Flexure models, which contain the mean crust (Tc) and elastic plate thickness (Te) parameters.The final values of Te and Tc were obtained through the selection of the lowest mean square error between the observed admittance curve and each of the theoretical curves.
The theoretical curves for the admittance of the Free-Air anomaly for Airy isostatic compensation (Z(k) Airy ) and for the flexure or plate isostatic compensation model (Z(k) Flex ) were calculated following the method given by McKenzie and Fairhead [47], where G is the gravitational constant, E = is Young's modulus, υ is Poisson's ratio, g is acceleration due to gravity, ρ c and ρ m are the average crustal density and density of material below the assumed flexed elastic plate, Tc is the effective depth of compensation, and Te is the effective elastic thickness.
For the admittance analysis, in the present work, three areas were selected: Mérida Andes (Zone A), the Cordillera de la Costa range (Zone B), and the Guyana shield (Zone C).For each area, a maximum of nine profiles were extracted from the Free-Air anomaly map and topography map, each of ∼500 km (Figure 9).Elevation data profiles were extracted from the Shuttle Radar Topography Mission (SRTM) with a grid spacing of 90 m.Each dataset along each profile was regularly spaced at intervals of 0.5 km.The coherence, phase of admittance, and energies (7) for each area were plotted with respect to wavenumber and are shown in Figure 10.
The observed admittance was computed from (5) for 2n samples, where n was taken as 9, which correspond to half of the longitude of the regularly spaced profiles.Theoretical admittance curves were computed (9) for Te values between 5, 10, 15, 20, 25, and 30 km.Te values were obtained after studying the square medium errors (RMS) of the correlation between observed and theoretical admittance curves.A best fit was considered whenever RMS dropped below a minimum value with the correlated theoretical Te value.
Finally, the observed and theoretical admittance curves were plotted with respect to wavenumber (Figure 11). 5 can be used to identify gravity anomalies associated with geological bodies approximately 50 km or greater in size.These maps are therefore suitable for correlating gravity with tectonic provinces having dimensions of several hundred kilometers.In addition, in the southern part of Venezuela, large areas have poor coverage of gravity observations, especially in the Amazon region (i.e., Guyana shield).

Regional and Residual Gravity Anomalies and Their Correlation with Major Tectonic Provinces. The 3 min resolution maps shown in Figure
Considering major topography and gravity changes, the study areas were separate in major topography/gravity provinces.The offshore provinces encompass the Venezuelan basin, South Caribbean accretionary prism, Aves Ridge, Grenada back arc basin, Lesser Antilles arc, Bonaire basin, and Tobago trough.The onshore provinces include the Guyana shield, Barinas-Apure basin, Maracaibo basin, Falc ón basin, East Venezuela basin, Perijá-Mérida, and coastal ranges.
Large areas of negative residual anomalies onshore (Figure 5(b)) are associated with Mesozoic and Cenozoic sedimentary basins (Figure 1).These intracratonic Mesozoic basins such as the Barinas-Apure, Maracaibo, Falc ón, and East Venezuela basins have distinct gravity signatures.For instance, the Barinas-Apure basin is characterized by long wavelength anomalies with low amplitude (−20 to +20 × 10 −5 m/s 2 ).Another negative residual anomaly, near inland, but not as large in magnitude, is a FAA low through Bonaire basin.This basin has a quasirectangular shape with low negative residual anomalies which are emphasized by the surrounding gravity highs associated with the coastal ranges and the Leeward Antilles arc.
The negative Bouguer anomaly observed in the Eastern Venezuelan basin (Figure 3) is characterized by a SW-NE  8. Red lines indicate cutoff values for coherence, phase, and coherent (thick line) and incoherent (thin line) energies.(a) For Mérida Andes, the coherence is high for k < 0.04, the estimated phase is close to zero for k < 0.08, and the coherent energy is higher than the incoherent energy for values of k < 0.08.trending of a negative residual anomaly (∼ −50 × 10 −5 m/s 2 ) (Figure 5(b)).This gravity low is associated with Cenozoic sediments, which accommodates the result of thrust-sheet loading (i.e., Serranía thrust belt) that forced the American continental lithosphere to flex downward between the Guyana shield and El Pilar fault [48].Schmitz el al. [49] used seismic refraction data to estimate a maximum of 10-13 km of sedimentary infill for this basin.According to Jácome et al. [14], the amount of sediments is the result of multiple processes: (a) flexural loading of the Cordillera de la Costa range, (b) large and continuous deposition of continental material from the Guyana shield as the main source, and (c) the subduction dynamics of oceanic South America underneath the Caribbean plate.The spectral analysis of the Bouguer anomaly (Figure 6(d)) estimated an anomaly source located 12 km deep that may correspond to the top of the basin basement.
In the Mérida Andes, the residual gravity map ( Maracaibo basin.These gravity highs are associated with exposed Archean to Lower Proterozoic high-grade metamorphic basements within the mountain range.Relative maximums follow dikes (diabase) that intruded the Precambrian and Paleozoic rocks of Mérida Andes, which formed as a consequence of the convergence of the Panama arc and western South America, which initiated the formation of this mountain range [16][17][18][19].Negative residual anomalies flanking the Mérida Andes showed similar shapes but very different amplitudes.These anomalies are mainly caused by strong flexes of the crust and the infill sediments of the Barinas-Apure basin in the south, and the Maracaibo basin in the north, in association with a chain-scale thrust and back thrust systems along the Mérida foothills.The northern flank has lower amplitude (< −40-40 × 10 −5 m/s 2 ) with a relative maximum located in the northwest.In contrast, the southern flank barely reaches −40-40 × 10 −5 m/s 2 in the most southeastern area of the anomaly.These anomalies were interpreted by F. E. Audemard and F. A. Audemard [19] to be the consequence of the rheological characteristics of two different continental crusts (the crust underneath the Maracaibo basin and the South American crust underneath the Barinas-Apure basin).In other words, the Maracaibo crust underwent a more recent tectonic and thermal event due to continental rifting during the Jurassic than the longcooled Precambrian crust of the SA craton in the south.In summary, the Andes of Venezuela is a floating orogen involving incipient, gently NW-dipping continental subduction that generates a shallow foreland basin on the Barinas-Apure basin side, while it strongly flexes the Maracaibo crust on the forearc-equivalent side, where a deep flexural basin develops in association with a chain-scale backthrust along the Maracaibo foothills on the northwest [19].
Other regions where positive residual anomalies were observed were the Lesser Antilles arc and the Aves ridge.They were characterized by shorter wavelengths (<100 km) but high amplitudes (>100 × 10 −5 m/s 2 ) in the residual anomalies map.These pronounced gravity highs and steep gradients were the combined effect of the steep slope of the bathymetric and high-density volcanic rock outcroppings in these localities, which displayed a wide range of compositions.
The regional gravity anomaly map provided a general view of the areal extent of the gravity response from the major geological units and reflects the gross crustal structure of the area.The calculated regional anomaly map displayed several well-defined gravity highs and lows of varying dimensions and relief.For instance, under the Guyana shield, the regional gravity field showed no evidence of significant crustal thickness variations (Figure 5(b)).This result is consistent with seismic studies that indicate a crust thickness of ∼45 km [38].In this region, regional and residual gravity maps suggest that the observed anomalies can be explained by upper crustal lateral density variations alone rather than by changes in the thicknesses of the lower crust and Moho.This is supported by the small range of ∼20 to ∼ −20 ×10 −5 m/s 2 for the regional field in this area.
A conspicuous gravity high (∼40 × 10 −5 m/s 2 ) in central Venezuela trending from NE to SW separates the Barinas-Apure basin and Eastern Venezuela basin.This anomaly extends from the El Ba úl high to granitic rocks of the Santa Rosalía complex in the Precambrian Guyana shield, which crop out in the southwestern margin of the Orinoco River.In addition, Viscarret et al. [50] determined U-Pb zircon ages and interpreted the El Ba úl massif to be a Paleozoic basement belt, which shows more affinity with the Mérida Andes than with rocks of the Guyana shield.Therefore, this long wavelength anomaly is very likely created by the overlapping of anomalies caused by the El Ba úl high and the Santa Rosalía complex.
After the separation of the regional component of the BA, the Eastern Venezuelan basin still has the most relevant anomaly.This pronounced negative Bouguer gravity anomaly has been studied by several authors [3,[51][52][53][54].The negative gravity anomaly is roughly parallel to the arc platform extension, but it does not extend west of Gulf of Paria.According to gravity and recent seismic refraction studies, this negative gravity anomaly is indicative of the very large load on the South American lithosphere here as well as the large amount of sediments in the basin.
Regional gravity anomalies of the South Caribbean accretionary prism and Leeward Antilles arc (i.e., ABC Islands) form a positive-negative gravity pair characteristic of subduction zones.The high-density rock in the Leeward Antilles islands are characterized by residuals with very high amplitudes of 100-120 × 10 −5 m/s 2 , but regional field anomalies also have very high amplitudes.Whereas, the gravity anomalies caused by strong density contrast in the masses that form these volcanic bodies, that are or are being accreted to the continent, produce local anomalies that are very emphasized by the steep bathymetry.This interpretation is similar to the one by Bonini et al. [51], based on gravimetric data.

Effective Elastic Thickness and Moho
Depth.The admittance and coherence results presented here for the Mérida Andes (Figure 11) best fit a flexural model with Te = 15 km.The Airy model fits the observed admittance when the mean crustal thickness is 70 km.This value is incompatible with the estimates of the typical continental crust thickness (35-40 km).Consequently, the Airy model cannot be accepted as a possible isostatic compensation mechanism.A simple plate flexure model with Te = 15 km is more reasonable.
Along the Cordillera de la Costa range, the admittance calculation agrees more closely with the theoretical values calculated for a plate thickness of 10 km.Te is around ∼10 km north of the Eastern Venezuelan basin (i.e., along the Cordillera de la Costa range), indicating that low rigidity amplified the subsidence of the basin.
Coherence, phase of admittance, and coherent and incoherent energies calculated for the Cordillera de la Costa range and Guyana shield showed poor quality and reliability for Te estimations.In other words, the phase of the admittance was not close to zero, the coherence was less than 0.5, and incoherent energy was higher than coherent energy for all wavenumbers in the selected profiles (Figure 9).The admittance estimates for the set of profiles crossing Mérida Andes showed a noticeable data scattering at wavelengths longer than 160 km.This scattering could be caused by the wide variability of the thickness of sediments and also by very short wavelength anomalies that could be considered as noise.
Tassara et al. [55] used a wavelet formulation of the classical spectral isostatic analysis to invert satellite-derived gravity and topography/bathymetry for Te over South America.According to their calculations, Mérida Andes, which is the region located between the coastline (i.e., Cordillera de la Costa range) and Colombian Eastern Cordillera, exhibits variable Te from 35 to 40 km with an uncertainty of ±10 km.Beneath the Cordillera de la Costa range, Te decreases to ∼35 ± 10 km.This high Te value in the Cordillera de la Costa range may reflect the combined effect of the strength of the upper continental and the partially subducted Caribbean oceanic lithospheres.Beneath the Eastern Venezuela basin, the elastic thickness ranges between 40 and 50 km with an uncertainty of ±10 km.The Guyana shield shows Te as low as 10 km, although with an uncertainty of up to 25 km.
The Te values estimated in this study did not match the elastic plate thicknesses estimated by Tassara et al. [55].These discrepancies between the admittance results are the most probable cause for the presence of a considerable amount of loading from the top, which corresponds to stacking of thrust sheets and uplifts of basement rocks in coastal ranges, and from the bottom, that is, the partially subducted Caribbean slab.Another factor to consider is the size of the window used to compute admittance.According to Pérez-Gussinyé et al. [56] windows that are too small introduce spurious spatial variations, and windows that are too large tend to average the spatially varying Te values and smooth the true structure.These discrepancies could be caused by the large and tectonically heterogeneous area required by the admittance technique, which would tend to bias their result towards a weaker plate [57].
The results obtained in the spectral analysis confirm the previously established value of 35 km as the mean reference Moho depth [38].In addition, the crust thickness is not homogeneous.For instance, the Moho topography shows a NE-SW depression beneath the Mérida Andes with a maximum depth of ∼55 ± 5 km; this value gradually decreases toward the south and north of the main strike of this feature.The maximum crustal thickness in Mérida Andes (Figure 7) was partially constrained by the receiver function analysis [37], which shows Moho depths of 50 and 47 km near this maximum.The Maracaibo basin seems to have a crustal thickness of 35-40 km.The estimated Moho geometry under the area of the Falc ón basin agreed well with the refraction seismic data modeled along the 70 • W profile, which indicates a crustal thickness of ∼27 km [58].The crustal thickness beneath the Eastern Venezuela basin ranges from ∼50 to ∼40 km.The Guyana shield is regionally underlain by a crust of ∼45 km.Northwards beneath the Guárico and Barinas-Apure basins, the crustal thickness reaches ∼40 km.Along the coastline, the Caribbean Mountain System crustal thickness oscillates around 30 km and decreases slightly toward the Leeward Antilles (Figures 7  and 8).
Figure 8 show that there are significant differences among Moho values estimated with spectral analysis and Moho estimated with receiver function.On the other hand, Moho depth values in Figure 7 and closely match the values estimated by Schmitz et al. [38] and also produces the best fits to the observed Bouguer gravity data [54].According to Schmitz et al. [38] the mismatch between seismic Moho and receiver function Moho is caused by the resolution of these technics.

Conclusions
A new Bouguer gravity map of Venezuela was developed based on an up-to-date dataset available in the country.All data were reprocessed and homogenized according to gravity processing standards, and the effects of the topography were corrected with digital topographic maps.The final dataset was comprised of more than ∼80,000 observation points that are now available for future detailed interpretations and future crustal investigations, such as 2D and 3D gravity modeling.
The results obtained in the Bouguer gravity map presented in this work can be greatly improved by adding more gravity data in areas where scarce and poor quality data are the only data available, such as in mountain ranges (e.g., Andes) and the Amazon forest (e.g., Guyana shield).
The Free-Air gravity and topography admittance analysis of the data windows over the Mérida Andes provided elastic thicknesses in the range of 30-35 km.
Finally, Te values presented in this work using the admittance method were found to be lower than Te calculations carried out in previous studies on the scale of South America and surrounding plates.These differences in the Te calculation could be caused by tectonical heterogeneities and problems associated with the admittance technique (window size, data coverage and presence of noise), but these arguments are still in discussion.It is possible that satisfactory results for this area that has a very complicated structure could only be obtained by a more sophisticated approach.

Figure 1 :
Figure 1: Study area and general tectonic and kinematic features of the northern margin of South America setting (modified from F. E. Audemard and F. A. Audemard [19]) on a digital elevation map.Arrows indicate relative regional motions.Texts indicate the name of geological or tectonic features.A, B, and C mean Aruba, Bonaire, and Curazao (ABC islands).

Figure 2 :
Figure 2: Station distribution map of gravity stations used for the new gravity anomaly map of Venezuela.
values range from -225 to 225 × 10 −5 m/s 2 .Offshore, in the Venezuela Basin (VB), magnitudes of the FAA barely reach low positive values.In contrast, the South Caribbean accretionary prism shows a broad gravity low with a WNW-ESE trend.A gathering of local anomaly highs separated by relative low gravity values appears as local highs along the Leeward Antilles (Aruba, Bonaire, and Curazao).The highest FAA values (more than 200 × 10 −5 m/s 2 ) are caused by the subduction in northeastern Venezuela (i.e., Lesser Antilles volcanic arc), which extends further eastward.In the area of the Grenada basin (GrB) and Tobago trough (TT), gravity values decrease down to −80 × 10 −5 and −45 × 10 −5 m/s 2 , respectively.Onshore, the most relevant BA is the prominent anomaly low observed in the Eastern Venezuelan basin (EVB) mostly caused by a large amount of infill sediments.The BA has NE-SW trending positive values in the area between the Guárico basin (GB) and Guyana shield (i.e., Precambrian rocks).

Figure 5 :
Figure5: (a) Regional gravity map obtained by applying a Butterworth filter defined by results of the spectral analysis (Figure4) of the gravity anomaly map.(b) Residual gravity map obtained by subtracting the regional gravity map from the observed gravity map in Figure3.Both maps: contour interval = 5 × 10 −5 m/s 2 .Graphical indications as described in Figure3.
)-6(d) show results for windows over the anomalies observed in the Barinas-Apure basin, Mérida Andes, Maracaibo basin, Guárico basin, Eastern basin, and Leeward Antilles.The final results are shown in Figures 7 and 8.

Figure 6 :
Figure 6: Results for the power spectrum analysis of the Bouguer gravity field using 2D fast Fourier transform.Black lines represent the trends of the linear regression separating the spectrum of mass points (blue points) located at a maximum estimated depth: (a) Barinas-Apure basin, (b) Mérida Andes, (c) Maracaibo basin, (d) Guárico basin, (e) Eastern basin, and (f) Leeward Antilles.Estimated depth (Sd).

Figure 7 :
Figure 7: Map of source-depth estimation (Moho) obtained by spectral analysis of gravity data.Triangles are the receiver function Moho depth estimations in kilometers [37].Gray lines are Moho depths estimated from seismic studies [38].Location of profiles are shown as black lines labeled by 70 W and 65 W.

Figure 10 :
Figure 10: Observed admittance values computed from FAA and topography for areas (a), (b), and (c) in Figure 8. Red lines indicate cutoff values for coherence, phase, and coherent (thick line) and incoherent (thin line) energies.(a)For Mérida Andes, the coherence is high for k < 0.04, the estimated phase is close to zero for k < 0.08, and the coherent energy is higher than the incoherent energy for values of k < 0.08.

Figure 11 :
Figure 11: Observed admittance values computed from the FAA and topography profiles in Figure 9 for low wavenumbers, where k = 0.0-0.04.Observed (solid circles) and theoretical (solid curves) admittance for the (a) Flexure model, where Te is the effective elastic thickness, and (b) Airy isostatic model, where Tc is the mean thickness of the crust as a varying parameter of the model.Admittance errors were computed from the coherence [46].The best fitting curve for Mérida Andes was the simple flexure model of the lithosphere with Te = 17 ± 3 km.