INTRODUCCIÓN

La erodabilidad (susceptibilidad a la erosión) del suelo es un fenómeno que se presenta porque las propiedades físicas, químicas y biológicas del suelo, facilitan la acción de los agentes erosivos externos, como lluvias y vientos. Para el estudio de estos procesos habitualmente se emplean técnicas estadísticas como el análisis de componentes principales, que permite establecer la posible correlación de las características del suelo sobre algún posible evento (^{Kooch, Jalilvand, Bahmanyar & Pormajidian, 2008}; ^{Ghaemi, Astaraei, Emami, Nassiri Mahalati & Sanaeinejad, 2014}).

El análisis del componente principal es una de las técnicas multivariadas que estructura un grupo de datos multivariados, creando nuevas variables denominadas componentes principales, que contienen la información de un conjunto original de variables y puede minimizar la dimensionalidad del evento a estudiar (^{Pla, 1986}; ^{Srivastava, 2002}; ^{Anderson, 2003}; Aldás & Jiménez, 2017).

Los objetivos de un análisis del componente principal son: generar nuevas variables que puedan expresar la información contenida en el conjunto original de datos, reducir la dimensionalidad del problema de estudio y eliminar, cuando sea posible, algunas de las variables originales si estás aportan poca información (^{Pla, 1986}; ^{Srivastava, 2002}; ^{Aldas Manzano & Jiménez, 2017}).

Aldas y Jiménez (2017), describen la primera componente principal como la combinación lineal de las p variables originales como:

(1)

Se pueden generar componentes principales, y a estos se les puede realizar un análisis geoestadístico y así poder interpretar de manera conjunta la variación espacial de las propiedades físicas de un territorio.

El objetivo del presente trabajo consistió en identificar zonas diferenciales de suelo en la hoya del río Suárez, con base en el análisis espacial de las propiedades físicas, para proponer en un futuro prácticas apropiadas de conservación, lo cual resulta muy importante dado que actualmente en dicha hoya se hace un uso intensivo de la tierra en agricultura y ganadería. También se busca evaluar si el análisis de componentes principales permite una reducción de las variables observadas que facilite la zonificación sin mayor pérdida de información, e inclusive, favorezca una mayor claridad en la interpretación.

MATERIALES Y MÉTODOS

El estudio se localizó entre los departamentos de Boyacá (municipios de Chitaraque, Moniquirá, Togüí, San José de Pare, Santana) y Santander (municipios de Barbosa, Chipatá, Güepsa, San Benito, Suaita y Vélez), en Colombia (Figura 1). En general, esta zona presenta topografías de montaña ondulada y quebrada, con pendientes catalogadas de medianas a altas (^{Rodríguez, García, Roa Díaz & Santacoloma, 2004}). En esta área se encuentran terrenos a alturas sobre el nivel del mar entre los 1200 y 1900 msnm con temperaturas promedio entre los 18 a 23 grados centígrados y precipitaciones anuales entre los 1200 y 2700 milímetros al año (^{Insuasty Burbano, Manrique & Palacio, 2003}). Los suelos de la región se dedican especialmente al cultivo de caña, cuya producción se destina a la elaboración de panela, un alimento dulce de alto consumo en Colombia (^{Vergara, Rodríguez & Alarcón, 2018}). En menor escala también se presenta la ganadería, que también tiene un impacto negativo en las propiedades del suelo y su conservación.

Se ubicaron en total 932 puntos georreferenciados en un área de 40.000 hectáreas en la zona de estudio, siguiendo una de grilla rectangular, con una separación entre puntos de muestreo de 700 metros, tanto a lo largo como a lo ancho.

En cada punto, en el perfil de los primeros 20 cm de profundidad, se tomaron muestras de suelo para el análisis de laboratorio donde se determinó el porcentaje de arena en la fracción sólida del suelo (A) por el método del hidrómetro de Bouyoucos, el diámetro medio ponderado (DMP) de los agregados por el método de Yoder, la capacidad de retención de agua disponible (CRAD) por el método de succión con cámaras de presión, la densidad aparente (DA) del suelo por el método del volumen conocido y la densidad real (DR) por el método del picnómetro. En los análisis descritos se utilizó la guía para métodos analíticos del laboratorio de suelos descritos en el Instituto Geográfico Agustín Codazzi (^{IGAC, 2006}).

Primero, las variables A, DMP, DA, DR, CRAD se cruzaron por pares para determinar los coeficientes de la correlación lineal (Pearson), mediante la función cor.test de la librería stats (^{R Core Team, 2016}). Luego los componentes principales (CP) se calcularon con la función princomp de la librería stats (^{R Core Team, 2016}). Para establecer la conveniencia de usar 2 o 3 CP (las nuevas variables a espacializar), se tuvo en presente que el porcentaje de varianza explicada acumulada fuera al menos del 60% y así disminuir la dimensionalidad de las variables a estudiar.

Con los CP, se realizó el análisis geoestadístico, tal como se detalla a continuación. En primer lugar, se observó la tendencia espacial de los CP. La estacionariedad de los CP con respecto a su localización geográfica se realizó mediante la función sm.regression de la librería sm (^{Bowman & Azzalini, 2014}). Con esto se estimó un modelo de regresión no paramétrico Kernel y un gráfico de dispersión en dos dimensiones, con ejes Y (variable evaluada), X (latitud o longitud) y una banda de confianza de no efecto de la superficie de regresión estimada (Bowman & Azzalini, 2014). Para evaluar la anisotropía de los CP, se construyeron semivariogramas direccionales en 4 orientaciones a 0, 45, 90 y 135 grados, con la función variog4 de la librería geoR (^{Ribeiro & Diggle, 2016}).

Posteriormente, se elaboraron los semivariogramas para las CP, calculando el semivariograma experimental para cada CP, mediante la función variogram de la librería gstat y se ajustó un modelo teórico con la función fit.variogram de la librería gstat, ambos según ^{Pebesma (2004}). También, se estimó el estadístico de ajuste (), para seleccionar el mejor modelo.

Si los gráficos de estacionariedad y los valores de pepita (), meseta () y rango ( de los semivariogramas indican que las CP son estacionarias se implementa predicción por kriging ordinario (KO) y si se observa tendencia espacial, se utiliza kriging universal (KU). Para esto se usó la función krige de la librería gstat (^{Pebesma, 2004}) que calcula los predictores de kriging ordinario y kriging universal, con sus respectivas varianzas. Con los valores estimados por KO o KU se elaboran los mapas con la función spplot de la librería sp (^{Pebesma & Bivand, 2005}) para cada componente principal.

RESULTADOS Y DISCUSIÓN

En la Tabla 1 se presentan los coeficientes de correlación de las propiedades de física de suelo estudiadas, los cuales en general, fueron significativamente diferentes de cero (valor p<0,05, Ho: Rho=0), aunque el grado de correlación fue variable. Dos excepciones a esto se observan entre el diámetro medio ponderado (DMP) de los agregados y la densidad aparente (DA) y entre DMP y la capacidad de retención de agua disponible (CRAD) (Tabla 1), las cuales no fueron significativas, indicando que posiblemente el DMP de los agregados no está afectado por las características físicas consideradas; por el contrario, es posible que las variaciones del DMP de los agregados se puedan explicar por alguna variable no analizada, como la materia orgánica o el contenido de arcillas, entre otras (^{Burbano Orjuela & Silva Mojica, 2010}; ^{Montenegro & Malagón, 1990}).

Como la mayoría de las correlaciones presentadas para las propiedades físicas de suelo son significativas, es viable un análisis conjunto de las variables de estudio. Es de resaltar que la significancia de los coeficientes de correlación mostrados en la Tabla 1 debe interpretarse con precaución, ya que las características evaluadas pueden estar autocorrelacionadas espacialmente, pero, considerando la gran distancia entre puntos (700 m), es una situación poco probable.

La correlación directa (positiva) se observa en las parejas de DA - DR y A - CRAD, mientras que en las otras asociaciones cuando son significativas se presentan de forma inversa (negativa). ^{Vargas, Serrato y Trujillo (2015}), establecen que el contenido de arena del suelo tiene una relación opuesta con su capacidad para almacenar agua, pero en este estudio no se presentó este comportamiento, debido posiblemente al efecto de los contenidos de las otras fracciones texturales (limos y arcillas) o de la materia orgánica que contienen los territorios de la Hoya del río Suárez.

La asociación entre DR - DA es mayor a 0,5, similar a lo encontrado por ^{Salamanca y Sadeghian (2006}), demostrando con ello la tendencia según la cual los lugares con valores altos de DR también presentan una DA alta.

La correlación entre DA - CRAD (inversa), muestra que a menor DA existe más CRAD en el suelo, lo cual indica que en los suelos con baja DA y, por consiguiente, mayor espacio poroso, pueden contener y retener más agua (^{Blanco Sepúlveda, 2009}; ^{Varón, Camacho & Ordoñez, 2014}). Este comportamiento inverso se observa también entre A - DA, ya que un suelo con mayor espacio poroso por su contenido de A tiene menor DA, pero si el A es alto la estabilidad del suelo baja y como resultado el DMP de los agregados presenta valores bajos (^{Vargas et al., 2015}; ^{Novillo et al., 2018}).

Se utilizó la ecuación (1) para calcular las nuevas variables con el análisis de componentes principales. En las Tablas 2 y 3 se observan los valores y vectores propios generados con el análisis multivariado. El primer componente principal explica el 40,5% de la variabilidad de las propiedades físicas de suelo medidas en la Hoya del río Suárez, el cual está asociado directamente con el A, DMP de los agregados y CRAD e inversamente con el DR y DA (Tabla 2). El efecto del aumento de la DA disminuye los valores de las propiedades físicas que son importantes (DMP de los agregados y CRAD) para el desarrollo fisiológico de las plantas. El componente principal 1 posiblemente esté asociado inversamente a los efectos de la compactación de los terrenos en la zona de estudio. Por lo tanto, valores altos del componente principal 1 se pueden asociar con suelos sueltos, y valores bajos o negativos con suelos más compactados.

Por otro lado, el componente principal 2 recoge cerca del 22% de la varianza de las características físicas estudiadas, y recibe aportes positivos del valor de DMP de los agregados, y negativos por el valor de A. Este componente evidencia la relación inversa entre la estabilidad del suelo y el contenido de A en los terrenos. Cuando un área presenta valores bajos de DMP en los agregados son más propensos a la erodabilidad, por lo cual, este componente resulta importante en la detección de posibles zonas susceptibles a la erosión. Valores altos del componente principal 2 se pueden asociar a suelos resistentes a la erosión, y valores bajos o negativos a suelos susceptibles. Los otros componentes principales no son considerados debido a que su aporte es bajo en la explicación de la variabilidad de las variables físicas consideradas.

La Figura 2, muestra la existencia de observaciones atípicas, en especial, valores altos del componente principal 1, los cuales pueden presentar condiciones favorables para el crecimiento radicular de las plantas debido a los valores bajos de DA y DR y altos en cuanto a CRAD, presentes en el suelo. Además, el dato negativo del componente principal 1, presentó valores altos en DA y DR, o sea, mayor compactación. Mientras que en el del componente principal 2, los valores negativos se pueden asociar a problemas de estabilidad del suelo debido a los valores bajos en el DMP de los agregados.

La Figura 3 muestra los modelos de regresión no paramétricos de los componentes 1 y 2, respecto a su ubicación espacial. Se observa que las curvas de componente principal 1 se encuentran ligeramente fuera de las bandas de confianza de no efecto, mientras que el componente principal 2 no muestra una clara tendencia espacial y puede asumirse estacionariedad para el análisis geoestadístico.

En la Figura 4 se observan que los semivariogramas son semejantes a distancias menores a 10 Km. Esto indica que los componentes principales evaluados no tienen una variación considerable en su estructura de correlación. Se presentan diferencias en los semivariogramas a distancias grandes, pero esto es de esperarse debido al número pequeño de parejas que se tiene en estos lugares para hacer la estimación del semivariograma (^{Warrick, 2002}).

Se aprecia (Figura 5) que los valores de (pepita) son mayores o cercanos al valor de (meseta), mostrando que el espaciamiento utilizado en el levantamiento de información (700 m) no permite estimar este parámetro cerca al origen, debido a que no hay datos a distancias menores a 700 m, aunque también puede indicar problemas al realizar las mediciones, como alguna variabilidad instrumental en la toma o análisis de las muestras. Además, los (rangos) son menores a 7000 m, sugiriendo baja o media dependencia espacial. Teniendo en cuenta lo anterior y la Figura 3 se utilizó predicción kriging universal para el componente principal 1 y ordinario para el componente principal 2.

La Figura 6 muestra los mapas generados por kriging para los componentes principales 1 y 2. Con respecto al componente principal 1, sobresalen dos zonas, una con valores negativos (en los municipios de Chipatá, Barbosa, Vélez, Güepsa, San Benito y parte de San José de Pare y Santana) que indica valores altos de DR, DA y bajos de CRAD. Estos sitios son menos productivos para el establecimiento de cultivos, debido a que al tener valores de DA altos, sugiere un riesgo de compactación, con porosidad reducida que afecta las condiciones de retención de humedad del suelo, limitando un desarrollo adecuado de las raíces de las plantas (^{Salamanca & Sadeghian, 2006}). Mientras que las áreas con valores bajos (valores positivos) de DA, DR y altos CRAD (en los municipios de Chitaraque, Togüi y parte de San José de Pare y Santana), estos lugares son influenciados por los valores atípicos mostrados en la Figura 2 y son regiones de mayor productividad agronómica, caracterizadas por tener estabilidad estructural de agregados de moderadamente estables a muy estables (^{IGAC, 2005}).

En la Figura 6, los valores con altos contenidos de A y menor DMP de los agregados (más negativos), se establecen al sur de San Benito y en la parte oriental de los municipios de San José de Pare y Chitaraque. Estas áreas pueden ser susceptibles a erosionarse, puesto que, al tener un mayor A y baja estabilidad, se fomenta el desprendimiento de los agregados del suelo (^{Rivera, Lal, Amézquita, Mesa & Chaves, 2010}). Mientras que los territorios hacia el sur de Chipatá y oriente de San Benito presentan valores por debajo del promedio en A con suelos más estables, sugiriendo que en estos lugares hay menor riesgo de procesos erosivos.

Los mapas de las varianzas para los componentes principales 1 y 2 se observan en la figura 7, donde se ilustra que la incertidumbre es homogénea en la mayor parte de la Hoya del río Suárez, y que la variabilidad más alta se ciñe al contorno de la región de estudio; esto se debe al efecto de borde por el menor número de sitios de muestreo que se puede aplicar para el análisis.

CONCLUSIONES

Con el uso de estadísticas multivariadas y geoestadística se logró definir dos componentes que resumen las propiedades físicas del suelo y permiten zonificar la región de la Hoya del río Suárez: la compactación y la susceptibilidad a la erosión (erodabilidad). La zonificación resultante permitirá concretar las acciones necesarias para favorecer la conservación del suelo y reducir el impacto de las actividades agrícolas y ganaderas de la región.

Se pudo identificar el caso del sur de San Benito y en la parte oriental de los municipios de San José de Pare y Chitaraque, donde se presentan áreas susceptibles a erosionarse por sus altos contenidos de arena y baja estabilidad estructural de los agregados.

Además, en los municipios de Chitaraque, Togüi y parte de San José de Pare y Santana se encuentran zonas donde la productividad agronómica es mayor y se constatan unas condiciones de física del suelo, que se reflejan en el componente 1 del modelo analizado, con una menor compactación, dados los valores bajos en cuanto a densidad aparente y real y además una mayor estabilidad estructural de los agregados.