1. Introducción
Hacia el flanco oriental de la Cordillera Oriental de Colombia, existen importantes vías de comunicación que unen la capital del país con zonas planas (Llanos Colombianos), caracterizadas por ser carreteras de piedemonte con frecuentes procesos de remoción en masa. Los deslizamientos pudieran causar otros fenómenos como inestabilidad de infraestructuras, pérdidas humanas y materiales, mientras modifican el relieve, generando un alto impacto sobre la evolución del paisaje (Campforts et al., 2020).
En virtud de que posibles interacciones entre procesos climáticos y tectónicos que ocurren en espacios geográficos humanizados podrían ocasionar daños materiales y pérdidas humanas, el estudio de zonas tectónicamente activas, y en particular del sistema montañoso andino, adquiere relevancia desde la perspectiva de la ingeniería geológica o geología aplicada debido a la complejidad geológica de este sistema. Mediante este enfoque se puede lograr un acercamiento a la comprensión de los procesos naturales que operan en la superficie terrestre y que determinan la evolución del paisaje, como el desplazamiento de fallas geológicas, el transporte de agua y sedimentos en las redes fluviales, y la ocurrencia de procesos erosivos; dinámicas que afectan el desarrollo y operación de infraestructura de diversos tipos y por lo tanto definen la gestión del riesgo de origen natural.
La presente investigación explora la forma en que operan las interacciones entre el patrón de precipitaciones actuales y la tectónica como procesos controladores del relieve actual a corto y largo plazo en la cuenca del río Guayuriba, donde en los últimos años, importantes procesos de inestabilidad asociados a movimientos en masa se han reportado (p.e., km 58) (Figura 1), Cordillera Oriental de Colombia, tomando como base el modelamiento numérico de datos termocronológicos existentes, en particular edades de trazas de fisión en apatito (por sus siglas en inglés: Apatite Fission-Track, AFT) y trazas de fisión en circón (por sus siglas en inglés: Zircon Fission-Track, ZFT), parámetros geomorfométricos, datos sismológicos y litología.
Para ello se siguió la metodología propuesta por Bermúdez et al. (2011a, 2011b, 2013), Bustos y Bermúdez (2015), quienes combinaron datos geológicos recolectados en campo con métodos de termocronología cuantitativa (Braun, 2003; Braun et al., 2006, 2012), con lo que se desarrolló el modelado de la historia de exhumación de al menos 35 superficies de erosión del altiplano antioqueño en Colombia, discriminando en cada una de esas superficies los posibles controles tectónicos, climáticos o la interacción de ambos como agentes rejuvenecedores del paisaje actual, lo cual constituyó la primera aplicación del software PECUBE® en Colombia.
2. Marco geológico
La Cordillera Oriental de Colombia hace parte de los Andes del Norte de Suramérica, de las tres cordilleras colombianas (Figura 1a), esta se exhumó como resultado de una interacción compleja entre las placas Caribe, Nazca-Cocos, Suramérica y acreciones de bloques (e.g.,Cediel et al., 2003). Según Mora et al. (2006), al igual que muchos otros autores, esta cordillera corresponde a un orógeno que se invirtió durante una fase compresiva cenozoica, valiéndose de estructuras reactivadas de edad mesozoica producidas durante una fase de rift. Entre esas importantes estructuras reactivadas que hoy en día producen la carga tectónica y el acortamiento de la montaña, están las fallas Servitá, San Juanito y Naranjal, situadas en el flanco oriental del orógeno (Figura 1b). Mora et al. (2008), gracias a la incorporación de datos sísmicos de la Agencia Nacional de Hidrocarburos (ANH), una estratigrafía detallada, estimaciones de paleoelevación calibradas con paleobotánica, reflectancia de vitrinita y nuevos datos de trazas de fisión en apatito (AFT), reconstruyen la retrodeformación del flanco Oriental de la Cordillera Oriental en el Macizo de Quetame explicando la asimetría de esta estructura regional debido a su acortamiento y migración del frente orogénico al oriente. El crecimiento topográfico máximo ocurrido entre 6 Ma a 3 Ma formó una barrera orográfica en el flanco oriental que posteriormente interceptó los vientos húmedos provenientes del Amazonas hasta la actualidad, lo cual ocasionó que la erosión por precipitaciones se concentrara en esta parte de la cadena. Mora et al. (2008) sugieren que esta erosión focalizada a lo largo del flanco oriental del orógeno aumentó la pérdida de masa (erosión) en este sector a valores comparables con tasas de exhumación controladas solamente por la tectónica. En contraste, debido a la barrera orográfica oriental, la precipitación se redujo en los sectores occidentales del orógeno durante el mismo período, y la cuenca se drenó internamente. En consecuencia, la eliminación de masa y el flujo tectónico se redujeron en el flanco occidental, por lo cual Mora et al. (2008) proponen que un comportamiento heterogéneo eventualmente produjo la asimetría tectónica regional. El flanco oriental de la cordillera destaca por sus contrastes topográficos y climáticos en comparación con la Sabana de Bogotá en la parte central y el flanco occidental, debido a que el flanco oriental tiene las elevaciones medias más altas y una topografía irregular con profundos cañones fluviales que han disectado hasta las rocas del basamento.
El área de estudio (Figura 1) se encuentra en el flanco oriental de la Cordillera Oriental, localmente aparece como margen oriental de la cordillera el sistema de fallas de Guaicáramo al suroriente de la zona de estudio, siendo éste de los sistemas de fallas más activos en Colombia (Cooper et al., 1995; Audemard, 1999; Chicangana et al., 2007). En la Orogenia Andina, estas estructuras basales mesozoicas de origen extensional como las fallas normales de Naranjal, Servitá y San Juanito jugaron un papel fundamental en la determinación del lugar de deformación durante la reactivación contractiva cenozoica durante la cual sufrieron inversión (Mora et al., 2006), concentrando el acortamiento y engrosamiento en el flanco oriental. Dicho engrosamiento ha exhumado las rocas del basamento a elevaciones más altas que en el flanco occidental, lo que ha provocado la exposición del basamento en elevaciones de 4000 msnm.
La mayor parte de la zona considerada en esta investigación hace parte de la cuenca del Río Guayuriba (Figura 1b), caracterizada por el desarrollo de un valle ancho con relieve ondulado y de clima frío; con topografía más abrupta hacia el sur de la cuenca en los dominios del Macizo de Quetame conformada por las rocas metamórficas y meta-sedimentarias de dicho macizo, con relieve escarpado, altas pendientes, valles profundos y topografía abrupta.
A lo largo del flanco oriental de la Cordillera Oriental, en el Macizo de Quetame (Figura 1a), el basamento exhumado constituye hasta 4 km de rocas metamórficas de grado bajo y medio, del Paleozoico medio (Segovia, 1963; Ulloa y Rodríguez, 1979). Las rocas superpuestas en disconformidad del Paleozoico Tardío comprenden del Devónico al Carbonífero el Grupo Farallones (Ulloa y Rodríguez, 1979). Luego ocurrió el episodio de rifting en el Mesozoico el cual generó un sistema de cuenca extensional amplio delimitado por el sistema de fallas La Salina hacia el oeste, y el sistema de fallas Servitá-Lengupá, hacia el este (Figura 1a, Mora et al., 2006; Sarmiento-Rojas et al., 2006). A lo largo del pie de estas fallas normales, la erosión condujo a la eliminación total o parcial de las secuencias del Paleozoico (Mora et al., 2006). La paleotopografía resultante se conservó debajo de las subsiguientes unidades marinas poco profundas del Cretácico Inferior, depositadas durante la etapa activa de rifting. Como resultado, cambios drásticos de facies y espesores (4.5-6.0 km) caracterizan los sedimentos del Cretácico Inferior (Mora et al., 2006, 2009) (ver Figura 1). Estas unidades están cubiertas por 1.5-2 km de unidades del Cretácico Superior (Formaciones Une y Chipaque, y el Grupo Guadalupe), depositadas durante una etapa de subsidencia termal post-rift (Sarmiento-Rojas et al., 2006). La sección superior de la secuencia, que va desde el Mioceno Temprano al Mioceno Tardío (Parra et al., 2006), está compuesta por estratos continentales cuyo ambiente depositacional registra la transición de condiciones marinas a terrestres que se caracterizaban por la existencia de ríos meandriformes y sistemas de abanicos aluviales proximales (Mora et al., 2009).
3. Bases de datos y métodos
La litología y las estructuras geológicas de la zona de estudio fueron extraídas del mapa geológico 1:500,000 del Servicio Geológico Colombiano (Gómez et al., 2015). Adicionalmente, los valores de espesor cortical hc para la Cordillera Oriental fueron obtenidos del modelo CRUST 1.0 (Bassin et al., 2000; Laske et al., 2013). Los principales rasgos geomorfométricos se seleccionaron a partir de varios mapas a escala 1:100,000 (Plancha 266 Villavicencio, Plancha 247 Cundinamarca) publicadas por INGEOMINAS (2011).
Información morfométrica adicional fue calculada a partir del empleo de un modelo de elevación digital (DEM por sus siglas en inglés) con una resolución de 30 m por píxel descargados del sensor ALOS (https://www.eorc.jaxa.jp/ALOS/en/aw3d30/index.htm). Además, fue recopilada información histórica sobre movimientos en masa ocurridos en el área de estudio a partir del Sistema de Información de Movimientos en Masa (SIMMA, 2020) del Servicio Geológico Colombiano. Las unidades litológicas de las cuales se tienen datos termocronológicos se resumen en la Tabla 1.
Símbolo | Grupos y Unidades |
Descripción litológica | Edad |
---|---|---|---|
K2cp | Formación Chipaque |
Lutitas, calizas, fosforitas, cherts y cuarzoarenitas. Predominio de facies finas al norte del Cocuy y facies más arenosas al sur (Guerrero, 2002). |
Cenomaniano- Aastrichtiano |
K1c | Formación Fomeque |
Alternancia de conjuntos de lutitas negras y conjuntos de calizas y areniscas calcáreas (Renzoni, 1962). |
Berriasiano- Cenomaniano |
Grupo Cáqueza |
Arcillolitas y limolitas negras con intercalaciones menores de arenitas y calizas. Segmentos de cuarzoarenitas de grano fino a grueso y conglomerados (Pulido y Gómez, 2001). |
||
Jsb | Formación Brechas de Buenavista |
Conglomerados y brechas; arenitas de grano fino a conglomeráticas y calizas estromatolíticas, e intercalaciones de lodolitas negras (Renzoni, 1968; Beltrán, 2005). |
Titoniano |
DCf | Grupo Farallones |
Cuarzoarenitas, arcillolitas, lodolitas grises y, ocasionalmente, calizas y conglomerados (Pulido y Gómez, 2001). |
Devónico- Carbonífero |
PEqgu | Filitas y Cuarcitas de Guayabetal |
Filitas verdosas y cuarcitas (Campbell y Bürgl, 1965). |
Ordovícico- Llandovery |
3.1. BASES DE DATOS TERMOCRONOLÓGICAS
En estas formaciones fueron tomados 25 datos termocronológicos de edades AFT y ZFT, por Mora et al. (2008) y Parra et al. (2009), los cuales fueron usados en este trabajo para la modelización numérica (Tabla 2).
Código | Latitud | Longitud | Nomenclatura | Unidad o Grupo | AFT | error AFT | ZFT | error ZFT |
---|---|---|---|---|---|---|---|---|
BV-MP3-F BV-MP7-F |
4.4427 4.4349 |
-74.0207 -73.8003 |
K2cp K2cp |
F. Chipaque | 13.8 14.3 |
2.6 3.8 |
N/A N/A |
0 0 |
SJ-2a | 4.4304 | -74.0094 | K2cp | G. Cáqueza | 1.2 | 0.4 | N/A | 0 |
BV 120-F | 4.3979 | -73.6802 | K2cp | 2 | 0.7 | N/A | 0 | |
BV-151 | 4.3705 | -73.9543 | K2cp | 3.8 | 0.7 | 18.5 | 1 | |
BV126 | 4.3772 | -73.9008 | K2cp | N/A | 0 | 13.3 | 1 | |
BV-65 | 4.3970 | -73.9087 | K2cp | N/A | 0 | 61.3 | 6.4 | |
SJ-1 | 4.4438 | -73.9266 | K2cp | N/A | 0 | 5.9 | 0.4 | |
BV-121 | 4.3968 | -73.6809 | K2cp | N/A | 0 | 136.5 | 13.4 | |
BV-427 | 4.1986 | -73.9533 | Jsb | Formación Brechas de Buenavista |
3 | 0.4 | N/A | 0 |
WPT | 4.1962 | -73.6438 | DCf | G. Farallones | 2.9 | 2.1 | N/A | 0 |
BV-283 | 4.2796 | -73.7124 | DCf | N/A | 0 | 145.2 | 17.3 | |
BV-192 | 4.2674 | -73.6091 | PEqgu | Filitas y Cuarcitas de Guayabetal |
2.3 | 0.6 | 9 | 0.5 |
BV-195 | 4.2933 | -73.8003 | PEqgu | 2.6 | 0.5 | N/A | 0 | |
BV-196 | 4.2925 | -73.7828 | PEqgu | 2.6 | 1.1 | 11.4 | 1.1 | |
BV-277 | 4.1614 | -73.7834 | PEqgu | 2.7 | 0.3 | N/A | 0 | |
BV-279 | 4.1624 | -73.6702 | PEqgu | 2.6 | 0.3 | 165.9 | 12.9 | |
FT-1 | 4.4889 | -73.6690 | PEqgu | N/A | 0 | 13.1 | 0.9 | |
SJ-5 | 4.4716 | -73.7312 | PEqgu | N/A | 0 | 11.5 | 0.6 | |
BV-423 | 4.2029 | -73.6931 | PEqgu | 1.6 | 0.4 | N/A | 0 | |
BV-90 | 4.2245 | -73.7951 | PEqgu | 2.8 | 0.4 | N/A | 0 | |
BV-86 | 4.2789 | -73.8941 | PEqgu | N/A | 0 | 8.3 | 0.5 | |
BV-194 | 4.2976 | -73.8346 | PEqgu | N/A | 0 | 9.8 | 0.6 | |
BV-90 | 4.2251 | -73.7771 | PEqgu | 2.8 | 0.4 | N/A | 0 | |
FT-6 | 4.4599 | -73.8938 | PEqgu | 2.4 | 0.7 | N/A | 0 |
Adicionalmente, se compilaron 7 edades y tasas de incisión fluvial mediante isótopos cosmogénicos (Berilio-10) realizadas en terrazas del Río Guayuriba por Struth et al. (2015), las cuales fueron interpoladas a lo largo de la zona de estudio usando el método IDW de ArcGIS 10.8.
3.1.1. ÍNDICE DE EROSIÓN A LARGO PLAZO
Las distribuciones de las edades termocronológicas correspondientes a edades AFT y ZFT se utilizaron para calcular las tasas de exhumación o de erosión a largo plazo, asumiendo que la cantidad de erosión es igual al levantamiento, esto equivale a decir que la topografía se encuentra en estado de equilibrio (Perne et al., 2017).Con este propósito se empleó la rutina Age2Edot desarrollada por Brandon et al., (1998). Los parámetros de entrada de esta rutina son: 1) el tipo de mineral fechado (apatito o circón), 2) el método utilizado (trazas de fisión o (U-Th)/He), 3) la difusividad termal, habiendo en este caso usado 25 km2My-1, 4) el espesor del modelo de 40 km basado en el espesor cortical promedio del modelo CRUST1.0 usado, y 5) una temperatura en la base del modelo fija de 400 °C, al igual que la considerada para los modelos termocinemáticos. Esto permitió convertir cada edad termocronológica en una tasa de erosión que es comparada con los otros métodos.
3.1.2. MODELAMIENTO TERMOCINEMÁTICO TRIDIMENSIONAL
Con el fin de establecer la historia de exhumación de este sector de la Cordillera Oriental, se utilizaron datos termocronológicos AFT y ZFT existentes (Tabla 2) tomados de Mora et al. (2006) y Parra et al. (2009), se utilizó el DEM, y los parámetros resumidos en la Tabla 2. Todos estos fueron incorporados al código termocinemático 3D PECUBE en su versión directa ("forward") desarrollado por Braun (Braun, 2003; Braun et al., 2006, 2012) el cual permite integrar distintos parámetros termales, flexurales, topográficos y estructurales (orientaciones y geometrías de fallas). El código resuelve la ecuación del calor en tres dimensiones, permitiendo predecir edades termocronológicas para distintos sistemas, y elevaciones. De estas predicciones, fue posible generar perfiles edad-elevación sintéticos, estimar tasas de erosión a largo plazo; además de estudiar la influencia de zonas de fallas en la región analizada.
Los parámetros texturales y flexurales restringen la respuesta termal del área modelada a cambios de exhumación y relieve. La historia termal de puntos en la superficie se usa para predecir las edades de AFT y ZFT tomando en cuenta la temperatura de cierre de cada sistema termocronológico (Wagner y Van den Haute, 1992).
Se realizó un total de 10000 modelos directos "forward" asumiendo tres y cuatro fases de exhumación. Se asume el periodo inicial del modelo y algunos parámetros flexurales y termales fijos como datos de entrada del código, los cuales son mostrados en la Tabla 3.
Parámetro | Valor | Referencias |
---|---|---|
Tiempo inicial del modelo | 40 (Ma) | |
Densidad de la corteza | 2700 (kg/m3) | Gómez et al. (2005) |
Densidad de la astenosfera | 3200 (kg/m3) | Gómez et al. (2005) |
Módulo de Young | 1.d11 | Gómez et al. (2005) |
Coeficiente de Poisson | 0.25 | Gómez et al. (2005) |
Espesor del modelo | 40 (km) | (Bassin et al., 2000; Laske et al., 2013) |
Difusividad termal | 25 (km2/Ma) | Braun y Robert (2005) |
Temperatura en la base del modelo | 400 (°C) | Bermúdez et al. (2011a) |
Temperatura a nivel del mar | 22 (°C) |
En cada iteración se busca la obtención del mejor ajuste entre edades observadas y predichas ("misfit"). Por esta razón, se computan diversos modelos forward para abarcar una amplia gama de posibles soluciones en diferentes escenarios, predominantemente bajo el control de datos estructurales, litológicos y geotermales existentes en la literatura. La selección del mejor modelo se basa en la minimización del error cuadrático medio o 'misfit' en la Ecuación 1.
Donde n es el número de observaciones, tobs y tcal corresponde a las edades observadas y predichas por PECUBE para cada una de las muestras i.
3.2. PRECIPITACIÓN
Los registros meteorológicos de precipitación, humedad y viento durante el período 1990 - 2019 fueron provistos por el Instituto de Hidrología, Meteorología y Estudios Ambientales de Colombia (IDEAM). En la Tabla 4, se presentan las estaciones meteorológicas utilizadas y los valores promedios en el intervalo de tiempo antes mencionado.
Estación | Gutiérrez | Servita | El Calvario | Monfort | Fómeque |
Buena- Vista |
Monter- Redondo |
Susu-Muco | Las Casas | Llano Largo |
---|---|---|---|---|---|---|---|---|---|---|
Código | 35020300 | 35030291 | 35030010 | 35030021 | 35020290 | 35030091 | 35020011 | 35020020 | 35030080 | 35025051 |
Latitud | 4.253917 | 4.188806 | 4.349997 | 4.310083 | 4.486528 | 4.274778 | 4.254444 | 4.196167 | 4.44116 | 4.482833 |
Longitud | -74.00269 | -73.6931 | -73.7100 | -73.648028 | -73.89041 | -73.61783 | -73.82333 | -73.77191 | -73.9363 | -74.03027 |
Altitud | 2300.0 | 191.0 | 1800.0 | 1100.0 | 1900.0 | 1280.0 | 1300.0 | 10.0 | 2100.0 | 2980.0 |
Enero | 70.0 | 393..9 | 91.0 | 229.2 | 38.4 | 508.6 | 97.6 | 217.4 | 34.0 | 67.2 |
Febrero | 69..9 | 109.9 | 30.0 | 94.6 | 16.2 | 164.5 | 41.1 | 81.2 | 12.0 | 20.7 |
Marzo | 84.4 | 132.0 | 61.0 | 167.9 | 29.0 | 164.7 | 69.7 | 86.9 | 27.1 | 56.4 |
Abril | 128.1 | 478.5 | 161.6 | 341.3 | 71.8 | 533.9 | 167.0 | 312.6 | 67.2 | 112.2 |
Mayo | 226.6 | 764.4 | 292.0 | 639.6 | 120.6 | 898.9 | 264.7 | 505.2 | 133.6 | 176.5 |
Junio | 316.1 | 922.6 | 378.3 | 927.2 | 159.1 | 1040.1 | 374.2 | 708.6 | 129.8 | 187.5 |
Julio | 358.3 | 811.5 | 398.1 | 768.3 | 123.3 | 933.1 | 354.7 | 638.9 | 99.4 | 203.2 |
Agosto | 259.0 | 742.7 | 444.2 | 726.4 | 124.9 | 953.0 | 414.1 | 647.3 | 114.4 | 209.5 |
Septiembre | 252.1 | 619.0 | 360.0 | 636.2 | 107.8 | 688.0 | 346.0 | 575.5 | 84.5 | 152.1 |
Octubre | 149.0 | 567.2 | 271.3 | 632.2 | 68.4 | 636.8 | 230.5 | 455.4 | 72.2 | 71.1 |
Noviembre | 131.6 | 598.5 | 198.8 | 590.1 | 100.6 | 789.5 | 191.0 | 326.8 | 66.6 | 114.3 |
Diciembre | 132.6 | 707.1 | 151.9 | 404.3 | 97.8 | 837.6 | 132.5 | 308.1 | 81.1 | 140.2 |
Prec Anual | 2177 | 6847 | 2839 | 6157 | 1058 | 8149 | 2683 | 4864 | 922 | 1511 |
Evapotrans piración |
1293 | 1389 | 1324 | 1422 | 889 | 1366 | 1002 | 1648 | 753 | 742 |
Humedad relativa |
84 | 89 | 90 | 90 | 74 | 79 | 80 | 89 | 62 | 61 |
La precipitación mínima media encontrada fue de 12 mm/mes en la estación de "Las Casas", mientras la máxima corresponde a 1040.1 mm/mes encontrada en la estación de Buenavista para el mes de junio. Estos valores fueron interpolados en ArcGIS usando el método de distancia inversa ponderada, en el cual se determinan los valores de celda a través de una combinación lineal ponderada de un conjunto de datos de precipitaciones (Watson y Philip, 1985), presentando el modelo de precipitaciones anual en la Figura 2. Según el análisis del clima en la cuenca del Río Guayuriba del "POMCA (Plan de Manejo y Ordenamiento de una Cuenca) Guayuriba" (UTGS, 2018) la dirección predominante del viento en 3 de las estaciones meteorológicas (Monteredondo, Servitá y Susumuco; Tabla 3) en todos los meses es desde el norte y el noreste, si bien en el segundo trimestre del año se tiene una mayor proporción de vientos que corren desde esta dirección. Se destaca también que en el trimestre de junio a agosto se incremente la proporción de viento que viene desde el oeste y el noroeste; esto puede estar asociado de cierta manera a la temporada de lluvias, pues los eventos de precipitación en el piedemonte pueden generar vientos a escala local, cerca de la superficie.
3.3. SISMICIDAD
A partir de la base de datos del Servicio Geológico de los Estados Unidos (USGS por sus siglas en inglés) y del catálogo de la Red Sismológica Nacional -RSN- del Servicio Geológico Colombiano, se compiló una base de datos de sismicidad instrumental para el período 1974 al 2020. A partir de esta última, se construyó el mapa mostrado en la Figura 3, y se calcularon distintos parámetros sísmicos como: deformación (DS), energía (ES) y levantamiento (U), con los parámetros flexurales de entrada resumidos en la Tabla 3. Cabe acotar que sólo se consideraron los sismos con profundidad menor o igual a 15 km (Bermúdez et al., 2013) con el fin de analizar el aporte de estos a la deformación vertical.
Los principales eventos sísmicos en la zona, ocurrieron en Fómeque el 18 de octubre de 1734 y en El Calvario el 24 de mayo del 2008, cuyas magnitudes fueron de 6.3 y 5.7 en la escala de Richter respectivamente. El evento más reciente, con Mw 6.0, localizado en las coordenadas 4.454° latitud Norte y 73.635° longitud Oeste, a una profundidad de 3.9 km, generó diversos deslizamientos en la vía nacional Bogotá-Villavicencio, que acabaron con la vida de seres humanos en el año 2008 durante el bloqueo de la vía (Clarín, 2008).
3.4. DIVISIÓN EN BLOQUES LITOTECTÓNICOS
En la Figura 4, se presenta la delimitación litotectónica realizada de acuerdo con las fallas presentes en la zona de estudio e información estructural obtenida del perfil geológico provisto por Parra et al. (2009). El perfil antes mencionado se construyó tomando como base las líneas sísmicas de la Agencia Nacional de Hidrocarburos (ANH), siendo posible comparar por sectores los distintos parámetros obtenidos en cada uno de los once dominios litotectónicos discriminados en esta investigación.
Posteriormente, empleando el modelo de elevación digital y la base de datos de precipitaciones de cada uno de los bloques litotectónicos, se calcularon los siguientes atributos primarios y secundarios del terreno: Relieve local a 1 y 2.3 km; integral hipsométrica (HI); capacidad de transporte de sedimentos (STI); índice topográfico de humedad (WI); índices de erosión: poder erosivo total (TSP, Total Stream Power), poder erosivo por cizalla SSP (Shear Stress Power) y poder erosivo unitario USP (Unit Stream Power); índice de gradiente longitudinal del canal (SL) e índice de empinamiento (ksn).
Las ecuaciones para cada uno de estos parámetros se resumen en la Tabla 5.Con la base de datos de sismicidad se calcularon: energía, deformación y levantamiento sísmico denotados como ES, DS y U, respectivamente (Tabla 6). Con la base de datos de isótopos cosmogénicos (Struth et al., 2015), se obtuvieron las tasas de erosión actual calculada por isotopos cosmogénicos.
Parámetro | Ecuación | Referencia | |
---|---|---|---|
Relieve local (R) | 𝑅= 𝐷𝐸𝑀−𝑚𝑖𝑛(𝐷𝐸𝑀) Calculado usando radios de 1 y 2,3 km en ArcGIS |
(2) | Montgomery y Brandon (2002) |
Integral hipsométrica (HI) |
|
(3) | Pike y Wilson (1971) |
Capacidad de transporte de sedimentos (STI) |
STI = (𝑚 + 1) × (𝐴𝑠 / 22.13)𝑚 × 𝑠𝑖𝑛(𝑆 / 0.0896)𝑛 𝐴𝑠 es el área de acumulación de flujo, y S es la pendiente. En esta investigación se asumió que m=0,4; n=1. |
(4) |
Moore y Burch (1986) Goulsbra et al. (2012) |
Índice topográfico de humedad (WI) |
𝐴𝑠 y 𝑆 se definen igual al caso anterior |
(5) | Beven y Kirkby (1979) |
Índice de erosión potencial a corto plazo (TSP, USP, SSP) |
Ap es el área de acumulación de flujo pesada por las precipitaciones (P) para cada píxel 𝑘 es el factor de erodabilidad, se asume en esta investigación k=1. Si m = n = 1, se obtiene la potencia total de flujo (TSP) Si m = 1/2 y n = 1, se obtiene la potencia por unidad del ancho del canal (USP) Si m = 1/3 y n = 2/3 se obtiene la potencia por poder de cizalla (SSP) |
(6) |
Tucker et al. (1997) Bermúdez et al. (2013) |
Índice de Gradiente longitudinal del canal (SL) |
ΔH es el cambio en la elevación de la sección, ΔL es la longitud de ésta, y L es la longitud total desde el punto medio de la sección de interés hasta el punto más alto del canal aguas arriba. |
(7) | Hack (1973) |
Índice de Empinamiento (ksn) |
ksn y Ө son el índice de empinamiento, y la concavidad, respectivamente. Este índice fue calculado usando la herramienta TopoToolBox de Matlab, para lo que es necesario determinar el área de acumulación de flujo tomando en cuenta el patrón de precipitaciones de la zona, y la pendiente a partir del modelo de elevación digital |
(8) |
Flint (1974) Wipple y Tucker, (1999) Schwanghart, y Kuhn (2010) Schwanghart and Scherler (2014) |
Parámetro | Ecuación | Referencia | |
---|---|---|---|
Energía sísmica (Es) |
𝑙𝑜𝑔10(𝐸𝑠)=𝑏𝑀𝑖+𝑎 Los parámetros a y b de esta relación se calculan mediante regresión lineal de la relación magnitud-frecuencia de los sismos. El área se dividió en 4 y 16 cuadrantes con el propósito de analizar su variabilidad espacial. M es la magnitud en escala de Richter |
(9) | Gutenberg y Richter (1954) |
Deformación Sísmica horizontal (DSH) |
a y b se derivan de la ecuación (9) Mmax es el máximo de la magnitud observada, μ es el módulo Young, ΔV representa el volumen de la corteza, Δt es el periodo de tiempo de la base de datos sismológicos |
(10) |
Braun et al. (2009) Bermúdez et al. (2013) |
Deformación sísmica vertical (𝐷𝑆v) |
𝐷𝑆𝑣=−𝐷𝑆𝐻×2.5𝑀𝑎 Se calcula durante los últimos 2.5 Ma para lograr extrapolar el levantamiento sísmico en el tiempo y hacerlos coincidir temporalmente con los datos termocronológicos |
(11) | Braun et al. (2009) |
Levantamiento sísmico (𝑢) |
hc es el espesor de la corteza tomado del modelo CRUST 1.0 mencionado anteriormente, ρc y ρm son los promedios de las densidades de la corteza y el manto, respectivamente. |
(12) |
Bassin et al. (2000); Braun et al. (2009) Laske et al. (2013), |
A partir de los datos termocronológicos se calcularon las tasas de erosión a largo plazo mediante el uso de los códigos: PECUBE (Braun, 2003; Braun et al., 2012) y Age2Edot (Brandon et al., 1998). Finalmente, para hallar relaciones significativas entre estos parámetros, y facilitar la discusión, se realizó un análisis de correlaciones empleando los coeficientes de Pearson con un valor p de 0.2 seguido de un análisis de regresión múltiple (Kendall y Stuart, 1973).
4. Resultados
4.1. TASA DE PRECIPITACIONES, ELEVACIÓN, Y RELIEVE LOCAL
El promedio anual de precipitación ha fluctuado en la zona de estudio en las últimas décadas desde 3 m/año hasta 7 m/año; localmente, los bloques surorientales tienen valores de precipitación más altos (5 ± 1 m/año) que el resto del área de estudio, siendo los bloques de Servitá y Guayuriba los que más reciben precipitación con un promedio anual de 6 m/año y 6.9 m/año, respectivamente, a diferencia de los bloques noroccidentales como Teusacá, Quetame y Chingaza que no superan el promedio de precipitación anual de 2.8 ± 0.6 m/año.
En cuanto a la distribución de elevaciones, éstas se encuentran en el rango comprendido entre 400 y 4000 msnm. Adicionalmente, el relieve de la zona se caracteriza por tener valores mayores a los 1000 m. La Figura 5 muestra los mapas de relieve calculados a 1 km y 2.3 km, respectivamente. En cuanto a los bloques litotectónicos discriminados, se observa que el bloque Guayabetal presenta el valor máximo de relieve (2342 m) usando un radio de 1 km, seguido del bloque Guatiquía. Por su parte, el promedio del relieve en toda la zona es de 1970 m, lo cual habla de una topografía abrupta en casi toda el área, donde el bloque de Servitá es el que presenta menor promedio de relieve con 1477 m.
4.2. INTEGRAL HIPSOMÉTRICA
La hipsometría del área corresponde a un paisaje muy erosionado en general, donde en la Figura 6 se observan áreas en estado de desequilibrio (Hi>0.6) en naranja y rojo.
4.3. ATRIBUTOS PRIMARIOS Y SECUNDARIOS DEL TERRENO
Los atributos del terreno se presentan en la Figura 7. En particular, el índice de capacidad de transporte de sedimentos (STI) se visualiza en la Figura 7b. El índice de humedad topográfica (WI) se presenta en la Figura 7c, en la cual los sitios muy planos coinciden con las partes más húmedas del paisaje, tal como el sureste del área. El valor máximo del índice topográfico de humedad (WI) se presenta en el bloque de Servitá con un valor de 3.17 debido a su bajo relieve, el promedio en la zona es de 2.2 para el WI; mientras el índice STI presenta valores negativos y positivos asociado con zonas de sedimentación y erosión respectivamente, con valores positivos mayores para el bloque de Guayuriba de 0.045, lo que sugiere que este bloque, el cual se encuentra en el sector Chirajara, tiene un alto potencial erosivo y coincide con importantes deslizamientos en la zona.
Los índices de erosión potencial a corto plazo (EI), como el poder erosivo por cizalla (SSP, Shear Stress Power), se visualizan en la Figura 7d, el poder erosivo unitario (USP, Unit Stream Power) en la Figura 7e y el poder erosivo total (TSP, Total Stream Power) en la Figura 7f. El índice TSP presenta un promedio muy alto de 2.23E7 W/m2 en el bloque Guayuriba, comparado con el promedio de todos los bloques de 3.6E6 W/m2. Este comportamiento persiste con el índice USP, donde el bloque Guayuriba presenta el valor promedio de 685 que sobresale entre el promedio de todos los bloques de 275 W/m2. También en el índice SSP el bloque Guayuriba exhibe un valor alto de 54 W/m2, siendo la mayor magnitud de este índice en los bloques donde se tiene una media de 32 W/m2, siendo el bloque Guayuriba el que más presenta potencial erosivo a corto plazo, lo cual se debe a la confluencia del Río Negro y el Río Blanco en este bloque, por lo cual se concentra el alto poder erosivo en el Río Guayuriba. Después del bloque Guayuriba, el bloque Guatiquía exhibe un alto potencial erosivo en el índice USP y TSP con valores de 290 y 2.85E6 W/m2, respectivamente.
4.4. ÍNDICE DEL GRADIENTE LONGITUDINAL DEL CANAL (SL) Y EMPINAMIENTO (KSN)
El índice SL se presenta en la Figura 8 para diferentes segmentos de perfiles longitudinales del Río Guayuriba. Se muestra además el intervalo de precipitación y la litología presente en cada segmento. El índice SL presenta su valor máximo de 1,321 en el bloque Guayuriba seguido del bloque de Río Blanco y Guayabetal, con 1.112 y 939 respectivamente, siendo el promedio del índice SL en todos los bloques de 709.
En la Figura 9, se presenta el mapa con valores del índice ksn, Se observa en esta figura los sitios de muestreo con edades termocronológicas existentes. El índice ksn presenta una media de 243 para toda la zona siendo un valor alto de incisión, mientras el valor mínimo lo presenta el bloque de Servitá con 66.
4.5. PARÁMETROS SISMOLÓGICOS (ENERGÍA, DEFORMACIÓN Y LEVANTAMIENTO)
En los bloques litotectónicos de Naranjal y Guayabetal se registra la mayor actividad sísmica de la zona. En la Figura 10, se muestran los eventos sísmicos desde 1974 hasta el 2020 y los parámetros sísmicos calculados a partir de esa base de datos. En términos de la subdivisión en 4 y 16 cuadrantes, respectivamente, no introduce diferencias significativas en los valores b calculados. Sin embargo, debido a que la subdivisión está también relacionada con la escala de los fenómenos tectónicos observados, se decidió a lo largo de la investigación mantener ambas subdivisiones. La mayor deformación sísmica y levantamiento sísmico ocurre en la zona delimitada entre las fallas Naranjal y Quebrada Grande (Figuras 10a y 10b). La energía sísmica liberada es mayor en el norte de la zona de estudio (Figuras 10c y 10d). La energía sísmica varía entre 1.18E5 J a 2.41E6 J para todos los bloques, y el valor promedio es 5E5 J. El valor máximo de 2.41E6 J se encuentra en el bloque de Guayabetal, concentrándose la energía sísmica en este sector. Le siguen los bloques Quetame y Servitá con energía sísmica promedio de 5.545 J y 4.21E5 J, respectivamente. Se observa que estos bloques presentan los valores máximos de energía sísmica debido a que sus márgenes corresponden con las principales estructuras del área como las fallas de Servitá, Quetame y Naranjal.
El resto de los bloques presentan energía sísmica baja de 2.5E5 ± 0.6E5 J, comparado con los bloques de Guayabetal, Quetame y Servitá.
4.6. ÍNDICE DE EROSIÓN ACTUAL POR DATACIÓN DE ISÓTOPOS COSMOGÉNICOS (10BE)
En las terrazas del Río Guayuriba, existen 7 estimaciones de tasas de incisión (01, 03-08 de la Figura 11) previamente determinadas por Struth et al. (2015) a partir de mediciones de Berilio-10, las cuales fueron interpoladas en esta investigación, obteniéndose tasas de ~ 75 mm/ka (muestras 01, 03, 04, 08; Figura 11), mientras que las tasas de erosión aumentan significativamente a valores de 100, 479 y 670 mm/ka (07, 05 y 06) hacia el sur de la zona de estudio, respectivamente.
4.7. TASA DE EROSIÓN A LARGO PLAZO
La tasa de erosión a largo plazo promedio para toda el área es de 2.37 mm/año, siendo una tasa muy alta de erosión, donde los valores más altos se concentran al oriente, en el bloque Guatiquía, con una tasa de 3.47 mm/año, seguido por los bloques Guayuriba además Río Blanco con una tasa de 2.97 mm/año y 2.63 mm/ año, respectivamente.
4.8. RESULTADOS DEL MODELO DIRECTO TERMOCINEMÁTICO
Las predicciones del modelo termocinemático directo de menor misfit (3.327) son mostradas en la Figura 12. Los parámetros de entrada como se mencionó anteriormente, fueron temperatura basal de 400°C y la producción de calor de 1.2 mW/m2. Se consideró además una topografía en equilibrio para simplificar el tiempo de cómputo.
La Figura 13 resume la variación de las tasas de exhumación predichas para el modelo de menor misfit. Considerando tres fases de exhumación, se obtuvo el primer pulso entre 40 Ma y 25 Ma con una tasa de exhumación de 0.5 km/Ma. Luego en el segundo pulso entre 25 Ma y 15 Ma, la tasa de exhumación decrece a 0.1 km/Ma, y a partir de los 15 Ma hasta el presente, la tasa de exhumación aumenta abruptamente a 2 km/Ma. Las edades predichas se convirtieron a tasas de exhumación a largo plazo, las cuales están resumidas en la tabla 6.
Se presenta en la Figura 14, la comparación de perfiles edad-elevación sintéticos versus observados. Para la Figura 14, se realizó una regresión lineal con el propósito de obtener estimaciones de tasas de exhumación de los dos sistemas termocronológicos utilizados, obteniéndose valores de coeficientes de determinación: R2=0.9194 y R2=0.9856, para las edades AFT y ZFT respectivamente. Adicionalmente, se obtuvieron las siguientes expresiones:
La tasa de erosión derivada de las relaciones edad AFT versus elevación es 2.6 km/Ma; mientras que la tasa de erosión estimada de las relaciones edad ZFT versus elevación es de 1 km/Ma. Tal diferencia entre ambas estimaciones sugiere que la velocidad de exhumación se ha incrementado desde los últimos 5 Ma. La Figura 15a muestra las edades de AFT predichas por el código PECUBE para toda la zona de estudio, y en más detalle las del Bloque Guayabetal en la Figura 15b, en el cual coinciden las mayores edades de AFT predichas con la deformación sísmica.
Finalmente, para facilitar la discusión, se presentan en la Tabla 7 los valores promedios de los parámetros geomorfológicos y sísmicos de cada uno de los bloques litotectónicos de la Figura 4.
Bloque | ES4 |
DS4 ×10-17 |
U4 | ES16 |
DS16 ×10-17 |
U16 | P | STI | WI | TSP | USP | SSP | SL | ksn | HI | R1 | erAFT | erCN |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Teusacá | 6084.1 | 0.9 | 7.9 | 2178.7 | 3.0 | 19.6 | 2.0 | 2.2E-02 | 1.9 | 286046.3 | 108.4 | 18.0 | 465.7 | 259.9 | 0.5 | 1589.3 | 1.5 | 103.9 |
Quetame | 6425.6 | 1.6 | 69.2 | 4688.1 | 2.8 | 18.1 | 2.9 | -2.2E-02 | 2.0 | 672530.1 | 165.7 | 24.1 | 503.7 | 211.1 | 0.5 | 1718.5 | 1.9 | 135.4 |
Q. San Martín | 3301.7 | 1.2 | 61.3 | 1590.9 | 0.01 | 0.2 | 4.1 | 3.7E-02 | 1.9 | 2347313.8 | 318.4 | 36.3 | 892.0 | 284.4 | 0.4 | 2212.3 | 2.2 | 346.1 |
R. Chiquito | 1894.7 | 1.2 | 10.6 | 746.6 | 0.03 | 0.2 | 3.1 | -2.2E-02 | 1.7 | 1076527.4 | 249.9 | 32.8 | 617.3 | 295.2 | 0.4 | 2224.8 | 2.0 | 263.5 |
Chingaza | 8814.3 | 4.0 | 7.9 | 7121.6 | 4.2 | 27.5 | 3.7 | 1.2E-02 | 2.0 | 305113.8 | 183.7 | 26.1 | 339.9 | 287.9 | 0.6 | 2159.4 | 2.4 | 145.3 |
Naranjal | 3438.1 | 4.7 | 7.8 | 2317.8 | 4.8 | 31.2 | 4.1 | 8.4E-04 | 2.1 | 2534552.1 | 266.6 | 32.7 | 794.8 | 281.3 | 0.4 | 2074.0 | 2.3 | 208.3 |
Guayabetal | 15527.2 | 1.1 | 13.1 | 6906.3 | 4.6 | 30.1 | 5.0 | -6.9E-03 | 2.2 | 1634755.9 | 278.8 | 34.3 | 939.4 | 273.6 | 0.5 | 2342.9 | 2.5 | 317.1 |
Guatiquía | 7025.4 | 9.5 | 75.4 | 3997.6 | 3.3 | 21.6 | 5.2 | -6.1E-02 | 2.4 | 2847799.5 | 290.7 | 33.8 | 426.3 | 258.3 | 0.4 | 2272.9 | 3.4 | 265.4 |
Servitá | 4032.7 | 0.1 | 30.1 | 1253.9 | 0.04 | 0.3 | 6.0 | -4.5E-02 | 3.2 | 2677215.6 | 177.2 | 22.3 | 387.1 | 66.3 | 0.3 | 1467.4 | 2.2 | 397.4 |
Guayuriba | 5691.6 | 0.1 | 26.2 | 1794.4 | 0.04 | 0.2 | 6.9 | 5.4E-02 | 2.5 | 22311850.4 | 685.1 | 54.3 | 1321.3 | 219.9 | 0.4 | 1733.6 | 3.0 | 567.4 |
R. Blanco | 2649.7 | 2.0 | 70.8 | 1239.4 | 0.04 | 0.2 | 5.0 | -2.3E-02 | 2.3 | 2965740.6 | 306.1 | 35.4 | 1112.7 | 244.0 | 0.4 | 1883.7 | 2.6 | 527.6 |
5. Discusión
Los Himalayas, Los Andes y los orógenos de Nueva Zelanda constituyen regiones claves donde las interacciones entre tectónica y aspectos climáticos, en particular el patrón de precipitaciones, han sido documentadas, las cuales pueden influenciar de manera fundamental la evolución individual del relieve (Koons et al., 2002; Sobel et al., 2003). Para el área de estudio estas relaciones son desconocidas. Con el propósito de facilitar la discusión, trataremos de dar respuestas a preguntas claves como, por ejemplo: ¿Cómo influye la sismicidad en los bloques con más deslizamientos (bloques surorientales)? ¿Qué bloques presentan correlaciones altas entre relieve y clima o relieve y tectónica? ¿Existirá acople entre los procesos climáticos y tectónicos?, si esto es así ¿En qué escala de tiempo se presentan tales interacciones?
Con los resultados de la Tabla 7, se realizó un análisis de correlaciones de Pearson, presentado bajo la forma de correlogramas y mostrado en la Figura 16. En la Figura 16a, se incorporan todos los bloques considerados, mientras que en la Figura 16b se consideraron sólo los bloques litotectónicos con mayor actividad sísmica cercanos a la zona de deslizamiento (Naranjal, Guayabetal, Guatiquía, Servitá, Guayuriba y R. Blanco). El incremento en los valores de la excentricidad (relaciones entre eje mayor y eje menor) de la elipse sugiere una mayor dispersión entre los pares de variables que están siendo comparados. En ambos correlogramas, las relaciones entre todos los parámetros sísmicos ES4, DS4, ES16, DS16, U4 y U16 son esperadas ya que son calculadas de la misma forma, por lo que son descartadas. Análogamente, las correlaciones entre las variables P, TSP, SSP, USP, WI, y STI tampoco deben ser consideradas porque todos los índices de erosión a corto plazo fueron pesados por las precipitaciones.
La Figura 16a muestra importantes relaciones (r = 0.7) entre la deformación sísmica calculada para toda el área dividida en cuatro cuadrantes DS4 y las tasas de erosión erAFT y erCN. Sin embargo, cuando se refina la escala incluyendo procesos locales, no se observan correlaciones entre DS16 y las tasas de erosión antes mencionadas. En contraste, correlaciones entre el patrón de precipitaciones P y las tasas de erosión erAFT y erCN soportaría los efectos de un control climático; la correlación entre P y erCN (r = 0.8) es ligeramente mayor que la de P y erAFT (r = 0.7). Todos los índices de erosión a corto plazo TSP, USP y SSP correlacionan con ambas tasas de erosión (r ≥ 0.5). Sin embargo, los mayores valores de r se observan en las relaciones entre TSP, USP, SSP y erCN (r ≥ 0.7). La correlación entre ksn y R1 (r = 0.7), sugiere que el relieve actual está siendo afectado por los procesos de incisión de ríos; pero a su vez, la pérdida de correlación entre R1 y las tasas de erosión erAFT y erCN (r ≤ 0.4), al igual que R1 y P (r <0), sugieren que al considerar todos los bloques no se puede discriminar controles tectónicos o climáticos sobre el relieve.
En contraste, la Figura 16b, para los bloques con mayor actividad sísmica, refleja una importante relación entre R1 y ES4 (r = 0.6), R1 y ES16 (r = 0.8), R1 y DS16 (r = 0.8) y U16 (r = 0.8). Para estas zonas, la correlación entre relieve R1 y P es negativa (r = -0.6). Con respecto a las tasas de erosión, se observa que localmente la tectónica no pareciera controlar las tasas de erosión obtenidas de las edades AFT ya que no existe correlación entre ES16, DS16 y U16 con erAFT. Sin embargo, la correlación entre U4 y erAFT (r = 0.6) sugiere que globalmente la tectónica tiene una influencia importante sobre las tasas de erosión derivada de las edades AFT y sobre el relieve. Observando las correlaciones entre P y erCN (r = 0.7), y los índices TSP con erCN (r=0.6) y USP versus erCN (r = 0.6) sugiere un desacople entre aspectos climáticos y tectónicos. A partir de cierto momento, la precipitación pareciera controlar las tasas de erosión más recientes.
El análisis realizado sugiere que la morfología actual del área de estudio no es afectada uniformemente por la tectónica y las precipitaciones. Se pueden distinguir dos áreas distintas de control, la primera localizada en al noroccidente (Bloques de Teusacá, Chingaza y Quetame; Figura 1b) con tasas bajas de erosión y actividad tectónica donde la evolución del paisaje es más pasiva, y el relieve correlaciona fuertemente con la precipitación, siendo ésta el principal factor que controla esta área noroccidental. Por otra parte, la relación de la precipitación con las tasas de erosión a corto y a largo plazo en toda la zona, sugieren el control del clima en el rejuvenecimiento de la superficie.
La segunda zona se encuentra en el macizo de Quetame y en el piedemonte llanero, representado por los bloques surorientales (Guatiquía, Guayabetal, Guayuriba, Naranjal, Servitá y Río Blanco). En estas zonas es donde ocurrió el deslizamiento que motivó el presente estudio y posiblemente exista una concentración de actividad sísmica a lo largo del tiempo donde se evidencia el aporte tectónico para la formación de estos bloques, sugerido por la relación del relieve con los parámetros sísmicos (DS16 y U16). Sin embargo, esta observación debe ser tomada con cautela ya que la ventana temporal de observación es menor a 50 años. Aunque el control tectónico es significativo, recientemente el clima ha rejuvenecido el paisaje, posiblemente desde el Pleistoceno. Una vez que se produce la rápida exhumación (~2 km/Ma) del macizo de Quetame (Parra et al., 2009), se genera una barrera orográfica en la cual los vientos húmedos provenientes del Amazonas se concentran, lo cual explicaría por qué los termocronómetros muy someros, como en el caso de los isótopos cosmogénicos por 10Be (erCN), correlacionan con la variable precipitación.
Los análisis anteriormente realizados demuestran que no es una tarea sencilla discriminar cuál de los dos procesos, clima, por un lado, en este caso representado por el patrón de precipitaciones actuales, y la tectónica por otro, controlarían el paisaje actual de la zona de estudio. Con el propósito de discriminar la importancia de estos procesos sobre el relieve actual y sobre las tasas de erosión, se realizó un análisis de regresión múltiple donde las variables respuesta o dependiente son: el relieve (R) y las tasas de erosión (erAFT y erCN), mientras que las variables predictoras o independientes serían el patrón de precipitaciones, los parámetros morfométricos y los parámetros sismológicos. Así, en la Figura 17, se muestra el aporte de cada variable predictora sobre las dos variables respuestas consideradas.
Estos análisis permiten discriminar a los procesos de incisión fluviales SSP como una variable predictora importante para el relieve, seguido del índice ksn y la deformación sísmica DS4, la cual involucra la tectónica global del área. Se observa, por ejemplo, que la deformación sísmica a una escala local no tiene un impacto significativo sobre el relieve considerado. En contraste, el comportamiento de las tasas de erosión erAFT y erCN parecieran ser explicadas por el patrón de precipitaciones global. Existen diferencias en el comportamiento de las otras variables para estas tasas de erosión, mientras TSP, ksn y DS16 tienen una influencia positiva sobre las tasas de erosión erAFT y sobre las tasas de erosión derivada de los isótopos cosmogénicos. Por su parte, la deformación sísmica no pareciera afectar a nivel regional ni local estas últimas tasas.
6. Conclusiones
El modelado directo termocinemático 3D permitió establecer la contrastante historia de exhumación de este sector de la Cordillera Oriental, al menos tres fases de exhumación son necesarias para reproducir el patrón de edades termocronológicas reportadas en el área. Durante estos pulsos se podrían producir diferencias en cuanto al mecanismo controlador de la exhumación, así los efectos del clima y la tectónica podrían acentuarse o disminuir en distintos intervalos de tiempo. Principalmente, en el primer pulso entre 40 Ma y 25 Ma se observó una tasa de exhumación lenta de 0.5 km/Ma posiblemente controlada por tectónica, seguido de una disminución en la tasa de exhumación entre 25 Ma a 15 Ma, lo que indicaría un intervalo de quietud tectónica, finalmente, y a partir de los 15 Ma al presente las tasas de exhumación se han venido incrementando a valores de 2 km/Ma, donde este incremento posiblemente se acentuó en el Pleistoceno, y tal vez tenga relación al incremento en la interacción del clima y la tectónica para el área de estudio.
El principal factor dominante del paisaje hoy en día es la incisión fluvial primariamente, impulsada por el incremento en la pendiente regional del orógeno por acumulación progresiva del acortamiento y engrosamiento cortical. A nivel general, la tectónica ejercería un control importante sobre el relieve, pero su efecto disminuiría sobre las tasas de erosión a largo y corto plazo, tomando un importante rol las precipitaciones para explicar las tasas de erosión.
La actividad tectónica pudiese estar creando topografía en zonas de alta sismicidad (Bloques Guayuriba, Río Blanco y Guayabetal), y en conjunción con la existencia de precipitaciones de valores mayores a 4 m/año han incrementado la incisión fluvial del Rio Blanco y el Río Negro, ambos tributarios del Río Guayuriba.
La morfología actual del área de estudio no es afectada de la misma forma por la tectónica y las precipitaciones, donde se pueden distinguir dos áreas distintas: una primera localizada en el noroccidente caracterizada por tasas bajas de erosión y de actividad tectónica, donde la evolución del paisaje es más pasiva y donde el relieve correlaciona fuertemente con la precipitación, siendo éste el principal factor que controla los bloques noroccidentales: Teusacá, Chingaza y Quetame. La segunda zona se encuentra desde el macizo de Quetame hasta el piedemonte oriental hacia el sureste. Esta zona muestra mayor actividad sísmica en tiempos instrumentales, donde están presentes las principales estructuras tectónicas, como las fallas de Naranjal, Quetame, San Juanito y Servitá con alto potencial sismogénico. Las fuertes relaciones de los parámetros sismológicos con el relieve local y el índice de incisión ksn sugieren que la tectónica controló la formación de esta zona por medio del acortamiento y engrosamiento del orógeno, el cual incrementó el relieve y la pendiente topográfica, aumentando el poder erosivo de incisión (ksn) hasta el Pleistoceno, cuando el clima rejuveneció los bloques surorientales.
La incisión fluvial se concentra además de los valles en forma de V (cuenca del río Guayuriba), localmente en las zonas de deslizamientos, donde existe una mayor tasa de precipitación, lo que aumenta la criticidad, aunado a una importante deformación sísmica (como la observada en la falla de Servitá).