Delineation of irrigation management zones in a Quartzipsamment of the Brazilian semiarid region

The objective of this work was to delineate irrigation management zones using geostatistics and multivariate analysis in different combinations of physical and hydraulic soil properties, as well as to determine the optimal number of management zones in order to avoid overlaping. A field experiment was carried out in a Quartzipsamment, for two years, in an irrigated orchard of table grape, in the Senador Nilo Coelho Irrigation Scheme, in the municipality of Petrolina, in the state of Pernanbuco, Brazil. Soil samples were collected for the determination of soil physico-hydraulic properties. A portable meter was used to measure soil apparent electrical conductivity. Spatial distribution maps were generated using ordinary kriging. Management zones for five different combinations of soil properties were defined using the fuzzy c-means clustering algorithm, and two indexes were applied to determine the optimal number of management zones. Two combinations of soil properties can be used in the management zone planning in order to monitor soil moisture.


Introduction
The knowledge of the spatial variability of soil physical and hydraulic properties which influence the soil water dynamics can improve irrigation management.Geostatistics can contribute to the understanding of this variability (Grego et al., 2014), but each soil property should be assessed separately.Kitchen et al. (2005) and Valente et al. (2012) suggest that management zones should be defined by evaluating more than one soil property.Therefore, multivariate analysis can help to delineate management zones, besides facilitating the study of various soil properties.Li et al. (2007) reported that multivariate analysis clustering effectively identifies management zones because it is simple, functional, and economically viable.Morari et al. (2009) claim that the combination of geostatistical interpolation and multivariate clustering analysis aids the precision viticulture by enabling the efficient division of a cultivation area into management zones.
Several studies have used a combination of clustering analysis with fuzzy c-means algorithm (FCM) (Bezdek et al., 1984) and geostatistical analysis to delineate management zones (Li et al., 2007;Morari et al., 2009;Moral et al., 2010;Wang et al., 2012).Odeh et al. (1992) found that this algorithm generates continuous clustering and it is, therefore, helpful in defining management zones with continuous soil properties.Other clustering analyses used for this purpose include Ward minimum variance (Fleming et al., 2000), fuzzy k-means (Song et al., 2009), and colony (Hou et al., 2016).
Soil apparent electrical conductivity (ECa) has also been used as an auxiliary parameter in delineating management zones (Kitchen et al., 2005;Li et al., 2007;Valente et al., 2012) because of its strong correlation with different soil properties (Moral et al., 2010;Molin & Rabello, 2011;Kweon et al., 2013).Nevertheless, there are very few studies using ECa to delineate management zones for sandy and sandy loam soils of largely irrigated agricultural areas, such as those in the Submédio São Francisco River Basin, located in the semiarid region of Brazil.
ECa indicates seasonal variations based on soil water content changes (Brevik et al., 2006;McCutcheon et al., 2006).These variations can affect the interaction between ECa and other stable soil properties.Cox (2005) reported that cluster analysis identifies patterns and characteristics via data partition.Jaynes et al. (2005) used cluster analysis to find temporal and spatial patterns related to soybean yield using data collected for several years.The application of cluster analysis in ECa data, when temporally and spatially represented, can assist in the identification of temporal patterns in the data.
The objective of this work was to delineate irrigation management zones by assessing different combinations of physical and hydraulic soil properties, using geostatistics and multivariate analysis.The output would be used to determine the optimal number of management areas, in order to avoid overlapping.

Materials and Methods
This study was carried out in a vineyard of 'Thompson Seedless' grapevines grafted on 'SO4' rootstock, planted in May 2004, and trellised in a multiwire, horizontal system.The 1.6 ha experimental area had 20 rows and 82 plants per row, in 2.5 m spacing between vines, and 4.0 m between rows.The area is located in the farm 180 of Senador Nilo Coelho Irrigation Scheme, Sector 5 (9°23'12.8"S,40°38'13.8"W,at 394 m altitude), in the Submédio São Francisco River Basin, in Petrolina county, in the state of Pernambuco, Brazil.The irrigation system was a microsprinkler with one emitter per plant.Its estimated flow was 30 L h -1 , and its wet area was 2.4x2.5 m.The soil in the area was classified as a Neossolo Quartzarênico (Santos et al., 2013), which corresponds to a Quartzipsamment in the soil taxonomy classification.
The variability of physico-hydraulic soil properties was evaluated by defining two transects in the study area, located in rows 5 and 15.Each transect had forty points separated every 5 m.There were eighty sampling points (Figure 1).Deformed soil samples were collected using an Edelman auger, and undeformed soil samples were collected using a cylindrical ring (0.05x0.05 m).Samples were taken in 2011 at 0.00-0.20 and 0.20-0.40m soil layers.The grower numbered the rows and plants, which facilitated georeferencing of each collection site.
In the deformed samples, sand, silt, and clay granulometric fractions were determined (kg kg -1 ), accordingly to Donagema et al. (2011).The centrifugation method (Silva & Azevedo, 2002) was used to determine soil moisture (m 3 m -3 ) at field capacity (θ fc , moisture retained at 0.006 MPa), and the permanent wilting point (θ pwp , moisture retained at 1.5 MPa).Available water (AW, m 3 m -3 ) was obtained by measuring the difference in the water content between θ fc and θ pwp .Soil density (Ds, kg dm -3 ) was determined in undeformed samples using the volumetric ring method.Representative values of the soil layer at 0.00-0.40m soil layer were obtained by calculating the arithmetic averages in each sampled layer.
ECa (dS m -1 ) was measured in 2012 at 0.00-0.40m soil layer, using a portable meter (Rabello et al., 2010) at 100 and 101 days after pruning (DAP) of the grapevines.Sampling all twenty rows, separated at every second plant, resulted in a 4.0×5.0m grid and 820 sampling points (Figure 1).Other measurements were performed in 2012 at 57 and 60 DAP, and in 2013 at 63, 78, and 91 DAP, using rows 1, 5, 10, 15, and 20, separated at every second vine, for a total of 205 sampling points.
The spatial distribution of variables was characterized by building experimental semivariograms, and adjusting theoretical models and their respective setting parameters (nugget effect, sill, and range of spatial dependence).The quality of the adjustments was assessed by calculating the spatial dependence index (SDRz), proposed by Zimback (2001), and SDR, by Seidel & Oliveira (2014).Cross-validation was performed for both the regression coefficients (RC), which indicate whether the kriging underestimates or overestimates extreme values, and the determination coefficients (R²) were evaluated.The adjusted models were then used in an ordinary kriging procedure to interpolate the sampled values in unsampled locations.Isoline maps were built for each property, employing the interpolated values in a single grid of 2.58 × 1.31 m for a total of 4,992 interpolation points.The geostatistical and map-building procedures were performed using GS + (Robertson, 1998).
For the cluster analysis, the interpolated values were standardized to obtain the same amplitude for all sets of properties, with a mean value equal to 0 and standard deviation equal to one.Guo et al. (2012) reported that this procedure is a prerequisite for multivariate cluster analysis.
Multivariate cluster analysis was performed using the FCM algorithm (Bezdek et al., 1984) and the e1071 package (Meyer et al., 2014) from R software (R Core Team, 2012).The area was divided into different management zones, according to the following interpolated attributes: C1, physical and hydraulic soil properties; C2, ECa for 205 points (57 and 60 DAP), and 820 points (100-101 DAP The FCM clustering analysis assigns degrees of association with each management zone (grouping) to the points of interpolated values (individuals), allowing the evaluation of the relationship between each individual and each group (Cox, 2005).This analysis is an iterative process that results in the optimal division of data into a number of groupings.
The initially assigned variables were: the number of interactions (K); the fuzzy weight exponent (m), ranging from 1 to ∞; and the number of management zones.K should meet the criterion of completing the iterative process (Ɛ).For this study, the criteria used by Odeh et al. (1992) were adopted with ɛ = 0.001.Therefore, K = 50 sufficed to meet the adopted Ɛ value.Cox (2005) indicated that m is usually used in the range of 1.25 to 2.0.Therefore, m was defined as 1.25 for all combinations of variables.
The optimal number of management zones was chosen by testing the cluster analysis for 2, 3, 4, Pesq.agropec.bras., Brasília, v.51, n.9, p.1283-1294, set.2016 DOI: 10.1590/S0100-204X2016000900028 5, and 6 zones.The optimal number of zones was validated using the fuzziness performance index (FPI) and the modified partition entropy (MPE).These parameters were used in several studies (Odeh et al., 1992;Fridgen et al., 2004;Sun et al., 2012).According to Fridgen et al. (2004), FPI measures the separation degree of individuals (classified by cluster analysis) from the generated clusters.Its value varies between 0 and 1.The closer the FPI is to zero, the smaller is the overlap between the management zones, and, consequently, the better defined is the partition (Odeh et al., 1992).The MPE measures the degree of disorganization of individuals between groups.Its value varies between 0 and 1, and values closer to zero indicate better organized groups (Fridgen et al., 2004).Therefore, the optimal number of management zones should give the lowest values for these functions.
After the optimal number of zones has been identified, for each combination of properties, each interpolated point will belong only to the area where it showed the highest degree of association.This process is known as defuzzification (Guastaferro et al., 2010).After the partitions and spatial locations were obtained for each interpolated point, management zone maps were made using GS + (Robertson, 1998).
Soil sampling points and ECa measurement sites were superimposed onto the management zone maps of each property combination to associate points with management zones.Therefore, the values of each physico-hydraulic soil property and each ECa measurement, representative of 0.00-0.40m soil layer, were separated into different data sets corresponding to the points of each management zone.Their averages were calculated.Significant differences in the average values of each data set from each management area were evaluated by analysis of variance using Tukey's HSD (honestly significant difference) test (Figure 2).

Results and Discussion
Data for sand, silt, clay, and AW (available water) contents were normally distributed whereas the Ds data distribution was not, according to the results of the Kolmogorov-Smirnov test.For ECa, only measurements at 57 DAP (2012) and 63 DAP (2013) showed normal distribution.
Although data normality is desirable, it is not absolutely necessary for kriging (Moral et al., 2010), as long as the distribution does not include many elongated tails (Cressie, 1993).In the present study, data were not normally distributed (Table 1); however, since the kurtosis coefficients were relatively close to zero, the tails were only slightly elongated and allowed the application of geostatistics.
The coefficient of variation (CV) for sand and Ds data was lower than 12%, indicating low variability.For the remaining properties, however, the variation was intermediate (12%< CV <60%), as proposed by Warrick & Nielsen (1980).Corwin et al. (2006) found similar CV values for clay, Ds, AW, and ECa at different depths in clay loam soils.Molin & Faulin (2013) classified CV in a similar way for clay and ECa data in two areas with medium soil texture.
The model, when adjusted to the experimental semivariogram, ranged between Gaussian and spherical for the physical and hydraulic soil properties (Table 2).  (2SDRz, spatial dependency ratio (Zimback, 2001).ECa, soil apparent electrical conductivity; DAP, days after pruning; and CR, coefficient of regression.
The adjusted models were highly accurate for most variables, except for silt, which showed lower R² values than other parameters (Table 2).The range values indicated that the sampling density was adequate to achieve spatial dependence for all studied variables.The range value of attribute AW was higher than those of the other variables, which indicates greater homogeneity of AW throughout the study area.For future samplings, the distance among samples for each physico-hydraulic soil property should be less than its range values to ensure spatial dependence among samples.
The SDR proposed by Seidel & Oliveira (2014) indicated that AW, Ds, and ECa data at 60 and 63 DAP (in 2012 and 2013, respectively) showed higher spatial dependence (Table 2).The SDR proposed by Zimback (2001) was strong (SDRz >75%) for the clay, AW, Ds, and ECa, after 91 DAP (2013); the remaining attributes and ECa measures had moderate spatial dependence (25%< SDRz ≤75%).Thus, spatial dependence was observed for all ECa data.AW and Ds had the highest spatial dependence, indicating higher quality of the spatial distribution maps for these variables.
CR and R 2 values of the cross-validation were adequate, except for the sand and silt fractions, and ECa at 57 DAP ( 2012) and 91 DAP ( 2013), which showed a lower quality of the model adjusted for the semivariogram when compared with the other data sets.Parfitt et al. (2009) obtained 0.88, 0.90, and 0.91 CR values, and 0.27, 0.68, and 0.78 R 2 values for clay, silt, and sand, respectively, in medium-textured soils.
The results of validation for the optimal number of management zones for each combination are presented in Figure 3.The lowest values of each index for an equal number of management zones were obtained for the combinations C1, C3, C4, and C5, with 3, 2, 2, and 5 management zones respectively.For C2, the lowest FPI was two management zones, and the lowest MPE was four zones.Nevertheless, with three management zones, these two indexes were simultaneously low, so this number of zones was adopted as optimal for this combination.Overall, the range of the ideal number of management zones determined in several studies is similar to that found in the present study.Li et al. (2007) and Morari et al. (2009) suggested that three zones were ideal, whereas Guastaferro et al. (2010) proposed two zones as the optimal number.Valente et al. ( 2012) used different combinations of variables, and found different optimal numbers of zones; however, they adopted three zones as the standard for all combinations because this number predominated.It was found that the clay attribute had the greatest heterogeneity in the study area (Figure 4).This observation is consistent with its lower range value (Table 2).
The spatial distribution was similar for the AW and Ds maps (Figure 4).The lowest values were concentrated in the upper right side of the maps.This result indicates that AW can only be estimated with Ds data for this soil type.This limitation may reduce testing costs.The similarity between the spatial distributions of AW and Ds may be due to the higher percentage of macropores in sandy and sandy loam soils than that in clay soils; when Ds increases, this percentage decreases while the water retention capacity of those soils increases.In those soils, the micropores -which retain moisture at > 1.5 MPa tensions owing to the increase in Ds -are in small proportion and have limited influence on AW, when compared with the decrease of the macropore proportion.Grego & Vieira (2005) reported similar results for Oxisols, in which AW was higher in regions with higher Ds.
Comparison of management zone maps with those for the spatial distributions of the physico-hydraulic soil properties (Figure 4) indicates that clustering separated the larger values of AW and Ds in zone 3, the higher percentage of clay and the lower percentage of sand in zone 2, and the higher percentage of clay and the lower percentage of silt in zone 1.
Higher values were observed on the left side of the ECa spatial distribution map and on the right side above the center (Figure 5).During the two-year study period, as the DAP increased, the regions with larger ECa values decreased.Since there was little variation in soil texture, ECa varied due to a decrease in soil moisture (θ).In 2012 and 2013, the grapevines were irrigated throughout their growing seasons, but the amount of irrigation water applied was lower in the final stages of each season, due to the lower water demand by plants.According to Nascimento (2013), in 2012, the average θ (0.0093 m 3 m -3 ) at 0.00-0.45m soil layer, at 100 and 101 DAP, was lower than it was at 57 DAP (0.101 m 3 m -3 ) and at 60 DAP (0.0097 m 3 m -3 ).In 2013, the average θ value at 91 DAP (0.0092 m 3 m -3 ) was lower than that at 63 DAP (0.109 m 3 m -3 ) and at 78 DAP (0.108 m 3 m -3 ).
Comparison of the spatial distribution maps with those for the management zones of combinations C2 and C3 (Figure 5) indicated higher ECa in zone 3 (from combination C2) and zone 2 (from combination C3).Although the distribution of the ECa was improved by using three management zones in 2012, and two in 2013, it is clear that this attribute changed little during the studied period, indicating temporal uniformity.The θ values did not vary significantly over the studied period.
The management zone map for combination C4 (Figure 6) summarizes ECa monitoring during grapevine growing season in 2012 and 2013.The coverage area of zone 2 was similar to those in zones 3 and 2 from combinations C2 and C3, respectively (Figure 5).They showed the highest ECa levels.
The low variation of θ due to irrigation in the area throughout the growing seasons allowed the determination of the spatial and temporal patterns of ECa.It represents various soil properties with instantaneous measurements, and enables greater detailing of the area without increasing costs.
Comparison of the management zone maps from combination C1 (Figure 4) with those from combination C5 (Figure 6) indicated that area 3 was divided into three other zones (zones 3, 4, and 5).This division is due to the addition of ECa to the cluster analysis, which increased detailing in the area.
The variance analysis applied to the average values of the physical and hydraulic soil properties, and ECa in each management zone obtained by FCM cluster analysis, indicated a significant difference between at least two management zones for all evaluated properties, with a confidence level of at least 95% for all combinations.This result shows that cluster FCM analysis, along with geostatistical analysis, was effective in dividing the area into management zones.Several other studies that used geostatistics and FCM cluster analysis also defined management zones satisfactorily (Li et al., 2007;Valente et al., 2012;Wang et al., 2012).
The average values of texture in the management zones formed by combinations C1 and C5 were similar in zones 1 and 2 (Table 3).In the remaining zones, which subdivided zone 3 (from combination C1) into three other zones (3, 4, and 5), higher amounts of silt were observed in zone 4, and higher amounts of clay were seen in zone 5.Both of these zones were formed from combination C5.
The ECa averages in combination C3 were the following (Table 3): high in zone 3, intermediate in zone 2, and low in zone 1.In the management zones formed from combinations C3 and C4, ECa were higher in zone 2 and lower in zone 1. Corwin et al. (2006) reported that several soil properties influence ECa levels.Morari et al. (2009) reported a close association between ECa and the physical properties of the soil.It was positive for the fine soil fractions (silt and clay), and negative for the coarse soil fractions (gravel and sand).In combination C5, the analysis of zone 5 indicated the influences of clay and Ds on the ECa levels.There was a strong association of high ECa values with high clay values and Ds, in comparison to the other zones.Corwin & Lesch (2005) observed that   Stadler et al. (2015) evaluated the correlation between soil texture and ECa, in three different areas with clay and clay loam soils.They found correlation between soil texture and ECa in only one area.These authors concluded that the absence of correlation was due to low sample variability for the other soil properties.For the present study, in sandy and sandy loam soils, ECa was influenced by Ds, when combined with the content of clay or silt fractions.Nevertheless, only the division of the area into zones with homogeneous soil properties, and the separate analysis for each zone, allowed the assessment of the influence of these properties on ECa.
The management zone maps from combination C1 or C5 (Figures 4 and 6, respectively) and the average values for the physico-hydraulic soil properties, in each management zone (Table 3), allowed the planning of the zones for monitoring θ values in them, and the adoption of a differentiated irrigation management.Nevertheless, combination C5, which used ECa data, enabled greater detailing of the fine soil fractions, and could, therefore, help improve water use management.

Conclusions
1.The delineated management zones using the physico-hydraulic soil properties (C1), and the measurement of the soil apparent electrical conductivity combined with those properties (C5) allows the delimitation of the management zones; Table 3. Mean values of soil attributes in the management zones arising from combinations of physical and hydraulic soil properties, at 0.00-0.40m soil depth in Quartzipsamment (1) .
2. The use of FCM clustering analysis associated with geostatistical analysis enables the standardization of the soil apparent electrical conductivity data from different time periods.

Aknowledgements
To Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (Capes), for granted scholarship to the first author; and to the Sasaki Farm, who provided the experimental area.

Figure 1 .
Figure 1.Representation of soil sampling sites, and readings of soil apparent electrical conductivity (ECa), in different days after pruning (DAP) of the grapevines.

Figure 2 .
Figure 2. Flowchart of the procedure stages.

Figure 3 .
Figure 3. Fuzziness performance index (FPI) and modified partition entropy (MPE), for different numbers of combinations in management areas, wherein: C1, combination of physical and hydraulic soil properties; C2, data of soil apparent electrical conductivity (ECa) measured in 2012; C3, ECa measured in 2013; C4, ECa measured in 2012 and 2013; and C5, ECa measured in 2012 and 2013, and physical and hydraulic soil properties, at 0.00-0.40m soil layer in Quartzipsamment.

Figure 4 .
Figure 4. Spatial distribution maps of physical and hydraulic soil properties, and management zones obtained with the combination of the physical and hydraulic soil properties (C1) at the 0.00-0.40m soil layer in Quartzipsamment.

Figure 5 .
Figure 5. Spatial distribution maps of soil apparent electrical conductivity (ECa), in different days after pruning (DAP) of grapevines in 2012 and 2013, and management zone maps obtained with ECa data measured in 2012 (C2) and in 2013 (C3), at 0.00-0.40m soil layer in Quartzipsamment.

Figure 6 .
Figure 6.Management zone maps obtained with soil apparent electrical conductivity (ECa) measured at different days after pruning (DAP) of grapevines, in 2012 and 2013 (C4), and the areas obtained by the combination of ECa data, in 2012 and 2013, and the physical and hydraulic soil properties (C5), at 0.00-0.40m soil depth in Quartzipsamment.

Table 1 .
Descriptive statistics and normality test of soil properties at the 0.00-0.40m soil layer in Quartzipsamment.

Table 2 .
Theoretical semivariogram models with their respective parameter settings, and cross-validation of soil properties at 0.00-0.40m soil layer in Quartzipsamment.

Table 3
) had the lowest average clay fraction and generally intermediate ECa values (owing to the high Ds values of that zone), and intermediate to high silt values in comparison to the other zones.In general, zone 2 had high ECa values, possibly because of the higher amounts of clay and silt in that area.Ekwue & Bartholomew (2011) analyzed ECa in clay, clay loam, and sandy loam soils, and concluded that the clay fraction influenced ECa in proportion to the θ values.ECa values, like those for Ds, were similar to AW values.Although ECa is theoretical, increases in its values do not necessarily indicate larger θ values.
Means followed by equal letters, among management zones, do not differ by Tukey's HSD test, at 5% probability.AW, available water; Ds, soil density; ECa, soil apparent electrical conductivity.DAP, days after prunning.