Delineation of management zones based on soil physical attributes in an irrigated guava field in the Semi-Arid region , Brazil

Delineating management zones based on soil physical attributes may be an important technique in order to optimize the use of irrigated water and fertilizers in the irrigated agriculture poles in the semiarid region of Brazil. Therefore, the objectives of this study were: (1) To determine the spatial dependence of some soil physical attributes; and (2) To identify management zones using fuzzy cmeans clustering analysis in a guava field in the semi-arid region. The experimental site was structured in a grid of 99 referenced points spaced 4 m, which correspond to the number of guava trees in an irrigated field in the semi-arid region of Brazil. Undisturbed and disturbed soil cores were taken from the 0 to 0.2 m layer in each point. Soil total porosity (TP), soil bulk density (BD), and soil texture were determined. Additionally, soil penetration resistance (PR) was measured using a dynamic hammer penetrometer. Statistical analyses consisted of data description and geostatistics. Management Zone Analyst software was used to carry out the fuzzy c-means clustering algorithm. The spatial dependence of soil texture and BD was confirmed. It was possible to delineate management zones based on soil physical attributes for site-specific irrigation and nutrient management and clay content and BD were the best attributes to delineate management zones in the study area.


INTRODUCTION
The fruit cultivation in the semi-arid region of northeastern Brazil has grown and been an attractive commercial activity.This is due, mainly, to the public investment in irrigated agriculture in this semi-arid region, which created what they called "irrigated agriculture poles".These poles favor the cultivation of various fruit species successfully and allow the production of fruit throughout the year, including during periods when international markets (Europe, Asia and USA) have reduced supplies.One the most important irrigated poles is in Petrolina, Pernambuco, located in the Sao Francisco Valley region of northeastern Brazil.Among the many fruit species, the guava cultivation (Psidium guajava L.) distinguishes itself, since it is the third most important fruit crop cultivated in this region, where about 41% of Brazilian production of guava is produced (Instituto Brasileiro de Geografia e Estatística -IBGE, 2013).Its production has been targeted for the domestic and international markets, as fresh or processed fruits (sweets and juices) (Instituto Brasileiro de Frutas -IBRAF, 2015).
Despite the productive potential of these irrigated poles of the semi-arid region, its singular characteristics such as uneven distribution of rainfall, high rates of evaporation, and frequently shallow and sandy soils increases the risk of soil compaction and salt accumulation (Montenegro and Montenegro, 2006;Santos et al., 2012).Therefore special care in soil management in the semi-arid irrigated areas should be required, mainly related to soil physical-hydric attributes.
Soil attributes vary across the landscape due to differences in soil formation and related processes, as well as in land management and agricultural practices.It means that a crop field is not homogeneous, thus soil management should be carried out taking into account soil variability (Akbas, 2014).Otherwise, important resources can be wasted such as water for irrigation and fertilizers.
Cropping intensity has increased significantly in the last half of century (Valipour, 2015a) and there are still great potential increasing of crop yield (that is, fruit cultivation) in South America (Valipour et al., 2014), including the semi-arid region of Brazil.However, in order to provide food for sustainable development, since water resources are limited, agricultural water management must be taken into account (Valipour, 2014;Valipour, 2015b).Therefore, studying agricultural water management is very important for this region.One approach for agricultural water management is site-specific management based on soil physical attributes.
Site-specific management of soil physical attributes is an attractive and intuitive approach for increasing input efficiency of agricultural systems by adjusting fertilizer and water rates based on the soil characteristics (Peralta et al., 2015).Thus, the concept of management zones can be introduced, which are sub-regions of a crop field that may differ in factors such as the soil type, topography, water or nutrient availability (Bullock et al., 2009).There are a lot of techniques used for the delineation of management zones, such as stochastic interpolators that use geostatistical analysis (that is, ordinary Kriging), deterministic interpolators (that is, inverse distance weight), Co-Kriging, when there are two spatial correlated variables, etc (Dalchiavon et al., 2011).
However, it is not easy to define how many unique zones a field should be divided into, because if the user defines a different number of classes, the results of the classification will be different, and consequently, different management zones will be produced (Fraisse et al., 2001).Fuzzy c-means cluster analysis is an alternative in order to overcome this problem, because this procedure groups similar individuals into distinct classes through an iterative process called clustering (Fridgen et al., 2004), and the number of management zones can be chosen based on the fuzziness performance index (Odeh et al., 1992) and/or normalized classification entropy (Bezdek, 1981).Many studies have used fuzzy c-means clusters analysis to delineate management zones in cereal fields (Jiang et al., 2011;Li et al., 2007;Rodrigues and Corá, 2015), but there are few studies using this procedure in crop fruit fields in the semi-arid region.
Therefore, delineating management zones based on soil physical attributes can be an important technique in order to optimize the use of irrigated water and fertilizers in the irrigated agriculture poles in the semi-arid region.The objectives of this study were: 1) to determine the spatial dependence of some soil physical attributes and; 2) to identify management zones using fuzzy c-means clustering analysis in a guava field in an irrigated pole in the semi-arid region.

Site description
This experiment was carried out in the Nilo Coelho irrigated perimeter, Petrolina county, Pernambuco state, northeastern Brazil, (9°19'10.47"S, 40°33'48.91"W,elev.400 m a.s.l.).The climate of the region, according to the Köppen classification, is of the 'hot semi-arid' (Bsh`) type, characterized by high temperatures (average 26°C), low humidity, high evaporation rates, and especially marked by the scarcity and irregularity in rainfall distribution (400 mm).The soil of the experimental area was classified as Yellow Argisol (Ultisol -American classification Soil Taxonomy).
The guava plants cv.'Paluma', were spaced 4 m × 4 m in planting holes that measured 80 x 80 x 80 cm.The starter fertilization was carried out based on soil analysis and consisted of 0.036 m 3 of bovine manure plus 400 g of single superphosphate (18% P2O5) 30 days before transplanting, and fertilization at each 15 days after transplanting through irrigation water, following the recommendations of Natale and Prado (2006).Soil samples were collected about 60 days after planting.Plants were daily irrigated with one micro sprinkler per plant for a flow of 45 L h -1 , following the crop coefficient (Kc) defined by Bassoi et al. (2001).

Soil sampling and laboratory analyses
The experimental site was structured in a grid of 99 referenced points, which correspond to the number of guava trees.The experimental scheme is depicted in Figure 1.Undisturbed soil cores (0.05 x 0.05 m) were taken from the 0 to 0.2 m layer in each point using a double-cylinder, hammerdriven core sampler.Soil total porosity (TP) and soil bulk density (BD) were determined from these samples using the methods proposed by Flint and Flint (2002).Additionally, soil samples of each point were taken from the layer 0 to 0.2 m using a dutch auger.For each sample, particle size was determined by the standard pipette method (Gee and Dani, 2002).
A dynamic hammer penetrometer with a 65 cm long guide rod, 1.28 cm basal cone diameter and 30° cone angle was used to measure soil penetration resistance (RP) at each grid point to a depth of 20 cm.The cone geometry conforms to American Society of Agricultural and Biological Engineers Standard S313.3 (ASABE, 2010).The mass (M) of the employed drop hammer was 3.9 kg and the mass (m) of the penetrometer system (that is, guide rod and cone) was 3.6 kg.The PR in response to a hammer drop is given as the following equation (Equation 1) (Stolf, 1991): . ( where h is the impact drop height in meters, g is the acceleration due to gravity (9.807 m s − 2 ), A is the cone basal area (1.287 × 10 −4 m 2 ) and Δx is the cone displacement depth in meters.Considering an impact height of 0.4 m and geometrical and physical penetrometer parameters discussed above Equation (1) reduces to Equation (2) (Stolf, 1991): with PR in MPa and Δx in meters.Soil moisture was determined in order to verify spatial correlation with PR.However, it was not verified spatial correlation between PR and soil moisture, therefore, soil moisture was not used in the analysis.

Statistical analyses and data processing
Descriptive statistical analyses (mean, standard deviation, minimum, maximum, coefficient of skewness, and coefficient of kurtosis) were calculated and the Shapiro-Wilk normality test was performed.Estimated semivariogram models were produced to estimate the spatial dependence of samples and to identify systematic and random variations (Isaaks and Srivastava, 1989).Models were selected based on the smallest residual sums of squares (RSS) and the highest determination coefficient (R 2 ).Kriging (stochastic interpolation method) was used to interpolate the variable that showed spatial dependence and Inverse Distance Weighed (IDW) (deterministic interpolation method) was used to interpolate data without spatial dependence (pure nugget effect).The interpolated maps were produced in raster format of 0.5 m pixels, afterwards data files were imported as a comma delimited text file with the first row as labels for the columns to Management Zone Analyst (MZA) software (Fridgen et al., 2004).A spatial correlation matrix using the interpolated map data was performed in order to select the variables that were correlated for clustering (Bazzi et al., 2013).
Fuzzy clustering algorithms were used to classify data into homogenous zones by MZA software.The options settings chosen were measure of similarity as Euclidean when only one variable was clustered (that is, clay content) and Mahalanobis when more than one variable was clustered (that is, clay + bulk density), fuzziness exponent = 1.3, maximum number of iterations = 300, convergence criterion = 0.0001, minimum number of zones = 2, and maximum number of zones = 6 (Fridgen et al., 2004).
A number of validity statistics were used to determine the best combinations of fuzziness index and number of clusters, as well as the overall clustering performance.The fuzziness performance index (FPI) is a measure of the degree to each different classes share membership (fuzziness) and values are constrained between 0 and 1 (Odeh et al., 1992).The normalized classification entropy (NCE) was used to decide how many clusters were the most appropriate for creating management zones, which can help the users obtain management zones simply and quickly (Bezdek, 1981).Minimizing the value for NCE and FPI, the optimum number of clusters can be found.

RESULTS AND DISCUSSION
By comparing average values of soil particle size, the soil was classified as sandy according to the United States Department of Agriculture (USDA) texture method; however, minimum and maximum values indicates that there is high variability in the study field, mainly for clay and silt content (Table 1).Based on the coefficient of variation classification (CV) suggested by Pimentel-Gomez and Garcia (2002), the variability of clay and silt content was high (CV = 20-30%), while the variability in sand content was low (CV < 10%) (Table 1).These results demonstrate that management of water and fertilizer by average values may be not efficient, even in small fields, and wrong decisions could be made about soil management.Based on the coefficients of skewness and kurtosis and the Shapiro-Wilk normality test, soil particle size fractions were considered normally distributed (Table 1).
Soil bulk density (BD) showed values varying from 1.28 to 1.79 g cm -3 (Table 1), which suggests that there are regions in the study field that BD values are higher than the critical interval (1.70 to 1.80 g cm -3 ) for sandy soils as proposed by Reichert et al. (2003).Soil total porosity (TP) values also showed a high range (0.13 to 0.69 cm 3 cm -3 ),  indicating that the guava plant suffered water and/or air stress in some part of the study field due to some low values (< 0.20 cm 3 cm -3 ) (Table 1) (Reichert et al., 2003).Based on CV values, BD and TP were classified as low (CV < 10%) and medium (CV = 10-20%) variability, respectively.Soil penetration resistance (PR) values varied from 0.99 to 11.73 Mpa (Table 1).PR showed the highest range among all study variables and it can be classified as a very high variability (CV > 30%).It can be expected since PR is related with many variables such as soil moisture, soil texture and soil bulk density (Reichert et al., 2007).According to the coefficients of skewness and kurtosis and normality tests, the distribution of BD, TP and PR are non-normal; however, this did not prevent geostatistical analyses, since normality is not a requirement.
Spatial dependence was observed in granulometric attributes and the exponential model was adjusted for the three fractions (Table 2).In an Oxisol, Bottega et al. (2014) adjusted a Gaussian model for clay, sand and silt content at the 0 to 0.2 m layer.The range of granulometric attributes, which is the parameter that indicates the distance at which the samples are spatially dependent in the semivariograms, varied from 10 to 15 m, which was greater than the grid sample spacing (Table 2).It means that the grid spacing was enough to allow modeling spatial dependence in this study area for granulometric attributes (Rodrigues et al., 2012).
Soil bulk density showed spatial dependence and an exponential model was adjusted (Table 2).The range value (13.2 m) was similar to those of granulometric attributes (Table 2).The spherical model was adjusted for BD by Tesfahunegn et al. (2011) in the Mai-Negus catchment of Tigray, northern Ethiopia.The pure nugget effect was observed for TP and PR, which means that the nugget effect (error) values were the same as sill values, that is, there is non-spatial-dependence for these variables in the study area (Table 2).Similar results were found by Rodrigues et al. (2012), who observed nonspatial dependence for TP in a latosol (Oxisol) under a no-tillage system in southwestern Brazil, while different results were found by Campos et al. (2014), who observed spatial dependence for PR at the 0 to 0.3 m layer in an Argisol (Ultisol) in the Amazon region.It appears that the grid spacing was not sufficient to modeling spatial variability for TP and PR in this study area and smaller spacing should be required.These results were expected since TP and PR are not only related to soil formation as soil texture, but also related to some soil managements, thus, their variability tends to be higher than those found for soil texture variables.In addition, PR is also related to soil moisture.
Ordinary Kriging was used to interpolate the variables that showed spatial dependence, which were clay, sand, silt and BD, while Inverse Distance Weighed (IDW) was used to interpolate the variables without spatial The clay content map showed a region with values lower than 10% in the left central region of the map (Figure 2A) and a small region in the lower right corner with values varying from 12 to 15%.This difference in clay content could be enough to interfere with other soil attributes.For instance, the increase of clay content can influence the water retention curves, increasing capillarity and water adsorption (Carducci et al., 2011).In addition, clay content can influence the cation exchange capacity and nutrient availability, which can modify the fertilization management (Power and Prasad, 2010).
The soil bulk density map showed a region with values lower than 1.54 g cm -3 , which are adequate for sandy soils (Reichert et al., 2003) (Figure 2D).The higher BD values were found at the edges of the map.BD are related to soil compaction and soil water storage, thus different values can change irrigation management (Reichert et al., 2003) in the study area.It was possible to interpolate TP and PR variables using IDW; however, due to the non-spatial-dependence their maps tend to create unrealistic island like concentric shapes (bull's eyes).This can be a problem for delineating management zones, possibly producing maps with low precision and accuracy.Although soil maps have been created by interpolation techniques, an important question remains: how many unique zones should a field be divided into?To answer this question fuzzy cluster analyses were performed.
For farmers to adopt site-specific management, the development of management zones must be simple, functional, and economically feasible (Li et al., 2007); thus, the variables more useful to the irrigation and fertilizer management were selected.First of all, clay content was selected because it is the most important fraction of the soil texture, since it is related to majority of soil chemical and physical attributes (Reichert et al., 2007).Second of all, a correlation matrix with interpolated data of the soil maps (pixel by pixel) was performed and the relationship between clay content and BD (r = 0.354) (Table 3) was verified as the best correlation, as suggested by Bazzi et al. (2013).BD and TP showed the highest correlation (r= -0.614); however, management zones with these variables were not performed, because they are directly related, thus, there is collinearity (Table 3).Last of all, the variables that showed spatial dependence were selected.This is the reason that TP and PR were not selected.Therefore, management zone maps were performed for clay content, BD and Clay+BD (Figure 3A to C).In addition, these variables are commonly used by farmers and they can be determined in most soil labs.
The minimum least membership sharing (FPI) and amount of organization (NCE) indicate that the optimum number of clustering zones, considering both indices, was three for clay content, four for BD and also four for Clay + Bulk Density (Figure 4A to C).Similar results were found by Jiang et al. (2011) in a dry field of 1 ha in the western semi-arid region of Chinathe study area was divided into two irrigation management zones by MZA,  and soil physical attributes had high uniformity in each subzone and significant difference between subzones which confirmed the partition; Li et al. (2007) argued that cluster analysis, based on the assumption that grouping data points into naturally occurring clusters reduces within-zone variability, provides an opportunity for identifying management zones in a site and, potentially, applying site-specific management to maximize crop production across the entire area.Fuzzy c-means clustering analysis divides a field in zones, but it does not show the variables values.However, the main objective this analysis is not to apply variable-rate of fertilizer or water, but reduce the variability within-field.For example, it is possible to reduce the number of samples, since it can be performed within each zone.Some practical implications can be suggested for this study area: (1) Nutrient management could be realized using the clay zone map, which could optimize the use of fertilizers as stated by Bazzi et al. (2013).For example, in this present study, the field should be divided in three zones based on clay content, where Zone 1 has the lowest clay content, Zone 2 has medium clay content and Zone 3 has the highest clay content in the field (Figures 2a and 3a).Therefore, fertilizer management could be different for each zone.(2) Irrigation management also could be optimized using the clay + BD zone map.For instance, in order to determine irrigation in a crop field it is necessary to determine the water retention curve (WRC), which is an expensive and time-consuming analysis.In addition, the WRC should be correctly determined and it should represent the water behavior of the area where the irrigation system will be implemented.Since clay and bulk density are the two most important variables related to water retention curve, the samples to perform determine the WRC could be collected within each delineated zone.Another way to increase irrigation efficiency, which always has been one of the main concerns of experts and farmers all of the world as stated by Valipour (2013), it is through simulation models.For this purpose, it is necessary input some data such as depth of flow, volume of infiltrated, field slope, soil texture, etc. (Valipour, 2012).Therefore, management zones based on soil attributes may be an efficient technique to improve the precision and accuracy of irrigation simulation models, since those models takes into account the spatial variability of soil physical attributes, such as soil texture.
Complex field assessments and data manipulation may not be justifiable in terms of time, benefit or economics (Li et al., 2007).Therefore, the results of the present study showed that the MZA software used in the study to delineate management zones reduced study requirements  of basic theories for operators and overcame the disadvantages of a high theoretical level and difficult to master of quantitative analysis.In addition, strong use of science, high precision and fast speed are the characteristics of this method which provides a simple and quick way for the delineation of site-specified irrigation and nutrient management zones as analysis, as argued by Jiang et al. (2011).Minimizing the value for NCE and FPI, the optimum number of clusters can be found.

Conclusions
The spatial dependence of soil texture and soil bulk density was confirmed.However, soil total porosity and soil penetration resistance did not show spatial dependence in the present study area.It was possible to delineate management zones based on soil physical attributes for site-specific irrigation and nutrient management in an irrigated guava field in the semi-arid region in Brazil.Clay content and soil bulk density were the best soil physical attributes to delineate management zones in the study area.

Figure 1 .
Figure 1.Sampling scheme for soil attributes in an irrigated guava field in the semi-arid region, Brazil.

Figure 2 .
Figure 2. Interpolated maps of clay (A), sand (B) and silt (C) content, soil bulk density (D), soil total porosity (E) and soil penetration resistance (F) from 0 to 0.2 m depth in an irrigated guava field in the semi-arid region, Brazil.

Figure 3 .
Figure 3. Management zone maps using fuzzy c-means clustering analyses of clay content (A), soil bulk density (B) and Clay+Bulk density (C) in an irrigated guava field in the semi-arid region, Brazil.
(A) Clay content zones (B) Bulk density zones (C) Clay + bulk density zones

Figure 4 .
Figure 4. Fuzziness performance index (FPI) and normalized classification entropy (NCE) as calculated by Management Zone Analyst software for clay content (a), soil bulk density (b) and Clay+Bulk density zones in an irrigated guava field in the semi-arid region, Brazil.

Table 1 .
Descriptive statistics of soil physical attributes from 0 to 0.2 m depth in an irrigated guava field in the semi-arid region, Brazil.

Table 2 .
Variogram model parameters and interpolation method of soil physical attributes from 0 to 0.2 m depth in an irrigated guava f ield in the semi-arid region, Brazil.

Table 3 .
Spatial correlation matrix of maps (pixel by pixel) of soil physical attributes from 0 to 0.2 m depth in an irrigated guava field in the semi-arid region, Brazil.