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

Figura 1
Localización del área de estudio. Izquierda: ubicación de la HDRS en Colombia. Derecha: esquema de sitios de muestreo.

Figure 1. Location of the study area. Left: Location of the HDRS in Colombia. Right: Scheme of sample sites.

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 [R²], los cuales se establecieron con la función criterio.cv de la librería geospt (Melo et al., 2015).

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.

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 t²(pepita) con un cincuenta por ciento más alto que el valor de s²(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 (t² = 108.88, s² = 81.64 & f= 2,639.03)for the sand content (A), spherical (t² = 0.34, s² = 0.31 & f= 1,922.91) for the weighted mean diameter (DMP), exponential (t² = 0.021, s² = 0.007 and f= 2,818.13) for bulk density, exponential (t² = 0.011, s² = 0.005 & f= 1,637.79) for real and spherical density (t² = 5.55, s² =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 m²,

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.

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.

Implicaciones éticas

Los autores declaran que no existen implicaciones éticas.

Conflicto de interés

Los autores declaran que no existen conflictos de interés financieros o no financieros que podrían haber influido en el trabajo presentado en este artículo.

Agradecimientos

Al Ministerio de Agricultura y Desarrollo Rural de Colombia (MADR) por la financiación y a la Corporación Colombiana de investigación Agropecuaria (Agrosavia) por la ejecución del proyecto “Recomendaciones técnicas preliminares de manejo de suelos en ladera para el sistema de producción de caña panelera en la Hoya del Río Suárez (HDRS)” el cual pertenece al macroproyecto “Generación de recomendaciones agronómicas y nuevas tecnologías en la producción de panela en Colombia” de la red de innovación de transitorios.

Referencias

Ahumada González, J. A. (2015). Efectos de la tecnología utilizada en la producción panelera en las transformaciones ambientales de la Hoya del Río Suárez. Universidad Nacional de Colombia. https://repositorio.unal.edu.co/handle/unal/56732

Barrios Maestre, R., y Florentino de Andreu, A.. (2009). Variabilidad espacial de las propiedades físicas de dos suelos cultivados con palma aceitera en el estado Monagas, Venezuela. Revista Cientifica UDO Agrícola, 9(4), 912-924. http://udoagricola.orgfree.com/V9N4UDOAg/V9N4Barrios912.htm

Bivand, R. S., Pebesma, E., y Gomez-Rubio, V. (2013). Disease Mapping. En Applied Spatial Data Analysis with R (pp. 319-361). Springer. https://doi.org/10.1007/978-1-4614-7618-4_10

Bosque Sendra, J. (1992). Sistema de Información Geográfica (I y II ed.). Ediciones Rialp.

Bowman, A., y Azzalini, A. (2014). R package 'sm': nonparametric smoothing methods (version 2.2-5.5). http://www.stats.gla.ac.uk/~adrian/sm

Bravo, C., Cabrera, J., Carvallo, M., Cánchica, H., Ramírez, E., y González, I. (2011). Variabilidad espacial de propiedades físicas y químicas del suelo y su relación con el rendimiento de frijol (Vigna unguiculata L.) en una parcela experimental de los llanos centrales venezolanos. En J. Martínez Fernández y N. Sánchez Martín (eds.), Estudios en la zona no saturada del suelo (pp. 79-84). Vol 10. ZNS’11. Universidad de Universidad de Salamanca. http://www.zonanosaturada.com/zns11/publications.html

Burbano Orjuela, H. y Silva Mojica, F. (2010). Ciencia del Suelo: Principios Básicos. Sociedad Colombiana de la Ciencia del Suelo.

Ceddia, M. B., Vieira, S. R., Villela, A. L. O., Mota, L. dos S., Anjos, L. H. C. dos, y Carvalho, D. F. de. (2009). Topography and spatial variability of soil physical properties. Scientia Agricola, 66(3), 338-352. https://doi.org/10.1590/S0103-90162009000300009

Guarín Abril, J. B. (2000). Evaluación de cuatro cebos atrayentes y cinco ingredientes activos utilizados para el control de la hormiga loca Paratrechina fulva, Mayr (Hymenoptera: formicidae) en caña panelera, en la hoya del rio Suarez.Universidad Nacional de Colombia.

Hengl, T. A. (2007). Practical Guide to Geostatistical Mapping of Environmental Variables. Office for Official Publications of the European Communities. https://publications.jrc.ec.europa.eu/repository/handle/JRC38153

Henríquez, C., Killorn, R., Bertsch, F., y Sancho, F. (2005). La geostadística en el estudio de la variación espacial de la fertilidad del suelo mediante el uso del interpolador kriging. Agronomía Costarricense, 29(2), 73-81. http://www.mag.go.cr/rev_agr/v29n02_073.pdf

Huang, L., Wang, C. Y., Tan, W. F., Hu, H. Q., Cai, C. F., y Wang, M. K. (2010). Distribution of organic matter in aggregates of eroded Ultisols, Central China. Soil and Tillage Research, 108(1-2), 59-67. https://doi.org/10.1016/j.still.2010.03.003

Instituto Geográfico Agustín Codazzi [IGAC]. (2003). Estudio General de Suelos y Zonificación de Tierras. IGAC.

Instituto Geográfico Agustín Codazzi [IGAC]. (2005). Estudio General de Suelos y Zonificación de Tierras del Departamento de Boyacá. IGAC.

Instituto Geográfico Agustín Codazzi [IGAC]. (2006). Métodos analíticos del laboratorio de suelos (6ª ed.). IGAC.

Jaramillo Jaramillo, D. F., Anaya Gómez, M. L., Restrepo Moná, C. A., González Sánchez, H. A., y Álvarez Mejía, F. (2011). Variables físicas que explican la variabilidad de suelo aluvial y su comportamiento espacial. Pesquisa Agropecuaria Brasileira, 46(12), 1707-1715. https://doi.org/10.1590/S0100-204X2011001200017

Lal, R., y Shukla, M. K. (2004). Principles of Soil Physics. CRC Press. https://doi.org/10.4324/9780203021231

Manrique Estupiñan, R., Insuasty Burbano, O., Mora Padilla, C. J., Rodrigez Borray, G., Blanco Suárez, R., Mejía Flores, L. A., Pinto, J., y Sandoval Sandoval, G. (2000). Manual de caña de azúcar para la producción de panela. Corporación colombiana de investigación agropecuaria [CORPOICA] y Fondo Nacional de la Panela [FEDEPANELA]. http://hdl.handle.net/20.500.12324/18504

Medina, C., Camacho-Tamayo, J. H., y Cortés, C. A. (2012). Soil penetration resistance analysis by multivariate and geostatistical methods. Engenharia Agrícola, 32(1), 91-101. https://doi.org/10.1590/S0100-69162012000100010

Melo, C., Santacruz, A., y Melo, O. (2015). geospt: Geostatistical Analysis and Design of Optimal Spatial Sampling Networks (version 1.0-2). https://cran.r-project.org/web/packages/geospt/index.html

Monroy-Rodríguez, F. L., Álvarez-Herrera, J. G., y Alvarado-Sanabria, Ó. H. (2017). Distribución espacial de algunas propiedades físicas del suelo en un transecto de la granja Tunguavita, Paipa. Revista U.D.C.A Actualidad & Divulgación Científica, 20(1), 91-100. https://doi.org/10.31910/rudca.v20.n1.2017.66

Montenegro, H., y Malagón, D. (1990). Propiedades Físicas de los Suelos. Instituto Geográfico Agustín Codazzi.

Moreno Araujo, M. (2011). Variabilidad espacial de las propiedades físicas, químicas y biológicas de dos suelos agrícolas de la provincia de Quebec, Canadá. Revista de la Facultad de Agronomía de la Universidad del Zulia, 28(2), 223-241. https://produccioncientificaluz.org/index.php/agronomia/article/view/26882

Mousavi, S. R., Sarmadian, F., Dehghani, S., Sadikhani, M. R., y Taati, A. (2017). Evaluating inverse distance weighting and kriging methods in estimation of some physical and chemical properties of soil in Qazvin Plain. Eurasian Journal of Soil Science, 6(4), 327-336. https://doi.org/10.18393/ejss.311210

Pebesma, E. J. (2004). Multivariable geostatistics in S: the gstat package. Computers & geosciences, 30(7), 683-691. https://doi.org/10.1016/j.cageo.2004.03.012

Pebesma, E. J., y Bivand, R. S. (2005). Classes and methods for spatial data in R. R News, 5(2), 9-13. http://cran.r-project.org/doc/Rnews/Rnews_2005-2.pdf

Petersen, G. W., Bell, J. C, McSweeney, K., Nielsen, G. A., y Robert, P. C. (1995). Geographic Information Systems in Agronomy. Advances in Agronomy, 55, 67-111. https://doi.org/10.1016/S0065-2113(08)60538-6

R Core Team. (2020). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. https://www.r-project.org/

Ribeiro Jr., P. J., y Diggle, P. J. (2016). geoR: Analysis of Geostatistical Data (versión 1.7-5.2). http://www.leg.ufpr.br/geoR/

Rivera-Posada, J. H., Lal, R., Amézquita-Collazos, É., Mesa-Sánchez, Ó., y Chaves-Córdoba, B. (2010). Predicción de la erodabilidad en algunos suelos de ladera de la zona cafetera colombiana. Cenicafé, 61(4), 344-357. http://hdl.handle.net/10778/508

Rodriguez, G., Garcia, H., Roa Díaz, Z., y Santacoloma, P. (2004) Producción de panela como estrategia de diversificación en la generación de ingresos en áreas rurales de América Latina. Food and Agriculture Organization of the United Nations.

Rueda Calier, F., Peñaranda Mallungo, L. A., Velásquez Vargas, W. L., y Díaz Báez, S. A. (2015). Aplicación de una metodología de análisis de datos obtenidos por percepción remota orientados a la estimación de la productividad de caña para panela al cuantificar el NDVI (índice de vegetación de diferencia normalizada). Ciencia y tecnologia agropecuaria, 16(1), 25-40. https://doi.org/10.21930/rcta.vol16_num1_art:377

Schabenberger, O., y Gotway, C. A. (2005). Statistical methods for spatial data analysis: Texts in statistical science. Chapman and Hall/CRC.

Soil Survey Staff. (2014). Keys to Soil Taxonomy (12^th ed.). USDA-Natural Resources Conservation Service. https://www.nrcs.usda.gov/resources/guides-and-instructions/keys-to-soil-taxonomy

Torres, D., Álvarez, J., Contreras, J., Henríquez, M., Hernández, W., Lorbes, J., y Mogollón, J. P. (2017). Identificación de potencialidades y limitaciones de suelos agrícolas del estado Lara, Venezuela. Bioagro, 29(3), 207-218. http://www.ucla.edu.ve/bioagro/Rev29(3)/6.%201639.pdf

Vargas Díaz, R. E. (2019). Zonificación de la Hoya del Río Suárez por propiedades físicas del suelo, para el cultivo de caña panelera. Universidad Nacional de Colombia. https://repositorio.unal.edu.co/handle/unal/76943

Vásquez Polo, J. R., Baena Garcia, D., y Menjivar Flores, J. C. (2010). Variabilidad espacial de propiedades físicas y químicas en suelos de la granja experimental de la Universidad del Magdalena (Santa Marta, Colombia). Acta Agronómica, 59(4), 449-456. https://revistas.unal.edu.co/index.php/acta_agronomica/article/view/20129

Vergara, R., Rodríguez, L., y Alarcón, Y. (2018). Estudio de la caracterización de la actividad productiva del sector de la caña panelera en la hoya del río Suárez. NOVUM, 1(8), 32-48. https://revistas.unal.edu.co/index.php/novum/article/view/69931

Villatoro, M., Henríquez, C., y Sancho, F. (2008). Comparación de los interpoladores IDW Y Kriging en la variación espacial de pH, Ca, CICE y P del suelo. Agronomía Costarricense, 32(1), 95-105. https://revistas.ucr.ac.cr/index.php/agrocost/article/view/6773

Waller, L. A., y Gotway, C. A. (2004). Applied spatial statistics for public health data. John Wiley & Sons. Inc.

Warrick, A. W. (2002). Soil physics companion. CRC press.

Notas de autor

rvargas@agrosavia.co

Información adicional

e-location: e4118