Se obtuvo un modelo tridimensional de la corteza mediante un proceso de inversión de datos gravimétricos para la región oriental de Cuba. Los datos y el modelo cubren un área rectangular de 64 600km2. El modelo inicial fue constreñido con la geología de superficie, la información sísmica y de perforación. Se aplicó un algoritmo de inversión que utiliza los datos de gravedad para estimar las topografías 3D a partir de las unidades geológicas principales. El modelo nos proporciona información cuantitativa sobre las profundidades y espesores de las formaciones geológicas más importantes. En el mismo se observan las secuencias alóctonas de diferente composición y origen sobre el basamento carbonatado de la Plataforma de Bahamas. La mayoría de los máximos en la anomalía de la gravedad se deben a la presencia de mantos más densos de ofiolitas poco profundas. Se destaca al suroeste el máximo gravimétrico provocado por la presencia de la corteza oceánica más densa generada en el Centro de Dispersión de Caimán.
A three-dimensional crustal model for Eastern Cuba, obtained through a process of gravity data inversion is presented. The study area cover a rectangular area of 64 600km2. The initial model for the inversion was constrained by surface geology, seismic and drilling data. The inversion algorithm uses gravity data to estimate 3-D topographies from the main geological units. The model provides quantitative information on the depths and thicknesses of the geological formations. The resulting model provides new information about the regional composition of the crust. Alien sequences are observed with different compositions and origin over the basement of Bahamas carbonate platform. Most of the maximum gravity anomalies are associated with presence of dense shallow ophiolite sheets. The most remarkable detail is the gravity “southwest” maximum, related to the presence of denser oceanic crust generated in the Cayman spreading center.
Models of genesis and evolution for the Caribbean-Cuba region show little agreement; in-situ (Giunta et al., 1997, James, K.H; 2003) and allochthonous (Pindell et al., 1990, 2009; Iturrade-Vinent, 1998, 2002; Cobiella, 2005; García-Casco et al., 2008 and Sommer et al., 2011). Allochthonous models indicate the area is formed by fragments of the ancient Caribbean plate that overrode the continental Bahamas margin. Over this margin are the Paleocene Volcanic Arc (PVA) and Neogene-Quaternary sediments forming basins. Other studies took into consideration seismic refraction, gravity, deep drilling boreholes and satellite images. Studies from 30 years ago showed that the island is on transitional crust, with 17 to 30km thick and characterized by three layers: an upper volcanic-sedimentary layer with P velocity of 4.0-4.8km/s, a lower layer with P of 5.8-6.4km/s and the deepest layer with P of 6.3-6.7km/s (Bovenko et al., 1982; Otero et al., 1998; Bush and Shcherbakova, 1986). Alternatively, Eastern Cuba is constituted to the south by an oceanic crust and to the north by continental crust (Tenreyro et al., 1994). Otero et al. (1998) based on re-interpreted seismic and gravity data proposed that the crust southward of Cauto-Nipe (Figure 1) is oceanic, ~20km thickness, and below the basin there is a fine-transitional crust with 20 to 30km thick and northward continental crust. Through teleseismic data, Palau et al. (2006) determined the presence of a 1km thick shallow layer with P velocity of 3.6km/s, underneath a 6km thick layer with P of 5.8km/s and a deeper 13km thick layer with 6.9km/s for P. Recently, Gonzalez et al. (2011) defined a 16 to 30km earth crust thickness for eastern Cuba, through joint inversion of Rayleigh waves dispersion and receptor functions. All these models are assumed one-dimensional or horizontal stratified.
The study area is located at southeastern Cuba Island. The area was divided in a regular grid of 646 prisms of 10kmx10km at surface. One grid represents one geological unit. Nine grids are located stratifie. Density contrast is kept constant for every unit. OFZ: Oriente Fault Zone. NHFZ: North Hispaniola Fault Zone.
Such varied results and explanations indicate that a further research is required. The main purpose of our research is to obtain the structure of the earth crust for Eastern Cuba using a dense gravity data. We extended a method that was previously used in a reduced area of Eastern Cuba with mining purposes and using only magnetic data (Batista et al., 2007). We used the Gallardo et al. (2003) algorithm, which minimizes the quadratic norm of differences between gravity data and the model response, constraining the solution or model with the surface geology, boreholes (~3.5km as maximum) and seismic reflection profiles. Every unit is simulated with a regular mesh of 10kmx10km prisms. Along the iterative inverse process, the depths to every prism is moved automatically in order to fit the gravity data and restricted to obey the constraints imposed.
The study area is a rectangle of 190km by 360km that covers the southeast of Cuba. A large land area and two pieces of oceanic crust are involved; the Atlantic Ocean and the Caribbean Sea where the Oriente Fault Zone (OFZ) is located (Figure 1). The area limits in “Southern Cuba coordinates system” are (100,000; 290,000) m North and (420,000; 760,000) m East. The rectangular region was divided into a grid of prisms with 10kmx10km on surface as shown in Figure 1.
Geologic contextIturralde (1998) recognizes two levels of the geological structure of Cuba: the substrate folding and the neo-autochthonous. Each consisting of different geological units.
The existing folded substrate is formed by pieces from the North-American and the ancient Caribbean and Pacific plates. Neoautochthonous units are sediments from the Neogene-Quaternary.
In our study area (Figure 2), the folded substrate is composed by: ophiolites, Cretaceous and Paleogene Volcanic Arcs. Ophiolites (West-East strips) located at the north, are over the Bahamas platform and under the Cretaceous Volcanic Arc (CVA). Sometimes, ophiolites and CVA are mixed forming an ophiolitic mélange. It is assumed that this ophilites were emplaced when the collision between the extinct CVA and the Bahamas Platform occurred.
Outline of the different geologic units and faults presented in the study area (Iturralde-Vinent., 1998).
Ophiolites eastward the area (Figure 2) are located in the Mayarí-Sagua-Baracoa massif. These ophiolites sheets are over olistostromes and over the CVA.
The Albian-Campanian volcanic arc is Cretaceous in age. It lies in tectonic contact with the northern ophiolites. Near the contact, the arc rocks are even more deformed, with fissured and foliated areas, and with chaotic masses that contain a mixture of ophiolites blocks, vulcanite and plutons.
For simplification, we refer to this volcanic arc as CVA. This is constituted by volcanogenic-sedimentary complexes (calc-alkaline and alkaline composition), plutons and the metamorphic complex.
The Paleocene volcanic arc (VPA) is a characteristic of the oriental south portion of the island (our area). Its age goes from the upper Daniense to the lower Eocene in the western part of Cuba, Jamaica, Hispaniola, Puerto Rico and Virgin Islands. The Paleocene arc was formed over the deformed remains of the CVA-Ophiolites units. VPA is constituted by volcano-sedimentary and plutonic rocks with different composition.
Neo-autochthonous units are represented by sedimentary rocks originated from the upper Eocene to recent. Three sedimentary cycles can be recognized; first, a stadium of the upper Eocene to the Oligocene, second, the lower to upper Miocene, and the Pliocene to the recent.
MethodGravity data consist of a rectangular mesh of 340km East direction and 190km North direction, with interpolated data every 3km. This design is optimal because we are looking for low spatial frequency structures. This produces at least nine observations over every 10kmx10km prism. Data was collected and processed by Instituto de Geología y Paleontología de Cuba (IGP) and represents the complete Bouguer anomaly (Blakely, 1996), using 2.3gr/cm3 for the earth crust density (Figure 3).
Complete Bouguer anomaly map for Southeastern Cuba island. Gravity anomaly highs are named; 1 Levingstone, 2 La Guira, 3 Piloto, 4 El Salvador, 5 La Perrera, 6 Eje Magmático Sur, 7 New oceanic crust coming from the Cayman dispersion center. Capital letters and lines indicate the six 2D cross-sections made to the 3D density model.
To derive the 3D density model from the complete Bouguer gravity anomaly the software by Gallardo et al. [2005] was used. The top and bottom depths for multiple rectangular prisms were determined using inequality or equality constraints for those depths. We assume that the ground consists of geological units with irregular bottom and top topography in contact with other units. We simulate every unit with a conglomerate of rectangular prisms as shown in Figures 1 and 4. The whole 3D model is constituted of separate geological units or set of prisms with different density contrasts. In Figure 4 we show an example of a 3D model with four geological units and their respective set of prisms. A constant horizontal cross-section area for all the prisms is assumed.
Example of a 3D model with only four geological units. Grids are stratified at the beginning but through the iterations the prisms depths vary. Some prisms’ thickness collapses to zero allowing the outcropping of the lower geological units. Numerically, those zero thickness prisms exist but geologically, they do not exist.
The inversion process moves the top and bottom depth for every single prism at every geological unit. Restrictions are imposed to not allow overlap or spaces between prisms. The quadratic norm of the differences is minimized between data (go) and model response (gr) plus a smoothing term (equation 1).
subject to
Where m is the unknown vector containing the depths from every prism. Matrix D is the horizontal (x, y) first derivatives of the depths. This term minimizes the top depths differences between adjacent prisms. Term β magnifies or dismisses this term. When it is zero the model shows very rough top topography for every unit; when large, every topographic unit looks very smooth, except where the data (first term in equation 1) requires larger jumps. This can happen where geological faults are located.
Depth determination is quoted by means of quadratic programming (Gill et al., 1986), using inequalities or equalities. This allows introduction of surface geology, wells and seismic data as constraints.
Surface geology is introduced as a priori information. Figure 4 shows that unit-1 prisms are displaced to allow units-2 and 3 prisms outcropping. This is performed in the algorithm by collapsing the prisms’ thickness to zero.
For the modeling, the horizontal prism area was fixed in 10kmx10km as shown in Figure 1, giving a set of 646 prisms for every geological unit.
In order to reduce the non-uniqueness, density contrasts are considered as known. Densities were obtained by direct sampling on the surface. Those densities have a variance range due to heterogeneities inside the geological unit. In the inversion process we adjusted the density contrasts along those ranges. The inversion is not completely automatic because we had to try different density contrasts as fine adjustments.
Geologic models allowed us to establish seven geological units (Iturralde-Vinent 1998, 2002; Cobiella, 2005; Sommer et al., 2011). We expanded from the simplest model (Occam’s razor) with only seven units to nine units (Table 1) for the inversion process. We considered only those units that exhibit a density change, plus the gravity response of the sea, and the mantle response which was subtracted when corrected by theoretical ellipsoid (Blakely, 1996). Twenty iterations were performed to arrive at the final model.
Listed are the nine geological units used to obtain the 3D density model. We worked with the shown density contrasts in gr/cm3.
| Units | Density (g/cm³) | Density Contrast | |
|---|---|---|---|
| 1 | Sea | 1.03 | -1.27 |
| 2 | Neogene-Cuaternary deposits | 2.25 | -0.05 |
| 3 | Paleogene Volcanic Arc | 2.29 | -0.01 |
| 4 | Ophiolite 1 (Maffic Mayari -Sagua - Baracoa) | 3.05 | 0.75 |
| 5 | Cretaceous Volcanic Arc (Tunas - Holguín y Sagua- Baracoa) | 2.95 | 0.65 |
| 6 | Ophiolite 2 (North Holguín) | 3.1 | 0.7 |
| 7 | Bahamas Platform | 2.2 | -0.1 |
| 8 | Oceanic crust | 3.15 | 0.85 |
| 9 | Mantle | 3.3 | 0 |
The most recent hypothesis about the Southeastern Cuba states the crustal structure consists of folded basement overriding Bahamas platform. The folded basement is constituted of ophiolites thin sheets intercalated over and under the Cretaceous Volcanic Arc (CVA). At the Southeast of our study area, the Paleogene Volcanic Arc (PVA) rocks are predominant and lay over the ophiolites flakes and CVA package (Iturralde-Vinent 1998; Iturralde-Vinent et al. 2002; Cobiella et al. 2011). These geological models were constructed from petrological data taken from surface rocks. The geological unit at depth is inferred with a large probability error that increases with depth. In contrast, our 3D model obtained from gravity data has a quantitative character. We can estimate depths, thickness and dip angles with some probability error that also increases with depth, but these errors are smaller than the geological ones. Thus the algorithm takes the surface geology and boreholes information as geologists would and then guides the 3D structures at depth obeying the physics behind the gravity data. Our procedure reduces the uncertainly at depth considerably.
Before presenting the 3D model obtained, it is important to show the model response (Figure 5B) and the differences between data and response (Figure 5C). Model response has a 6% data misfit. We reproduced the gravity highs (from 140 to 214mGal), and also the Cauto-Nipe basin gravity low at the NW with values from 0 to 10mGal. The differences map has a minimum and maximum of -4 to 4mGal and shows almost a random behavior. The main differences concentrate at the south border of the island because the 10kmx10km prisms do not fit exactly the steep coastline (Fig 1). There is a surface excess or deficiency of mass producing a misfit. Excluding that border misfit, general misfit must be lower than 6%.
In Figure 6 is shown the 3D density model obtained. We show the bottom topography for every geological unit. Figure 6A shows the surface topography and bathymetry. This surface is known and therefore constrained. Sea effect was taken in consideration. According with the units shown in Table 1, sea bottom topography corresponds with sediments top topography, sediments bottom corresponds with PVA top and so on. Top mantle topography is not shown because it corresponds with the bottom oceanic crust depth. In the bottom depth maps, positive value means above and negative value below sea level are shown. For example: PVA outcrops at Sierra Maestra, therefore showing positive elevation levels at sediments bottom topography (Figure 6A). While these contour plots are not visually informative, they can be digitized and used for future research.
3D density model obtained. Maps represent the bottom topography for every grid or geological unit with respect to the sea level. Figure 6. (A) Terrain topography (positive levels) and bathymetry (negative levels).
For Central and Eastern Cuba (see oval in Figure 7), Otero et al. (1998) argued that there is an oceanic crust transitioning to continental northward. In order to test that hypothesis, we obliged the 3D model to put oceanic crust at the very bottom of the crust (denser than continental). The gravity response was so high that the algorithm never got a good convergence as we can see in Figure 7B. If we put a continental crust instead (Bahamas platform), the fitness is optimal as we can see in Figure 7A. We therefore postulate that Bahamas platform must be there because it is less dense than oceanic crust.
Figure 8A shows the surface geology declared at every surface prism. It also shows the location of four cross-section of the 3D model. Figure 8B show those cross-sections. At cross-section AA’ and BB’, we can see that the Bahamas platform dips inside the mantle with a slope close to 45o. However at the East (cross-section CC’), the same platform arises forming the Mayarí-Moa-Baracoa massif. Cuevas (1998) believes that this massif was raised by isostatic compensation movements. Chang (2003) argues that the Mohorovicic boundary is moving upward. Based on our model, the latter hypothesis is more probable, because we do not see isostatic roots in cross-section DD’.
The Bahamas platform belongs to the North-America continental plate (Figure 1). Cross-section AA’ (Figure 8) shows clearly how this platform penetrates the mantle in a kind of slab with an approximate thickness of 10 to 12km, dipping southward. It begins with a low angle at North, where it almost outcrops, increasing the angle southward to 45o below the Cauto-Nipe basin. Cross-section BB’ is similar but the PVA is less overlapped. Cross-section CC’ shows a rebound or vertical uplift of the Bahamas platform. It seems that the PVA is distributed more at the West of this part of the island. Iturralde-Vinent (1998), Cobiella (2005), García-Casco et al. (2008) and Sommer et al. (2011) have suggested this slab before, but this is the first geophysical evidence. The Bahamas platform density is very close to 2.3g/cm³, meaning that density contrast is almost zero. The shallow presence of Bahamas platform at NW (~2km) justifies the low gravity values.
At cross-section AA’, the Cauto-Nipe fault is signed as a high depth gradient. This fault cuts the northern side of the Sierra Maestra massif, dipping northward.
At the SW corner of the study area, the model needed a high density body (cross-section AA’; Figure 8 B) in order to fit the gravity high (number 6 at Figure 3). Chang (2003) has suggested a pluton below the CVA rocks.
Southward in the Cuba island, the Oriente fault is present (Figure 1). The fault signals the change from continental to oceanic crust. Even southward of cross-section AA’ (Figure 8B), the 3D model cuts this fault. The high gravity anomaly requires a denser body which could correspond to oceanic crust. This new oceanic crust has been shifted from West to East by the Caiman trench. The gravity anomaly decreases southward in cross-sections BB’ and CC’ (Figure 8B). The 3D model justifies this with the less dense Gonave continental micro plate (Calais et al. 2002, 2006).
Cross-section BB’, at the center cuts an ophiolite body of 4km thick. This high density body causes the Levingstone high gravity anomaly (number 1 at Figure 3). Levingstone had been interpreted as caused by an ultrabasic mantle intrusion (Otero et al. 1998). We sought to understand the origin of such gravity anomaly highs (Figure 3). We made two additional cross-sections over the 3D density model. Cross-section EE’ (Figure 8B) crossed Levingstone, La Guira and Piloto gravity highs (number 1, 2 and 3 at Figure 3). Figure 8A shows a great correspondence between the ophiolites bodies and the gravity anomaly highs. This is explained because the density contrast is +0.75gr/cm3 (Table 1) and also because those bodies are shallower producing three high-frequency features over the gravity anomaly. Cross-section CC’ passes over the El Salvador gravity high (number 4 at Figure 3). Despite the smoothing (Figure 8), it is clear that ophiolites at km 80 are producing the gravity highs.
Cross-section FF’ (Figure 9) passes over La Perrera gravity high (number 5 at Figure 3). Figure 9B shows that in this case the ancient ocean crust keeps strong correlation with this anomaly.
The Eastern Island has a more complex geology, including the Mayarí-Moa-Baracoa ophiolitic massif. Iturralde-Vinent (1998), Cobiella (2005) and Sommer et al. (2011) suggest that the evolution of this area is different. Cross-section CC’ (Figure 8B) shows that the Bahamas carbonate basement raised up, forming a horst bellow the Mayarí-Moa-Baracoa massif. Here, crustal thickness is 15km.
The southern corner of cross-section CC’ (Figure 8B) reached the deformed Santiago de Cuba belt, formed by the rising of the ancient ocean crust. Tectonically this was due to a transpressive process produced by the oblique contact between the Caribbean plate and Gonave microplate. It is important to emphasize that here the Gonave microplate is constituted by the CVA and ophiolites over the ancient oceanic crust that have been migrating from the Yucatan basin and is displaced by the Oriente fault to the present position. Similar crust composition is observed at La Española Island (Case et al. 1990) where ophiolites and CVA outcrop.
Current discussion seeks to explain whether the Bahamas platform and Caribbean plate collisioned or subducted each other. Iturrade-Vinent (1998), Cobiella (2005) and Sommer et al. (2011) argue that the Bahamas platform subducted under the Caribbean plate (see cross-section AA’ at Figure 8B). García-Casco et al. (2008) argues for a collision. From our results, we propose a tectonic emplacement. At Manaibon elevation, the ophiolites can be seen over the Bahamas platform. At Manaibon and Cupey Sierra, olistolites appear inside the ophiolites, indicating that the ophiolites formation dragged over the top of the Bahamas platform when the overriding process occurred. Thus a low density Bahamas platform was overridden by a denser Caribbean plate. Since normal subduction phenomena occurs when a denser layer is overridden by a less dense layer, in this case, collision occurred rather than subduction.
When we conceived the conceptual model, we assumed that might be in the bottom, remains of an ancient crust and an oceanic crust younger (CSC) from the area of generation (Cayman Trench). Both were modeled as a single unit because they have approximately the same density, but in the same cross-section AA’ (figure 8) they are shown as two types oceanic crusts.
A limitation of this research is that we used 10kmx10km prisms area, therefore, we cannot resolve high spatial frequencies. Nevertheless, the gravity anomaly is very smooth (Figure 3). There is no high frequency information in such a gravity anomaly. Another limitation is that we could not differentiate the Cuban metamorphic complex from the Cretaceous volcanic arc because densities are similar, and also between sedimentary volcanogenic sequences (Paleogene age).
ConclusionsThree-dimensional gravity inversion can be done using the same size grid of prisms and assuming densities as unknowns implying a lineal process with a single iteration solution, but the non-uniqueness is quite high by the huge number of prisms. In this research, our inversion procedure considers depths as unknown and densities as known, bringing a non-lineal problem that must be linearized by an iterative process. In this inversion, the non-uniqueness diminishes considerably because the number of prisms is much smaller. However, we also applied geological and geophysical constraints to diminish such non-uniqueness even further. We used the Occam’s razor criteria that warrants a most probable 3D model.
The uncertainty inside the model is not constant. It increases with depth but also at the 3D model edges because of the lack of data. High frequency details may also have high uncertainty because the gravity anomaly is very smooth. Consequently, we looked for the simplest but the most probable 3D model. For scale purposes, 10kmx10km area prisms were optimal.
The smoothing parameter or roughness penalization avoids large depth steps between contiguous prisms. However, if the resolution of gravity data is optimal, the fitness term becomes more important than the smoothing term. In our 3D density model, large depth steps appear when data require them, because we are using an intermediate value for the smoothing parameter.
Parallel 3D topographies mean linear correlations between the unknowns. Problem increases when increasing the prism number and also with deeper formations. However, our model does not exhibit high linear correlation at the deepest part, because we did not use a large number of prisms.
The gravity response from the 3D density model obtained reproduces very well the general shape of the data, particularly the gravity lows at Cauto basin (NW of the map), with values between 0 and 10mGal. It also fits the named gravity maximums. We determined that those maximums are caused by the presence of shallow ophiolite sheets mainly.
Crustal thickness obtained from the 3D model is 10 to 15km southward and thicker northward. This result is agree with previous qualitative geological models.
The density model further shows the complex 3D topographies of the Bahamas platform, the Cretaceous volcanic arc, the Paleogene volcanic arc, and the ophiolites sheets, and how they intrude each other.
The 3D density model shows the new oceanic crust from the Gonave microplate (SE of the area) caused by the pull of the Cayman spreading center.
The response of the ophiolites sheets and Cretaceous volcanic arc overriding the Bahamas platform provides evidence for collision rather than subduction between the ancient Caribbean plate and the Bahamas platform (present North-American plate).
We thank the Institute of Geology and Paleontology of Cuba and especially Eric Escobar and Fernando Mondelo for giving us the data of the gravity anomaly used in this study.