Naturaleza del problema investigado

Los movimientos en masa investigados se localizan en la cuenca del río Negro, entre las poblaciones de San Bartolo y La Dolorosa, en la provincia de Morona Santiago, Ecuador (Figura 1).

El área de estudio está conformada por depósitos de origen aluvial, de origen de roca metamórfica con fragmentos de diverso tamaño, transportados por gravedad y la acción del agua. Su origen es local, producto de la alteración in situ de las rocas y posterior transporte como derrubios de ladera o depósitos de solifluxión (^{de Vallejo, 2004}). Frecuentemente están asociados a masas inestables. Los depósitos tomados en esta investigación se caracterizan por rocas metamórficas (mayormente pizarras y filitas grafíticas), englobados en una matriz (suelo areno-limoso), esta se distribuye irregularmente en las vertientes del territorio montañoso formándose por alteración y desintegración in situ de las rocas ubicadas en las laderas superiores adyacentes y por la acción de la gravedad. Presentan inestabilidad con grietas y diaclasamientos en la superficie.

Las laderas están afectadas por procesos gravitacionales, que consisten en el movimiento pendiente abajo de suelo, roca y regolito, agravados por el agua: Lluvias y caudales erosivos. El agua afecta a la estabilidad de laderas de tres formas básicas:

La infiltración en la masa del suelo proviene del agua de vertientes en la ladera y del agua lluvia, y se produce a través de los sistemas de grietas de tracción existentes que evidencian la descompresión del basamento rocoso, provocando el incremento de la presión de poros y del campo de las fuerzas de filtración, cuya magnitud depende de la permeabilidad de la matriz del coluvión. Cuando los poros del sedimento se llenan de agua, la cohesión de la matriz de limo-arena disminuye, permitiendo que los fragmentos o bloques de material se deslicen unos sobre otros con relativa facilidad. La saturación reduce la resistencia interna de los materiales, los cuales son puestos fácilmente en movimiento por la fuerza de la gravedad (^{Tarbuck, 2005}). La socavación lateral que producen los cuerpos de agua de ríos, los cuales fluyen en el pie de las laderas, disminuye su capacidad de sostenimiento.

Los procesos gravitacionales en la zona son movimientos en masa de tipo deslizamiento superficial, que es un movimiento muy lento, prácticamente imperceptible, el cual afecta suelos y materiales alterados provocando deformaciones continuas que se manifiestan al cabo del tiempo en la inclinación o falta de alineación de árboles, vallas, muros, postes, etc., en las laderas. Esta es una de las formas más lentas de procesos gravitacionales (^{De Vallejo, 2004}). Por otro lado, ^{Cruden & Varnes (1996)} la definen como movimientos graduales de rocas, escombros o suelo a nivel superficial y sub superficial pendiente abajo, sobre una o varias superficies de ruptura. Esta definición se apega más a la naturaleza de los movimientos de coluviales estudiados, ya que por su heterogeneidad, las superficies de ruptura suelen alcanzar diversas profundidades. Otros autores como ^{Highland & Bobrowsky (2008)} la definen como corrientes de tierra lenta (arrastre), causada por una tensión cortante intensa; pero insuficiente para causar una falla.

El desencadenamiento de los movimientos de ladera por causas meteorológicas y climáticas está relacionado fundamentalmente con el volumen, intensidad y distribución de las precipitaciones, así como con la erosión fluvial. (^{De Vallejo, 2004}).

La descripción del tipo de movimientos se realizó con base en la escala propuesta por ^{Cruden & Varnes (1996)}, (Figura 2). Es aplicable a movimientos de reptación, excluyendo hundimientos, colapsos o avalanchas que suelen producirse de modo violento.

^{Cruden & Varnes (1996)}, consideran que la velocidad de movimiento puede ser el mejor parámetro para describir el potencial de daño de un deslizamiento; sin embargo, es muy difícil de predecir e incluso medir. El modelo que se presenta permitirá superar estas dificultades, y se espera contribuir al avance en posteriores investigaciones de esta importante área de la geología.

Los movimientos son de tipo traslacional, ya que los movimientos rotacionales se dan mayormente en laderas conformadas por suelos homogéneos (^{De Vallejo, 2004}).

Se han establecido umbrales de intensidad y duración de la lluvia para el desencadenamiento de movimientos en masa en diferentes lugares. El principal inconveniente para su aplicación como medida preventiva es que no son extrapolables a otras zonas (^{De Vallejo, 2004}).

Por lo señalado, es de interés deducir una ley matemática que permita predecir la velocidad de reptación que se presentaría bajo la acción de determinadas condiciones hidrometeorológicas propias de una región (precipitación media y caudal medio del río) aplicable a coluviones. En estos materiales, la hidrogeología puede ser muy compleja por su heterogeneidad; y como consecuencia muestran una respuesta errática ante los efectos de las lluvias (^{Díaz, 2009}); por lo que el análisis estadístico de la información reveló una correlación moderada entre valores de velocidad observados y calculados.

La ley de reptación es una ecuación que expresa la dependencia entre la velocidad que adopta una masa de coluviones, con un conjunto de ocho variables, de las cuales se eligen como fundamentales a la precipitación media y el caudal medio del río Negro. Estas variables fundamentales se les conoce como predictivas, y a la variable dependiente (velocidad) se le conoce como predecible.

La información de las características hidrometeorológicas de la zona, así como de los estudios geofísicos y geotécnicos, fue proporcionada por HIDROSANBATOLO, que aprovecha las aguas del río Negro en esta zona para la generación de energía eléctrica.

El período de la investigación comprendió desde el 02 agosto de 2016 hasta el 19 noviembre de 2019. Se investigaron seis zonas inestables, las cuales se designaron mediante el abscisado de la tubería de conducción de la central HIDROSANBARTOLO, cuyo trazado se desarrolla en este sector (Tabla 1).

Marco Geológico Regional

En la Figura 3a se muestran los materiales constitutivos de las laderas. La roca de la formación Napo deformada (K_ND) subyace a los depósitos coluviales (Q_C).

La zona se ubica en el frente Subandino, la cual se encuentra afectada por fallas inversas que han causado fracturamientos y cizallamientos en el basamento rocoso. La Figura 3b muestra la presencia de grandes extensiones agrietadas que atraviesan quebradas y cursos de agua naturales y donde se producen infiltraciones.

Existen rasgos morfológicos indicativos de que esta región ha sido afectada por movimientos en masa de tipo reptación de diversa magnitud y edad (Figuras 4 y 5).

Explicación del método. Fundamento matemático

La base matemática de esta investigación es el análisis dimensional, basado en el teorema Π “pi” de Buckingham. El procedimiento aquí expuesto se basa en la obra de ^{Gibbinngs (2011)}. Su aplicación se expone a continuación:

La variable predecible considerada es la velocidad media de reptación V, que adopta una masa de coluviones durante un período de tiempo.

La Tabla 3 expone el simbolismo para las dimensiones básicas adoptado por el Sistema Internacional de Unidades (^{Gibbinngs, 2011}).

Las dimensiones de las cantidades físicas que intervienen en el problema se expresan a partir de las básicas (Tabla 4).

Elección de variables

Las variables se eligen de manera que sean las más convenientes experimentalmente (^{Gibbinngs, 2011}). Cada variable y sus dimensiones se detallan como sigue (para denotar las dimensiones de una variable se utilizó el corchete “[ ]”):

V: Velocidad de reptación [V]=LT^-1

P: Precipitación media [P] = LT⁰

Q: Caudal medio fluvial [Q]=L³T^-1

g: Aceleración de la gravedad [g]=LT^-2

A: Área del deslizamiento [A]=L²T⁰

h_s : Altura de la masa afectada por la presión de poros [h_s]=LT⁰

k: Coeficiente de permeabilidad del suelo [k]=LT^-1

j: Pendiente media natural de la ladera [j]=L⁰T⁰

tanφ: Tangente del ángulo de fricción interna de suelo [tanφ]=L⁰T⁰

Cuando las dimensiones están elevadas al exponente cero significa que la variable es adimensional o que carece de esa dimensión. Por ejemplo, la precipitación tiene dimensión de longitud (L), y carece de la dimensión tiempo T, por lo que [P]=LT⁰. Esta notación será de utilidad en la deducción de los parámetros adimensionales.

Se adoptaron como variables fundamentales P y Q, que contienen las dimensiones L y T, pues engloban de mejor manera las características cinemáticas, geométricas, geotécnicas e hidrológicas del problema. A estas variables se las conoce como predictivas, y para las mismas, se verificó que el determinante de los exponentes de sus dimensiones es diferente de cero. Así, en la variable P: los exponentes son 1 y 0, y en la variable Q son 3 y -1, luego:

Cálculo del número de parámetros:

Número de variables = 9

Número de dimensiones fundamentales = 2

Número de parámetros adimensionales = 9 - 2 = 7.

Las variables restantes se combinan con el producto de potencias de las variables predictivas para formar los parámetros adimensionales pi, cuyas expresiones se obtienen como sigue:

Velocidad de reptación V:

Aceleración de la gravedad g:

Área del deslizamiento en planta A:

Espesor de la masa afectada por la presión de poros h_s:

Coeficiente de permeabilidad k:

Las variables j y tanφ son adimensionales, por lo que las ecuaciones de sus parámetros toman la misma expresión:

El teorema de Buckingham establece que existe alguna función f₁ tal que:

Dado que se considera que estas variables explican en sí mismas el fenómeno ^{Gibbinngs (2011)}.

Es posible mediante el producto de los parámetros pi combinar estos constituyendo nuevos parámetros, ya que son adimensionales.

Combinando los parámetros Π1, Π3 y Π5, formaremos el parámetro Π8 que será igual a:

Combinando los parámetros Π2, Π4 y Π6 y el inverso del parámetro Π7, formamos el parámetro Π9 que es igual a:

Podemos formar la función f₂ tal que:

Finalmente, si se asumen ciertas condiciones de regularidad y diferenciabilidad sobre la función anterior, podrá usarse el teorema de la función implícita para expresar el parámetro Π8 que contiene la variable predecible V, mediante una función f₃ del parámetro Π9:

La función f₃ es desconocida, y se determinó calculando Π8 y Π9 con los datos de campo para cada una de las variables.

Escarpes. Implantación de mojoneras en el terreno

Los escarpes se graficaron en un mapa topográfico, para determinar el área de cada deslizamiento, representada por la variable A. En campo se procedió a la colocación manual de 77 mojoneras cilíndricas de hormigón de resistencia 180 kg/cm² a 28 días, de diámetro 15 cm y altura 30 cm.

Para que los datos obtenidos reflejen mejor el resultado de la investigación, se colocaron las mojoneras en el terreno, siguiendo los siguientes criterios:

Las Figuras 6 y 7 muestran el escarpe principal para cada zona, los mojones instalados y rasgos morfológicos.

Monitoreos topográficos. Cálculo del desplazamiento y de la velocidad

Una vez implantadas las mojoneras se realizó la toma de sus coordenadas iniciales con base en el sistema WGS84, ZONA 17S y se ubicaron en los mapas geomorfológicos.

El valor para la variable V se determinó mediante toma semanal a quincenal de coordenadas con estación total de las 77 mojoneras. Los hitos geodésicos del polígono de referencia han sido fijados en sitios de afloramientos de roca, en donde es muy poco probable que existan movimientos, y así evitar que el error de medición sea excesivo.

Si el vector posición del mojón i en la fecha 1 es r₁ y en la fecha 2 es r₂, el vector desplazamiento Δr del mojón i entre estas dos fechas es ∆r=r2-r1. Si se conocen las coordenadas topográficas para estas fechas (E1, N1, Z1) y (E2, N2, Z2), la magnitud del desplazamiento en cm se calculó como:

Si Δt es el tiempo en días que transcurrió entre las fechas 1 y 2, la magnitud de la velocidad del mojón i en cm/d se calculó como:

La velocidad promedio observada Vo entre dos fechas para un deslizamiento se calculó como la media aritmética de las velocidades individuales de los mojones.

con l = Número de mojones (del inglés landmarks).

La velocidad promedio observada Vo, representa la medida de la variable V. Esta velocidad es causada por una precipitación media P y un caudal fluvial medio Q en el período de tiempo correspondiente a Vo.

A modo de ejemplo explicativo, en las Tablas 5 y 6 se expone el cálculo del desplazamiento y la velocidad para los mojones 5 y 7 del deslizamiento en la Zona 1, en el periodo junio-julio 2018.

Las coordenadas se tomaron con el Datum WGS84.

La velocidad promedio observada Vo entre dos fechas, se calculó como la media aritmética de las velocidades individuales de los mojones.

siendo l = número de mojones (del inglés landmarks). El valor Vo es la medida para la variable predictiva V del modelo, en un periodo de tiempo determinado.

Geotecnia e hidrología

Los deslizamientos investigados están afectados por la presión de poros, en un espesor h_s que comprende la altura desde el nivel de la capa con Vs = 3 000 m/s (roca impermeable) hasta el respectivo nivel freático. El modelo geológico proporciona información útil para estimar h_s: Dado que la capa con velocidad de onda 2000 a 2500 m/s se conforma de roca fracturada, posee permeabilidad secundaria debido al fracturamiento provocado por las fallas inversas, características del frente subandino, y estará afectada por la presión de poros. Como la capa con Vs = 3000 m/s corresponde a roca impermeable, el espesor h_s se determina a partir del nivel superior de la roca impermeable, hasta el límite superior dado por la posición del nivel freático.

En la Figura 10, si Hm es el espesor total promedio de la masa deslizante, y NF la profundidad del nivel freático:

La profundidad del nivel freático para cada zona inestable se investigó mediante observaciones hechas en calicatas y sondeos eléctricos verticales (SEV) (Figura 11). Los valores son estimados, debido a la heterogeneidad del suelo que constituye la matriz del coluvión, lo que implica la existencia de acuíferos confinados reflejando profundidades freáticas irregulares.

Para cada zona, en la Tabla 7 se reporta el promedio de las profundidades observadas.

En la Tabla 8 se presentan los espesores promedios Hm y el cálculo de h_s para cada zona.

Los parámetros de resistencia al corte de la matriz del coluvión (cohesión C y ángulo de fricción interna ϕ), se determinaron mediante ensayos de corte triaxiales no drenados (Tabla 9).

La permeabilidad del coluvión se determinó insitu en cuatro puntos de cada ladera en estudio, mediante el método de ensayo Lefranc (^{Marinas, 2002}). La medida de la permeabilidad de un suelo expresa la facilidad con que los flujos subterráneos pueden circular en la masa (^{Marinas, 2002}). En suelos con permeabilidad muy pobre, el flujo subterráneo de agua presentará mayor resistencia, lo que incidirá en una mayor velocidad (Tabla 10).

A partir del mapa topográfico, se delimitó el área interior del escarpe para cada zona. El límite inferior es la orilla del río Negro. El límite superior o la “corona” del deslizamiento lo constituyen las zonas altas ubicadas en las cercanías del camino que une las poblaciones de San Bartolo y La Dolorosa. La delimitación lateral del área se realiza con un trazo que esté acorde con el direccionamiento de las grietas y rupturas de pendiente, y que abarque hondonadas y pantanos, los cuales son zonas de potencial infiltración. El área de cada deslizamiento se expone en la Tabla 11.

La pendiente se calculó mediante la media ponderada en la longitud total de los rumbos o direcciones que sigue la descompresión de las masas, los cuales se estiman mediante la morfología de las grietas de tracción. Los valores se presentan en la Tabla 12.

Para la zona 1, por ejemplo, se pueden observar en la Figura 6 dos direcciones de movimientos: La primera de longitud 648.96 m entre las elevaciones 865 y 760 msnm, con pendiente j1=105/648.96=0.162. La segunda, de longitud 630.74 m entre las elevaciones 760 a 605 msnm, con pendiente j2=155/630.74=0.256 . La pendiente media ponderada en la longitud total es:

La precipitación media se calculó a partir de los registros diarios de lluvias de un pluviómetro ubicado en la cuenca del río Negro a 2.2 kms de la zona 1. Se procedió de la siguiente manera: En cada intervalo de tiempo de medición de la velocidad, el área bajo la curva de precipitación, dividida para el intervalo de tiempo representa la precipitación media. Por ejemplo, entre el 11 y 19 de junio 2018 el área bajo la gráfica es igual a 108.40 mm x d (Figura 12), que dividida para 8 días muestra 13.55 mm.

De manera similar, el caudal medio en este periodo de tiempo es igual al área bajo la curva de caudales, 2103.52 (m³/s) x d, dividida para 8 días, que da 262.94 m³/s (Figura 13).

Al valor del caudal medio se restó el caudal de 30 m³/s que capta la central HIDROSANBARTOLO para la generación de electricidad.

Cálculo de parámetros adimensionales Π

Relacionando los valores de velocidad observados Vo de cada zona con los valores de sus variables correspondientes y con las variables predictivas P y Q, se calcularon los parámetros Π8 y Π9 mediante (3) y (4). Para la aceleración de la gravedad se ha tomado el valor estándar g= 9.81 m/s².

Las Figuras 14 a 16 muestran el cálculo de los parámetros adimensionales Π para cada zona inestable.

Con los valores obtenidos se realizaron gráficos de dispersión entre estos parámetros para determinar la función f₃ que tenga el mayor coeficiente de determinación R². (Figura 17).

Sustituyendo (3) y (4) en Y y X, respectivamente, se determinaron las siguientes expresiones:

Zona 1:

Zona 2:

Zona 3:

Zona 4:

Zona 5:

Zona 6:

Ley general de reptación para coluviones

Las ecuaciones halladas difieren para cada deslizamiento, pero presentan la misma estructura algebraica, lo que condujo a deducir la ley general de reptación de coluviones, la cual se puede expresar mediante la ecuación:

con P≥0.54mm;n>0.500; k˄A>0.

De (10) a (15), tomando para cada deslizamiento el coeficiente m, se verificó una buena correlación con el área correspondiente de cada zona inestable (Figura 18).

La función potencial mostrada, muestra el mayor poder predictivo para el coeficiente m:

Para definir los exponentes, se calculó la media aritmética de los valores obtenidos en (10) a (15):

Reemplazando (17) y (18) en (16), se obtuvo la Ley General de reptación de coluviones:

con P en mm, Q en m³/s, A en m², k en cm/s, g = 9.81 m/s², h_s en m, j y tanφ adimensionales.

La ley expresada mediante (19) es válida para deslizamientos de laderas en cuencas hidrográficas constituidas por coluviones, caracterizados por rocas metamórficas en matriz limo-arenosa, y fue determinada para deslizamientos con extensiones desde 228500 hasta 518100 m², en una región afectada por fallas inversas, que han causado fracturamientos y cizallamientos en el basamento rocoso.

Los caudales medios han sido calculados tomando en cuenta la regulación de la obra de captación de la central HIDROSANBARTOLO, que utiliza el caudal de 30 m³/s para la generación de electricidad.

La velocidad de reptación es inversamente proporcional a la permeabilidad k del suelo de la matriz, siendo mayor el riesgo en laderas con permeabilidad muy baja.

Correlación lineal entre valores observados (Vo) y valores del modelo (V)

Los gráficos de correlación nos indican un ajuste muy débil para valores de velocidad menores a 1 cm/d.

Se obtuvieron coeficientes de correlación lineal “r” mayores a 0.50. El análisis reveló que este se obtiene para caudales del río mayores a 26 m³/s. Para caudales menores a 26 m³/s, el modelo no es confiable. El análisis realizado con datos del modelo para Q < 26 m³/s dio como resultado coeficientes de correlación muy bajos, por lo que estos caudales fueron omitidos. En estadística es permisible eliminar ciertos datos adecuadamente, siempre que exista una justificación apropiada (^{Navidi, 2006}).

Estadísticamente, un coeficiente de correlación mayor o igual a 0.50 representa un ajuste moderado, mientras que valores inferiores a este representan un ajuste débil (Navidi, 2006). La Figura 19 muestra los ajustes lineales y los coeficientes de correlación respectivos para cada zona inestable.

El poder predictivo del modelo está en el rango 25 al 26 %.

La restricción para Q > 26 m³/s se fundamenta en la protección que confieren las terrazas aluviales en el pie de las laderas, por lo que caudales menores no constituyen la causa para acelerar los movimientos por socavación lateral.

La heterogeneidad del suelo también ocasiona que el poder predictivo del modelo tome el máximo de 26 %.

Descripción de los movimientos en masa

La Figura 2 muestra la escala de velocidades, en donde se describen los movimientos como: Muy lento, Lento, Moderado, Rápido, Muy rápido y Extremadamente Rápido para 6 valores típicos de la velocidad.

En un período de siete días, se registraron:

P máxima = 28.15 mm

Q máximo = 259.21 m³/s

Bajo estas condiciones hidrometeorológicas, aplicando (19) se determinó la siguiente descripción de los movimientos con base en la escala de Cruden y varnes (Figura 2, Tabla 13).

La velocidad máxima es 2.90 cm/d y se produce en la zona 2. Está en el intervalo 1.5 m/año - 1.5 m/mes (0.411 - 5.00 cm/d), para el cual, el movimiento se describe como “Lento” según la escala de Cruden y Varnes (Figura 2).

El mes de mayo es el más lluvioso en la región. En mayo 2020 (31 días), se registraron para las variables predicitvas los siguientes valores medios mensuales (Tabla 14):

P máxima = 18.02 mm

Q máximo = 183.14 m³/s

En la zona 2 se presentaría la máxima velocidad, de 1.99 cm/d, movimiento del tipo “Lento”, con desplazamiento de 61.69 cm en 31 días.

Umbrales de precipitación y caudal para movimiento moderado

Aplicando (19) se calculó la precipitación media y caudal medio, para el cual se desencadenó el movimiento moderado en las laderas, el cual se dará cuando V supere los 5.0 cm/d. Los valores fueron: P = 18 mm y Q = 600 m³/s. Bajo estas condiciones, en los sectores 1 y 2 se desencadenará un movimiento “Moderado” (Tabla 15).

Se registraron durante la investigación, en un período de 24 horas, precipitaciones de 39.5 mm y caudales medios de 761.58 m³/s, que son mayores a los umbrales calculados, y desencadenarán movimientos moderados en zonas 1, 2, 4 y 5 (Tabla 16).