Introducción

La conservación del suelo es esencial para la seguridad alimentaria y un futuro sostenible. Sin embargo, cerca de 33% de los suelos a nivel mundial presentan degradación de moderada a grave, siendo la erosión hídrica una de las principales causas. En México este tipo de erosión ocupa 66% de la superficie nacional (^{INEGI, 2015}) y 1.97% de esta superficie, es producto de la acción del flujo concentrado que causa incisión y socavamiento sobre la superficie del suelo, formando canales de 0.5 a 30 m (Soil Science Society of America, 2001 citado por ^{Rahmati et al., 2017}).

El desarrollo y crecimiento de cárcavas sobre el terreno, depende de la combinación de diversos factores, como: precipitación, topografía, patrón de drenaje, litología, características del suelo, uso y manejo del mismo (^{Poesen et al., 2003}; ^{Conoscenti et al., 2014}; ^{Rahmati et al., 2017}). Comúnmente, el crecimiento lateral de las cárcavas se ha estimado como un incremento de su área o longitud a través del tiempo, mediante clavos de erosión, fotografías aéreas, imágenes LIDAR y reconstrucción de fotografías - 3D (^{Castillo y Gómez, 2016}; ^{Rahmati et al., 2017}), con especial atención en la cabecera de la cárcava. Sin embargo, la erosión en los taludes de cárcavas es la fuente principal de sedimento durante la etapa final del desarrollo de la cárcava (^{Qin et al., 2018}) y el proceso principal de su evolución (^{Chaplot et al., 2011}).

En el estado de Hidalgo (México), se reportaron afectaciones por cárcavas en 1.44% de su superficie (^{INEGI, 2014}), particularmente en el ejido de Santo Domingo (Huasca de Ocampo), donde la erosión en cárcavas ha provocado la pérdida de fertilidad del suelo y reducido la superficie agrícola. Ante esta problemática, la presente investigación tuvo como objetivo evaluar las características topográficas y propiedades físico-mecánicas en taludes de cárcavas para identificar los principales factores que inciden en su estabilidad, mediante análisis de componentes principales, y generar un modelo de regresión lineal múltiple que permita predecir su avance lateral.

Materiales y Métodos

Área de estudio. Comprende dos cárcavas ubicadas en una cuenca al sur del municipio de Huasca de Ocampo, pertenecientes al ejido de Santo Domingo (Figura 1). La cárcava 1 se localiza en las coordenadas geográficas 20° 09’ 30.16” y 20° 09’ 26.07” N; 98° 29’ 09.96” y 98° 29’ 19.99” O, y la cárcava 2 en las coordenadas 20° 08’ 55.57” y 20° 8’ 53.26” N; 98° 29’ 32.31” y 98° 29’ 32.95” O. Las cárcavas presentaron taludes con horizonte A visible y presas de control de azolves, que han favorecido la estabilización del lecho. El área de estudio se encuentra a 2360 m de altitud, presenta un clima templado sub-húmedo, con precipitación media anual de 862 mm y temperatura promedio de 17 °C. El suelo dominante es vertisol crómico, asociado a agricultura de temporal, bosque de encino y vegetación secundaria arbustiva de bosque de encino.

Figura 1: Localización del área de estudio. 

Levantamiento topográfico. Se realizó un levantamiento topográfico con la estación total SOKIA SET 630RK3 de 20 secciones trasversales. La información se procesó con QGIS, y se obtuvieron las covariables mostradas en la Figura 2, así como la pendiente longitudinal y transversal del terreno a la altura de la cresta, la curvatura vertical (CV) y horizontal (CH), y orientación de los taludes.

Figura 2: Variables determinadas en sección transversal, longitud (L) y altura del talud (A_t). 

Determinación de propiedades físico-mecánicas. Se realizó un muestreo de suelo en 30 sitios sobre los taludes de las secciones transversales. La densidad aparente (ρ _a ) se determinó con el método del cilindro de volumen conocido de 0-5 y 5-10 cm. La resistencia a la penetración (RP) se obtuvo con un penetrómetro de ultrasonido (FIELDSCOUT SC 900, SPECTRUM) a 2.5, 5 y 10 cm de profundidad. La textura se obtuvo con el método de Bouyoucos modificado (previa separación de arenas); la erodabilidad (K) se calculó de acuerdo a la ecuación de ^{Wischmeier y Smith (1978)}; los límites de Atterberg: limite líquido (LL), plástico (LP) e índice de plasticidad (IP) con la copa de Casagrande (^{ASTM 4318, 2005}) y el índice de contracción (IC) de acuerdo al manual de la ^{SRH (1962)}, modificando el uso de Hg por parafina. La distribución del tamaño de partículas (DTP), se realizó con seis tamices con aberturas de 2.0, 1.0, 0.6, 0.25, 0.15 y 0.075 mm, la conductividad hidráulica (Ks) se estimó con el permeámetro de carga constante (^{Klute y Dirksen, 1986}) en muestras inalteradas y la cohesión mediante un ensayo de compresión simple no drenado (^{ASTM 2166, 2000}).

Estimación de la velocidad de avance lateral. El avance lateral (AV) se obtuvo en m año^-1, con la comparación de una ortofoto de 1995 (^{INEGI, 1995}) y el levantamiento topográfico realizado en el año 2017.

Análisis estadístico. Se utilizó un ACP para identificar las variables de mayor importancia en el avance lateral de las cárcavas (AV), seleccionando los componentes (CP) con valores propios >1 (^{Kaiser, 1960}) y en cada componente las tres variables con los mayores vectores propios. Posteriormente, a través de RLM, se determinó una ecuación de predicción para el AV, utilizando como variables independientes aquellas seleccionadas en cada CP. La RLM utilizó un método escalonado, y en el modelo obtenido se evaluaron los siguientes criterios: a) factor de inflación de varianza (VIF), para descartar la multicolinealidad, VIF >2 (^{Lin, 2008}); b) valor de S, para evaluar qué tan bien el modelo describe a la variable dependiente; c) valor de R²; d) valor de R² pronóstico, para determinar la predicción en nuevas observaciones; e) valor P ≤ α de las variables seleccionadas en el modelo, f) valor P ≤ α de la regresión y g) PBIAS (porcentaje de sesgo) de acuerdo a ^{Moriasi et al., 2007}. El valor de α fue de 0.10. El análisis estadístico se realizó con los softwares R Studio y Minitab Inc.

Resultados y Discusión

Los terrenos adyacentes a las secciones transversales tuvieron una pendiente longitudinal desde 3.62 a 6.8%, y una pendiente transversal promedio de 2.45%, consideradas suaves, con un potencial débil a la presencia de erosión y de tipo laminar (^{Ferrando, 1993}). La altura y longitud de los taludes varió de 2.78 a 5.23 m, con una media de 3.88 m y de 4.62 a 10.41 m, con una media de 7.27 m (Figura 3a). La media del ángulo de la pendiente de los taludes fue de 26.13° con variaciones de 15 a 40° (Figura 3b) que, desde el punto de vista geomorfológico, se consideraron moderadas o moderadamente escarpadas (^{Ferrando, 1993}).

Figura 3: Altura, longitud (a) y pendiente de taludes (b). 

El espesor promedio de la cresta fue de 42.20 cm, variando de 30 a 62 cm (Figura 4), ubicándose en los horizontes húmico de textura franco migajonosa y de textura franco arcillosa, con presencia de caras de deslizamiento (HA- HBSs del suelo), siendo la porción del terreno aledaño que forma parte de la sección transversal de la cárcava y donde se inicia el talud de la cárcava. En esta zona del terreno se presentan grietas producto de los procesos de humedecimiento y secado del suelo, que provocan su desprendimiento en forma de bloques que caen sobre la pendiente de los taludes.

Figura 4: Espesor de la cresta en secciones transversales. 

La curvatura vertical de los taludes fue convexa (valor negativo) en 46.67%, indicando que el flujo desacelerará, 46.67% fue cóncava (valor positivo), lo que aumenta la velocidad del flujo y riesgo de erosión y 6.67% restante fue lineal (valor de 0) (Figura 5a), sin presentar cambios en las condiciones de erosión o deposición (^{Toledo, 2017}). En la curvatura horizontal 50% fue cóncava, y el flujo se concentrará sobre la ladera, incrementando el desarrollo de la cárcava (valor negativo), 43.33% fue convexa (valor positivo), por lo que el flujo se dispersará y no presentará riesgo sobre el talud y el resto fue lineal (valor de 0), donde el flujo se conducirá sobre el talud de forma paralela (Figura 5b). Todos los taludes tuvieron como común denominador la orientación norte, correspondiendo 66.67% una orientación noreste, 23.33% norte y 10% noroeste. De forma general, los taludes con orientaciones al norte suelen recibir menor luz y, por consiguiente, tener mayores contenidos de humedad con respecto a la orientación sur (^{Nadal et al., 2007}).

Figura 5: Curvatura vertical (a) y horizontal (b) sobre taludes.  

Clasificación mecánica del regolito. En términos de la carta de plasticidad y la clasificación del Sistema Unificado de Clasificación de Suelos, los valores medios de limite líquido (LL) e índice de plasticidad (IP), indicaron que el regolito se comporta como un limo inorgánico con ligera plasticidad (ML), pero próximo al límite de las arcillas minerales de plasticidad media (Figura 6).

Figura 6: Limite líquido (LL) e índice de plasticidad (IP) del regolito, representados sobre la carta de plasticidad de Casagrande. Arcillas inorgánicas de plasticidad baja a media (CL), arcillas inorgánicas de alta plasticidad (CH), limos inorgánicos de ligera plasticidad (ML), limos orgánicos o arcillas limosas orgánicas de baja plasticidad (OL), Limos inorgánicos (MH), arcillas orgánicas de plasticidad media a alta o limos orgánicos (OH). 

Los limos inorgánicos con ligera plasticidad son inestables por su propia naturaleza, particularmente cuando aumenta la humedad, con tendencia a fluir cuando hay saturación, difíciles de compactar, muy susceptibles a erosionarse y sujetos a tubificación. Asimismo, se caracterizan por tener regular resistencia al corte en estado compacto y de saturación, nula a ligera resistencia a la disgregación y riesgo medio a deslizamientos (^{Lambe y Whitman, 1974}).

Descripción de propiedades físico-mecánicas del regolito. De acuerdo al contenido medio de arenas, limos y arcillas, el regolito presentó una textura arcillosa. La distribución del tamaño de partículas (Figura 7), muestra que el regolito presentó mayor porcentaje de microagregados (<0.25 mm) con respecto a macroagregados (0.25 a > 2 mm); representando una degradación de la estabilidad estructural del suelo y una mayor susceptibilidad a los procesos erosivos (^{Volverás y Amézquita, 2009}).

Figura 7: Porcentaje de la distribución de tamaño de partículas (DTP).  

La media de la densidad aparente (ρ _a ) de 0-5 cm, fue de 1.06, variando de 0.72 a 1.52 gr cm^-3, y de 5 a 10 cm de 0.73 a 1.59 gr cm^-3 con una media de 1.11 gr cm^-3, observando un aumento de la ρ _a con la profundidad. La resistencia a la penetración (RP) (Figura 8a, b), tuvo un comportamiento similar a la ρ _a , con valores promedio a los 2.5, 5 y 10 cm de 603.97, 972.31 y 1682.32 kPa, respectivamente. Las variaciones a los 2.5 cm fueron de 140 a 1544 kPa, a los 5 cm de 103.21 a 1404 kPa y los 10 cm de 412.78 a 3195 kPa. De forma general, a menores valores de ρ _a , RP y compactación, mayor espacio poroso, capacidad de almacenamiento de agua, y viceversa.

Figura 8: Densidad aparente (a) y resistencia a la penetración (b) sobre taludes. 

Los valores medios en el índice de plasticidad (IP) (13.90%) y limite líquido (LL) (42.81%), clasificaron a su vez al regolito como mediano y altamente plástico (Figura 9a), con un potencial medio a alto de hinchamiento (^{Murthy, 2002}; Muñoz, 2018¹). El índice de contracción (IC), representó el cambio de volumen del regolito ante procesos de humedecimiento y secado, con un valor medio de 31.57% y variaciones de 27.13% a 38.84% (Figura 9b).

Figura 9: Límites de Atterberg (a) e índice de contracción (b )sobre taludes. 

La conductividad hidráulica (Ks) fue moderada con una media de 7.95 y variaciones de 1 a 16 cm h^-1 (Figura 10). Aunque se esperaría una Ks lenta (1-2 cm h^‑1) debido a la textura del regolito, los valores de Ks se favorecieron por la presencia de grietas a lo largo del perfil, causando un incremento en el movimiento del agua y altas tasas de infiltración.

Figura 10: Conductividad hidráulica sobre taludes. 

La erodabilidad del suelo (Factor K) tuvo un valor promedio de 0.024 Mg ha^-1 año^-1 /MJ h^-1 mm h^-1 y variaciones de 0.018 a 0.034 Mg ha^-1 año^-1 /MJ ha^-1 mm h^‑1 (Figura 11a), lo que indica una resistencia media del regolito a la erosión (^{Figueroa et al., 1991}). La cohesión (Figura 11b) tuvo una media de 3.07 Mg m^-2 y variaciones de 0.64 a 5.38 Mg m^-2, clasificando al regolito como cohesivo rígido (^{GEO5, 2018}), donde la aireación es escasa y las raíces penetran con mucha dificultad y suelen crecer a través de las grietas; caracterizándose por retener gran cantidad de agua y drenaje escaso (^{Murthy, 2002}).

Figura 11: Erodabilidad (a) y cohesión (b) sobre taludes.  

Estimación de la velocidad del avance lateral. El mayor avance lateral (Figura 12) se tuvo en los sitios muestreados sobre los taludes 1, 2 y 3, con valor de 0.17 m año^-1, considerando que los sitios 2 y 3 pertenecen a la misma sección transversal. El avance en estos sitios, se asoció al tipo de la curvatura vertical que favoreció la aceleración del flujo, el arrastre de partículas y el riesgo de erosión sobre el talud, así como la curvatura horizontal para los sitios 2 y 3, que, al ser de tipo cóncava, también aumentó el desarrollo de la cárcava. El porcentaje de limos para los sitios 1 y 2, fue alto con 39% y 35%, incrementando la erodabilidad del suelo. Mientras que, el sitio 3 tuvo uno de los valores más bajos de cohesión (1.66 Mg m^-2). Los tres sitios reportan valores mínimos en el porcentaje de partículas con diámetros de 0.015 y 0.075 mm.

Figura 12: Avance lateral (1995- 2017) sobre taludes. 

El menor avance lateral se tuvo en los sitios de muestreo 7 y 9 (0.04 m año^-1), correspondiendo a la misma sección transversal, donde la curvatura vertical y horizontal disminuyeron la aceración del flujo y su dispersión sobre el talud. Ambos sitios de muestreo, observaron los valores más altos en el porcentaje de partículas con diámetros de 0.015 y 0.075 mm, y valores bajos en el espesor de la cresta y la conductividad hidráulica.

Ahora bien, al considerar el mayor avance lateral (0.17 m año^-1) para una proyección a 50 y 100 años, se tiene que, estos taludes presentarán un avance lateral de 8.5 y 17 m, respectivamente, lo que implica una gran pérdida en la superficie agrícola en un futuro.

En general, en el área de estudio, el comportamiento del flujo de escurrimiento se ve influenciado por el tipo de curvatura en el talud, favoreciendo o no su aceleración y con ello la erosión. Asimismo, el material fino (limos y arcillas) que compone al suelo tiene un comportamiento inestable al formar grietas, que causan el desprendimiento del suelo en forma de bloques y que caen sobre el talud. Este material, ante procesos de humedecimiento y secado, se disgrega fácilmente y cuando ocurre un evento de lluvia es transportado, favoreciendo el crecimiento lateral de las cárcavas.

La presencia de grietas favorece la disgregación, disminuye la cohesión y la resistencia al esfuerzo de corte (^{Bradford y Piest, 1980}), y cuando estas se llenan de agua, la presión de poros aumenta, provocando el estallamiento al interior del suelo y por consiguiente la falla y la caída del material (^{Collison, 2001}; ^{Martínez-Casanovas et al., 2004}). Asimismo, la nula o escasa cobertura vegetal en los taludes de las cárcavas, favorece el crecimiento lateral e inestabilidad de los taludes, al no haber una relación entre las raíces y las partículas, que mejore la cohesión y la resistencia del suelo.

Esta información sugiere que el control de cárcavas, no solo consiste en construir presas de control de azolves, sino que, se debe propiciar la revegetación natural sobre los taludes como ha sido recomendado por ^{Heede (1976)} y ^{Bastola et al. (2018)}, así como algunas técnicas que favorezcan la cementación de las partículas del suelo, a fin de disminuir la energía de impacto del escurrimiento, mejorar las fuerzas cohesivas y la resistencia al esfuerzo de corte. Complementariamente, se debe considerar el cabeceo de las cárcavas, la construcción de canales de desviación del flujo superficial, la revegetación o la construcción de sistemas de terrazas en los terrenos aledaños de las cárcavas.

Análisis de componentes principales. Previo al ACP, se realizó la selección especifica de covariables a fin de evitar información redundante y que de acuerdo a ^{Ferrero et al. (2002)}, cuando algunas covariables no aportan información adicional, o es recogida por otra covariable, es posible eliminarlas sin afectar significativamente los resultados. El resultado del ACP, fue la selección de ocho componentes, con auto valores de 3.47 a 1.07 y contribución total de 76.40% de la variación total (Cuadro 1).

Los vectores propios obtenidos para cada CP y las covariables seleccionadas se muestran en el Cuadro 2. El CP1 y el CP2, concentraron información de la distribución de partículas del regolito; el CP3 y CP6, representaron las características topográficas del talud y espesor de la cresta; el CP4 y CP7, mostraron información de la resistencia del regolito a la compactación y penetración; el CP5, reflejó información sobre la consistencia y la conductividad hidráulica del regolito y el CP8 mostró información sobre la consistencia y resistencia del regolito a la erosión.

Regresión lineal múltiple para pendiente del talud. Las covariables seleccionadas de cada CP se utilizaron para generar un modelo de regresión múltiple para predecir el avance lateral (Cuadro 3). El mejor modelo obtuvo ajuste de 84.02% y R² para la predicción de nuevas observaciones de 76.99%, expresándose de la siguiente forma:

AV = 0.2795 - 0.0102 * G0.15 ‑ 0.001617 * G0.075 + 0.000999 * M_talud + 0.001622 * Espesor

donde:

AV = avance lateral (m año^-1);

G0.15 y G0.075 = partículas de 0.15 y 0.075 mm, respectivamente (%);

M_talud = pendiente del talud (°);

Espesor = espesor de la cresta (cm).

De acuerdo al modelo de regresión, el avance lateral, disminuirá 0.0102 unidades por cada incremento en el tamaño de partícula de 0.15 mm, y 0.001617 con partículas de 0.075 mm, y aumentará 0.000999 y 0.001622 por cada unidad de incremento en la pendiente del talud y espesor de la cresta.

Finalmente, la Figura 13 muestra los valores observados del avance lateral y la predicción de los mismos. En esta se observa que los puntos por encima de la línea son sobrestimados por el modelo, mientras que los puntos por debajo son subestimados. De acuerdo al PBIAS (porcentaje de sesgo), el modelo obtenido tiende a sobrestimar 0.07%.

Figura 13: Valores originales contra valores predichos para el avance lateral del talud.  

Conclusiones