Español (pdf)
Articulo en XML
Referencias del artículo
Traducción automática
Enviar articulo por email
Citado por SciELO 
Similares en
SciELO 
Permalink
1. Introducción
Los suelos colombianos presentan gran variabilidad en sus propiedades físicas, químicas y biológicas, lo cual amerita la aplicación del análisis espacial, con el fin de realizar un manejo racional y eficiente de las áreas productivas. Cuando una propiedad del suelo presenta alta heterogeneidad, es importante determinar si existen tendencias o zonas que permitan describir espacialmente dicho comportamiento (Jaramillo Jaramillo et al., 2011). Para tal fin, el análisis de variabilidad espacial y los métodos de estimación o predicción de valores en puntos no evaluados a través del uso de métodos de interpolación y predicción geoestadística juegan un papel clave (Henríquez et al., 2005). Técnicas como la distancia inversa ponderada (interpolación) y kriging (predicción) permiten la elaboración de mapas con diferentes temáticas y finalidades (Mousavi et al., 2017).
La interpolación es una técnica que calcula el valor de una variable para un sitio no evaluado en el espacio, como una estimación a partir de los valores conocidos de esa variable en otras posiciones cercanas donde si hay registros (Bosque Sendra, 1992). La metodología de interpolación más empleada se conoce como el método de la distancia inversa ponderada (IDW, por sus siglas en inglés) (Bivand et al., 2013) y se define como un promedio ponderado del valor de colindantes al lugar de predicción (Hengl, 2007), aunque puede tener variaciones. También se puede recurrir a técnicas de predicción espacial, como el kriging, que además de predecir valores determina la varianza de la estimación realizada. Esta metodología se encuentra dentro del área de la estadística espacial, particularmente en la geoestadística. Se deben cumplir tres fases en un estudio geoestadístico para aplicar este método de predicción: primero, un análisis de la estacionariedad, en segundo lugar, una evaluación de la correlación espacial y la isotropía y, por último, la predicción correspondiente por kriging (Schabenberger y Gotway, 2005; Waller y Gotway, 2004).
En la agricultura, estas metodologías se han aplicado para la delimitación de áreas en el control de enfermedades, el análisis de la contaminación de suelos, la evaluación de aspectos en entomología, nematología y fertilidad de suelos, entre otras (Petersen et al., 1995). A nivel mundial existen varios trabajos relacionados con la aplicación de métodos de interpolación y predicción para el estudio de la variabilidad espacial usando propiedades tanto físicas como químicas del suelo (Henríquez et al., 2005; Moreno Araujo, 2011; Mousavi et al., 2017; Villatoro et al., 2008). Estos autores indican que en cuanto a la predicción espacial de características presentes en el suelo en ocasiones un método presenta una mejor estimación que el otro. En Colombia, los estudios como los de Vásquez Polo et al. (2010) y Jaramillo Jaramillo et al. (2011), demuestran que no existe un patrón único en la distribución espacial de dicha variabilidad y la misma puede ser influenciada por el manejo agrícola y la escala de muestreo (Bravo et al., 2011). Por lo tanto, es importante establecer cual método (interpolación o predicción) puede modelar de forma más adecuada las características físicas del suelo para un lugar determinado y así conocer cómo se distribuyen estas propiedades del suelo, puesto que inciden de manera importante sobre el crecimiento de las plantas y la disponibilidad de nutrientes; aspectos fundamentales para la sostenibilidad y actividades asociadas al manejo agrícola (Burbano Orjuela y Silva Mojica, 2010).
2. Materiales y Métodos
2.1. Localización
El estudio se realizó en la hoya del río Suárez [HDRS], localizada en los departamentos de Boyacá y Santander en los municipios de Chitaraque, Moniquirá, Togüí, San José de Pare, Santana, Barbosa, Chipatá, Güepsa, San Benito, Suaita y Vélez (Figura 1). El área de estudio se ubica entre los 1.000 y 2.000 m s.n.m., con variación de la precipitación anual entre los 1.200 y 2.700 mm, y humedad relativa del ambiente entre 65 y 85 %. El brillo solar anual promedio se encuentra entre las 1.900 y 2.000 horas (luminosidad de 5,7 horas día-1) y las temperaturas medias de 20 °C (Vergara et al., 2018).
2.2. Muestreo y caracterización del suelo
Los suelos predominantes en la HDRS se encuentran en paisajes de montaña y lomerío, en relieves moderados a fuertemente quebrados, con pendientes entre 7 y 50 %; algunos sectores están afectados, principalmente, por movimientos en masa y erosión, hay ligera presencia de pedregosidad en la superficie. Son suelos cuya taxonomía corresponden a Typic Dystrudepts, Oxic Dystrudepts, ChromicHapluderts, LithicUdorthents, entre otros, según la clasificación USDA (Soil Survey Staff, 2014), con texturas finas y moderadamente finas, ligera y moderadamente ácidos, profundos y superficiales, limitados por alta saturación de aluminio y presencia de abundantes fragmentos de roca. Tienen su origen en la formación geológica Arcabuco y Villeta. La formación Villeta presenta suelos con afloraciones de roca caliza, conformadas por pizarras arcillosas con buen contenido de materia orgánica (Guarín Abril, 2000). La actividad agrícola más importante que se desarrolla en esta región son la producción de panela a partir del cultivo de caña de azúcar y el procesamiento de la guayaba (Rodríguez et al., 2004).
En la zona de estudio se ubicaron 932 puntos georreferenciados siguiendo una forma de rejilla cuadrada, con un distanciamiento entre puntos de muestreo de 700 m, tanto a lo largo como a lo ancho (un punto cada 49 hectáreas). En cada sitio se tomaron muestras de suelo a una profundidad de 20 cm para análisis en laboratorio, donde se determinaron: 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 de las partículas de suelo [DMP] 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 del suelo [DA] por el método del volumen conocido y la densidad real [DR] por el método del picnómetro (Instituto Geográfico Agustín Codazzi [IGAC], 2006).
2.3. Análisis geoestadísticos
Los datos de las propiedades físicas del suelo fueron analizados con el software estadístico libre R® 3.6.3 (R Core Team, 2020). Las interpolaciones se generaron con el método de la distancia inversa ponderada [IDW], usando la función IDW de la librería gstat (Pebesma, 2004). Para el análisis geoestadístico se visualizó la estacionariedad de las propiedades físicas del suelo de la HDRS con respecto a los diferentes ejes espaciales, utilizando la función sm.regression de la librería sm (Bowman y Azzalini, 2014); se determinó si las variables estudiadas tienen un comportamiento anisotrópico, mediante semivariogramas direccionales (a 0°, 45°, 90° y 135°) con la función variog4 de la librería geoR (Ribeiro y Diggle, 2016); la construcción de los semivariogramas de cada variable, se calculó mediante la función variogram de la librería gstat (Pebesma, 2004) y con fit.variogram de la librería gstat (Pebesma, 2004) ajustando al modelo teórico de forma visual. Cuando el proceso fue estacionario, se usó el predictor KO (kriging ordinario) y en caso de presentarse tendencia espacial, se utilizó el predictor KU (kriging universal). En ambos casos se calcularon los predictores con la función krige de la librería gstat (Pebesma, 2004).
Con el resultado de las interpolaciones y predicciones de las propiedades físicas del suelo se elaboraron mapas temáticos por variable analizada, con ayuda de la función spplot de la librería sp (Pebesma y Bivand, 2005), utilizando los criterios de la Tabla 1 y, por último, se compararon para seleccionar la mejor alternativa con los criterios de validación cruzada usando la función krige.cv de la librería gstat (Pebesma, 2004). Los estadísticos de comparación usados fueron el error medio de predicción [EP], la raíz del error cuadrado medio de predicción [RCP] y el coeficiente de determinación [R2 ], los cuales se establecieron con la función criterio.cv de la librería geospt (Melo et al., 2015).
Tabla 1 Valores utilizados para definir mapas de zonificación de las siguientes propiedades: contenido de arenas (A), diámetro medio ponderado (DMP), densidad aparente (DA), densidad real (DR) y capacidad de retención de agua disponible (CRAD).Table 1. Values to define zoning maps of the following properties: sand content (A), weighted mean diameter (WMD), bulk density (DA), real density (DR) and available water retention capacity (CRAD).
3. Resultados y Discusión
Las medidas de tendencia central y variabilidad para A, DMP, DA, DR y CRAD se muestran en la Tabla 2. Los valores medios de las propiedades físicas del suelo evaluadas son semejantes a los mencionados en evaluaciones localizadas en los municipios de Santana y Güepsa (Ahumada González, 2015; Rueda Calier et al., 2015). La variabilidad de A y CRAD fue sobresaliente, como se muestra en los coeficientes de variación [CV] mayores al 35 %, lo cual puede deberse a la gran variedad de topografías y suelos que presenta la zona de estudio (IGAC, 2003, 2005). El contenido de arenas puede afectar la capacidad de retención de agua, puesto que a mayor contenido de A disminuye la CRAD (Montenegro y Malagón, 1990), mostrando que, si una de estas dos características es heterogénea, la otra presentó un comportamiento similar. Para estas variables se mencionan CV más bajos en el trabajo de Vásquez Polo et al. (2010), pero los muestreos usados fueron a longitudes inferiores a los empleados en este estudio.
Tabla 2 Medidas de tendencia central y variabilidad para contenido de arenas (A), diámetro medio ponderado (DMP), densidad aparente (DA), densidad real (DR) y capacidad de retención de agua disponible (CRAD). Calculados con datos obtenidos en la HDRS.Table 2. Measurements of central tendency and variability for sand content (A), weighted mean diameter (WMD), bulk density (DA), true density (DR) and available water retention capacity (AWR). Calculated with data obtained in the HDRS.
La DA y el DMP presentaron CV entre 16-35 %; estas propiedades son afectadas primordialmente por las labores agrícolas de preparación del suelo para el establecimiento de los cultivos, puesto que al intervenir el suelo con el uso de implementos como arado o rastrillo se busca que las plantas tengan un medio adecuado para su desarrollo (Manrique Estupiñan et al., 2000) dando como resultado la homogenización de la DA y DMP en los primeros centímetros de profundidad de la capa arable. La DR presentó valores de CV menores al 15 % los cuales son similares a los presentados por Monroy-Rodríguez (2017).
En la Figura 2, se observan las regresiones no paramétricas (con las bandas de referencia para la hipótesis de no efecto) para las características de estudio. Estos gráficos indican que para las variables no se presenta una tendencia espacial del valor de la media, semejante a lo expuesto por Ceddia et al. (2009) y Jaramillo Jaramillo (2011) que asumieron estacionariedad con distancias de muestreo inferiores a 20 metros (en suelo perteneciente al orden gleysol (frigid Aeric Haplaquept) y aluvial). La DA presenta un comportamiento opuesto al de las otras variables, indicando una tendencia marcada respecto a la ubicación geográfica (longitud). Mayores valores de DA pueden ser atribuidos entre otras, a la implementación de labores de preparación del suelo que permiten una mayor aireación transitoriamente durante las primeras fases del cultivo, y al paso del tiempo estas condiciones pueden disminuir por la intervención de precipitaciones, riego y uso de maquinaria agrícola que puede compactar el suelo y, por ende, incrementar la DA (Montenegro y Malagón, 1990). Otro factor determinante para la presencia de altos valores de DA, es la característica de las operaciones de carga y transporte de la caña realizadas en la región, en donde son empleados animales (mulares) que retiran el producto cosechado del área del cultivo, llevándolo hasta las unidades de procesamiento (trapiche). Esta actividad resulta en intenso pisoteo de la pezuña de los animales en el suelo (Medina et al., 2012), transmitiendo las cargas al suelo y resultando finalmente en la compactación de este.
Figura 2 Dispersión de los valores de contenido de arenas (A), diámetro medio ponderado (DMP), densidad aparente (DA), densidad real (DR) y capacidad de retención de agua disponible (CRAD), respecto a coordenadas planas (puntos negros). Curva de regresión ajustada por métodos Kernel (línea continua negra). Banda de confianza de no efecto (banda azul).Figure 2. Dispersion of the values of sand content (A), weighted mean diameter (WMD), bulk density (DA), real density (DR) and available water retention capacity (CRAD), with respect to plane coordinates (points blacks). Regression curve adjusted by Kernel methods (continuous black line). Confidence band of no effect (blue band).
En la Figura 3, los semivariogramas presentaron un comportamiento semejante en distancias por debajo de los diez kilómetros partiendo del origen, indicando que no tiene una estructura de correlación notable, igual a lo expuesto por el estudio realizado por Ceddia et al. (2009). Además, se evidencia que a distancias muy grandes los semivariogramas de las características evaluadas son diferentes. Sin embargo, este comportamiento es de esperarse, debido a que en esos sitios se encuentra una cantidad reducida de parejas de datos para realizar la estimación (Warrick, 2002).
Figura 3 Semivariogramas a diferentes direcciones (0, 45, 90 y 135 grados) para las variables contenido de arenas (A), diámetro medio ponderado (DMP), densidad aparente (DA), densidad real (DR) y capacidad de retención de agua disponible (CRAD).Figure 3. Semivariograms at different directions (0, 45, 90 and 135 degrees) for the variables sand content (A), weighted mean diameter (WMD), bulk density (BD), real density (DR) and retention capacity of water available (CRAD).
En la Figura 4, se presentan los modelos adaptados a los semivariogramas experimentales. Para todas las variables de estudio se presentan valores t2(pepita) con un cincuenta por ciento más alto que el valor de s2(meseta), eso puede indicar situaciones que inducen error, como el caso de presencia de incertidumbre por el uso de diferentes instrumentos de medición al realizar el muestreo, o también puede sugerir que las distancias utilizadas en el levantamiento de información no posibilitan un cálculo adecuado del semivariograma a distancias cercanas al origen. Además, el f(rango) en todas las propiedades físicas del suelo son inferiores a cinco mil metros, mostrando una dependencia espacial de baja a media.
Figura 4 Los puntos representan el semivariograma (experimental) y la línea continua es el ajuste de los modelos: exponencial (t2 = 108,88,s2 = 81,64 y f= 2.639,03) para el contenido de arenas (A), esférico (t2 = 0,34, s2 = 0,31 y f= 1.922,91) para el diámetro medio ponderado (DMP), exponencial (t2 = 0,021, s2 = 0,007 y f= 2.818,13) para densidad aparente, exponencial (t2 = 0,011, s2 = 0,005 y f= 1637,79) para densidad real y esférico (t2 = 5,55, s2 = 1,98 y f= 4.638,14) para la capacidad de retención de agua disponible (CRAD). Figure 4. The points represent the semivariogram (experimental) and the continuous line is the fit of the models: exponential (t2 = 108.88, s2 = 81.64 & f= 2,639.03)for the sand content (A), spherical (t2 = 0.34, s2 = 0.31 & f= 1,922.91) for the weighted mean diameter (DMP), exponential (t2 = 0.021, s2 = 0.007 and f= 2,818.13) for bulk density, exponential (t2 = 0.011, s2 = 0.005 & f= 1,637.79) for real and spherical density (t2 = 5.55, s2 =1.98 & f= 4,638.14) for the available water holding capacity (CRAD).
El modelo esférico se ajustó al DMP y CRAD, mientras que el exponencial se adaptó a A, DA y DR, esto se observa también en los estudios de Bravo et al. (2011), Barrios Maestre y Florentino de Andreu (2018), y Vásquez Polo et al. (2010). Sin embargo, estos autores muestran mayor dependencia espacial comparada con este estudio, posiblemente debido a las menores distancias de muestreo que utilizaron (inferiores a treinta metros). De acuerdo con todo lo anterior se utilizó la predicción kriging universal para la DA y kriging ordinario para el resto de las características físicas de suelo.
Según los gráficos de estimación y predicción mostrados en la Figura 5, los valores de A menores al 30 % se observan al norte de San Benito y en Chipatá, mientras que los mayores porcentajes se muestran hacia las zonas del occidente de Chitaraque, oriente de San José de Pare, el norte de Togüí y sur de Santana. Estas áreas con altos contenidos de arena pueden tener algunos problemas en cuanto a mayores probabilidades de que se presente lixiviación de sustancias como los agroquímicos, o limitaciones en el almacenamiento de agua que puede ser aprovechable para las plantas (Torres et al., 2017).
Figura 5 Zonificación para el contenido de arenas (A), diámetro medio ponderado (DMP), densidad aparente (DA), densidad real (DR) y capacidad de retención de agua disponible (CRAD), con la técnica de la distancia inversa ponderada (IDW) y usando kriging ordinario (KO) y universal (KU). Cada píxel corresponde a 6.561 m2.Figure 5. Zoning for sand content (A), weighted mean diameter (DMP), bulk density (DA), true density (DR) and available water retention capacity (CRAD), with the weighted inverse distance technique (IDW) and using ordinary (KO) and universal (KU) kriging. Each pixel corresponds to 6.561 m2,
La Figura 5 muestra la zonificación del DMP en la HDRS, donde las áreas moderadamente estables se encuentran en los alrededores de Chitaraque, San José de Pare, Santana y Togüí en el departamento de Boyacá. La estabilidad de agregados del suelo es mayor en las diferentes zonas del departamento de Santander. La buena condición de estabilidad de agregados del suelo en la HDRS a pesar de los problemas de compactación se puede atribuir al constante depósito de materia orgánica en la superficie del suelo como resultado de la actividad agrícola en los ciclos semipermanentes de la caña de azúcar. El aporte de materia orgánica contribuye a la agregación del suelo, y consecuentemente a la formación de macroagregados (Huang et al., 2010).
La DA del suelo puede relacionarse con la susceptibilidad a erosión, porque entre más alto sea el valor de esta propiedad los terrenos serán más propensos a erosionarse (Rivera-Posada et al., 2010). Los valores mayores de DA en la HDRS se ubican en los territorios de Chipatá, Güepsa y al sur de San Benito, zonas que tienden a ser más erosionables en comparación con el resto del área de estudio.
La DR tiene una relación con el contenido de materia orgánica del suelo, puesto que entre más bajo el valor de DR los terrenos tienden a tener más partículas orgánicas en su composición (Lal y Shukla, 2004). En la HDRS se observa que los sitios con menores valores de DR se localizan en Barbosa, Chitaraque, San José de Pare, Suaita y Togüí. En estos lugares, al establecer los cultivos, se puede requerir menor cantidad de aplicación de nutrientes a las plantas, puesto que el contenido orgánico es de gran importancia para el desarrollo fisiológico de las plantas.
La CRAD en la HDRS presenta baja capacidad de almacenamiento de agua con valores menores al quince por ciento. Se observa una mejor CRAD en más área del territorio del departamento de Boyacá, comparado con el departamento de Santander (Figura 5). Esto puede indicar que los cultivos tengan dificultades en su desarrollo y el suelo puede ser más propenso a erosionarse (Montenegro y Malagón, 1990) en algunas zonas de Santander como el oriente de Chitaraque.
El uso de las técnicas geoestadísticas permitió obtener zonificaciones con mejor definición, calidad visual de los sitios, al mismo tiempo que permiten observar la incertidumbre de las predicciones realizadas (Figura 6), en comparación con el método determinístico [IDW] para el caso de las propiedades de física de suelo medidas en la HDRS.
Figura 6 Gráfico de varianzas de la predicción por kriging de las variables: contenido de arenas (A), diámetro medio ponderado (DMP), densidad aparente (DA), densidad real (DR) y capacidad de retención de agua disponible (CRAD).Figure 6. Graph of variances of the kriging prediction of the variables: sand content (A), weighted mean diameter (DMP), bulk density (DA), real density (DR) and available water holding capacity (CRAD).
El uso de las predicciones por kriging presentaron valores más bajos en los estadísticos de la validación cruzada (error medio de predicción y la raíz del error cuadrado medio) y mayores valores en el coeficiente de determinación comparados con las estimaciones calculadas por IDW (Tabla 3), mostrando mejores resultados con el uso del kriging. Además, los gráficos de varianza de los métodos de predicción indican que la incertidumbre fue uniforme en la HDRS (Figura 6). La mayor variabilidad de predicción se presentó en el perímetro de la zona de estudio, causado por la menor cantidad de lugares de muestreo en estos sitios.
Tabla 3 Error medio de predicción (EP), raíz del error cuadrado medio de predicción (RCP) y coeficiente de determinación (R2 ) del análisis de validación cruzada para el kriging ordinario (KO), universal (KU) e interpolación por distancia inversa ponderada (IDW), realizados con las variables contenido de arenas (A), diámetro medio ponderado (DMP), densidad aparente (DA), densidad real (DR) y capacidad de retención de agua disponible (CRAD).Table 3. Mean prediction error (EP), root mean square error of prediction (RCP) and coefficient of determination (R2) of the cross-validation analysis for ordinary (KO), universal (KU) kriging and interpolation by inverse distance weighted (IDW), carried out with the variables sand content (A), weighted mean diameter (DMP), bulk density (DA), real density (DR) and available water retention capacity (CRAD).
4. Conclusiones
El uso de herramientas estadísticas como las metodologías de interpolación y predicción pueden modelar el comportamiento de las propiedades físicas del suelo en la HDRS y ayudan para la interpretación de posibles fenómenos relacionados con estas propiedades.
Las dos metodologías de predicción llevan a mapas similares en todos los casos, sin embargo, los valores de las predicciones por kriging son más homogéneos en los pixeles de las figuras, mostrando zonas con mayor uniformidad para las propiedades físicas del suelo evaluadas, lo cual fue validado con mejores resultados en los estadísticos de la validación cruzada para las estimaciones realizadas usando kriging.
Contribuciones de los autores
Ruy Edeymar Vargas Diaz: redacción – borrador original, investigación, curación de datos, análisis formal, visualización, redacción – revisión y edición.
Julio Ricardo Galindo Pacheco: conceptualización, redacción – revisión y edición.
Ramón Giraldo Henao: conceptualización, redacción – revisión y edición.
Viviana Marcela Varón Ramírez: redacción – revisión y edición.
Wilmar Alexander Wilches Ortiz: redacción – revisión y edición
Clara Viviana Franco Florez: redacción – revisión y edición.
