INTRODUCTION

In conventional agriculture, a large percentage of crop yield variability is due to the differentials that may exist in the soil properties (^{Bramley, 2009}). This variability can be analyzed from different perspectives, scales based on both the intervention objectives and the level of decision required. It also depends on natural disturbances and human actions, especially in agricultural practices (^{Fraterrigo & Rusak, 2008}).

Mapping, the use of cartographic tools, is a valuable approach to adapt work processes to landscape changes. Mapping is referred to as a tool because it is an important means to strengthen decision making, in this case, to improve agricultural production through nutrient use efficiency (NUE) (^{Murillo, 2006}).

The recognition of the role of fertilizers in increasing agricultural production, and consequently, in the production of food, fiber, and even energy, contrasts severely with the negative nature of the information that is currently being released on fertilizers use (^{Murillo, 2006}).

For agricultural sciences, information at different scales is often used to support decision making with respect to diagnosis, monitoring, and predictions related to the management of productive systems (^{Pachepsky & Hill, 2017}). At the research level, the spatial variability of physical and chemical properties of soils at the field scale has been more widely considered than the patterns of such distribution (^{Reza et al., 2017}). In this sense, the modeling of these properties will allow, among others, to estimate the distance in which different soil samples are independent as well as to improve the sampling process either for planning activities in a traditional way or under a precision agriculture scheme (^{Peralta & Costa, 2013}; ^{Peralta et al., 2015}).

In terms of production and academic contributions, various activities are carried out in most production units. Regarding Botana Experimental Farm of Universidad de Nariño - Colombia, is a site where different activities take place, including growth, development, and productivity assessment of certain crops from the Andean region. In each of these, physical and chemical studies of the soils are generally conducted to determine their fertility level, which benefits the crop or system under study.

However, in most of the agricultural production units in the region, the results of soil analysis are not very useful over time. One of the reasons that account for this is the fact that the analyses are not properly stored or processed in databases, which hinders production planning in the long run. If properly used, these databases or mapping, which depict the spatial distribution of the physical and chemical properties of soils, would make it possible to estimate the productive potential of the farm. This would result in appropriate agricultural, forestry, and livestock planning, without the need to repeat the already executed processes.

According to ^{Álvarez et al. (2015)}, today, several research centers are using geographic information systems (GIS) to elaborate soil maps, which are of vital importance for planning and correct soil management. In turn, appropriate management increases the competitive capacity of soils and the possibilities to better channel agricultural and agroforestry programs and projects. Using GIS, soil fertility zonation studies facilitate adequate planning in terms of agricultural use of soil and its optimal management (^{Álvarez et al., 2015}). Similarly, these tools allow making timely recommendations that would permit the optimization of resources and the definition of productive systems that become a very important input for farm management and thus take measures that promote the implementation of Precision Agriculture in these research centers.

Site-specific fertilization involves the assessment of the spatial relationships between soil properties and crop yield. In this sense, the spatial distribution of physical, chemical, and biological parameters of soil has been observed to affect crop yield (^{Machado et al., 2000}). However, the spatial dependence of soil biological parameters has been little explored (^{Moreno, 2011}). One of the most important uses of GIS is spatial analysis, especially the use of interpolations of different types of variables. In agriculture, the use of this tool makes it possible to analyze the variability of different characteristics in the landscape such as soil, diseases, and pests, among others (^{Clay et al., 2007}). This undoubtedly helps to quantify the impact of this variation on production and the possible management guidelines required to optimize farm yields (^{Bertsch et al., 2002}).

Considering the above, this research aimed to assess the implementation of soil fertility zoning at Botana Experimental Farm of Universidad de Nariño as a tool to contribute to productive processes, research, and academic practices of the agriculture programs of the University and productive systems of the region.

MATERIALS AND METHODS

Location. The present study was carried out at Botana Experimental Farm of Universidad de Nariño - Department of Nariño, at 2820 m asl, 12.4ºC and an average annual rainfall of 694mm. It is located in a low mountain dry forest life zone (bs-MB) at 1°09'30.8"N; 77°16'31.8"W.

Fieldwork and laboratory methodology. Through the use of a Geographic Information System, a sampling grid of 0.5 hectares was established, for a total of 106 soil samples in the productive areas of Botana Experimental Farm. That is, the samplings were drawn from an overall area of 55 hectares (productive area). This aspect guaranteed the significance of the study area about the established interpolation method and the areas dedicated to production within the system.

Once the sampling grid was defined, the sampling points were identified and 15 soil subsamples were taken in each grid using an Edelman auger. These subsamples were mixed into a single sample, placed in plastic bags and labeled. The samples were sent to the laboratory of Universidad de Nariño to determine their physical-chemical characteristics (Table 1 and Table 2).

Statistical and geostatistical analysis. The data were descriptively analyzed and a principal component analysis was performed using the SPAD software. The geostatistical analysis was carried out with ArcGIS 10 software (^{ArcGIS, 2010}). The semivariogram was obtained for each characteristic and through an iterative process where the active lag and the step were modified; the theoretical model of best fit was established, taking into account as parameters of decision the coefficient of determination (R²) and the sum of squares of the residuals (RSS). Subsequently, a cross-validation analysis was carried out, of importance in the estimation of values in unsampled sites, through the point Kriging interpolation method, because the samples come from points and not from combinations or mixtures (^{Balzarini, 2014}). This interpolation was the basis for the construction of thematic maps that allow visualizing the spatial variability of each of the physical and chemical properties analyzed. In the mentioned software, spatial distribution maps were also made for each variable, subsequently reclassifications of the results were performed to generate other fertility distribution maps for each variable in the three ranges, according to the values for the Andean region (Table 3).

RESULTS AND DISCUSSION

General statistics. When analyzing the soil fertility traits of the total area samples, a high percentage of variation was found, which suggested that there was a high intrazonal variability (Table 4).

Taking into account some classifications for this type of results in similar studies, such as the one carried out by ^{Vásquez (2009)}, we can classify the variability of the data as follows:

Relatively low or homogeneous variations are only presented for pH and bulk density (DBH), as shown in Table 5, whose similar results were obtained in a study of the variability of soil fertility of the experimental farm of the University of Magdalena (^{Vásquez, 2009}). There is also agreement with the results found by ^{Henríquez (2013)}, ^{Londoño & Moreno (2014)} and ^{Ibarra et al. (2009)}.

The traits that presented extremely heterogeneous variability were P, K, Cu, Zn, N and S. These results also coincide with what was found by ^{Londoño & Moreno (2014)}. Likewise, ^{Gutiérrez et al. (2010)}, ^{Moreno (2011)} and ^{Alesso et al. (2017)} found in a soil fertility distribution study that available phosphorus and potassium showed high coefficients of variation. ^{Moreno (2011)} argues that although pH is a factor that influences the availability of phosphorus in the soil, he did not find a relationship between the high variability shown by this element and the low coefficients of variation observed in pH.

The high value of the coefficient of variation for P (158%) agrees with several works reported in the literature. In this sense, ^{Molina et al. (2005)} cite similar results. While the following: ^{Ovalles (1991)}; ^{Paz et al. (1996)}; ^{Melchiori & Echeverría (2000)}, ^{Ponce et al. (1999}); ^{Sadeghian et al. (2001}); ^{Jaramillo (2002)} and ^{Silva et al. (2003)}, say it may be due to residuality due to fertilization with P.; and that other factors influencing phosphorus retention by the soil, such as: a) Organic matter. In general, it has been found that the effect that the organic phase of the soil P-S diminishes the fixation of phosphorus. b) Presence of hydrated iron and aluminum oxides. In acid soils, these abound and form insoluble compounds with phosphorus. e) Amount and type of clay.

Geostatistical analysis of soil fertility. An exploratory data analysis was performed using the Infostat software, applying a basic statistic to identify the mean, variance, standard deviation, kurtosis and skewness of the data. The SPAD software was applied for multivariate analysis to determine the explanatory traits of the main soil components. In the ArcGIS 10 program, the geostatistical analysis was performed (^{ArcGIS, 2010}). For each trait, the semivariogram was obtained, and through an iterative process in which the active lag and the step were modified. A theoretical model of best fit was established, taking into account as decision parameters the coefficient of determination (R2) and the sum of squares of the residuals (RSS), for which the first must be the closest to 100% and the second the smallest within the situations raised. The parameters of the best-fit theoretical semivariograms obtained for each of the properties studied are recorded in Table 6.

The models of the selected semivariograms correspond to those with the greatest spatial adjustment that explain the distribution of the evaluated variables. Table 6 shows that the semivariograms with the spherical model are the best fits for CIC, Fe, Zn, B, N and C. For MO, P, Ca, K and S, it was the Gaussian, and for Mg and texture, the Exponential model corresponded.

The range indicates the distance from which the samples are spatially independent of each other, and represents the size of grain or stain that the variable represents; in this sense, the traits with the greatest scope of spatial dependence were pH, P, Zn and B. On the other hand, the traits with the minor spatial dependence were Ca, Mg, N, C and texture. It is also observed in Table 6 that in all cases, the effective range exceeded the minimum sampling distance used in this study, which was 50m, suggesting an adequate sampling distance (^{Van et al., 2000}).

The spatial dependence was analyzed through the relative Nugget effect [%C0 / (C0 + C)], cited by ^{Trujillo (2011)}. That revealed a strong dependence for pH, Ca, K, Fe and C (<25%), while MO, P, Mg, Cu, B, N and texture have a moderate dependence (25-75%) and weak spatial dependence (> 75%) for CIC, Zn and S.

Spatial distribution maps. According to ^{Balzarini (2014)}, Kriging interpolations allow the generation of spatial distribution maps for each of the physical and chemical traits studied. As an example, the spatial distribution maps for soil organic matter (SOM) and Mg are indicated below in Figure 1.

Figure 1 Interpolation of the spatial distribution according to the Kriging method for SOM (a) and Mg (b). 

Considering the samplings carried out and the interpolation of the evaluated variables, Figures 2 and Figure 3 show the values reclassified according to the fertility ranges for the Andean zone mentioned in Table 3 of this document. First, Figure 2 indicates the spatial distribution for Boron B, Calcium Ca, Carbon C, Cation Exchange Capacity CIC, Copper Cu, Potassium K and Magnesium Mg, while Figure 3 shows the corresponding values for Manganese Mn, Organic Matter SOM, Nitrogen N, Phosphorus P, Zinc Zn, pH and Texture. In this way, maps with high, medium and low ranges were generated for the said traits. From the data observed in the maps, it is inferred that variables such as Mg, Mn, Ca and K occupy most of the surface of Botana Experimental Farm in a high range, both distributed in the upper, middle and lower part of the farm. The Cu and P presented a certain range, especially towards the middle and lower zones. The variables B, C, N and MO, did not present high ranges, only medium and low.

According to the texture data obtained (Figure 3), the vast majority of the soils of Botana Experimental Farm (12ha) have a sandy clay loam to clay loam texture, and only 5ha are of a sandy loam to loamy texture. The thematic map for pH (Figure 3) shows areas ranging from red (neutral) to green (very strongly acidic). In general, most of the farm has moderately acidic soils (16has).

Figure 2 Classification of spatial distribution of B (a), Ca (b), C (c), CIC (d), Cu (e), K (f) and Mg. The colors red, blue and yellow represent high, medium and low contents, respectively. 

Figure 3 Classification and spatial distribution of Mn (a), SOM (b), N (c), P (d), and Zn (e). Where the red, blue, and yellow colors represent high, medium and low content of the evaluated variables, respectively; For pH (f), where red for pH represents values ​​greater than 6.5 and green less than 5.0; For Texture (g), red represents sandy loam to loam soils, dark green loam to clay loam, orange clay loam, and light green sandy clay loam. 

Classification analysis. To define relatively homogeneous zones based on the principal components, the first three components were selected, explaining 80% of the variability of the data. This analysis allowed obtaining three clusters (Figure 4).

Figure 4 Location of the samples in each Cluster obtained and IDW interpolation of the Clusters. 

Cluster 1. Making an interpolation of the sampled points that correspond to this group, an area of 20.7has was obtained, located towards the western side of Botana Experimental Farm, which corresponds to the lots: 1, 2, 3, 4, 5, 6, 7, 11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 22, 32, 24, 35, 26 and 27. The traits of these batches are: high content of Mg, K, Zn, Mn and Ca; low C and SOM content; medium to low P and Cu content; medium CEC; moderately acidic pH, and sandy clay loam to loam texture.

Cluster 2. This group corresponds to an area of ​​13.9has, located towards the eastern part of the farm, which corresponds to lots: 28, 29, 30, and 31, part of 32, 7 and 11. The traits of this cluster are strongly acidic pH; medium SOM content; low content of P, B, Cu; high content of K, C, Ca, Mg, Mn; high CIC; texture Sandy Clay Loam to Sandy Loam.

Cluster 3. This group corresponds to an area of ​​three hectares, located towards the southeastern part of Botana Experimental Farm. The traits of this zone are neutral pH; medium content of B; high content of Zn, K, P, Ca, Cu, Mg, Mn; high CEC; low C content; texture Sandy Clay Loam to Clay Loam.

According to the previous characteristics, we can assure them that the third cluster presents the best conditions; therefore, it corresponds to an area of high fertility. It should be noted that the administrative infrastructure of the farm is located in this area; however, there are many agricultural vocations with good fertility traits. The second Cluster presents quite favorable conditions for the development of cold-climate crops; therefore, it corresponds to a zone of medium fertility. The first Cluster presents quite limiting traits since it presents the lowest content in some elements, such as SOM, Cu, C and P, with a loamy-to-loamy sandy clay texture (Figure 5).

Figure 5 Zoning of soil fertility in Botana Experimental Farm (the numbers indicate the numbering of the lots managed by the Farm administration). 

CONCLUSIONS

In general, the properties evaluated, except for pH and bulk density, showed significant variability and spatial distribution throughout the study area. Mg, Mn, Ca and K occupy most of the surface, while organic matter and phosphorus content are the most limiting properties. In this case, the distribution maps made it possible to identify spatial trends that, independently of signifying a specific fertility level, could be projected as references for differential nutrition actions in farm management.

Based on the restriction levels for the evaluated traits and possibly due to the soil mineralogy dominated by non-crystalline aluminosilicates, the low cation exchange capacity and soil management factors, three fertility zones could be identified on the farm: high, medium and low, the latter being the most predominant.