Skip to main content

Science – Society – Technology

Curie point depth, thermal gradient, and heat flow in the Colombian Caribbean (northwestern South America)

A Publisher Correction to this article was published on 03 August 2019

This article has been updated

Abstract

Geothermal gradients were estimated at the points of a grid with cell size of 50 km by 50 km in northwestern Colombia. Depths to the bottom of magnetic sources were assumed to represent Curie depths, and were estimated by means of statistical-spectral analysis of aeromagnetic data contained in square windows of 100 km by 100 km, 200 km by 200 km, and 300 km by 300 km. The centroid method and a variation of it that considers the fractal distribution of magnetization were applied. The modified centroid method provided better estimations of Curie depths, which range between 13 km and 47 km. Obtained Curie depths are comparable to those reported in other studies of regional character. The corresponding estimated geothermal gradients correlate quite well with estimations based on the bottom simulating reflector and bottom hole temperatures. Observed differences are small and can be accounted for by local heat flow variations due to shallow ground water flow, or recent sedimentation. Elaborated geothermal gradient and heat flow maps for northwestern Colombia are accurate and consistent with estimates of the thickness of the oceanic crust and continental cortical thicknesses, reported by previous authors. The maps presented in this study represent a contribution to the heat flow studies in northwestern South America.

Introduction

Contributions to the heat flow mapping of Colombia were reported in regional studies (IHFC 2008; Cardoso et al. 2010; Li et al. 2013; Davies 2013; Salazar et al. 2017). Previous studies were of a lower resolution and depicted long wavelength heat flow anomalies (IHFC 2008; Cardoso et al. 2010). Heat flow studies were reported for the Colombian-Caribbean oceanic zone with no detailed work on Colombian territory (Salazar et al. 2017; Li et al. 2013).

However, a higher resolution heat flow map of the area comprising the Colombian-Caribbean oceanic domain, and northern continental Colombia (Fig. 1) would be a contribution to the study of the interaction of the Caribbean, and South American tectonic plates.

Fig. 1
figure 1

Present configuration of the Caribbean plate (modified from Giunta and Orioli 2011). CHC: Choco; CHR: Chorotega; MAY: Maya; CHT: Chortis; SCDB: South Caribbean Deformed Belt; SB: Sinu Belt; SJFB: San Jacinto Folded Belt; COBA: Eastern Mountain Range of the Andean Block. The red and green dots correspond respectively to the Plato Basin, and to the San Jorge basin. The red inlay indicates the study zone

In particular, knowledge of the depth to the Curie isotherm (assumed to coincide with the depth to the bottom of magnetic sources, DBMS) will help to characterize the thermal structure of the Colombian-Caribbean lithosphere, and in particular its rheological strength (brittle, or fluent crust). Furthermore, a more detailed heat flow mapping of the above-mentioned region will help to evaluate the potential of geothermal resources, as well as those of hydrocarbures, and to better understand the geologic and tectonic evolution of this region.

The study area comprises the north-western part of South America (Fig. 1). The area covers parts of the basins of Colombia and Venezuela as well as areas of Panama, and continental Colombian and Caribbean domains (Fig. 1). The area includes the convergence zone of the Caribbean, South America, Cocos, and Nazca tectonic plates (Granja Bruña 2005). According to Salazar et al. (2017), the Caribbean plate subducts under the South American plate, while the Cocos plate subducts beneath the Caribbean plate.

The oceanic zone of the Caribbean Plate is approximately 10- to 15-km thick (Diebold et al. 1981; Mann 1999). The average global oceanic crustal thickness is 6 km (Mann 1999). An oceanic plateau is present in its interior, which is 12- to 15-km thick, of a transitional character, and with an approximately 2-km thick sedimentary cover. The continental zone presents thicknesses between 20 and 60 km (Salazar et al. 2017).

This study reports, for the described study zone, thermal gradient estimated from (1) surface temperatures, (2) Curie point depths (CDP), and (3) a Curie temperature of 580° corresponding to magnetite. The applied methods are described. Results are presented and analyzed.

Background

The geothermal gradient is a measure of how the temperature varies with depth. In a 1D Earth the geothermal gradient (Nondorf 2016; Sigismondi and Ramos 2008) and heat flow are related by the Fourier’s law (Abraham et al. 2015; Kasidi and Nur 2013):

$$q = \lambda \;{\text{d}}t/{\text{d}}z$$
(1)

where \(q\) is heat flow, \(\lambda\) is thermal conductivity, and dt/dz is the geothermal gradient. Temperature measurements at different depths along a borehole constitute the most common way to establish the geothermal gradients. At sea, the bottom simulator reflector (BSR) method enables one to estimate geothermal gradients at sea environments.

This method, proposed by Yamano et al. (1982), has been amply applied worldwide. The BSR is defined as the anomalous reflector that coincides with the stable zone of gas hydrate (López and Ojeda 2006; Ganguly et al. 2000; Shankar and Sain 2009). According to Liao et al. (2014), in low latitudes, methane gas hydrate is present at depths between 700 and 1000 m.

Depth to the Curie isotherm constitutes another tool to establish the geothermal gradient and the heat flow coming out of the mantle. According to Tanaka et al. (1999), Curie temperature (\(\theta_{\text{c}}\)) can be defined as:

$$\theta = \left[ {{\text{d}}t/{\text{d}}z} \right]z,$$
(2)

From Eqs. 1 and 2 heat flow (\(q\)) can be estimated in terms of the Curie temperature (\(\theta\)) and depth to the Curie Isotherm, (Zb):

$$q = \lambda \left( {{\raise0.7ex\hbox{$\theta $} \!\mathord{\left/ {\vphantom {\theta {Z_{\text{b}} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${Z_{\text{b}} }$}}} \right)$$
(3)

The Curie temperature of magnetite, 580 °C, is assumed for this study as the temperature at which crustal rocks lose their magnetic properties. In several studies (Abraham et al. 2015; Anakwuba and Chinwuko 2015; Kasidi and Nur 2013) this temperature was assumed. The Curie point depth (CPD) corresponds to the depth at which the rocks reach such a temperature (Tissot and Welte 1984). To estimate the CPD several spectral statistical magnetic techniques have been developed to determine the DMBS, which are assumed to represent the Curie isotherm. Among them, the centroid and fractal methods that have been applied in Australia, the United States, China, Japan, and Taiwan (Bansal et al. 2011; Bouligand et al. 2009; Espinosa-Cardeña et al. 2016). Following a description is given of the BSR method and of the spectral statistical methods used to estimate the depth to the Curie isotherm in this study.

Methods

Spectral statistical methods

Statistical-spectral analysis of aeromagnetic data is based on the work of Spector and Grant (1970), Bhattacharyya and Leu (1975), Okubo et al. (1985). The developed methods enable one to determine the DBMS, which can be associated to the Curie point depth or CPD (Bansal et al. 2011).

The centroid method of Tanaka et al. (1999) is the most currently used method. It has been applied in Turkey (Dolmaz et al. 2005), Nigeria (Kasidi and Nur 2012), South Africa (Nyabeze and Gwavava 2016) and East and South East Asia (Tanaka et al. 1999). The Fractal method has been applied in South Africa (Nyabeze and Gwavava 2018), Central India (Bansal et al. 2013), Western USA (Bouligand et al. 2009), Germany (Bansal et al. 2011), Nigeria (Abraham et al. 2015) and Mexico (Espinosa-Cardeña et al. 2016).

Centroid method

The centroid method is based on mechanical statistical principles (Spector and Grant 1970), and on the assumption that the magnetic anomalies constitute an uncorrelated random distribution of magnetic sources (Abraham et al. 2015; Bansal et al. 2010, 2011). The power spectrum of the reduced to the pole, magnetic anomalies is defined as (Blakely 1996; Spector and Grant 1970):

$$\left\langle {E_{\text{RP}} \left( k \right)} \right\rangle = \left\{ {\left| {F\left( {\Delta T\left( k \right)_{\text{RP}} } \right)} \right|^{2} } \right\} = {\text{Cm}}^{2} 4\pi^{2} M^{2} \left\langle {e^{ - 2hk} } \right\rangle \langle \left( {1 - e^{ - tk} } \right)^{2}\rangle \left\langle {S^{2} \left( {k,a,b} \right)} \right\rangle$$
(4)

where h is the mean depth to the top of the magnetic source ensemble, t is the thickness of the magnetic ensemble, M is the magnetic moment, \(S^{2} \left( {k,a,b} \right)\) is a function of the magnetic ensemble mean horizontal dimensions (a, b), \(k\) corresponds to the radial wave number in cycles/km, and \({\text{Cm}}\) = 10−7 is the proportionality constant to transform Eq. (4) to the SI unit system (Blakely 1996). Zhou and Thybo (1998) determined that the radial power spectrum with or without reduction to the pole are virtually identical. Then, the radial power spectrum \(P\left( k \right)\) can be rewritten as:

$$P\left( k \right) = A_{1} e^{{ - 2\left| k \right|Z_{k} }} \left( {1 - e^{{ - \left( {Z_{\text{b}} - Z_{\text{t}} } \right)}} } \right)^{2}$$
(5)

where \(A_{1}\) is a constant, \(Z_{\text{b}}\) and \(Z_{\text{t}}\) are respectively the depths to the bottom and to the top of the magnetic body ensemble. By simplifying Eq. 5, the centroid depth of the magnetic source can be calculated from the large wavelengths spectrum zone (Bhattacharyya and Leu 1975, 1977; Okubo et al. 1985) by the following expression:

$${ \ln }\left( {\frac{{P\left( k \right)^{1/2} }}{k}} \right) = A_{2} - \left| k \right|Z_{\text{o}}$$
(6)

where \(A_{2}\) is a constant, \({ \ln }\) is the natural logarithm, and \(Z_{\text{o}}\) is the centroid depth. If it is assumed that the top signals of the magnetic ensemble dominate the high-medium frequency portion of the spectrum, \(Z_{t}\) can be obtained similarly from Eq. 5 (Bhattacharyya and Leu 1975; Okubo et al. 1985; Spector and Grant 1970) by the following equation:

$${ \ln }\left( {P\left( k \right)^{1/2} } \right) = A_{3} - \left| k \right|Z_{\text{t}}$$
(7)

where \(A_{3}\) is a constant, and \(Z_{\text{t}}\) is the magnetic source top depth.

Fractal method

Gravity and magnetic fields follow a fractal behavior (Dimri and Ganguli 2019). Magnetic susceptibility (Fedi 2003) and crust magnetization, in particular, have fractal behaviors (Pilkington et al. 1994). Bouligand et al. (2009) used this magnetization model to estimate depths to the bottom of magnetic sources in southwestern USA.

The radial average power spectrum for a fractal distribution of the magnetic anomaly sources follows the relation:

$$E_{\text{c}} \left( k \right) = \frac{{\ln E_{\text{RP}} \left( k \right)}}{{k^{ - \alpha } }}$$
(8)

where \(E_{\text{c}} \left( k \right)\) is the corrected spectrum, \(E_{\text{RP}} \left( k \right)\) is the spectrum reduced to the pole, \(k^{ - \propto }\) is the correction factor (Fedi et al. 1997), and \(\propto\) is the decay factor (Pilkington et al. 1994). Hence, the centroid depth can be obtained from the spectrum corrected for its fractal behavior by the following equation:

$${ \ln }\left( {k^{\alpha } * \frac{P\left( k \right)}{{k^{2} }}} \right) = A_{4} - 2\left| k \right|Z_{\text{o}}$$
(9)

Similarly to the centroid method, one can calculate the depth to the mean top of the magnetic ensemble by combining Eqs. 7 and 8:

$${ \ln }\left( {k^{\alpha } * P\left( k \right)} \right) = A_{5} - 2\left| k \right|Z_{\text{t}}$$
(10)

Fractal scaling factors between 2 and 4 have been used by Pilkington et al. (1994), Fedi et al. (1997), Maus et al. (1997), Bouligand et al. (2009), Salem et al. (2014) and Nyabeze and Gwavava (2018). Values of 3 and 1.5 over correct the spectrum according to Ravat et al. (2007). According to Bansal et al. (2011) a value of \(\alpha\) = 2 works well, but Bansal et al. (2013) consider that the value that best corrects the power spectrum is \(\alpha\) = 1.

In both methods,\(Z_{\text{b}}\) results from a simple relationship (Bhattacharyya and Leu 1975, 1977; Okubo et al. 1985) between (\(Z_{\text{t}}\)) and (\(Z_{\text{o}}\)):

$$Z_{\text{b}} = 2Z_{\text{o}} - Z_{\text{t}}$$
(11)

The BSR method

On marine seismic lines, the BSR results from the large acoustic impedance contrast between the hydrate gas layer above the BSR (high seismic velocity zone) and a underlying layer of free gas that constitutes a low-seismic velocity zone (Dong et al. 2018). To calculate the geothermal gradient, the following steps are followed: (1) identification on the seismic lines of the reflectors corresponding to the sea floor and to the BSR, (2) depth conversion of the two selected seismic horizons, (3) conversion of depth to pressure, (4) use of phase diagram to obtain the temperature corresponding to the previous estimated pressure, and (5) calculation of the geothermal gradient that implies dividing the temperature difference by the corresponding depth difference (Minshull 2011).

In practice, the geothermal gradient is calculated assuming a sea bottom temperature of 4 °C (Vohat et al. 2003), and making use of the equation developed by Yamano et al. (1982) and reviewed by Grevemeyer and Villinger (2001):

$$q = k \Delta T$$
(12)

where

$$\Delta T = \frac{{\left( {T_{\text{BSR}} - T_{\text{sea}} } \right)}}{{\left( {Z_{\text{BSR}} - Z_{\text{sea}} } \right)}}$$
(13)

\(q\) is heat flow, \(k\) is thermal conductivity, \(\Delta T\) is geothermal gradient, \(T_{\text{BSR}}\) is BSR temperature, \(T_{\text{sea}}\) is sea bottom temperature, \(Z_{\text{BSR}}\) is BSR depth, and \(Z_{\text{sea}}\) is sea bottom depth.

In this study, 250 time-migrated seismic lines, located mainly in the continental slope, were examined; of these, in only 20, the BSR could be identified. The sea bottom and the BSR were converted to depth using velocities of 1450 m/s for the sea water, and 1850 m/s for the sediments (average velocity used in other converging zones, as for example Brown et al. 1996).

For this work, the pressure to temperature conversion, was based on the phase-stability diagram of Sloan (1998) because the obtained geothermal gradients are comparable to those reported by López and Ojeda (2006) for some locations of the Colombian-Caribbean offshore zone.

Aeromagnetic data and processing

The magnetic data used in this study were taken from the World Digital Magnetic Anomaly Map V2 (WDMAM, http://www.wdmam.org/download.php), approved by the International Association of Geomagnetism and Aeronomy (IAGA) in 2015 (Lesur et al. 2016). In the first version of the map (Korhonen et al. 2007), for the marine part of the EMAG2 grid (Earth Magnetic Anomaly Grid), for those areas without data (Maus et al. 2009), extrapolations were made based on the sea floor age map (Müller et al. 1997) and on a sea bottom expansion model (Dyment et al. 2015). Version 2 (Li et al. 2013) involves tectonic plate movements, and voids in continental information were covered with synthetic data arising from the lithospheric field model GRIMM_L120 (Lesur et al. 2016). The internal magnetic field was eliminated through an extensive model (Sabaka et al. 2004), that includes spherical harmonics up to the 13th degree to avoid magnetic contributions from the core (Ravat et al. 2007). In this way, wavelengths larger than 500 km were excluded, in particular, from the Colombian-Caribbean magnetic data (Fig. 2). This assures that no magnetic sources that may not correspond to CPD are present in the data. Similar procedure was applied to the data used by Bouligand et al. (2009) to determine the CPD for southwestern USA. Models elaborated in a similar way have provided good results in determining Curie isotherm (Bouligand et al. 2009; Manea and Manea 2011; Espinosa-Cardeña et al. 2015).

Fig. 2
figure 2

Total field magnetic anomalies of the study area. Database taken from the World Magnetic Anomaly Map V2 (Lesur et al. 2016). Short gray lines indicate location of seismic lines used to calculate the geothermal gradient by means of the Bottom Simulator Reflector (BSR) method. Symbols indicate location of used Bottomhole Temperature (BHT) and the attached letters to them refers to respective differences between the calculated gradient by the fractal method and the gradient estimated from BHT, presented in Fig. 9. The upper left inset shows the location of the study zone within South America

The optimal window size to estimate Zb should be carefully selected. The maximum CPD is limited by 2π/L, L being the window size length (Campos-Enriquez et al. 1990; Campos-Enríquez et al. 2019). Window sizes between 4 and 10 times the depth to the source have been used by Ross et al. (2006), Bouligand et al. (2009), Chopping and Kennett (2015), Vargas et al. (2015). However, windows of more than 200 km could include contributions of different tectonic and geological environments (Ravat et al. 2007) and because of this Bansal et al. (2011), (2013) and Saibi et al. (2015) recommend windows of 200 × 200 km. Windows of 100 × 100 km provide good resolution but shallow Zb as reported in several studies (Ravat et al. 2007), windows of 300 × 300 km provide lower resolution and emphasizes long wavelength trends (Li et al. 2013). Windows of 200 × 200 km enabled to obtain intermediate resolution, avoiding regional tendencies.

For the (Fig. 2) the square windows of 100 × 100 km, 200 × 200 km, and 300 × 300 km were used for Colombian-Caribbean studies. The respective radial anomaly power spectrums were obtained from the magnetic anomalies contained in such windows. Then the corresponding DBMS were estimated. The centers of these windows form a grid of 50 km in east–west direction, and 50 km in the north–south direction. Geosoft Oasis Montaje was used to obtain the average radial power spectrum of each window. Regional trends were removed, the grids were expanded 10% by the maximum entropy method to eliminate edge effects. FFT was used to obtain the power spectra.

The CPD depths were estimated following the centroid method (Eqs. 6 and 7), and the fractal method (Eqs. 9 and 10) with a fractal parameter \(\alpha\) = 1. With this value, the different spectra zones were well defined, in particular the low wavelength part of the spectrum that provides information on the bottom depth, Zb. At posteriori, results confirmed it was a good choice.

Selection of the two spectrum domains to estimate the top or bottom depends on the interpreter judgment based on study area geology (Bansal et al. 2011). In this work, three slopes (maximum, minimum, and intermediate) were selected, in each domain, and used to estimate the respective CPDs. The respective average CPD and the corresponding difference (in %) were calculated. Generally, the obtained spectra can be classified in two categories. In one group, spectra present clearly defined zones or domains, in the other group, the long wavelength spectra domain is not clearly defined; when analyzing these latter spectra, various estimations were made until a realistic estimation aroused (Bansal et al. 2011). It was assured that the estimated depths were consistent with bathymetry (oceanic zone), topography (continental zone), and crustal thickness, and available geological and geophysical information.

CPD estimations obtained from each method (centroid and fractal) were checked for outliers along W-E and N-S profiles. If outlier values were detected then the corresponding spectra was re-analyzed and depths re-estimated. In this way all outliers were eliminated. Subsequently, values were interpolated by means of a minimum curvature method, and again the resulting values were checked for outliers.

Results and discussion

Figure 3 shows estimated CPD (Zb) examples for ocean and continental zones (respectively parts 1 and 2 of the figure) obtained using the centroid and the fractal methods. Points and lines in green correspond to Zt, those in red are related to Zo. The spectra are presented just up to a wave number of 0.1 because it is the interest zone where the CPD is defined (De Ritis et al. 2013; Salem et al. 2014). Although, as these examples show, Zb estimations by both methods are consistent, centroid method provided underestimated depths. Similar result has been reported by De Ritis et al. (2013), Hussein et al. (2013), and Li et al. (2013).

Fig. 3
figure 3

Examples of spectra and respective CPD estimations for oceanic (Part 1) and continental (Part 2) zones. Only the long wavelength domain of the spectrum is presented. Points and lines in green, respectively, indicate the spectrum zone used to calculate the depth to the top (Zt) and the corresponding adjusted straight line. Points and lines in red, respectively, show the zone used to determine (Zo) and respective adjusted line. Equations of the adjusted lines and their respective correlation coefficients (R2) are given. At the left side of Part 1, representative CPD estimates (estimated values common for that area) shown for an oceanic area (estimations by both method) and for a maximum value (right side). Part 2 presents, for the continental area, CPD estimations, by both methods, for a minimal (left side) and for a maximal value (right side)

Curie point depth map

The Colombian-Caribbean Curie Point Depth maps (Fig. 4) show that the Curie isotherm lies between 18 km and 47 km in the continent. Depths in the oceanic area were between 13 km and 27 km. These estimations are within the range of reported crustal thicknesses, between 20 and 60 km for the Colombian Caribbean and northwestern Venezuela as reported by Salazar et al. (2017). The depths were consistent with those reported in the Global Reference Model of Curie Point Depths (Li et al. 2017). The Curie point isotherm is located in the lower crust.

Fig. 4
figure 4

Curie point depth maps. a Based on the Fractal method. b Based on the Centroid method. White line represents the country limit. Maximum and minimum CPD are 47 and 13 km, respectively

Differences in percentage were calculated at each grid point for the CPD maps obtained by the fractal method (Fig. 5a), and by the Centroid method (Fig. 5b). The largest differences correspond to the centroid method (between 20 and 30%); while for the fractal method the dominant differences are in a range between 4 and 8%. Because the values of differences in percentage by the fractal method are lower, they were considered more reliable, and proceed to elaborate the Curie point depth map based on these.

Fig. 5
figure 5

Differences in percentage between Curie Point Depth estimations by the fractal and centroid methods. a Differences corresponding to the fractal method range between 4 and 6%. b Distribution corresponding to the centroid method. It presents several maximums with differences between 20 and 30%

Geothermal gradient map

The geothermal gradient anomaly map (Fig. 6) indicates that geothermal gradients in the offshore zone of Colombian-Caribbean region vary between 44 and 20 °C/km, while in the continent, values range between 24 and 14 °C/km; an atypical maximum of 39.5 °C/km stands out in the area amid Barranquilla, Sabanagrande and Cartagena (indicated in Fig. 6 by 1, 2, and 3, respectively). Intermediate gradients are observed associated with the offshore Guajira basin, in the central part of the Colombia basin, in the Magdalena river delta, and outstanding low gradients in the Plato sub-basin (in the Lower Magdalena Valley) and in the Uraba basin (indicated in Fig. 6 by numerals 6, 7, 8, 9, and 10, respectively). For the Colombian Guajira and the Gulf of Maracaibo in Venezuela (4 and 5 respectively in Fig. 6), Arnaiz-Rodríguez and Orihuela (2013) reported larger CPD values that are consistent with these low gradients. The map of this study (Fig. 6) is featured by similar values and patterns as the geothermal gradient map of Salazar et al. (2017). As has been discussed before, the geothermal gradient map reliability for the Colombian-Caribbean area is supported by the close agreement, in the offshore zone, with geothermal gradients estimated by the BSR method and, in the continental zone with 1140 records of Bottom Hole Temperature (BHT). Thus, CPD-based gradient estimations here obtained are a reasonable and accurate estimation of the geothermal gradients in the study area.

Fig. 6
figure 6

Colombian-Caribbean Geothermal Gradient Map. White line indicates the littoral and border lines of Colombia Venezuela and Panama. The extreme values of the color bar result from the interpolation, 1: Barranquilla; 2: Sabanagrande; 3: Cartagena; 4: Guajira; 5: Gulf of Maracaibo; 6: Offshore Guajira Basin; 7: Colombia Basin; 8: Magdalena Delta; 9: Plato Sub-Basin (in the Lower Magdalena Valley); 10: Uraba basin

Estimation of geothermal gradient by the BRS method were based on the phase-stability diagram of Sloan (1998) of Fig. 7. In Fig. 8 are compared geothermal gradients calculated from Curie point depths obtained by the centroid method corrected for the fractal behavior (red lines) with geothermal gradients obtained by following the BSR method (blue lines). Their mean difference is annotated, the broken black line represents the geothermal CPD-based gradient corrected for this difference; the corresponding seismic line number is indicated in parentheses (see its location in Fig. 2). The mean differences vary between − 1.5 and 10 °C/km. In two cases the difference is zero.

Fig. 7
figure 7

Phase curves of Brown et al. (1996), Dickens and Quinby-Hunt (1994), Lu and Sultan (2008), Miles (1995), Østergaard et al. (2000), Sloan (1998), Sloan and Koh (2007). These phase curves are the result of P, T point pairs obtained in laboratory and their respective adjustments (Gómez et al. 2016)

Fig. 8
figure 8

Geothermal gradient comparison for the offshore part of the Colombian-Caribbean domain. Red lines represent the geothermal gradient calculated by the centroid-fractal method, the blue lines correspond to the BSR-based geothermal gradients. The broken black lines represent centroid-fractal method-based gradients corrected by the mean difference between both types of data. Figures next to right of the lines are the difference between the two types of gradients. Where no broken black line is present, the difference between the gradients is equal to zero. Numbers in parentheses, at lower left corners, identify seismic line (see location in Fig. 2)

Vargas et al. (2009) calculated geothermal gradients for Colombia based on bottom hole temperature data set. The objective of the map was to assess geothermal potential zones in continental Colombia. This data set was used here to assess the reliability of the CPD-based geothermal gradients. Figure 9 shows the geothermal gradients estimations derived from depths obtained by the centroid-fractal method (red lines) and those calculated from BHT values (blue lines). Differences range from 8 to − 10 °C/km. Broken black lines represent gradients calculated by the centroid-fractal method, but corrected for these differences. The magnitude difference varies at the different zones (zones a to h indicated in Fig. 2). In zones b, c, and f, both types of gradients, present no statistically meaningful differences, which means that there is a good relationship between the gradients in each of these zones (b, c, and f in Fig. 9); in zones a, d, and e, the BHT gradient is larger than the fractal method based gradient (8, 5, and 6 °C/km, respectively), which could be explained by the thick sedimentary cover in these zones; in zones g and h, the BHT gradient is about 10 °C/km less than the estimated gradient based on the fractal method, a situation that could correspond to a shallower crust in the coast line, and also possibly to effects of shallow water circulation. However, in general, it is noteworthy that these gradient differences are within the residual error ranging between 5 and 15 °C commonly reported in detailed studies based on the BHT at depths close to 3000 m (Sigismondi and Ramos 2008).

Fig. 9
figure 9

Comparison of gradients estimated from Bottom Hole Temperature (BHT) data and CPD-based gradients. ah in parentheses in lower left corners of the different panels correspond to zones depicted in Fig. 2. The blue line represents the BHT data based gradient estimations; the red line represents the CPD-based gradients. i, j represent the depth BHT gradient. Black lines correspond to linear adjustment of the gradient points

Note that the gradient can increase or decrease with respect to depth (respectively Figs. 9i, j), which could be explained by various mechanisms, such as water circulation. The above comparison analysis indicates that the CPD-based geothermal gradients are consistent within acceptable errors with gradients obtained by the BSR method (in the sea zones) and from BTH (in the continental areas). The very good correlation lends reliability to the Curie depth estimated, and derived geothermal gradient and heat flow estimations.

Heat flow map

The heat flow values used to elaborate the map of Fig. 10 were obtained using the Curie point depths of Fig. 4, and a mean average thermal conductivity of 3.0 W/m °C. According to this map (Fig. 10), the Colombian-Caribbean area is characterized by higher values in the offshore zone than on the continent. In the sea zone, the values vary between 80 and 100 mW/m2; with an atypical relative low of 60 to 80 mW/m2 located off the center of the Colombian Basin. In the transition zone between the oceanic and continental zones, there is a conspicuous change from 100 to 60 mW/m2, the last value corresponds to the average for the continental Colombian–Caribbean domain.

Fig. 10
figure 10

Heat Flow Map for the Colombian Caribbean. The cold tones correspond to low heat flow (onshore zone) and the warm tones correspond to high heat flow (offshore zone). Tectonic frame as that of Fig. 1

The International Heat Flow Commission (IHFC 2008) heat flow map, characterizes the Panama region with values from 85 to 120 mW/m. A very similar pattern is observed in the heat flow map presented here, however, with a smaller range of heat flow values (70 to 100 mW/m2). According to the IHFC map, the entire Colombian Caribbean is featured by values varying smoothly around 85 mW/m2, a value close to the onshore and offshore areas average. The heat flow map for the South American sub-continent (Cardoso et al. 2010) is characterized by values ranging from 60 to 80 mW/m2. Both values and regional trends of this continental map can be also observed in the continental part of the map elaborated in this study. The heat flow map proposed in the framework of North Atlantic thermal evolution study (Li et al. 2013) reports for the Colombian Caribbean offshore zone, values between 15 and 60 mW/m2, which are consistent with the results from this study. Davis’s heat flow map (Davies 2013) proposes higher heat flows in the continental crust with respect to the oceanic crust, which contrasts with most of the published maps. This analysis of available heat flow maps indicates that the elaborated map is a valuable contribution to the heat flow studies in northeastern South America by adding reliable estimations and more resolution to the mapping of the heat flow anomalies in the continent as well as in the oceanic part of the study area.

Magnetic anomalies in the Colombian-Caribbean region present two textures. In the oceanic zone, the high frequency character of the anomalies convey a rugged texture that contrasts with the smoother magnetic anomalies (longer wavelength anomalies) of the continent. These two contrasting zones find a good correspondence in the Curie isotherm map. The shallower CPD are located in the oceanic domain. The continent is characterized by larger CPD. The geothermal gradient map is observed to hold an inverse relationship between depth to the Curie isotherm (CPD) and the geothermal gradient magnitude. Geothermal gradients in the Colombian oceanic zone are between 32 °C/km and 43 °C/km. However, low values of 27 °C/km in the Colombia basin and the Magdalena Delta have been found. These values might be due to a thickening of the oceanic crust, or to large sedimentary cover in those areas. In the ocean zone, a high heat flow is observed at the Beata Ridge that separates the basins of Colombia and Venezuela.

The basins of Plato and San Jorge in the lower Magdalena Valley, and in Maracaibo (Venezuela) are featured by geothermal gradients of 16 and 19 °C/km, respectively. Lower values probably associated with the large sedimentary infilling of these basins. Patterns on the heat flow map correlate very well with major tectonic features of the Colombian-Caribbean region. The folded and deformed Belt of Southern Caribbean (SCDB), possibly constituting the contact zone between the Caribbean and South American plates, is characterized by large heat flow values around 100 mW/m2 and by a geothermal gradient that decreases landward. This gradient attains 70 mW/m2 on the litoral of Colombia and Venezuela and could be interpreted as due to the Caribbean plate subduction beneath the South American plate.

On the continent, the lowest heat flow values were observed between the Santa Marta-Bucaramanga Fault and the Oca Fault. Heat flow could be due to the Caribbean plate subduction beneath the South American plate. Low heat flow values, 40 mW/m2, are associated with the Venezuelan Andes.

Conclusions

The centroid method provided CPD estimations with larger percentage differences than its fractal modification. Centroid-fractal based CPD were considered more reliable estimations. The CPDs for the Colombian-Caribbean offshore zone range between 13 and 29 km. The maximum CPD is located off the Colombia basin center with a depth of approximately 23 km. This depth could be associated to a large sedimentary cover. In continental Colombia CPD depths vary between 20 and 44 km. The deepest CPD is located in continental Venezuelan, at a depth of 47 km. In the oceanic domain, CPDs are located at the bottom of the oceanic crust or immediately below. On the continent, the depths are consistent with reported thicknesses of the complex crust beneath Colombia. The respective and more detailed Curie point depth map reproduces with more detail trends already observed in previously reported Curie depth maps for continental Colombia. These depths to the Curie isotherm were used to estimate the geothermal gradient map presented here.

Comparison analysis with geothermal gradients obtained by the BSR method for oceanic zones, and from bottom hole temperature data for continental areas confirms the reliability of CPD-based geothermal gradients. Differences between geothermal gradients estimated by the BSR method and by BHT measurements on one hand, and gradients estimated from the fractal method based CPDs on the other hand, were between 4 and 10 °C/km in Caribbean Sea areas, and between 2 and 10 °C/km on land. These low differences indicate that gradients are reliable, and within error limits world zones of less than 20%. These differences might be due to local effects of heat flow and recent sedimentation.

There was a good correlation between the obtained geothermal gradient and known values for the Colombian-Caribbean zone. Application of the centroid-fractal method in the Caribbean-Colombia domain, as well as in other areas of Colombia might contribute to the study of the thermal and geologic structure of the crust and upper mantle, as well as helping to assess areas with geothermal and oil potentials.

Available data

The authors used published data and authorized data by the Colombian Geological Service.

Change history

  • 03 August 2019

    An error occurred during the publication of a number of articles in Geothermal Energy. Several articles were published in volume 7 with a duplicate citation number.

Abbreviations

DBMS:

depth to the bottom of magnetic sources

BHT:

bottom hole temperature

CPD:

Curie point depth

BSR:

bottom-simulating reflector

WDMAM:

World Digital Magnetic Anomaly Map

IAGA:

International Association of Geomagnetism and Aeronomy

EMAG:

Earth Magnetic Anomaly Grid

References

Download references

Acknowledgements

The authors express their gratitude to Dr. Alberto Ochoa Yarza, Technical Director of Basic Geosciences of the Colombian Geological Service, for supporting geoscience research and, in particular, authorizing the use of geoscientific data. The support from the “Colciencias” doctoral scholarship program is acknowledged.

Funding

Not applicable.

Author information

Authors and Affiliations

Authors

Contributions

Data processing and writing of the manuscript (WQ). Consulting, supervision, manuscript correction (OC). Discussion and review of data (OH). All authors read and approved the final manuscript.

Corresponding author

Correspondence to Wilson Quintero.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Quintero, W., Campos-Enríquez, O. & Hernández, O. Curie point depth, thermal gradient, and heat flow in the Colombian Caribbean (northwestern South America). Geotherm Energy 7, 16 (2019). https://doi.org/10.1186/s40517-019-0132-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40517-019-0132-9

Keywords