Boletín de la Sociedad Geológica Mexicana Volumen 73, núm. 2, A120121, 2021 http://dx.doi.org/10.18268/BSGM2021v73n2a120121

Discriminación de la relación precipitación-tectónica como agentes modeladores del paisaje en los alrededores del Municipio Guayabetal, Cordillera Oriental de Colombia

Discrimination of the precipitation-tectonic relationship as landscape shapers in the vicinities of the Guayabetal municipality, Eastern Cordillera of Colombia

Daniel R. Hernández-Chaparro^1,*, Mauricio A. Bermúdez¹, Greg Hoke², Helbert García-Delgado³,

Silvia Machuca⁴

¹ Escuela de Ingeniería Geológica, Universidad Pedagógica y Tecnológica de Colombia, Avenida Central del Norte 39-115, Sede Central Tunja,Boyacá, Colombia.

² Department of Earth Sciences, Syracuse University, Syracuse, 13244, Nueva York, Estados Unidos de América.

³ Servicio Geológico Colombiano, Diagonal 53 No. 34 – 53, Bogotá, Colombia.

⁴ Pontificia Universidad Javeriana, Cra 7 No. 40-62, Bogotá, Colombia.

* Autor para correspondencia: (D.R. Hernández-Chaparro) This email address is being protected from spambots. You need JavaScript enabled to view it.

Cómo citar este artículo:

Hernández-Chaparro, D. R., Bermúdez, M. A., Hoke, G., García-Delgado, H. y Machuca, S., 2021, Discriminación de la relación precipitación-tectónica como agentes modeladores del paisaje en los alrededores del Municipio Guayabetal, Cordillera Oriental de Colombia: Boletín de la Sociedad Geológica Mexicana, 73 (2), A120121. http://dx.doi.org/10.18268/BSGM2021v73n2a120121

RESUMEN

En el Municipio Guayabetal, flanco oriental de la Cordillera Oriental de Colombia se han presentado importantes fenómenos de remoción en masa que afectan las vías de comunicación (km 58 vía Bogotá-Villavicencio). En esta zona, significativas tasas de precipitación y de exhumación son observadas, las cuales pudieron desencadenar estos procesos. En la presente investigación, se analizan las interacciones entre precipitación y tectónica como procesos controladores del relieve actual. Se compilaron edades termocronológicas existentes, y se generó un modelo termocinemático directo 3D para estimar los pulsos y las tasas de exhumación. Los resultados de este último sugieren 3 pulsos de exhumación: el primero entre 40 Ma - 25 Ma a una tasa de exhumación de 0.5 km/Ma, seguido de un pulso entre 25 Ma - 15 Ma con una tasa de 0.1 km/Ma, y finalmente, desde 15 Ma al presente tasas de ~2 km/Ma. Adicionalmente, modelos de elevación digital, datos de precipitación y sismológicos fueron utilizados para estimar atributos primarios y secundarios del terreno, así como: deformación sísmica, energía sísmica y levantamiento vertical por deformación sísmica. Las distintas variables fueron comparadas estadísticamente. Se concluye que el paisaje actual de la zona no es afectado uniformemente por la tectónica y las precipitaciones. En el noroccidente del área existen bajas tasas de erosión y actividad tectónica, así la evolución del paisaje es más pasiva, y el relieve es controlado por el patrón de precipitaciones. En contraste, hacia el suroriente específicamente entre el macizo de Quetame y el piedemonte llanero, el relieve es controlado por la tectónica presente en el área. El efecto antrópico con las modificaciones al paisaje introducidas por la creación de una vía nacional, aunque es importante, no fue analizado en la presente investigación.

Palabras clave: Termocronología, modelamiento numérico, evolución del paisaje, Cordillera Oriental de Colombia, estabilidad de taludes, remoción en masa.

ABSTRACT

At the Guayabetal Municipality, on the eastern flank of Colombia’s Eastern Cordillera, there have been significant mass removals affecting communication routes (km 58 Bogotá-Villavicencio road). In this zone, significant rates of precipitation and exhumation are observed, which could trigger these phenomena. In the present investigation, the feedbacks between precipitation and tectonics are investigated as controlling processes of the current relief. Existing thermochronological ages were compiled, and a direct 3D thermokinematic model was generated to estimate pulses and exhumation rates. The results of this last one suggests 3 exhumation pulses: the first one between 40 Ma - 25 Ma at an exhumation rate of 0.5 km/Ma, followed by a pulse between 25 Ma - 15 Ma with a rate of 0.1 km/Ma, and finally, from 15 Ma to the present rates of ~2 km/Ma. Additionally, rainfall and seismological data were used to estimate primary and secondary attributes of the terrain, as well as: seismic deformation, seismic energy and seismic uplift. The different variables were statistically compared. It is concluded that the current landscape of the area is not uniformly affected by tectonics and rainfall patterns. In the northwest of the area there are low rates of erosion and tectonic activity, so the evolution of the landscape is more passive, and the relief is controlled by the pattern of precipitation. In contrast, to the southeast, specifically between the Quetame massif and the piedmont plains, the relief is controlled by the tectonics present in the area. The anthropic effect on modifications to the landscape introduced by the creation of a national road, although it is important, was not analyzed in the present research.

Keywords: Thermochronology, thermal-kinematic modeling, landscape evolution, Eastern Cordillera of Colombia, slope stability, Mass- wasting processes.

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.

Figura 1. A) Entorno geodinámico y topografía simplificada de los Andes del Norte, Modificado de Taboada et al. (2000), Unidades litológicas descritas en la Tabla 1. B) Mapa geológico (Gómez-Tapias et al., 2015) que comprende la vía Bogotá-Villavicencio y distribución de edades termocronológicas existentes (Mora et al., 2008; Parra et al., 2009) además del deslizamiento del Km 58 objetivo marcado por triángulo rojo, los puntos azules y morados corresponden con la ubicación de muestras previamente fechadas por los métodos de trazas de fisión en apatitos (AFT) y circones (ZFT), respectivamente; C) Perfil transversal de la zona de estudio (Modificado de Mora et al., 2008).

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 concentrará 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 h_c 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.

Tabla 1. Principales Unidades Litológicas aflorantes en el área de estudio.

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).

Tabla 2. Edades termocronológicas existentes en la zona de estudio: Edades por trazas de fisión en apatito (AFT) y circón (ZFT) tomadas de Mora et al. (2008) y Parra et al. (2009). En la En la Figura 1, se ubica cada muestra por su código en el área de estudio.

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 km²My^-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.

Tabla 3. Parámetros de entrada fijos considerados en los modelos.

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, t_obs y t_cal 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. 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.

Figura 2.Mapa de precipitaciones promedios (m/año) generado mediante interpolación tipo krigging de Arcmap 10.8 a partir de los datos del IDEAM para el periodo 1990 a 2019. En esta figura, los cuadros en azul claro corresponden a las estaciones climáticas de donde se extrajeron los datos, sobre este mismo se han colocado las principales fallas geológicas de las zonas de estudio (Tomadas de Gómez-Tapias et al., 2015). La ubicación de las estaciones puede ser detallada en la Tabla 4.

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.

Figura 3.Mapa de sismicidad del sector estudiado, basado en los catálogos del Servicio Geológico de los Estados Unidos en círculos verdes (USGS por sus siglas en inglés), de sismicidad histórica, e instrumental de la red nacional de sismos del Servicio Geológico Colombiano entre 1974 y 2020 en círculos huecos, reportados sobre mapa estructural del SGC (Gómez-Tapias et al., 2015). Igualmente reporta la probabilidad sísmica en escala de colores, realizada por cokriging interpolando magnitud y profundidad en Arcmap 10.8.

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.

Figura 4.Delimitación de bloques litotectónicos propuestos entre las principales fallas, tomadas del mapa geológico-estructural del SGC (Gómez-Tapias et al., 2015), donde la línea roja corresponde a la Vía Bogotá a Villavicencio S-E y los triángulos corresponden a los principales fenómenos de remoción en masa tomados del catálogo electrónico del Sistema de Información de Movimientos en Masa (SIMMA, 2020).

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.

Tabla 5. Variables morfométricas, atributos primarios y secundarios del terreno calculados para la zona de estudio, a partir del modelo de elevación digital y la base de datos de precipitaciones.

Tabla 6. Parámetros sísmicos calculados para la zona de estudio a partir de la base de datos de eventos sismológicos reportados entre 1974 y 2020.

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.

Figura 5.Mapas de relieve relativo calculado para A) un radio de 1 km y para B) con un radio de 2.3 km, respectivamente. Se reportan sobre éstos los deslizamientos o fenómenos de remoción en masa descargados del catálogo electrónico del Sistema de Información de Movimientos en Masa (SIMMA, 2020). Sobre estos mapas se han reportado las principales fallas de la zona de estudio (Tomadas de Gómez-Tapias et al., 2015).

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.

Figura 6. Mapa de valores de la integral hipsométrica en la zona de estudio, realizados con la herramienta TAK (Forte et al., 2019). En este mapa los símbolos morados corresponden con las estaciones meteorológicas del IDEAM. Sobre este mapa se han reportado las principales fallas de la zona de estudio (Tomadas de Gómez-Tapias et al., 2015).

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.

Figura 7. (A) Patrón de precipitación promedio (m/año) para el período 1990-2019. Los cuadros azul claro corresponden a estaciones climáticas del IDEAM; (B) Índice de transporte de sedimentos (STI); (C) Índice Topográfico de Humedad (WI); (D) Índice de erosión a corto plazo por poder de cizalla (SSP) (Wm–2); (E) Índice de erosión a corto plazo unitario por ancho del canal (USP) (Wm–2); (F) Índice de erosión a corto plazo total (TSP) (Wm–2). Sobre cada uno de los mapas se muestran las principales fallas de la zona de estudio las cuales fueron tomadas de Gómez-Tapias et al. (2015).

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/m² en el bloque Guayuriba, comparado con el promedio de todos los bloques de 3.6E6 W/m². 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/m². También en el índice SSP, el bloque Guayuriba exhibe un valor alto de 54 W/m², siendo la mayor magnitud de este índice en los bloques donde se tiene una media de 32 W/m², 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/m², 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.

Figura 8.Perfil del flujo para el río Guayuriba con valores del índice SL para cada segmento entre dos knickpoints, en el cual se interpretó visualmente colocando en yuxtaposición litologías, presencia de fallas y valores de precipitación media anual. Para la elaboración de esta figura se consideró el cambio de las tasas de precipitación con la altura y la distribución de fallas existentes en la zona, las cuales fueron tomadas de Gómez-Tapias et al. (2015).

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.

Figura 9. Mapa de índice de empinamiento (ksn). En este mapa se muestran las principales fallas de la zona de estudio las cuales fueron tomadas de Gómez-Tapias et al. (2015).

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.18E⁵ J a 2.41E⁶ J para todos los bloques, y el valor promedio es 5E⁵ J. El valor máximo de 2.41E⁶ 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.54⁵ J y 4.21E⁵ 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.

Figura 10.Parámetros sísmicos obtenidos al dividir la zona de estudio en 4 y 16 cuadrantes, respectivamente: (A) Deformación sísmica (DS4) y (B) Deformación sísmica (DS16). (C) Energía sísmica (ES4) y (D) Energía sísmica (ES16). (E) Levantamiento sísmico (LS4) y (F) Levantamiento sísmico (LS16), respectivamente. Se muestran las principales fallas geológicas de las zonas de estudio, las cuales fueron tomadas de Gómez-Tapias et al. (2015).

El resto de los bloques presentan energía sísmica baja de 2.5E⁵ ± 0.6E⁵ 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 (¹⁰BE)

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.

Figura 11.Distribución de tasas de erosión derivadas a partir de isótopos cosmogénicos (10Be), donde el triángulo rojo-amarillo corresponde al deslizamiento del km 58. Los números enteros corresponden a muestras en terrazas fluviales de distintas cuencas de la cordillera oriental (Modificado de Struth et al., 2015).

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/m². Se consideró además una topografía en equilibrio para simplificar el tiempo de cómputo.

Figura 12. Edades AFT predichas por el modelo directo termocinemático 3D de menor misfit (3.327).

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.

Figura 13.El gráfico ilustra las variaciones de las tasas de exhumación en el tiempo para el modelo de menor misfit (3.327).

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: R²=0.9194 y R²=0.9856, para las edades AFT y ZFT, respectivamente. Adicionalmente, se obtuvieron las siguientes expresiones:

Figura 14. Comparación de relaciones edad-elevación observados (azul: AFT y naranja: ZFT) y predicho (cuadrados: AFT y círculos: ZFT) para el modelo de menor misfit (3.327).

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.

Figura 15. A) Edades AFT predichas por PECUBE para la zona de estudio. Las líneas de color son las fallas en el área: FSJ= Falla de San Juanito, FN= Falla de Naranjal, FQ= Falla de Quetame, FRB= Falla de Río Blanco, FS= Falla de Servitá, tomadas de Gómez-Tapias et al. (2015). B) Edades AFT predichas en el Bloque de Guayabetal; C) Perfil transversal de la zona de Estudio (Modificado de Mora et al., 2008).

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.

Tabla 7. Valores promedios de control potencial para cada bloque litotectónico propuesto: ES = Energía Sísmica (Joules); DS = Tasa de deformación sísmica (s-1); U = Levantamiento sísmico (m); P = Precipitación (m/año) para el período 1990-2019; STI = Índice de transporte de sedimentos; WI = Índice de humedad; TSP = Poder de flujo total (Wm–2); USP = Poder de flujo unitario (Wm–2); SSP = Poder de flujo por cizalla (Wm–2); SL = Índice del gradiente longitudinal del canal (m); ksn = Índice de inclinación; HI = Integral Hipsométrica; R1 = Relieve radio 1-km (m); erAFT = Erosión a largo plazo observada para datos AFT, er-CN = Erosión por datación de isótopos cosmogénicos (mm/ka).

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.

Figura 16.Correlogramas (A) para todos los bloques litotectónicos considerados, (B) bloques con mayor actividad sísmica cercanos a la zona de deslizamientos.

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 ¹⁰Be (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.

Figura 17. Resultados del análisis de regresión múltiple para todos los bloques litotectónicos considerados. En este análisis se toma como variable respuesta: (A) Relieve local, (B) tasas de erosión a largo plazo derivadas de los datos AFT (erAFT) y (C) tasas de erosión derivadas de los isótopos cosmogénicos (erCN), mientras que las variables predictoras (independientes) son las que se sitúan en el eje “y”.

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, 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á).

Agradecimientos

A la Universidad Pedagógica y Tecnológica de Colombia (UPTC), Sogamoso Colombia, por el apoyo brindado durante el desarrollo de la presente investigación. Especialmente por el financiamiento proporcionado a través del proyecto DIN-SGI-3104 de la Universidad Pedagógica y Tecnológica de Colombia (UPTC), y por el financiamiento proporcionado a través de las convocatorias de movilidad internacional que permitieron la estancia del Dr. Greg Hoke en la UPTC durante las campañas de muestreo de terrazas cercanas a la zona de estudio. Agradecemos a Andrea M.J. Coronato del CADIC-CONICET Argentina, Edier Aristizábal de la Universidad Nacional de Colombia, y a los editores Laura Perucca y Franck Audemard por la revisión detallada del manuscrito, comentarios y sugerencias que ayudaron a mejorar la calidad del mismo.

Referencias

Audemard, F.A., 1999, Morpho-structural expression of active thrust fault systems in humid tropical foothills of Colombia and Venezuela: Zeitschrift fuer Geomorphologie 118, 1 – 18.

Bassin, C., Laske, G., Masters, G., 2000, The Current Limits of Resolution for Surface Wave Tomography in North America: EOS Trans American Geophysical Union, 81, F897.

Beltrán, A., 2005, Cartografía geológica y análisis estructural de las unidades pre y eo-cretáceas del flanco occidental del macizo de Quetame, área del municipio de Gutiérrez, Cundinamarca, Flanco Oriental de la Cordillera Oriental de los Andes colombianos: Trabajo de grado. Universidad Nacional de Colombia.

Bermúdez, M., van der Beek, P., Bernet, M., 2011a, Asynchronous Miocene–Pliocene exhumation of the central Venezuelan Andes: Geology, 39, 139-142. http://dx.doi.org/10. 1130/G31582.1

Bermúdez, M., Van der Beek, P., Bernet, M, 2011b, Exhumación Diacrónica durante el Mioceno-Plioceno de Los Andes de Mérida: Revista Tribuna del Investigador, Universidad Central de Venezuela, 12, 20 - 23.

Bermúdez, M., van der Beek, P., Bernet, M., 2013, Strong tectonic and weak climatic control on exhumation rates in the Venezuelan Andes: Lithosphere, 5(1), 3-16. https://doi.org/10.1130/L212.1

Beven, K.J. Kirkby, M.J., 1979, Un modelo de hidrología de cuenca de contribución variable con base física: Boletín de Ciencias Hidrológicas, 24 (1), 43–69. http://dx.doi.org/10.1080 / 02626667909491834

Brandon, M.T., Roden-Tice, M.K., Garver, J.I., 1998, Late Cenozoic exhumation of the Cascadia wedge in the Olympic Mountains northwest Washington state: Geological Society of America Bulletin 110, 985 – 1009. https://doi.org/10.1130/0016-7606(1998)110%3C0985:lceotc%3E2.3.co;2

Braun, J. 2003, PECUBE: A new finite-element code to solve the 3D heat transport equation including the effects of a time-varying, finite amplitude surface topography: Computers & Geosciences, 29, 787–794. https://doi.org/10.1016/s0098-3004(03)00052-9

Braun, J., Robert, X., 2005, Constraints on the rate of post-orogenic erosional decay from low-temperature thermochronological data: application to the Dabie Shan, China: Earth Surface Processes and Landforms, 30, 1203-1225. https://doi.org/10.1002/esp.1271

Braun, J., van der Beek, P., Batt, G., 2006, Quantitative Termochronology. Numerical Methods for the Interpretation of Thermochronological Data: Cambridge University. 272 p.

Braun, J., Burbidge, D.R., Gesto, F.N., Sandiford, M., Gleadow, A.J.W., Kohn, B.P., Cummins, P.R., 2009, Constraints on the current rate of deformation and surface uplift of the Australian continent from a new seismic database and low-T thermochronological data: Australian Journal of Earth Sciences, 56, 99–110. http://dx.doi.org/10.1080/08120090802546977.

Braun, J., van der Beek, P., Valla, P., Robert, X., Herman, F., Glotzbach, C., Pedersen, V., Perry, C., Simon-Labric, T., Prigent, C., 2012, Quantifying rates of landscape evolution and tectonic processes by thermochronology and numerical modeling of crustal heat transport using PECUBE: Tectonophysics, 524-525, 1-28. https://doi.org/10.1016/j.tecto.2011.12.035

Bustos, X., Bermúdez, M.A., 2015, Caracterización de las superficies de erosión en Los Andes del Norte a través de métodos numéricos termocinemáticos y geomáticos: Cuadernos Cenamb, Caracas, Universidad Central de Venezuela, 58 p.

Campbell, C.J., Bürgl, H., 1965, Section Through the Eastern Cordillera of Colombia, South America: GSA Bulletin, 76(5), 567-590. https://doi.org/10.1130/0016-7606(1965)76[567:stteco]2.0.co;2

Campforts, B., Shobe, C. M., Steer, P., Vanmaercke, M., Lague, D., Braun, J., 2020, HyLands 1.0: a hybrid landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution: Geoscientific Model Development, 13(9), 3863–3886, https://doi.org/10.5194/gmd-13-3863-2020

Cediel, F., Shaw, R. P., Cáceres, C., 2003, Tectonic assembly of the Northern Andean block: 815–848. Retrieved from http://archives.datapages.com/data/specpubs/memoir79/CHAPTER37/CHAPTER37.HTM

Chicangana, G., Vargas-Jiménez, C.A., Kammer, A., Hernández-Hernández, T.A., Ochoa-Gutiérrez, L.H., 2007, Caracterización sismotectónica regional preliminar de un sector del piedemonte llanero colombiano: corredor San Juan de Arama – Cumaral, Meta: Boletín de Geología, 29, 61-74.

Clarín. 2008, El epicentro fue en El Calvario, 50 km. al sudeste De Bogotá, 25 de mayo del 2008, Recuperado de: https://www.clarin.com/ediciones-anteriores/terremoto-colombia-dejo-muertos_0_SJ4GTz6RaFl.html

Cooper, M., Addison, F., Alvarez, R., Coral, M., Graham, R., Hayward, A., Howe, S., Martínez, J., Naar, J., Peñas, R., Pulham, J., Taborda, A., 1995, Basin development and tectonic history of the Llanos Basin, Eastern Cordillera, and Middle Magdalena Valley, Colombia: AAPG Bulletin, 79 (10), 1421-1443. https://doi.org/10.1306/7834d9f4-1721-11d7-8645000102c1865d

Flint, J., 1974, Stream gradient as a function of order magnitude and discharge: Water Resources Research, 10, 969–973. https://doi.org/10.1029/wr010i005p00969

Forte, A.M., Whipple, K.X., 2019, The topographic analysis kit (TAK) for TopoToolbox: Earth Surface Dynamics 7, 87-95. https://doi.org/10.5194/esurf-7-87-2019

Gómez, E., Jordan, T.E., Allmendinger, R.W., Cardozo, N., 2005, Development of the Colombian foreland-basin system as a consequence of diachronous exhumation of the northern Andes: GSA Bulletin, 117, 1272-1292. https://doi.org/10.1130/b25456.1

Gómez-Tapias, J., Montes, N.E., Nivia, Á., Diederix, H., 2015, Mapa Geológico de Colombia 2015: Escala 1:1 000 000, Servicio Geológico Colombiano, Bogotá, 2 hojas.

Goulsbra, C., Evans, M., Lindsay, J., 2012, Temporary streams in a peatland catchment: the pattern and timing of stream network expansion and contraction and controls on these: Geophysical Research Abstract, 14, EGU2012-11014-1.

Guerrero, J., 2002, A proposal on the classification of System Tracts: Application to the Allostratigraphy and Sequence Stratigraphy of the Cretaceous Colombian Basin Part 1: Berrriasian to Hauterivian: Geología Colombiana, 27.

Gutenberg, B., Richter, C.F., 1954, Seismicity of the Earth and Associated Phenomena: Princeton, New Jersey, Princeton University Press, 245 p.

Hack, J., 1973, Stream-profile analysis and stream-gradient index: Journal of Research of the U.S. Geological Survey, 1(4), 421-429.

INGEOMINAS, 2011, Cartografía Geológica De La Plancha 247 Cáqueza Escala: 1:100.000 levantada por Patiño A, Fuquen J, Ramos J, Pedraza A, Ceballo L y otros, Bogotá, Colombia.

Kendall, M.G., Stuart, A., 1973, The Advanced Theory of Statistics, 2: Inference and Relationship: Griffin.

Koons, P. O., Zeitler, P. K., Chamberlain, C. P., Craw, D., and Melzer, A. S., 2002, Mechanical links between river erosion and metamorphism in Nanga Parbat, Pakistan Himalaya: American Journal Science, 302 (9), 749– 773, https://doi.org/10.2475/ajs.302.9.749

Laske, G., Masters., G., Ma, Z., Pasyanos, M., 2013, Update on CRUST1.0 - A 1-degree Global Model of Earth’s Crust: Geophys. Res. Abstracts, 15, Abstract EGU2013-2658

Montgomery, D. R., Brandon, M. T., 2002, Topographic controls on erosion rates in tectonically active mountain ranges: Earth and Planetary Science Letters, 201(3), 481-489. https://doi.org/10.1016/S0012-821X (02)00725-2

Moore, I.D., Burch, G.J., 1986, Physical Basis of the Length Slope Factor in the Universal Soil Loss Equation: Soil Science Society of America, 50, 1294-1298. https://doi.org/10.2136/sssaj1986.03615995005000050042x

Mora, A., Parra, M., Strecker, M. R., Kammer A., Dimaté, C., Rodriguez, F., 2006, Cenozoic contractional reactivation of Mesozoic extensional structures in the Eastern Cordillera of Colombia: Tectonics, 25, TC2010. http://dx.doi.org/10.1029/ 2005TC001854.

Mora, A., Parra, M., Strecker, M. R., Sobel, E. R., Hooghiemstra, H., Torres, V., Vallejo-Jaramillo, J., 2008, Climatic forcing of asymmetric orogenic evolution in the Eastern Cordillera of Colombia: GSA Bulletin, 120 (7-8), 930 – 949. http://dx.doi.org/10.1130/ B26186.1

Mora, A., Gaona T., Kley J., Montoya D., Parra M., Quiroz L. I., Reyes G., Strecker M. R., 2009, The role of inherited extensional fault segmentation and linkage in contractional orogenesis: A reconstruction of Lower Cretaceous inverted rift basins in the Eastern Cordillera of Colombia: Basin Research, 21(1), 111 – 137. https://doi.org/10.1111/j.1365-2117.2008.00367.x

Parra, M., Mora, A., Romero, O., Moron, S., 2006, Secciones estratigráficas detalladas de las Formaciones Carbonera (C1), León y el Grupo Guayabo en los Sinclinales de Nazareth y Medina, Piedemonte Llanero, Cordillera Oriental de Colombia: Internal Report, Smithsonian Tropical Research Institute.

Parra, M., Mora, A., Sobel, E.R., Strecket, M., González, R., 2009, Episodic orogenic front migration in the northern Andes: Constraints from low-temperature thermochronology in the Eastern Cordillera, Colombia: Tectonics, 28, TC4004, http://dx.doi.org/10.1029/2008TC002423.

Perne, M., Covington, M. D., Thaler, E. A., Myre, J. M., 2017, Steady state, erosional continuity, and the topography of landscapes developed in layered rocks: Earth Surface Dynamics, 5(1), 85-100. http://dx.doi.org/10.5194/esurf-5-85-2017

Pike, R.J., Wilson, S.E., 1971, Elevation-Relief Ratio, Hypsometric Integral and Geomorphic Area—Altitude Analysis: Geological Society of America Bulletin, 82, 1079-1084. http://dx.doi.org/10.1130/0016-7606(1971)82 [1079: ERHIAG]2.0.CO;2

Pulido, O., Gómez, L., 2001, Geología de la plancha 266 – Villavicencio (Informe y mapa a escala 1:100.000): INGEOMINAS, Bogotá.

Renzoni, G., 1962, Apuntes acerca de la litología y tectónica de la zona al Este y Sureste de Bogotá: Servicio Geológico Nacional, Boletín Geológico, 10(1-3), 59-79.

Renzoni, G., 1968, Geología del Macizo de Quetame. Geología Colombiana, Universidad Nacional de Colombia, Facultad de Ciencias, Departamento de Geología, Bogotá, 5, 75-127.

Sarmiento-Rojas, L. F., Van Wess, J. D., Cloetingh, S., 2006, Mesozoic transtensional basin history of the Eastern Cordillera, Colombian Andes: Inferences from tectonic models: Journal of South American Earth Sciences, 21(4), 383 – 411. https://doi.org/10.1016/j.jsames.2006.07.003

Schwanghart, W., Kuhn, N. J., 2010, TopoToolbox: A set of Matlab functions for topographic analysis: Environmental Modelling & Software, 25(6), 770–781. https://doi.org/10.1016/j.envsoft.2009.12.002

Schwanghart, W., Scherler, D., 2014, Short Communication: TopoToolbox 2 – MATLAB-based software for topographic analysis and modeling in Earth surface sciences: Earth Surface Dynamics, 2(1–7), https://doi.org/10.5194/esurf-2-1-2014, 2014

Segovia, A, 1963, Mapa Geológico de la plancha L-12 (Medina), de la República de Colombia: Servicio Geológico Nacional Bogotá, Colombia.

SIMMA, 2020, Sistema de Información de Movimientos en Masa: Bogotá, Dirección de Geoamenazas del Servicio Geológico Colombiano, catálogo de movimientos en masa (http://simma.sgc.gov.co/#/, consultada en junio de 2020).

Small, E.E., Anderson, R.S., 1998, Pleistocene relief production in Laramide mountain ranges: western United States: Geology, 26, 123-126. https://doi.org/10.1130/0091-7613(1998)026%3C0123:prpilm%3E2.3.co;2

Sobel, E.R., Hilley, G.E., Strecker, M.R., 2003, Formation of internally drained contractional basins by aridity-limited bedrock incision: Journal of Geophysical Research: Solid Earth, 108(B7), 2344, https://doi.org/10.1029/2002JB001883

Strahler A.N., 1952, Hypsometric area-altitude analysis of erosional topography: GSA Bulletin, 63, 1117-1142. https://doi.org/10.1130/0016-7606(1952)63[1117:haaoet]2.0.co;2

Struth, L., Babault, J., Teixell, A., 2015, Drainage reorganization during mountain building in the river system of the Eastern Cordillera of the Colombian Andes: Geomorphology, 250, 370-383, https://doi.org/10.1016/j.geomorph.2015.09.012

Taboada, A., Rivera, L. A., Fuenzalida, A., Cisternas, A., Philip, H., Bijwaard, H., Olaya, J., Rivera, C., 2000, Geodynamics of the northern Andes: Subductions and intracontinental deformation (Colombia): Tectonics, 19, 787 – 813. http://dx.doi.org/10.1029/ 2000TC900004

Tucker, G., Slingerland, R., 1997, Drainage basin responses to climate change: Water Resources Research, 33(8), 2031-2047. https://doi.org/10.1029/97wr00409

Ulloa, C., Rodríguez, E., 1979, Geología del Cuadrángulo K12, Guateque: Boletín Geológico, 22(1), 3 – 55.

UTGS, 2018, Plan de Manejo y Ordenamiento de la Cuenca del Río Guayuriba, Unión temporal Guayuriba sostenible, Bogotá D. C.

Wagner G., Van den Haute P., 1992, Fission-Track Dating: Dordrecht: Mineralogical Magazine, 57(389), 760-761. http://dx.doi.org/10.1180/minmag.1993.057.389.28

Watson, D. F., Philip, G. M., 1985, A Refinement of Inverse Distance Weighted Interpolation: Geoprocessing 2,315–327.

Whipple, K.X., Tucker, G., 1999, Dynamics of the stream–power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs: Journal of Geophysical Research, 104, 17661–17674.

Manuscrito recibido: 31 de Agosto de 2020
Manuscrito corregido: 22 de Octubre de 2020
Manuscrito aceptado: 21 de Diciembre de 2020