Introducción
La historia termal de provincias ígneo-metamórficas tiene la particularidad de involucrar una serie de procesos que resultan de la relación magma-roca encajante (Ehlers, 2005; Fu et al., 2010; Murray et al., 2018; González et al., 2022; Biralvand et al., 2023; Abbey et al., 2024), en adición a eventos mayores vinculados con la interacción entre las placas tectónicas y/o entre el clima y el relieve (Reiners y Brandon, 2006; Prenzel et al., 2013; Restrepo-Moreno et al., 2019; Noriega-Londoño et al., 2020; Stalder et al., 2020; García-Delgado et al., 2021; Zapata et al., 2021). Los primeros incluyen, específicamente: 1) el emplazamiento de masas plutónicas; 2) el enfriamiento de estas masas por conducción; y 3) el enfriamiento del cuerpo intrusivo y el calentamiento de la roca encajante por conducción y por circulación hidrotermal. Para los segundos se tiene: 1) el enfriamiento por exhumación inducida por la actividad de estructuras en configuraciones tanto compresivas como extensionales; 2) el enfriamiento por erosión del sustrato suprayacente; 3) el calentamiento por soterramiento; y 4) el calentamiento por el desarrollo posterior de arcos volcánicos.
Lo anterior muestra lo complejo que puede llegar a ser dilucidar estas historias termales (Harrison y Clarke, 1979; Murray et al., 2018; Biralvand et al., 2023), más aún si se tiene en cuenta que estos eventos y procesos pueden superponerse entre sí (Fu et al., 2010), lo cual dificulta su evaluación de forma independiente (Cetina et al., 2020). Una posible solución a lo anterior es la aplicación de modelado numérico inverso (Fu et al., 2010; Murray et al., 2018), ya que, con este, la estimación del momento en que los minerales alcanzan la temperatura de cierre no se limita a comportamientos lineales, sino que además incluye aquellos escenarios donde las tasas de enfriamiento y de calentamiento varían en función de los procesos geológicos (Prenzel et al., 2013; Ketcham, 2024).
La Cordillera Central colombiana registra un entorno como el expuesto en el párrafo anterior, donde complejos metamórficos, compuestos principalmente por esquistos, anfibolitas y gneises, están afectados por intrusiones intermedias a ácidas dominadas por granodioritas (Vinasco et al., 2006; Villagómez y Spikings, 2013; Blanco-Quintero et al., 2014; Bustamante et al., 2016). Uno de los procesos que más ha marcado la historia de este ramal de los Andes es precisamente el magmatismo (Chicangana, 2005; Restrepo-Moreno et al., 2009; Bustamante et al., 2016; Duque-Trujillo et al., 2019; López-Isaza y Zuluaga, 2020; Rodríguez-García et al., 2022), como lo soportan la amplia ocurrencia espacial y temporal de cuerpos ígneos intrusivos (Leal-Mejía et al., 2019), el registro de rocas y depósitos volcanoclásticos (Bayona et al., 1994) y la presencia del actual arco volcánico que hace parte del Cinturón del Fuego del Pacífico (Marín-Cerón et al., 2019). Lo anterior debido a que, por la proximidad a la zona de convergencia entre la margen NW de Suramérica y las placas oceánicas ancestrales (Farallón) y actuales (Caribe, Nazca) (Figura 1), esta franja ha sido el tramo cortical a través del cual ha ascendido, cristalizado y emanado el material fundido, durante los distintos eventos de subducción, colisión y acreción a lo largo del Fanerozoico (Villagómez y Spikings, 2013; Spikings et al., 2015; Leal-Mejía et al., 2019; Marín-Cerón et al., 2019).

Figura 1 A. Mapa geológico (modificado de Gómez et al., 2023). El análisis se llevó a cabo en dos sectores: el primero, al norte, abarca las muestras de las secciones A-A' y B-B'; y el segundo, al sur, las muestras de las secciones C-C' y D-D' (ver Figuras 11 y 12). Las secciones y la ubicación de las muestras se presentan en la Figura 2B. B. Mapa tectónico del NW de Suramérica. Los Andes del Norte son el resultado de la extensa y compleja interacción entre placas oceánicas (al W) y continental (al E), y constan de tres ramales en Colombia: Cordillera Occidental (CO), Cordillera Central (CC) y Cordillera Oriental (CE), a las cuales separan los valles de los ríos Cauca (VC) y Magdalena (VM).
Estudios previos han estimado tasas de exhumación y enfriamiento utilizando termocronología de media y baja temperatura en la Cordillera Central colombiana (Restrepo-Moreno et al., 2009; Villagómez y Spikings, 2013; Noriega-Londoño et al., 2020; Duque-Palacio et al., 2021; Zapata et al., 2021) y profundidades de emplazamiento, mediante geoquímica y termobarometría (Bustamante et al., 2017; Cardona et al., 2018; Duque-Palacio et al., 2021; Chavarría et al., 2022), pero no se habían abordado las historias térmicas desde el emplazamiento de los cuerpos ígneos hasta su exposición, esto es, considerando además de la exhumación y la erosión el enfriamiento magmático. Como se mencionó antes, esto es posible con la obtención de modelos numéricos, a partir de la inversión de las edades geotermocronológicas en programas computacionales como 4DTherm (Fu, 2006; Fu et al., 2010), el cual no se ha empleado en Colombia. Esta herramienta deriva las historias termales y de exhumación de intrusiones ígneas, e integra datos de geocronómetros y termocronómetros (Fu, 2006; Fu et al., 2010), a la vez que considera las propiedades térmicas del cuerpo intrusivo y de la roca encajante, y por tanto la interacción termal entre ellos (Ehlers, 2005; Fu et al., 2010).
Teniendo en cuenta lo anterior, con este estudio se emplea modelamiento termal a partir de datos geotermocronológicos existentes para entender las relaciones entre intrusión, exhumación y erosión; contrastar estos resultados con aquellos de estudios previos, con el fin de proporcionar una referencia metodológica para futuros estudios. Para ello, se analizan cinco cuerpos ígneos intrusivos que registran los principales eventos magmáticos de la Cordillera Central (Chicangana, 2005; Figura 1): 1) del magmatismo Permo-Triásico, que al parecer ocurrió a lo largo del borde occidental de Pangea, un intrusivo granítico indenominado y el Intrusivo Gnéisico de La Línea (Cochrane, 2013 en Restrepo y Toussaint, 2020; Figura 2C); 2) del magmatismo del Jurásico superior, el cual se concentra en el flanco oriental, el batolito de Ibagué (Villagómez y Spikings, 2013; Bustamante et al., 2016; Rodríguez-García et al., 2022; Figuras 2B y 2C); 3) del magmatismo del Cretácico superior, el cual se concentra en el flanco occidental, el stock de Córdoba (Castrillón, 2003; Villagómez et al., 2011; Figura 2D), y 4) del magmatismo del Paleoceno-Eoceno, el cual se extendió ampliamente a lo largo de la Cordillera, el stock el Hatillo (Bayona et al., 2012; Bustamante et al., 2017; Figura 2A). Las localidades evaluadas se distribuyen en dos sectores: al norte-noreste, entre las fallas de Palestina e Ibagué (Figuras 1A, 2A y 2B); y al sur-suroeste, abarcando el sistema de fallas de Romeral y las fallas de Palestina, Chapetón-Pericos e Ibagué (Figuras 1A, 2C y 2D).

Figura 2 Secciones derivadas del mapa geológico. A. En el sector norte se evalúan muestras del stock El Hatillo (DV163) y B. En el batolito de Ibagué (DV04); en el sector sur de los intrusivos Permo-Triásicos C. Gnéisico de La Línea: DV18; Granítico: DV82 y del batolito de Ibagué: DV09, DV06 y DV05. D. En el stock de Córdoba (DV26). La falla de Armenia se desprende de la falla Cauca-Almaguer, y a diferencia de esta, se encuentra activa, afectando incluso los depósitos cuaternarios, como lo ilustran las secciones C-C' y D-D' (Guzmán et al., 1998 e Ingeominas, 2000 en Aguilar y Mendoza, 2002). El sistema de fallas de Romeral está constituido, de este a oeste (г. e., en orden de acreción), por las fallas San Jerónimo, Silvia-Pijao y Cauca-Almaguer (Chicangana, 2005). Detállese la posición de las muestras con respecto a las fallas. Note que las relaciones estratigráficas con la roca caja y estructurales con otras unidades, del batolito de Ibagué y del stock de Córdoba en B, C, y D, involucran una elevada exhumación de la roca caja debido al acople de los procesos de fallamiento y erosión. La escala horizontal de la sección C-C' se redujo a % para su completa visualización. Ver nomenclatura estratigráfica en la Figura 1. El buzamiento del plano de fallas es aproximado.
Este estudio también proveerá una alternativa para cubrir el vacío existente en los mapas geotérmicos de Colombia (Lagardère y Vargas, 2021), ya que, si bien la Cordillera Central tiene un alto potencial geotérmico, no cuenta con datos de temperatura a profundidad, como sí es el caso de las cuencas petrolíferas, donde se han colectado del fondo de pozos de gas o petróleo datos de temperatura en la corteza (BHT o Bore Hole Temperature). Para ello, se producirán, como una primera aproximación, mapas de flujo de calor y de gradiente geotérmico, a partir de los resultados del modelo obtenido para el tiempo actual.
Marco geológico
La evolución de la Cordillera Central se ha descrito en relación a la acreción y subducción de bloques alóctonos (Rodríguez-García et al., 2022) o a una compleja interacción de bloques o elementos para-autóctonos (Vinasco, 2019). Una primera interacción ocurre durante el Paleozoico tardío, un bloque denominado Chibcha colisiona con la zona como resultado de la interacción entre Laurentia y Amazonia (Restrepo y Toussaint, 2020). Dos eventos, uno en el Carbonífero y otro en el Pérmico, marcarían la primera subducción en el occidente colombiano, con un arco magmático que quedaría registrado por la intrusión de granitoides al basamento (Rodríguez-García et al., 2022).
Durante el Triásico se produjo la separación de Pangea (Wilson, 1965, 1966; Scotese, 1976), lo que implicó el desprendimiento, la anatexis y la dispersión de parte de la litosfera que contenía los gabroides (Villagómez y Spikings, 2013; Spikings et al.,, 2015; Bustamante et al., 2017); y la instauración de un rift hacia la margen continental, donde se acumularía una característica sucesión sedimentaria continental-marina, como la que registran las Formaciones Luisa y Payandé en el Valle Superior del Magdalena (Mojica y Herrera, 1986) (Figuras 1 y 2C). Inducido por la subducción de la placa Farallón por debajo del margen NW de Suramérica (Bustamante et al.,, 2016), en el lapso Jurásico temprano-Jurásico medio (Rodríguez-García et al., 2020), se desarrolló un arco volcánico que suministró la componente ígnea de las volcanosedimentitas de la suprayacente Formación Saldaña (Bayona et al., 1994; Figuras 1 y 2C).
El basamento de este bloque está constituido por las metamorfitas de Tierradentro (anfibolitas y gneises) (Villagómez y Spikings, 2013) y por cuerpos ígneos como el intrusivo granítico indenominado analizado en este estudio (Figuras 1 y 2), unidades que podrían corresponder a los fragmentos ígneo-metamórficos mencionados arriba (Bustamante et al., 2017). Sobre el basamento yacen discordantemente las metasedimentitas del Complejo Cajamarca, i.e., esquistos cloríticos, micáceos, grafiticos; anfibolitas; gneises; y mármoles (Villagómez y Spikings, 2013; Figuras 1 y 2), cuyo protolito y metamorfismo se consideran del Triásico (Villagómez et al., 2011; Villagómez y Spikings, 2013). No obstante, es posible que, por la interacción de los distintos bloques para-autóctonos, estas unidades hayan sido afectadas por un metamorfismo más reciente (Blanco-Quintero et al., 2014; Bustamante et al., 2017; Restrepo y Toussaint, 2020; Rodríguez-García et al., 2020). El Intrusivo Gnéisico de La Línea está emplazado dentro de este complejo (Cochrane, 2013 en Restrepo y Toussaint, 2020; Figura 2C). La falla Palestina (Figuras 1A, 2A y 2C), considerada como una falla de rumbo dextro-lateral, disecta al Bloque Tahamí y su edad es paleocena (Álvarez, 1983; Restrepo y Toussaint, 2020).
En asociación con la llegada del Bloque Tahamí, se formó el arco magmático que dio lugar al batolito de Ibagué (Rodríguez-García et al., 2022). Esta unidad se encuentra en el flanco oriental de la Cordillera Central (Figura 1), es el cuerpo intrusivo jurásico de mayor longitud y es disectado por la falla de Ibagué a la altura de la ciudad del mismo nombre (Bustamante et al., 2016; Rodríguez-García et al., 2020; Figura 1A). En el bloque sur dominan granodioritas a monzodioritas, y se reporta un rango de cristalización de ca. 192-152 Ma (Bustamante et al., 2016; López-Isaza y Zuluaga, 2020); en el bloque norte dominan tonalitas a granodioritas, con edades de 164 a 138 Ma (Villagómez y Spikings, 2013; López-Isaza y Zuluaga, 2020; Rodríguez-García et al., 2020).
Con el cese del magmatismo de este intrusivo hacia 138 Ma, la margen continental migró hacia el occidente como resultado de un proceso de slab roll-back (Villagómez et al., 2011; Villagómez y Spikings, 2013; Spikings et al., 2015; Bustamante et al., 2016), lo que involucró el desarrollo de una configuración extensiva donde se acumularía la sucesión continental-marina del Cretácico inferior (Sarmiento-Rojas, 2018) en el continente, y hacia la margen, el desarrollo de un arco de isla con su cuenca retroarco, donde se acumularían las volcanosedimentitas del Complejo Quebradagrande (Nivia et al., 2006). Villagómez y Spikings (2013) reportan un evento de enfriamiento por exhumación hacia 138-130 Ma en el bloque sur del batolito de Ibagué, y lo relacionan con esta configuración tectónica.
Posteriormente, se establece un régimen compresivo, donde la cuenca retroarco se cierra producto de la colisión y acreción del arco, dando lugar al metamorfismo de las sedimentitas y formando la falla de San Jerónimo y el canal de subducción que constituiría la masa del Complejo Arquía, cuya obducción-acreción sería inmediata, causaría su metamorfismo y marcaría la falla Silvia-Pijao (Villagómez y Spikings, 2013) (Figuras 2C y 2D). Con este evento también se reactivaría la falla de Ibagué, como lo sustentan estos autores con las altas tasas de exhumación hacia 117-107 Ma que se restringen al bloque norte.
La siguiente etapa comprendió la subducción de la placa oceánica y la formación de la Gran Provincia Ígnea del Caribe (meseta oceánica con su arco suprayacente; Figura 2D) (Kerr et al., 1997; Vallejo et al., 2006; Spikings et al., 2015). El material subducido fue la fuente del magmatismo que dio lugar al batolito Antioqueño (70-80 Ma; Leal-Mejía et al., 2019; Restrepo-Moreno et al., 2019) y al stock de Córdoba (~77 Ma; Castrillón, 2003; Villagómez et al., 2011; Duque-Trujillo et al., 2019; Restrepo y Toussaint, 2020). Este último es el resultado de un único pulso magmático (Castrillón, 2003), y su emplazamiento ocurrió a lo largo de la falla Silvia-Pijao (post-cinemático) (Figuras 2D), y entre los complejos Arquía (denominado Complejo Rosario al sur; Figuras 2C y 2D) y Quebradagrande (Castrillón, 2003; Villagómez et al., 2011; Espinosa-Baquero, 2020; Mora et al., 2020). Está constituido principalmente por granodioritas (Castrillón, 2003; Villagómez et al., 2011) y más localmente por xenolitos de anfibolitas del Complejo Quebradagrande (Castrillón, 2003).
Una vez consumida la porción de corteza oceánica entre Suramérica y la Gran Provincia Ígnea del Caribe, esta última colisiona y se acreciona a la margen hacia 70-80 Ma, forma la falla Cauca-Almaguer (Figura 2D), cesa con la actividad magmática continental e induce la exhumación de la Cordillera Central (Kerr et al., 1997; Villagómez et al., 2011; Spikings et al., 2015; Duque-Trujillo et al., 2019). Después de un cese de ~10 Ma (Duque-Trujillo et al., 2019), un nuevo evento magmático tuvo lugar en respuesta a la subducción somera de la boyante placa Caribe (Bayona et al., 2012). Lo anterior, sumado al movimiento transcurrente de las dos placas, explica la amplia distribución de los cuerpos ígneos a lo largo de la Cordillera y una duración del magmatismo entre 20 a 30 Ma (Bayona et al., 2012; Leal-Mejía et al., 2019; Restrepo y Toussaint, 2020). Entre estos cuerpos se encuentra el stock el Hatillo, unidad que intruye a las rocas metamórficas del Complejo Cajamarca (Bustamante et al., 2017) (Figura 2A) y para la cual se documenta una composición granodiorítica (Bustamante et al., 2017) y una edad U/Pb de ~54,6 Ma (Bayona et al., 2012).
En el lapso Eoceno medio-Eoceno tardío (~45-30 Ma), las tasas de enfriamiento y de exhumación se incrementan alcanzando su pico hacia los ~45-40 Ma (Restrepo-Moreno et al., 2009; Villagómez y Spikings, 2013); este evento cubrió toda la cadena andina y fue producto de un aumento en la convergencia entre Suramérica y la placa Farallón (Villagómez y Spikings, 2013; Restrepo-Moreno et al., 2019). Hacia el Oligoceno medio (~26 Ma), esta última se fragmenta y da lugar a las placas Nazca y Cocos (Marín-Cerón et al., 2019; Figura 1B). Hacia 23 Ma se da otro incremento en las tasas de exhumación en la porción norte de la Cordillera Central (Restrepo-Moreno et al., 2009; Duque-Palacio et al., 2021).
El arco Panamá-Chocó es el último bloque oceánico que se acreciona a la margen (Pardo-Trujillo et al., 2020) , entre el Oligoceno tardío y el Mioceno medio (Duque-Caro, 1990; Mann y Corrigan, 1990; Farris et al., 2011; Marín-Cerón et al., 2019), y ocasiona un cambio en la convergencia de la placa Nazca (Villagómez y Spikings, 2013; Duque-Palacio et al., 2021) . Lo anterior llevaría a la colisión de la dorsal de Carnegie, y esta, a su vez, a altas tasas de exhumación y enfriamiento en los Andes septentrionales, hacia ~15 Ma (Chicangana, 2005; Villagómez y Spikings, 2013).
Hacia el límite Mioceno-Plioceno, inicia una nueva actividad magmática producto de la subducción de la placa Nazca (Bayona et al., 2012; Marín-Cerón et al., 2019). Este evento incluyó la formación de la actual cadena volcánica (Marín-Cerón et al., 2019) y la intrusión de diques hipoabisales a lo largo de la zona (Restrepo y Toussaint, 2020), cuya emisión y erosión aportarían material a los depósitos sedimentarios del Neógeno (Núñez y Murillo, 1982). Finalmente, en el Pleistoceno, ocurre la sedimentación del abanico de Ibagué, alimentado por el río Combeima (Núñez y Murillo, 1982); este depósito cuaternario aflora el este y sureste de la zona de estudio (Figura 1A), y cubre parcialmente al batolito de Ibagué (Figura 2B), al intrusivo granítico del Pérmico y a las unidades sedimentarias del Triásico, Jurásico y Cretácico (Figura 1A).
Metodología
Inicialmente, se seleccionó la zona de estudio, con base en la disponibilidad de edades geotermocronológicas y teniendo en cuenta que su posición fuera próxima a las principales fallas y distal de las estructuras volcánicas, ya que de esta forma se tendría una mayor resolución para la exhumación, relacionada con el movimiento de las primeras, y una menor probabilidad de reseteo por recalentamiento inducido por las segundas; las edades consideradas se presentan en la Tabla 1.
Seguidamente, se introducen las edades, la posición de las muestras, datos de temperaturas de cierre, geometría y propiedades termales de la intrusión y de la roca encajante, en el programa 4DTherm (Fu, 2006; Fu et al., 2010), cuyas características serán detalladas en las secciones siguientes. Los últimos dos parámetros se utilizan para discretizar el enfriamiento magmático, lo cual es fundamental para una correcta discriminación de la exhumación por tectonismo y por erosión (Murray et al., 2018). Aunque la cantidad de edades utilizadas dentro del código está limitada a solo 6 muestras por área, y, por ende, modelar un sector se podría considerar insuficiente, además de que la selección de las muestras (tanto en el campo como en las referencias) no se hizo siguiendo una filosofía de muestreo como la sugerida por Schildgen y van der Beek (2019), la ventaja del código utilizado es que combina edades termocronológicas con edades U-Pb, esto no está incluido hasta ahora en ninguna otra herramienta computacional de las existentes en el área de la geotermocronología. Así, el modelamiento de la exhumación magmática y postmagmática, a través del software 4DTherm, tiene varias ventajas, especialmente por su capacidad para integrar la dimensión temporal, simular procesos complejos como la diferenciación magmática y entender las condiciones de temperatura y presión durante la exhumación. La consistencia termodinámica y la flexibilidad del modelo permiten simular una variedad de escenarios geológicos, lo que facilita la reconstrucción precisa de las historias térmicas y geológicas de las rocas magmáticas y su posterior exhumación, con aplicaciones clave en la investigación de procesos petrológicos, tectónicos y de formación de depósitos minerales.
Una vez que se obtiene la curva de enfriamiento para el conjunto de muestras que comprenden ambos sectores, se seleccionan los mejores modelos termales basados en varios criterios, entre los cuales destacan: 1) la cercanía de las edades existentes a la curva de enfriamiento (criterio de Fu); 2) edades predichas en el rango de temperatura entre 300 a 500°C, las cuales se comparan con las edades geotermocronológicas 40Ar/39Ar (criterio ventana 300-500°C); 3) mediante comparación entre las edades predichas versus las observadas; 4) se analiza el gradiente geotérmico obtenido y las medidas de flujo de calor que fuesen consistentes con la distribución de rocas presentes en superficie; finalmente, 5) la tasa de exhumación obtenida es comparada con otras existentes en la literatura. Como los criterios 1 y 2 denotan la precisión del modelo para estimar la temporalidad del emplazamiento (650-1000°C) y del enfriamiento a isotermas medias a altas (300-500°C), se valida su eficacia para el cálculo del enfriamiento magmático y más aún por el hecho de que para ello se consideran las propiedades térmicas del cuerpo intrusivo y de la roca encajante. La Figura 3 resume en forma de diagrama de flujo cada una de las etapas realizadas durante la etapa de modelamiento numérico de historias termales, considerando desde el enfriamiento hasta la exhumación.
Una vez escogido el mejor modelo numérico de las historias termales, se determinan el flujo de calor basal (Q) y la conductividad termal (λ) de la roca encajante (Figura 4), así se estima el gradiente geotérmico superficial g G mediante la expresión:
El valor de flujo de calor inicial asignado al programa se elige teniendo en cuenta que la zona es tectónicamente activa (1 W/m2), y se consideran los valores iniciales de g G y Q asumiendo equilibrio, por ello deben combinarse en los cálculos de las ecuaciones. El flujo de calor se presume constante desde la base del modelo, sin embargo, este varía al finalizar el modelamiento.
Posteriormente, los valores de flujo de calor superficial y las estimaciones del gradiente geotérmico superficial se interpolan mediante el algoritmo Natural Neighborhood de ArcMAP 10.8® para obtener los mapas respectivos. En la Figura 4, se resume la metodología empleada para estimar los gradientes geotérmicos de la zona de estudio, una vez que se ha seleccionado el mejor modelo de historias termales en 4DTherm. El paleogradiente geotérmico puede ser estimado junto con las tasas de enfriamiento, la historia termal, la paleoprofundidad de la muestra en la cámara magmática y el tiempo geológico transcurrido desde que la muestra se encontraba en la cámara hasta que asciende a superficie por las corrientes convectivas, y la exhumación por erosión.
El cálculo del g G proporciona una valiosa herramienta para identificar zonas de falla a partir de anomalías geotérmicas, las cuales pueden ser generadas por la emanación de fluidos hidrotermales, por la perturbación de la estructura termal "normal" de la corteza al superponer segmentos corticales de temperaturas diferenciales (Ring et al., 1999; Stüwe, 2007; Ring y Brandon, 2008) o por el emplazamiento de cuerpos magmáticos, como es el caso del stock de Córdoba; la discriminación de estos escenarios quedaría definida por la intensidad de la anomalía. Así mismo, si se hace una evaluación temporal del g G , se podría establecer para dichas fallas un rango de formación y uno de enterramiento (si aplica). Por otro lado, dado que los mapas generados muestran la distribución del gradiente geotérmico y el comportamiento de la temperatura, es posible estimar de forma aproximada la extensión de fallas y magmatismo intrusivo para un tiempo determinado.
Distribución de las muestras empleadas
Se modelaron dos sectores en el área comprendida entre las fallas de Armenia y de Ibagué (Figura 1). Para el sector norte se emplearon las muestras DV04 y DV163 (Figuras 2A y 2B), mientras que para el sector sur se emplearon las muestras DV18 a DV82 (Figuras 2C y 2D), resaltadas en color gris en la Tabla 1. Todas estas edades fueron obtenidas del trabajo de Villagómez y Spikings (2013).
Tabla 1 Edades utilizadas en el presente estudio (tomadas de Villagómez y Spikings, 2013). Los intrusivos: gnéisico de La Línea y el granítico indenominado se mencionan conjuntamente en el texto como "Intrusivos Permo-Triásicos", ya que corresponden al mismo evento tectono-magmático.
| Muestra | Lon. (°W) | Lat. (°N) | Elev. (m) | Edad Est. | AHe ± 1s | AFT ± 1s | ZHe ± 1s | ZFT ± 1s | Unidad geológica | Litología -Intrusivo | Unidad geológica | Roca encajante |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| DV18 | 75,5550 | 4,4719 | 3292 | 236 | 19,6 ± 4,2 | 54,2 ± 6,0 | Intrusivo gnéisico de La Línea | Gneis granodiorítico | Grupo Cajamarca | Esquistos | ||
| DV26 | 75,6900 | 4,4085 | 1389 | 80 | 6,1 ± 0,4 | 13,9 ± 2,4 | 45,6 ± 7,6 | Stock de Córdoba | Granodiorita | Complejo Rosario, Grupo Bugalagrande, Formación Quebradagrande | Anfibolitas, esquistos, diabasas y andesitas | |
| DV05 | 75,2681 | 4,4076 | 1064 | 170 ± 5,0 | 7,5 ± 0,5 | 59,8 ± 16,8 | 85,3 ± 18,2 | Batolito de Ibagué | Granodiorita | Grupo Cajamarca y Gneises y Anfibolitas de Tierradentro | Esquistos y anfibolitas | |
| DV06 | 75,2945 | 4,4024 | 1277 | 170 ± 5,0 | 35,8 ± 4,6 | 47,7 ± 2,7 | 81,3 ± 12,4 | Batolito de Ibagué | Granito | Grupo Cajamarca y Gneises y Anfibolitas de Tierradentro | Esquistos, filitas, anfibolitas y gneises | |
| DV09 | 75,3032 | 4,4082 | 1354 | 170 ± 5,0 | 30,6 ± 5,8 | 43,0 ± 3,0 | 74,6 ± 13,4 | Batolito de Ibagué | Granito | Grupo Cajamarca y Gneises y Anfibolitas de Tierradentro | Esquistos, filitas, anfibolitas y gneises | |
| DV82 | 75,2331 | 4,2876 | 919 | 272 | 11,1 ± 0,7 | 68,3 ± 8,2 | 77,6 ± 10,8 | Intrusivo granítico indenomi-nado | Granito | Gneises y Anfibolitas de Tierradentro | Gneises y anfibolitas | |
| DV04 | 74,9753 | 4,7833 | 933 | 170 ± 5,0 | 10,6 ± 0,8 | 59,2 ± 20,2 | Batolito de Ibagué | Gabrodiorita | Rocas metamorfizadas regionalmente, afectadas por metamorfismo dinámico y térmico | Filitas y esquistos | ||
| DV163 | 74,9782 | 5,1787 | 1100 | 53 | 39,6 ± 7,0 | 55,3 ± 5,4 | Plutón de El Hatillo | Granito | Rocas con metamorfismo regional posteriormente afectadas por metamorfismo dinámico y metamorfismo térmico | Filitas y esquistos |
Modelado numérico de historias termales mediante 4DTherm
En este estudio se utilizó el código 4Dtherm (Fu, 2006; Fu et al., 2010), el cual permite estudiar el enfriamiento magmático e inferir la edad y profundidad de emplazamiento de rocas ígneas intrusivas. El código resuelve un sistema de ecuaciones conservativas, momentos y energía que, según Cathles (1977), Parmentier y Spooner (1978) y Turcotte y Schubert (2002), para un fluido incomprensible en un medio poroso, puede ser expresado de la siguiente forma (la cual puede ser resuelta por esquemas de diferencias finitas explícitas bidimensionales):
Donde T es la temperatura; t, tiempo; ρ, la densidad de las rocas; C, calor específico; λ, conductividad termal; U, velocidad de Darcy; Q, la producción de calor interno; p, presión; k, permeabilidad; μ, la viscosidad dinámica del fluido, y g es la aceleración de gravedad. Los subíndices f denotan la propiedad del fluido, y b es la propiedad del fluido y roca combinados o bulk.
La ecuación 2 permite incorporar fenómenos como: 1) calor latente de cristalización/fusión asumiendo que la cristalización de las rocas ocurre entre 650 a 1000°C, 2) la convección termal en cuerpos magmáticos, 3) la convección hidrotermal en medios porosos tomando en cuenta el número de Rayleigh, 4) el enfriamiento por exhumación/erosión y 5) el flujo de fluidos.
De esta forma se combinaron los datos termocronológicos (edades AHe, AFT, ZHe y ZFT) y geocronológicos (edades 40Ar/39Ar) de 6 muestras (DV18, DV26, DV05, DV06, DV09 y DV82) (ver Tabla 1) para el sector sur y 2 muestras (DV04 y DV163) para el sector norte (Figuras 1 y 2) . Los termocronómetros en condiciones ideales documentan los episodios de enfriamiento a isotermas bajas (~60-260°C), y los geocronómetros, aquellos que ocurren a isotermas altas a medias (~300-500°C), y están vinculados a los procesos de exhumación y erosión en ambos casos; el modelo acopla dichas edades junto con la de cristalización (U/Pb), y predice la temporalidad del proceso de enfriamiento magmático, teniendo en cuenta las propiedades de los cuerpos rocosos. En este proceso se diferencian dos etapas: en la primera, el enfriamiento está marcado por la interacción con la roca caja y en la segunda, por la relajación térmica (Ring et al., 1999); el paso de una etapa a la otra se define como tiempo en tránsito (Fu et al., 2010). De esta manera, es posible recrear la historia térmica de las masas plutónicas incluyendo el enfriamiento magmático y el enfriamiento por exhumación. La Tabla 2 resume las temperaturas de cierre asumida para cada método termocronológico, mientras que en la Tabla 3 se presentan los distintos parámetros considerados en el modelamiento numérico.
El modelo 4DTherm es una aproximación de elementos finitos (Fu, 2006; Fu et al., 2010), esto permite tener en cuenta otros parámetros relacionados con el enfriamiento de los cuerpos intrusivos (por ejemplo, enfriamiento por conducción, calor latente de cristalización y fusión, convección térmica dentro de los cuerpos magmáticos y circulación hidrotermal inducida por la intrusión magmática), además de la exhumación y erosión, y facilita la cuantificación de la historia térmica y de exhumación de un cuerpo ígneo dado, ya que por su brevedad con respecto a la segunda, la primera queda delimitada por el cambio abrupto de la pendiente en la curva de enfriamiento.
Tales historias se evalúan en términos de tasas de enfriamiento y de exhumación, ya que estas permiten diagnosticar las variaciones en la magnitud de estos procesos y asociarlas a eventos magmáticos, tectónicos o morfodinámicos, según el contexto geológico. La primera es la derivada de la curva de enfriamiento (°C/Ma); la segunda es la relación entre la tasa de enfriamiento y el gradiente geotérmico (°C/ km) (corregido en el modelamiento) correspondientes al tramo analizado (km/Ma).
Tabla 2 Temperaturas de cierre para cada método de datación utilizado en la modelización mediante 4DTherm (modificado de Braun et al., 2006). Diferencias entre litologías, parámetros de conducción de calor, profundidades de emplazamiento, entre otros, fueron incorporados en el modelo.
| Método de datación | Mineral | Temperatura de cierre (°C) |
|---|---|---|
| (U-Th)/He | Apatito | 60-75 |
| Trazas de fisión | 80-110 | |
| (U-Th)/He | Circón | 150-180 |
| Trazas de fisión | 180-250 | |
| U-Pb | 700-850 |
Tabla 3 Parámetros empleados en la modelización de historias termales mediante 4DTherm.
| Parámetro | Valor | Unidad de medida |
|---|---|---|
| Gravedad | 9,8 | m/s2 |
| Calor específico (magma) | Variable | J/kg °C |
| Calor específico (agua) | 1000 | J/kg °C |
| Calor específico (roca encajante) | Variable | J/kg °C |
| Calor latente (magma) | 418000 | J/kg |
| Conductividad térmica (roca encajante) | Variable | W/m °C |
| Conductividad térmica (magma) | Variable | W/m °C |
| Coeficiente volumétrico de expansión térmica del agua | 10-3 | 1/°C |
| Coeficiente volumétrico de expansión térmica del magma | 2,0 x 10-5 | 1/°C |
| Densidad del agua | 1000 | kg/m3 |
| Densidad del magma | Variable | kg/m3 |
| Densidad de la roca encajante | Variable | kg/m3 |
| Flujo de calor basal | 0,1 | mW/m2 |
| Gradiente geotérmico inicial | Variable | °C/km |
| Permeabilidad de la roca | 1,0 x 10-15 | m2 |
| Porosidad de la roca | 1% | - |
| Profundidad de emplazamiento del magma | Variable | m |
| Temperatura inicial del magma | 1000 | °C |
| Temperatura superficial | 10-20 | °C |
| Viscosidad dinámica del magma | 1,0 x 1012 | Pa s |
| Viscosidad dinámica del agua | 2,82 x 10-4 | Pa s |
Resultados
Análisis de historias termales
Con la finalidad de facilitar la visualización de los resultados, solo se presenta la modelación mediante el código 4DTherm de la intrusión para las muestras DV06 y DV09, pertenecientes al batolito de Ibagué. Para este cuerpo se estableció como roca encajante el Complejo Cajamarca y las metamorfitas de Tierradentro (Figuras 1A, 2B y 2C), conformados por esquistos, filitas, anfibolitas y gneises. Se consideró un modelo de 20 km x 30 km para la roca encajante, desde una profundidad de 4 km se asumió una geometría rectangular para el cuerpo intrusivo de 15 km x 12 km. Así, ambas muestras se ubicaron a una profundidad cercana a los 16 km (Feininger y Botero, 1982; Figura 5A). La Tabla 4 resume los parámetros empleados para el modelamiento numérico de la historia termal, una temperatura superficial de 20°C fue asumida y se proporcionó un gradiente geotérmico inicial de 3,33°C/km, sin embargo, con cada iteración el modelo corrige este gradiente geotérmico. La historia termal completa se ilustra en las Figuras 5, 6 y 7, y se detalla a continuación.

Figura 5 Primeros estadios de la historia termal por enfriamiento magmático. A. Emplazamiento del batolito de Ibagué a los 160,11 Ma. B. Inicio inmediato del enfriamiento.
Tabla 4 Parámetros empleados en el modelamiento 4DTherm para las muestras DV09 y DV06.
El emplazamiento del magma ocurrió a 1000°C hacia 160,11 Ma (Figura 5A). La primera etapa del enfriamiento magmático, que comprende desde este evento hasta el tiempo en tránsito, tuvo una duración de 290.900 años, y una tasa promedio de 1535°C/Ma (Figuras 5 y 6). Tanto la convección dentro de la cámara magmática como la circulación hidrotermal inducida por la intrusión tuvieron una duración de 260.800 años, iniciando desde los 160,11 Ma hasta finalizar en 159,85 Ma, aproximadamente (Figuras 5B y 6A-6C). Por otro lado, con el enfriamiento desde su emplazamiento, el cuerpo ígneo alcanza una temperatura de 751,80°C hacia los 159,84 Ma, y a partir de ese momento es controlado termalmente por procesos de exhumación por erosión y fallamiento. La segunda etapa, que abarca desde el tiempo en tránsito hasta que la intrusión se enfría, esto es, cuando el equilibrio termal se alcanza y tanto el cuerpo ígneo como la roca caja se enfrían a una misma tasa aproximadamente: cooled point (Fu et al., 2010), tuvo una duración de 30.000 años y una tasa de enfriamiento promedio de 6607°C/Ma. Bajo tal escenario, se tiene que: 1) la tasa de enfriamiento magmática es ~952°C/Ma; 2) durante la primera etapa, esta tasa involucró el enfriamiento del cuerpo ígneo y el calentamiento de la roca caja, y fue incrementada por exhumación posiblemente derivada de fallamiento; y 3) durante la última etapa, el acople entre la exhumación por fallamiento y aquella por erosión (Figuras 2C y 2D) llevó a una tasa de enfriamiento mucho más alta y, por tanto, a que la estabilidad termal se alcanzara más rápido de lo esperado. En las Figuras 5 a 7, las líneas rojas, azules y verdes corresponden a las distintas isotermas. La leyenda de temperaturas consiste en 40 rectángulos coloreados gradualmente desde el verde, para temperaturas bajas, al rojo o morado para temperaturas altas. El rango de temperatura de las leyendas está fijado entre 0 y 1000°C. Por lo tanto, cada rectángulo o franja representa una isoterma de 25°C. Las temperaturas negativas o de 0°C se representan en color blanco, mientras que las temperaturas mayores o iguales a 1000°C se representan en color rojo.
El enfriamiento magmático ocurre entre 160-140 Ma (Figuras 6D y 7A-7D). De acuerdo con el modelo, para este periodo de enfriamiento se tiene que: 1) el batolito queda expuesto en superficie a ~132 Ma, con posteriores eventos de meteorización y erosión hasta la actualidad; 2) el cuerpo yace a una profundidad entre 8 y 14 km entre 90 y 148 Ma, a una tasa de ~4°C/Ma; y 3) las rocas pasan por las isotermas comprendidas entre 300 a 500°C a una edad de ~50-60 Ma. La tasa de enfriamiento de exhumación promedio fue de 4,70°C/Ma, y en cuanto a la tasa de exhumación, para el batolito de Ibagué fue de 0,09 km/Ma, y para la roca encajante se obtuvo 0,15 km/Ma; la tasa de exhumación promedio fue de 0,10 km/Ma. Estos resultados sugieren una exhumación relativamente lenta para la parte centro-sureste de la Cordillera Central (Figura 1B).

Figura 6 A-C. Historia termal para el rango entre 160,09 y 159,85 Ma, correspondiente a la primera etapa del enfriamiento magmático. En esta etapa, la interacción termal consiste en el enfriamiento del intrusivo y el calentamiento de la roca huésped, por conducción de calor y circulación hidrotermal; al interior del intrusivo tiene lugar el proceso de enfriamiento por convección termal. La combinación de estos procesos es reflejada por la cúpula de isotermas en la roca encajante, la cual indica que el gradiente geotérmico aumenta hacia las zonas de contacto, y que su comportamiento allí también depende de la distribución del calor dentro del intrusivo. Nótese la rapidez del enfriamiento en la cámara magmática: A. Dominan temperaturas altas porque gran parte del magma aún se encuentra fundido; B. Dominan temperaturas altas a medias; C. Dominan temperaturas medias a altas, pues ya se encuentra en la etapa final de la cristalización. D. Relajación de isotermas en el cuerpo intrusivo y roca encajante a los 109,00 Ma; episodio correspondiente a la historia de exhumación.

Figura 7 Episodios de erosión del cuerpo ígneo y la roca encajante desde 71,92 Ma al presente, correspondiente a la historia de exhumación. A. Enfriamiento a la temperatura de cierre del termocronómetro ZFT (240°C). B. Enfriamiento a la temperatura de cierre del termocronómetro ZHe (180°C). C. Enfriamiento a la temperatura de cierre del termocronómetro AFT (120°C). D. Enfriamiento a la temperatura de cierre del termocronómetro AHe (60°C).
Contexto de las historias termales
La historia termal del enfriamiento por la interacción cuerpo intrusivo-roca caja, se enmarca dentro el evento tectono-magmático del Jurásico superior. Este evento comenzó con la acreción de un bloque denominado Tahami, constituido por las metamorfitas de Tierradentro y del Complejo Cajamarca (roca caja) y por los intrusivos Permo-Triásicos (p. ej.: intrusivo granítico sin diferenciar), y prosiguió con la subducción de la placa Farallón y el crecimiento del arco magmático continental asociado, dentro del cual tendría lugar el plutonismo del batolito (Villagómez y Spikings, 2013; Bustamante et al., 2016; Rodríguez-García et al., 2022).
La historia termal del enfriamiento por exhumación comenzó con la superposición de esta sobre el enfriamiento magmático. Inicialmente estuvo relacionada con fallamiento inducido generado por la colisión del bloque para-autóctono Tahami y luego involucró procesos de erosión. Villagómez y Spikings (2013) reportan un evento de exhumación en el bloque sur del batolito de Ibagué para el lapso 138-130 Ma (posterior al magmatismo del Jurásico superior), asociado con la actividad de las fallas bajo el régimen extensivo resultante del proceso de slab roll-back. Es posible que este evento se haya extendido al bloque norte, ya que los tramos de sección sobre los que se ubican las muestras analizadas corresponden a estructuras horst (Figuras 2B y 2C). Sin embargo, la exhumación por erosión también fue de gran relevancia, como lo sugiere la rápida exposición del batolito (132 Ma) aquí reportada. En los siguientes estadios, el enfriamiento por exhumación fue un poco más lento y respondió tanto a erosión como al fallamiento inducido por los eventos del Cretácico superior, del Eoceno y del Mioceno (ver marco geológico).
Relación entre las historias termales
El resultado de los mejores modelos numéricos termales para el resto de las muestras es presentado en forma de diagramas tiempo-temperatura en las Figuras 8, 9 y 10. De estas figuras es posible identificar tiempo y temperatura de emplazamiento, edad y temperatura de solidificación del magma, tiempo y temperatura de exposición y las velocidades de enfriamiento magmático y post-magmático. Todas las diferentes estimaciones de tasas de exhumación y enfriamiento para la roca encajante y la roca ígnea son resumidas en la Tabla 5.

Figura 8 Historias termales obtenidas para el sector sur, usando el código 4DTherm: A. Muestra DV18 - Intrusivo Gnéisico de La Línea. B. Muestra DV26 - stock de Córdoba. C. Muestra DV05 - batolito de Ibagué. Nótese el contraste de la duración de la segunda etapa del enfriamiento magmático entre B y A-C. Ver Figuras 1A, 2C y 2D para localización de las muestras.

Figura 9 Historias termales obtenidas para el sector sur, usando el código 4DTherm: A. Muestras DV06 y DV09 - batolito de Ibagué. B. Muestra DV82 - Intrusivo Gnéisico de La Línea. Ver Figuras 1A y 2C para localización de las muestras.

Figura 10 Historias termales obtenidas para el sector norte, usando el código 4DTherm: A. Muestra DV04 - batolito de Ibagué; B. Muestra DV163 - stock El Hatillo. Nótese el contraste de la duración de la segunda etapa del enfriamiento magmático entre las dos muestras. Ver Figuras 1A, 2A y 2B para localización de las muestras.
Tabla 5 Tasas de enfriamiento y de exhumación promedio obtenidas.
| Muestra | Intrusivo | Tasa de enfriamiento promedio (°C/Ma) | Tasa de exhumación promedio (km/Ma) roca caja | Tasa de exhumación promedio (km/Ma) roca ígnea |
|---|---|---|---|---|
| DV18 | Intrusivo Gnéisico de la Línea | 3,37 | 0,09 | 0,11 |
| DV26 | Stock de Córdoba | 5,67 | 0,19 | 0,13 |
| DV05 | Batolito de Ibagué | 4,53 | 0,03 | 0,08 |
| DV09 y DV06 | Batolito de Ibagué | 4,70 | 0,15 | 0,10 |
| DV82 | Intrusivo granítico indenominado | 1,68 | 0,03 | 0,07 |
| DV04 | Batolito de Ibagué | 4,44 | 0,04 | 0,08 |
| DV163 | Plutón de El Hatillo | 4,25 | 0,25 | 0,08 |
Los modelos en el sector sur corresponden a los intrusivos Permo-Triásicos (Tabla 1), al stock de Córdoba y al batolito de Ibagué (Figuras 1A, 2C y 2D). Para el stock de Córdoba se registran tasas cercanas a 0,2 km/Ma (muestra DV26; Figura 8B), ya que se encuentra dentro del sistema de fallas de Romeral (específicamente, en la terminación norte de la falla Silvia-Pijao; Figuras 1 y 2D). Esta área está caracterizada por un tectonismo activo desde su formación (Vinasco et al., 2006; Marín-Cerón y Vinasco-Vallejo, 2012; Vinasco, 2019). Para el batolito de Ibagué, las tasas de exhumación varían de 0,03 a 0,15 km/Ma, y se observan diferencias significativas en los modelos de las muestras DV06 y DV09 (Figura 9A) que en el de la muestra DV05 (Figura 8C), posiblemente esta diferencia se deba a la cercanía de las fallas (Figuras 1A y 2C), que posiblemente haya alcanzado su mayor exhumación por adelgazamiento cortical durante el Cretácico inferior (Villagómez y Spikings, 2013; Figuras 2B y 2C).
Los modelos de las muestras DV163 y DV04 se encuentran en el sector norte y corresponden al stock El Hatillo (norte) y al batolito de Ibagué (sur), respectivamente (Figuras 1A, 2A, 2B y 10). Los resultados obtenidos muestran en promedio tasas similares de exhumación para la roca ígnea de antes de ~0,1 km/Ma, lo cual sugiere un control regional para este bloque. No obstante, las tasas de exhumación de la roca caja son notablemente más altas en la parte norte (Figura 10B) (~0,3 km/ Ma), a pesar de que en ambos casos corresponde a las metamorfitas de Tierradentro y del Complejo Cajamarca (Figuras 2A y 2B), la diferencia se debe a la presencia de las fallas (Figura 2A). Esto también aplica para el caso de las muestras DV06-DV09 y DV26 del sector sur (Figuras 2C y 2D), donde es evidente la elevada exhumación de la roca caja por la acción combinada del fallamiento y la erosión.
En general, se puede establecer que las tasas de enfriamiento y de exhumación más altas (0,1-0,2 km/Ma) se encuentran en el sector sur, dado que la actividad tectónica es mayor allí, como lo sugieren el estrechamiento de la corteza y el respectivo incremento de fallamiento (Figura 1A). También, es posible precisar un comportamiento uniforme en la historia termal del batolito de Ibagué. Por otro lado, se observa que, en contraste con este cuerpo, la segunda etapa del enfriamiento magmático fue más lenta que la primera en los demás intrusivos, lo cual sugiere que la exhumación (tanto por fallamiento como por erosión) jugó un rol más importante en el evento tectono-magmático del Jurásico superior, que en los del Permo-Triásico, Cretácico superior y Paleoceno-Eoceno.
Discusión
Comparación de tasas de exhumación y enfriamiento
Los resultados obtenidos sugieren una tasa de exhumación promedio de 0,1 km/Ma para los dos sectores analizados; las tasas más bajas oscilan entre 0,03 y 0,08 km/Ma, y las más altas, entre 0,09 y 0,25 km/Ma, que corresponden a zonas distantes y cercanas a fallas, respectivamente. La tasa de enfriamiento por exhumación promedio es de ~4,2°C/km.
Las tasas de exhumación son consistentes con los valores documentados por estudios previos para el altiplano Antioqueño (al norte de la zona de estudio), como el de Restrepo-Moreno et al. (2009, 2019), quienes reportaron un rango de tasas de erosión de 0,2-0,4 km/Ma para los eventos de exhumación del Cretácico superior, Eoceno y Mioceno, y una tasa promedio de ~0,04 km/Ma para los lapsos de quiescencia tectónica; y el de Noriega-Londoño et al. (2020), quienes obtuvieron resultados de 0,6-0,9 km/ Ma para el Paleoceno-Eoceno, y de 0,01-0,02 km/ Ma después del Eoceno. Los resultados obtenidos de este modelamiento son similares a los obtenidos por Duque-Palacio et al. (2021), quienes reportan un incremento gradual en las tasas de exhumación, desde ~50 Ma esta velocidad es de 0,1 km/Ma, posteriormente desde 23 Ma la tasa de exhumación se incrementa a 0,2 km/Ma, y desde 10 Ma aumenta a 0,3 km/Ma.
Sin embargo, todos los estudios antes mencionados y el presente difieren de los resultados obtenidos por Zapata et al. (2021), entre 72 a 60 Ma, quienes reportaron altas tasas para el batolito Antioqueño, de 0,9-2,5 km/Ma (según el modelo 1) y de 0,71,5 km/Ma (según el modelo 2), entre 72 a 60 Ma. Posiblemente, estas tasas de exhumación derivadas por esos autores resultan en una sobreestimación al considerar fijas temperaturas altas en la base del modelo, ya que en el código usado por esos autores (Pecube, Braun, 2003; Braun et al., 2006, 2012) no fue incorporado fenómenos de magmatismo. Para esto es necesario modificar el código, con el fin de incorporar anomalías térmicas, tal como lo muestra Brichau et al. (2023) para la zona de Chiapas. En contraste, las tasas obtenidas por Zapata et al. (2021) de 0-0,3 km/Ma entre 60 y 20 Ma y 0-0,2 de 20 Ma son más acordes a las obtenidas mediante modelamiento 4DTherm. En términos generales, podría decirse que las tasas de exhumación en la región son relativamente consistentes, lo que podría indicar un control tectónico regional más uniforme, por lo menos al norte de la falla de Ibagué y desde finales del Mesozoico.
Por otro lado, Villagómez y Spikings (2013) determinaron una tasa de exhumación de ~0,3 km/Ma, y una de enfriamiento de 9°C/Ma para el batolito de Ibagué y su roca caja en zonas cercanas a las fallas de Ibagué, Palestina y Otú-Pericos, durante el lapso 4530 Ma; se encontró concordancia y discrepancia con nuestros resultados, respectivamente. La diferencia entre las tasas de enfriamiento se debe a que estos autores abordan las historias termales a partir de temperaturas mucho más bajas (~350-300°C), sin considerar que los efectos de la exhumación son significativos desde los primeros estadios del enfriamiento magmático, como se evidencia en este estudio. Esto resalta la importancia de estimar las historias termales desde el emplazamiento del magma, y, en la misma medida, la utilidad del modelado numérico inverso en 4DTherm para tal fin.
En términos generales, los autores antes mencionados convergen en valores similares de tasas de exhumación a partir del Paleógeno, sin embargo, durante el Jurásico-Cretácico es importante considerar el flujo de fluidos, la profundidad de emplazamiento, el hecho de que las rocas se transporten a niveles más someros de la corteza a través de procesos magmáticos y como fluidos. En este caso, es aconsejable que la comunidad científica combine la mayor cantidad de geotermocronómetros apoyados con modelamiento de geotermobarimetría, o modelado termomecánico, métodos considerados en el modelamiento numérico computacional (Gerya, 2019).
Aplicaciones para estimaciones de mapas de gradiente geotérmico y flujo de calor
La distribución del flujo de calor en la región estudiada muestra valores uniformes entre 99,99 y 100 mW/m2 (Figuras 11 y 12), lo cual es sorprendentemente alto, pero se encuentra dentro del margen de error aceptable, según Davies (2013). Este patrón indica que, aunque el flujo de calor es alto, es relativamente constante a lo largo del área estudiada. Sin embargo, cuando se compara con el mapa de referencia (p. ej.: Davies, 2013), se observa una diferencia significativa en los rangos de valores, ya que este último reporta valores que varían entre ~18 y 120 mW/m2 para la región. Esta discrepancia podría deberse a variaciones en la resolución de los modelos o a diferencias en los datos geotermales utilizados en los estudios.
En cuanto a los gradientes geotérmicos superficiales, el rango de valores obtenido (31,25 a 40°C/km) es más alto que el rango reportado por Alfaro et al. (2009), que se encuentra entre 20,4 y 32,2°C/km para la misma región. Esta diferencia podría reflejar un mayor control tectónico y una mayor influencia del vulcanismo actual en los gradientes térmicos en la región, especialmente en las áreas cercanas a fallas activas. Las similitudes en la distribución espacial de los gradientes geotérmicos entre ambos estudios sugieren que, a pesar de las diferencias en los valores absolutos, los patrones generales de variación térmica están bien representados en ambos modelos.

Figura 11 Comparación de estimaciones de flujos de calor mínimos, máximos y promedios derivados del modelo de Davies (2013) y el obtenido mediante el modelado de la historia termal de los distintos sectores considerados mediante el código 4DTherm.

Figura 12 Comparación de mapas de gradientes geotérmicos: A. Mapa preliminar de gradientes geotérmicos de Colombia (Alfaro et al., 2009). B. Mapa obtenido a partir del modelamiento numérico de datos geotermocronológicos.
Implicaciones geológicas, climáticas y futuras investigaciones
Los resultados de este estudio aportan una comprensión más detallada de las tasas de exhumación y enfriamiento en la Cordillera Central, así como de los patrones de flujo de calor y gradientes geotérmicos. Los modelos acá utilizados permiten discriminar entre enfriamiento magmático y enfriamiento post-magmático por exhumación (erosión y/o fallamiento). La velocidad de enfriamiento magmático depende de varios factores, entre los cuales destacan: la temperatura inicial del magma, la profundidad a la que se encuentra, el tipo de roca y el ambiente en el cual se encuentre el magma. Generalmente, la roca caja que rodea una intrusión magmática, y la roca encajante que contiene esa intrusión, se exhuman simultáneamente, pero existe la posibilidad de que ambas experimenten exhumación diferencial debido a factores como composición mineralógica, estructura tectónica, profundidad de donde provienen, la velocidad de enfriamiento del magma, el flujo de fluidos a niveles someros y las condiciones geológicas en las que se encuentran.
En esta investigación se observan discrepancias respecto a las tasas de exhumación calculadas; tales discrepancias se pueden deber a la filosofía como fueron tomadas las muestras. Por ejemplo, las muestras no son tomadas a lo largo de perfiles verticales de elevación, que permitan estimar tasas de exhumación (Braun, 2002) o distinguir si las muestras hoy en día se sitúan en antiguas zonas de reseteo parcial exhumadas. Sin embargo, tales discrepancias subrayan la necesidad de realizar más investigaciones para comprender mejor las variaciones espaciales y temporales en la evolución térmica y tectónica de la región y considerar otros fenómenos como magmatismo, flujo de fluidos, fallamiento y cambio de relieve. Autores como Bustos et al. (2013) y Bustos y Bermúdez (2015) mencionan un comportamiento diferencial en los patrones de exhumación al norte y al sur del batolito Antioqueño, esto se explica en términos de los procesos de interacción entre clima, relieve y tectónica; estos últimos procesos junto con el magmatismo deben ser analizados con detalle a futuro.
Futuras investigaciones podrían beneficiarse de la integración de datos geofísicos y geológicos adicionales, así como de la utilización de modelos termocinemáticos más avanzados que puedan capturar de manera más precisa las complejidades de la dinámica tectónica y térmica en la Cordillera Central. Esto permitiría una mejor evaluación de los procesos geológicos que han moldeado la región y proporcionaría una base más sólida para la comparación con otras regiones andinas.
Conclusiones
Una particularidad de la historia termal de cuerpos ígneos intrusivos es que no solo está definida por el enfriamiento por exhumación o post-magmático, sino también por el enfriamiento magmático. Su estimación respecto a este último es compleja, ya que la exhumación generalmente se superpone, y dada su corta duración (del orden de cientos de miles de años), se requiere una muy alta resolución temporal. Los resultados de este estudio evidencian lo anterior y documentan cómo estas dos historias termales se pueden discriminar desde un enfoque cuantitativo, discretizándolas para cinco cuerpos ígneos que registran los principales eventos magmáticos de la Cordillera Central, mediante modelamiento numérico inverso en 4DTherm, código que integra datos geotermocronológicos y considera las propiedades térmicas de las rocas.
La historia termal por enfriamiento magmático comprende dos etapas. En la primera, el cuerpo intrusivo se enfría y la roca caja se calienta, por los procesos de convección del magma, conducción del calor y circulación hidrotermal; en la segunda, los dos se enfrían hasta la relajación de las isotermas. Los modelos obtenidos indican que la segunda etapa fue más lenta que la primera en todos los cuerpos evaluados, con excepción del batolito de Ibagué, para el cual se reporta lo contrario. Esto sugiere que la exhumación (por fallamiento y por erosión), tuvo una mayor implicación en el evento tectono-magmático del Jurásico superior que en los del Permo-Triásico, Cretácico superior y Paleoceno-Eoceno.
Las mayores tasas de enfriamiento y de exhumación se registran en las inmediaciones de las zonas de falla. Así pues, los valores más altos corresponden al stock de Córdoba (5,7°C/km; 0,13 km/Ma), el cual se encuentra dentro de la zona del sistema de fallas de Romeral (Figuras 1A y 2D), caracterizada por un alto tectonismo desde su conformación. Este patrón también se observa dentro de los mismos cuerpos, pues los modelos obtenidos de muestras ubicadas en cercanía a estas estructuras arrojan valores más altos que aquellos de muestras más retiradas. La magnitud de esta diferencia varía en función de la distancia: en el caso del batolito de Ibagué, en el sector sur, donde la dispersión es relativamente corta, se encuentra que las tasas varían de 4,70°C/km, 0,15 km/Ma (roca caja) y 0,10 km/Ma (roca ígnea) a 4,53°C/km, 0,04 km/Ma (roca caja) y 0,08 km/Ma (roca ígnea), al alejarse de la zona de falla (Figura 2C); en el caso de los intrusivos Permo-Triásicos, cuyas muestras se encuentran notablemente dispersas, se observa un marcado contraste al pasar de tasas de 3,37°C/km, 0,09 km/Ma (roca caja) y 0,11 km/Ma (roca ígnea) en el área de la falla de San Jerónimo, a tasas de 1,68°C/km, 0,03 km/Ma (roca caja) y 0,07 km/Ma (roca ígnea) mucho más al este de la falla de Ibagué (Figura 2C). De la misma forma, en el sector norte, se reporta una alta tasa de exhumación para la roca encajante del stock El Hatillo (0,25 km/Ma), cuyos contactos son fallados (Figura 2A)
Las tasas de enfriamiento y exhumación en los dos sectores estudiados varían de 1,68 a 5,67°C/Ma y de 0,03 a 0,25 km/Ma, respectivamente. Los resultados sugieren que las tasas de exhumación/erosión a lo largo de la Cordillera Central no pueden ser asumidas como constantes desde el inicio del emplazamiento hasta la actualidad, ya que pueden variar por afectaciones en la tectónica y demás procesos geológicos a lo largo del tiempo. Las tasas de exhumación están en línea con las reportadas en estudios anteriores para la zona de estudio y el área Antioquia, lo cual sugiere un posible control tectónico regional para la región al norte de la falla de Ibagué, a partir del Cretácico superior. Por el contrario, las tasas de enfriamiento difieren de las propuestas previamente, y ello se debe a que no se habían estimado desde el emplazamiento del magma. Todo lo anterior denota la importancia de evaluar completamente las historias termales de cuerpos ígneos intrusivos, lo que implica la discriminación entre enfriamiento magmático y enfriamiento por exhumación.
Con respecto a los mapas de flujo de calor, se determinan valores cercanos a los 100 mW/m2, mientras que el gradiente geotérmico superficial estaría entre 31,25 y 40°C/km; estos valores tienden a ser más representativos para zonas cercanas a la actual cadena volcánica. Este análisis espacial se podría utilizar para identificar zonas de fallamiento, plutonismo y vulcanismo en un momento dado, a partir de las anomalías geotérmicas que estos generan, y sí la temporalidad de dicha evaluación se amplía, se podría también detallar su evolución y, por tanto, su relación con el contexto tectónico.


















