Integration of PCA and Fuzzy Clustering for Delineation of Soil Management Zones and Cost-Efficiency Analysis in a Citrus Plantation

Citrus spp. are one of the most important commercial crops with global marketing potential in the world, as in Iran. A soil management zone (MZ) as an appropriate approach is necessary to achieve sustainable production, along with improving soil management and increasing economic benefits in the commercial citrus plantations of northern Iran. As the first report, the biological and terrain attributes along with the physicochemical properties (57 soil samples, 0–30 cm) were used for MZ delineation using the integration of principal component analysis (PCA) and the fuzzy c-means clustering methods. An economic analysis based on the MZ results was also performed to determine the changes in each MZ using a relative cost (RC) value. The high correlation between soil properties and terrain attributes and the considerable spatial variation of these factors in the study area call for site-specific nutrient management. The optimal number of MZs was six and there was a significant heterogeneity variation among different MZs. The ranking of the MZs were MZ5 > MZ2 > MZ6 > MZ1 > MZ3 > MZ4 based on higher soil quality and lower costs per tree. The MZ4, MZ3, MZ1, MZ6, and MZ2 required 34.4, 30.6, 29.4, 9.77, and 9.44% more costs than MZ5 (as reference MZ) for achieving similar productivity, respectively. Therefore, this simple and cost-effective approach could be an initial step to utilize fertilizers site-specifically for data-scarce areas and reduce the soil property variability within the delineated MZs, which is fundamental for precision agriculture management.


Introduction
Citrus is one of the most important commercial crops with global marketing potential, which is grown across the world, particularly in tropical and subtropical countries [1,2]. It is important due to its role in human diets, as citrus fruits contain a high level of vitamin C and other essential nutrients, such as folate, potassium, and dietary fiber [3]. Annual global production of fresh citrus fruits has been growing in the last several decades, from 30 thousand tons in the 1960s to 130,947 thousand tons in 2015 [4]. Iran, the world's seventh largest citrus-producing country, followed by Brazil, China, United States, Mexico, India, and Spain, has an annual production of 3739 thousand tons of fresh citrus fruits from an area of~230,000 ha [4]. The growing demands for crop products in the 21st century have increased the request for new methods to evaluate the soil fertility status in an economical and environmentally friendly manner [5,6]. Moreover, the knowledge of spatial soil fertility status is important to intensive agricultural practices for sustainable agriculture [6,7].
Defining the homogeneous management zones (MZs) is vital for sustainable use of resources and economical aspects for gardeners and farmers [8]. In this regard, the use of soil management zones is necessary to achieve sustainable agriculture and a cost-effective approach for improving the management of citrus plantations.
Considerable advances in statistical techniques and modeling approaches have allowed different agricultural zones to be accurately delineated with new methods like machine learning, GIS, remote sensing, geostatistics, and kriging at a regional and farm scale [5,9]. Several methods and approaches, such as delineating by an active canopy sensor [10], using soil properties by fuzzy k-means cluster analysis [11,12], fuzzy c-means clustering algorithm, principal component analysis (PCA) [12][13][14][15][16][17], and soil fertility analysis [13] have been used to delineate different MZs. However, the selection of the techniques or approaches depends on the purposes of the determination of MZ and the complexity of each approach.
Different types of environmental and soil data could be used for delineating the MZs [10]. For instance, pedo-geomorphological factors [18], soil nutrients [19,20], intrinsic soil properties (e.g., electrical conductivity (EC), pH, particle size distribution, soil organic matter, cation exchange capacity, nitrogen, phosphorous, and potassium) [9,10,12,15,21] can be used for this purpose. Different types of information can be integrated to define MZs by using a clustering algorithm. This effective approach is capable of identifying zones or regions in the fields or farms that internally have similar soil properties, conditions, and productivity potentials [5]. In the application of site-specific management zones, several studies reported successful application of MZ delineation for agricultural purposes, fertilizer application, and farm management [12,22,23]. Aggelopooulou et al. [24], in an apple orchard in Greece, delineated four zones using soil and tree properties, which required different management of irrigation, fertilizer, and cultural activities. In coffee fields, by using physiochemical soil properties and the advantage of geostatistical analysis with fuzzy c-means clustering, Valente et al. [25] reported that it is required to use more than one soil property for better MZ delineation. Furthermore, Peralta et al. [26] delineated MZs in wheat fields in Argentina for the increment of nitrogen efficiency.
In most countries, such as Iran, the fertilizer and nutrient recommendations for agricultural lands are usually uniform on a regional scale [15] or farm scale [10,26], with the spatial heterogeneity of nutrient content in soils being often ignored [9]. However, the source of chemical fertilizers in Iran commonly include nitrogen (urea, 46%), phosphorus (triple superphosphate, 46%), potassium (potassium sulfate, 46%), and also cow manure [27]. Therefore, the similar fertilizer recommendations in the farm or regional scales could lead to over-application of fertilizers in areas with high nutrient levels and vice versa [9,15]. Herein, Peralta et al. [26] concluded that MZ delineation could increase the fertilizer efficiency of nitrogen (or other nutrients) in the commercial fields, and the risk of pollution will decrease along with a decrease in the application of resources, resulting in an increase of economic benefits.
Mazandaran Province, the northern seaside lands of the Caspian Sea, is the main citrus producing area in Iran. The yield of citrus in northern Iran is crucially dependent on soil properties and agricultural practices [27]. Therefore, this study, as the first assessment in a citrus plantation in Iran, aimed to evaluate the potential tools and methods to delineate MZs in the commercial citrus plantations of northern Iran for increasing economic benefits and identifying homogenous zones in this region. The hypothesis for this research is that the MZ could decrease the production costs during the growing season by producing different zones to which to apply various amounts of fertilizers.

Study Area
One hundred and sixty hectares of citrus plantation were investigated in northern Iran (Sari region, Mazandaran Province). The study area was located between 52 • 57 54 and 52 • 59 13 E longitudes and 36 • 29 05 and 36 • 29 50 N latitudes, as shown in Figure 1. The mean annual temperature and precipitation were 17.9 • C and 789 mm, respectively. The soil moisture and temperature regime classes are udic and mesic in this region, respectively [28].

Study Area
One hundred and sixty hectares of citrus plantation were investigated in northern Iran (Sari region, Mazandaran Province). The study area was located between 52°57′54′′ and 52°59′13′′ E longitudes and 36°29′05′′ and 36°29′50′′ N latitudes, as shown in Figure 1. The mean annual temperature and precipitation were 17.9 °C and 789 mm, respectively. The soil moisture and temperature regime classes are udic and mesic in this region, respectively [28].

Field Sampling and Data Collection
The soil sampling of MZ delineation tends to decrease the number and cost of the sampling procedure, while it still prepares precise information on soil properties and their spatial variability [29]. Totally, 57 composite soil surface samples (at a depth of 0-30 cm) in the commercial citrus plantation (8)(9)(10)(11)(12) year old trees) were taken at randomly chosen sites across the study area, as shown in Figure 1. During soil sampling, we used a portable Global Positioning System (GPS) (Etrex Vista Garmin, Olathe, KS, USA) to accurately mark the sampling points.

Soil Physicochemical Properties Analysis
The soil samples were air-dried at room temperature (~25 • C), crushed, and passed through a 2 mm sieve for chemical and physical analyses. The wet combustion method [30] was used to determine soil organic carbon (SOC). The soil bulk density (BD) was estimated using the cylindrical core method after 24 h at 105 • C [31]. The soil reaction (pH) [32] and soil electrical conductivity (EC) [33] (both measured on a 1:2.5 suspension) were determined using a multi-parameter apparatus. The Kjeldahl method [34] was used to measure the total nitrogen (TN) of the soil samples. The cation exchange capacity (CEC) was determined using sodium acetate (NaOAc) at a pH of 8.2 [35]. The available soil phosphorus (P av ) [36] and soil available potassium (K av ) [37] were also measured for all soil samples.
Accordingly, 25 g of each air-dried soil sample was placed in a closed glass jar, then moistened (~85% of field capacity) and put in a beaker containing a NaOH (10%) solution therein. The samples were incubated at 30 • C for four days, then the microbial soil respiration (Mres4, mg CO2-C kg −1 soil d −1 ) was measured using the method suggested by Bakhshandeh et al. [38]. The chloroform fumigation-extraction method was used to determine the microbial carbon biomass (Cmic, mg C kg −1 soil) in the soil samples [39].

Descriptive Statistics Analysis
The descriptive statistics, including mean, maximum (Max), minimum (Min), median, standard deviation (STD), coefficient of variation (CV), kurtosis coefficients, and skewness, were calculated for all data using the Statistical Analysis System (version 9.4, SAS Institute, Cary, NC, USA, [40]) software. The normality (d) of all soil properties was tested by the Kolmogorov-Smirnov (K-S) test at a significance level of 5%. Based on Warrick and Nielsen [41], the mean values of all attributes around the mean could be classified as data dispersion to three categories as low (CV < 12%), moderate (12 ≤ CV < 60%), and high (CV ≥ 60%).

Spatial Variability Analysis
Considering the sparse sampling points due to the applicability of MZ delineation as an economically and timely approach in the study area, the dataset cannot be directly used for preparing the spatial distribution maps of variables, because many hectares include no data or only a single one, as shown in Figure 1. Although the kriging method is the most common method for interpolation in soil MZ [20,26,42], the semivariograms of the soil variables show a weak to moderate spatial dependence. Therefore, we used the inverse distance weighting (IDW) method to generate the spatial variability maps.
The leave-one-out cross-validation was implemented for testing the performance of the IDW method for predicting the soil and terrain variables in unsampled locations in the study area [43,44]. Three evaluation criteria, including the coefficient of determination (R 2 ), normalized root mean square error (nRMSE, indicating the relative difference (%) of predicted versus observed values), and coefficient of residual mass (CRM) were used for testing the reliability on the interpolation method. The equations of the three evaluation criteria can be written as where, n is the number of samples, Oi and Pi are the predicted and observed variables, and O and P are the means for the predicted and observed variables, respectively.Ň is the mean of observed values.
According to Jamieson et al. [45], when the nRMSE is <10%, the results are considered an excellent interpolation; when 10% < nRMSE < 20%, it is a good interpolation; when 20% < nRMSE < 30%, a fair interpolation; and when the nRMSE > 30%, it is considered a poor interpolation. For the CRM criterion, the negative and positive values indicate that the models overestimate and underestimate the observed data, respectively.

Delineation of Management Zone
Briefly, the methodology used for the delineation of MZs consists of the principal component analysis (PCA) on soil and terrain variables and then performing the fuzzy c-means clustering fed with the selected principal components (PCs) to delineate MZs. The PCA, as a dimension reduction approach, uses correlated variables to produce new orthogonal variables by the orthogonal transformation called PCs. The PCA aggregates the principal sources of variability in the data. The PCs, having eigenvalues of equal to or greater than one, were selected to develop MZ classes in this study.
A fuzzy c-means clustering, an unsupervised continuous classification procedure, was carried out using the Management Zone Analyst (MZA) (Version 1.0, University of Missouri, Columbia, MO, USA) software to produce homogeneous MZs (see Fridgen et al. [46] and Jiang et al. [20] for details). This algorithm takes the spatial PCs of related variables obtained through the PCA as inputs for fuzzy c-means cluster analysis. This clustering algorithm is adequate for grouping soil and terrain variables in the soil continuum because it produces a continuous grouping of objects by assigning partial class membership. The studied area was divided into two to seven clusters using MZA software with the maximum number of iterations of 300, and the minimum and the maximum number of zones were 2 and 7, respectively. The stopping criterion of 0.0001 and fuzziness exponent of 1.5 were used; a Euclidean measure of similarity for univariate clustering and Mahalanobis measure of similarity for multivariate clustering [16,47].
Two cluster validity functions, including fuzzy performance index (FPI, indicating the degree of fuzziness) and normalized classification entropy (NCE, indicating the amount of disorganization) were used as indicators to determine the best cluster number [42]: where c is the number of clusters, n is the number of observations, µik is the fuzzy membership, and log a is the natural logarithm. The best number of clusters achieved when the FPI and NCE were at the minimum level demonstrated the least membership sharing and the greatest amount of organization. A one-way ANOVA was performed to determine the significant heterogenic variations (p ≤ 0.05) between MZs by the least significant difference (LSD) test. ILWIS 3.3 software was also used for spatial analysis and mapping.

Economic Analysis
The cost-efficiency of MZ delineation was one of the objectives of this study and also a comparison of production costs for different MZs. Thus, an economic analysis based on the MZ results was conducted due to changes in different soil properties and levels of chemical fertilizer application per tree in each MZ. This analysis was performed to obtain the cost-saving of the MZs as a percentage in comparison to the lowest cost zone. Therefore, the necessary economic data and prevalent costs (e.g., price of manure and chemical fertilizers) declared and approved by the Soil and Water Research Institute of Iran [48] were used. Since the local currency and price levels were used in the present study, we have calculated a relative cost (RC) to make it more comparable in space and time. For this purpose, we used the best zone (lowest cost zone, MZ5) in MZ analysis as a reference zone and the other zones were compared with this zone based on the relative cost of changes suggested by Bakhshandeh et al. [38]. RC % of change was calculated by where X is the cost of chemical fertilizer consumption per tree in each zone and Y is the cost of chemical fertilizer consumption per tree in the best zone (MZ5).

Exploratory Data Analysis
The descriptive statistics of measured soil physicochemical attributes, elevation, and slope are given in Table 1. Bulk density, as the index of the soil compaction, showed a low coefficient of variation (CV) (8.92%), and no compaction was observed. Hosseini et al. [49] reported that surface disturbance increased the soil pore distribution, which could lead to a decrease in BD. Soil reaction was predominantly neutral to low alkaline with a pH value ranging from 5.11 to 8.35, as shown in Table 1. The mean of EC, SOC, clay, CEC, Mres4, and Cmic content was 251.37 (µS cm −1 ), 1.30%, 50.70%, 15.71 (meq 100 g −1 soil), 59.85 (mg CO2-C kg −1 soil d −1 ), and 395.63 (mg C kg −1 soil), respectively, as shown in Table 1. The range values of macronutrients including TN, P av , and K av contents were 0.06 to 0.45%, 10.17 to 84.0 (mg kg −1 ), and 49.87 to 277.76 (mg kg −1 ), respectively, as shown in Table 1. According to the optimum ranges for the measured soil properties and topography attributes proposed by Kangarshahi and Akhlaghi Amiri [27], based on mean values, no limitation was observed for BD, EC, CEC, K av , P av , and elevation. The mean soil pH was slightly higher than the optimum level, while the mean SOC and TN values were lower than the optimum levels, as shown in Table 1. The optimal range of nutrients and soil properties could lead to optimum plant growth [5]. On the other hand, the optimal soil nutrient level without concerning the other environmental factors will not ensure soil productivity [50].
Based on data dispersion proposed by Warrick and Nielsen [41], all data except EC, Mres4, Cmic, P av , elevation, and slope showed the distributions close to mean and median values, as shown in Table 1. The BD, pH, CEC, and clay showed low variability (CV < 12%), EC, SOC, TN, Mres4, Cmic, K av , and elevation were considered moderate variability (12 ≤ CV < 60%), while P av and slope showed high variability (CV ≥ 60%), as shown in Table 1 [41]. Several studies reported that soil phosphorus generally had more variability than many other soil properties [14,51,52], which is in line with our findings, as shown in Table 1. These CVs showed that there are considerable spatial variations in soil and topography attributes, as shown in Table 1. Therefore, MZ delineation may be a cost-effective method to apply fertilizers for decreasing the production cost in the study area. Table 2 shows Pearson's correlation coefficients between the measured soil properties and topography attributes in the study area. Significant correlations between almost all attributes except a few of them can be observed at the 5% probability level, as shown in Table 2 Table 2. Soil properties and topography attributes including Mres4 with pH and clay; CEC with pH and EC; elevation with TN and Cmic; and slope with pH and EC did not present significant correlations at the 5% probability level. When most of the properties showed a significant correlation, PCA could summarize the principal sources of variability in the data [14].  BD, soil bulk density (g cm −3 ); EC, electrical conductivity (µS cm −1 ); SOC, total soil organic carbon (%); Cl, clay (in the 0-25 cm soil surface layer (%)); TN, total nitrogen (%); K av , available potassium (mg kg −1 ); Mres4, soil respiration at the 4th day incubated at 30 • C (mg CO 2 -C kg −1 soil d −1 ); Cmic, microbial biomass carbon (mg C kg −1 soil); CEC, cation exchange capacity (meq 100 g −1 soil); P av , available phosphorus (mg kg −1 ). ** and * significant at 0.01 and 0.05 probability levels, respectively, and n.s is non-significant.

Spatial Variability Analysis
The accuracy of the interpolation method (e.g., IDW) for soil and terrain properties is presented in Table 3. Based on validation criteria, the highest and lowest accuracy were obtained for the elevation (R 2 = 0.92) and pH (R 2 = 0.73), respectively. The interpolation results for all terrain properties except for CEC (nRMSE = 13.8) interpolated as excellent prediction [45]. Furthermore, all properties interpolated with overestimation (positive CRM) in comparison to the observed data, as shown in Table 3. Table 3. Accuracy of the inverse distance weighting (IDW) method for interpolation of soil and terrain properties.

Variable
Weighting BD, soil bulk density (g cm −3 ); EC, electrical conductivity (µS cm −1 ); SOC, total soil organic carbon (%); Cl, clay (in the 0-30 cm soil surface layer (%)); TN, total nitrogen (%); K av , available potassium (mg kg −1 ); Mres4, soil respiration at the 4th day incubated at 30 • C (mg CO 2 -C kg −1 soil d −1 ); Cmic, microbial biomass carbon (mg C kg −1 soil); CEC, cation exchange capacity (meq 100 g −1 soil); P av , available phosphorus (mg kg −1 ); R 2 , coefficient of determination; nRMSE, normalized root mean square error; and CRM, coefficient of residual mass. Figure 2 shows the spatial distribution maps of all soil properties and topographic variables. BD shows no apparent spatial pattern with only small pitches in the northwest of the study area, as shown in Figure 2a. The spatial distributions of soil pH and EC, as shown in Figure 2b,c, and elevation, as shown in Figure 2l, were almost similar, with greater values of pH in the north where lower elevation occurred in the study area (r = −0.693, p < 0.01). The lateral downward movement of decalcified solution within the soils may lead to the enrichment of calcium carbonated in low-laying areas. Soil organic carbon, as shown in Figure 2d, and clay, as shown in Figure 2c, showed opposite spatial patterns (r = −0.654, p < 0.01), whereas SOC showed a similar pattern with TN, as shown in Figure 2f, and CEC, as shown in Figure 2g (r = 0.610, p < 0.01 and r = 0.590, p < 0.01, respectively). The upland soils, which are located in the margin of the study area and adjacent to existing forests, have higher SOC and CEC values, indicating the recent conversion of land use change from forests to citrus plantations. The K av values displayed a homogeneous spatial distribution pattern, as shown in Figure 2j. It could be attributed to the illite clay minerals that are the most clay-dominant minerals in the study area [53]. The P av showed low content in the west to northwest and almost homogeneous spatial patterns for the rest of the study area, as shown in Figure 2k.

Principal Component Analysis (PCA)
The PCA was used to reduce the dimensionality and summarize the variability of soil and terrain variables. Due to the high correlation between soil properties and topographic attributes, as shown in Table 2, PCA is an appropriate method to aggregate the effect of the parameters [14,44,54]. The PCs with eigenvalues greater than 1.0 could explain 78.66% of the dataset variability in the study area, as shown in Table 4. The PC1 explained the majority of the total data variance (35.52%) with a positive loading factor for EC, clay, TN, Mres4, and Cmic, and a negative loading factor for SOC and CEC. The PC2 explained 20.44% of the total variance and included BD, pH, and slope (negative loading factors), and Pav (positive loading factor). The PC3 included Kav, with a positive loading factor, which explained an additional 13.46% of the total variance in the dataset. The last PC (PC4) explained 9.25% of the data variability with a negative loading factor for elevation, as shown

Principal Component Analysis (PCA)
The PCA was used to reduce the dimensionality and summarize the variability of soil and terrain variables. Due to the high correlation between soil properties and topographic attributes, as shown in Table 2, PCA is an appropriate method to aggregate the effect of the parameters [14,44,54]. The PCs with eigenvalues greater than 1.0 could explain 78.66% of the dataset variability in the study area, as shown in Table 4. The PC1 explained the majority of the total data variance (35.52%) with a positive loading factor for EC, clay, TN, Mres4, and Cmic, and a negative loading factor for SOC and CEC. The PC2 explained 20.44% of the total variance and included BD, pH, and slope (negative loading factors), and P av (positive loading factor). The PC3 included K av , with a positive loading factor, which explained an additional 13.46% of the total variance in the dataset. The last PC (PC4) explained 9.25% of the data variability with a negative loading factor for elevation, as shown in Table 4. Therefore, the first four components (PC1, PC2, PC3, and PC4) were used to delineate the MZs for the commercial citrus plantation here. Behera et al. [42] delineated soil MZs of oil palm plantations by using PCA and fuzzy c-means in India and found that three PCs could explain 60.31% of soil property variation. Moreover, several studies have shown that PCA is an effective method to discriminate and delineate soils [44,55]. Therefore, the PCA analysis combined twelve input variables into four PCs, accounting for the majority of spatial variability in these properties. ; EC, electrical conductivity (µS cm −1 ); SOC, total soil organic carbon (%); Cl, clay (%); TN, total nitrogen (%); K av , available potassium (mg kg −1 ); Mres4, soil respiration at the 4th day incubated at 30 • C (mg CO 2 -C kg −1 soil d −1 ); Cmic, microbial biomass carbon (mg C kg −1 soil); CEC, cation exchange capacity (meq 100 g −1 soil); P av , available phosphorus (mg kg −1 ).

Clustering Analysis
A cluster analysis was applied to classify the four PCs into MZs. The MZA software was used to perform the fuzzy c-means cluster algorithm on the scores of the four PCs to define the optimum number of MZs. This method allows us to differentiate different zones with a similar value of properties and higher differences between them. To obtain the optimum number of MZs, the NCE and FPI parameters were calculated and plotted, as shown in Figure 3, against the number of clusters (or MZs) [46]. Therefore, when the FPI and NCE values were minimal, the optimal number of MZs can be identified, as shown in Figure 3. In this study, six different management zones were finally identified as the optimum number of MZs, as shown in Figure 3. Figure 4 shows the spatial distribution map of the six delineated MZs in the study area. In each delineated MZ, the measured soil properties and topographic attributes present the lowest variance and highest degree of membership [9]. Thus, in each zone, different management practices, such as nutrient application, can be carried out to increase crop production while decreasing the costs.
Several studies stated that the analysis of variance is an effective method to assess the differences between the delineated zones [9,14,15,17,42]. Therefore, the one-way ANOVA was carried out to evaluate the effectiveness of PCA and fuzzy c-means cluster algorithm combination to delineate MZs and also its spatial variability.  Figure 4 shows the spatial distribution map of the six delineated MZs in the study area. In each delineated MZ, the measured soil properties and topographic attributes present the lowest variance and highest degree of membership [9]. Thus, in each zone, different management practices, such as nutrient application, can be carried out to increase crop production while decreasing the costs. Several studies stated that the analysis of variance is an effective method to assess the differences between the delineated zones [9,14,15,17,42]. Therefore, the one-way ANOVA was carried out to evaluate the effectiveness of PCA and fuzzy c-means cluster algorithm combination to delineate MZs and also its spatial variability.
There were significant differences (p < 0.05) between all soil and terrain attributes except for BD among the six MZs, as shown in Table 4. Since the agricultural practices and machinery movement are similar in all MZ zones, the soil pore distribution (e.g., macropore and mesopore volumes) remains constant and did not show significant differences. Hosseini et al. [49,56] proved that agricultural practices and soil compaction were the main reasons for the BD variation and differences in the farms. The slope percentage only showed a significant difference in MZ4 in comparison to other MZs. The MZ4 with the lowest clay, TN, Kav, Cmic, CEC, and lower SOC and Pav   Figure 4 shows the spatial distribution map of the six delineated MZs in the study area. In each delineated MZ, the measured soil properties and topographic attributes present the lowest variance and highest degree of membership [9]. Thus, in each zone, different management practices, such as nutrient application, can be carried out to increase crop production while decreasing the costs. Several studies stated that the analysis of variance is an effective method to assess the differences between the delineated zones [9,14,15,17,42]. Therefore, the one-way ANOVA was carried out to evaluate the effectiveness of PCA and fuzzy c-means cluster algorithm combination to delineate MZs and also its spatial variability.
There were significant differences (p < 0.05) between all soil and terrain attributes except for BD among the six MZs, as shown in Table 4. Since the agricultural practices and machinery movement are similar in all MZ zones, the soil pore distribution (e.g., macropore and mesopore volumes) remains constant and did not show significant differences. Hosseini et al. [49,56] proved that agricultural practices and soil compaction were the main reasons for the BD variation and differences in the farms. The slope percentage only showed a significant difference in MZ4 in comparison to other MZs. The MZ4 with the lowest clay, TN, Kav, Cmic, CEC, and lower SOC and Pav There were significant differences (p < 0.05) between all soil and terrain attributes except for BD among the six MZs, as shown in Table 4. Since the agricultural practices and machinery movement are similar in all MZ zones, the soil pore distribution (e.g., macropore and mesopore volumes) remains constant and did not show significant differences. Hosseini et al. [49,56] proved that agricultural practices and soil compaction were the main reasons for the BD variation and differences in the farms. The slope percentage only showed a significant difference in MZ4 in comparison to other MZs. The MZ4 with the lowest clay, TN, K av , Cmic, CEC, and lower SOC and P av values could have lower potential soil fertility than others, as shown in Table 4, suggesting more chemical fertilizers are needed in this zone (particularly nitrogen). Thus, this MZ has a great potential for environmental risk via nitrogen leaching through the soil profiles and also nitrogen load in soil surface run-off with highest slope among the other zones. The MZ1, MZ2, MZ4, and MZ6 had a mean value of TN less than the limiting critical value for the citrus tree. Moreover, all of the delineated MZs had K av deficiency, while only MZ2 had P av deficiency, as shown in Tables 1 and 4 [27]. In general, the ranking was MZ5 > MZ2 > MZ6 > MZ1 > MZ3 > MZ4 in terms of soil fertility according to all the soil properties examined, as shown in Table 5. This finding is important for deciding the amount of fertilizers to be added to the soils in different zones. It emphasizes the needed site-specific management of citrus garden soils for sustainable citrus productivity and, in turn, for lowering the environmental hazards by ununiformed application of chemical fertilizers [26,57]. BD, soil bulk density (g cm −3 ); EC, electrical conductivity (µS cm −1 ); SOC, total soil organic carbon (%); Cl, clay (in the 0-25 cm soil surface layer (%)); TN, total nitrogen (%); K av , available potassium (mg kg −1 ); Mres4, soil respiration at the 4th day incubated at 30 • C (mg CO 2 -C kg −1 soil d −1 ); Cmic, microbial biomass carbon (mg C kg −1 soil); CEC, cation exchange capacity (meq 100 g −1 soil); P av , available phosphorus (mg kg −1 ). Means with the same letter are not significantly different at the 0.05 probability level based on the least significant difference test.

Economic Analysis of Delineated Zones
The economic analysis of fertilizer savings was performed in the MZs in comparison with uniform fertilization management. Based on the results, MZ5 and MZ4 had the lowest and highest total cost of chemical fertilizers per tree and animal manure, respectively, and the ranking was MZ5 < MZ2 < MZ6 < MZ1 < MZ3 < MZ4, as shown in Table 6. On the other hand, if MZ5 is defined as a reference zone due to the lowest cost needed, other zones, therefore, can be compared to this zone using the RC values. In fact, MZ4, MZ3, MZ1, MZ6, and MZ2 required 34.4, 30.6, 29.4, 9.77, and 9.44% more costs than MZ5 for achieving similar productivity, respectively, as shown in Table 6. The MZ delineation allows for the utilization of different doses of fertilizers, and cultural and irrigation practices as well that could lead to a decline in operational costs. Thus, this result highlighted that the combination of PCA and fuzzy c-means analysis could be used to effectively delineate MZs for soil nutrient management to maximize crop production and minimize its costs [42]. Our finding is also supported by other studies that the application of fertilizers can be based on MZ delineation [15,17,24,58]. Table 6. Economic analysis of the fertilizer savings in the management zones by comparing the difference between site-specific fertilization and uniform management. Although the MZ delineation based on soil and terrain attributes (topographic attributes) had good accuracy and was economically and environmentally sound, it is time-consuming and to some extent costly because of soil analysis [10]. Therefore, when soil data with reasonable distribution is not available, and also a time-efficient and cost-effective approach are needed, using a combination of remote sensing data and topographic attributes (environmental covariates with high resolution) can potentially be an alternative approach to delineate MZs in agricultural farms [59], but it needs to be verified in this region to decrease the methodology costs as much as possible.

Conclusions
The current study shows that the spatial variations of some selected soil properties and terrain attributes could be a possible approach to delineate soil MZs in commercial citrus plantations. The high correlation of soil properties and terrain attributes indicated considerable spatial variability and that site-specific nutrient management is needed in this region. The optimal number of MZs by PCA and fuzzy c-means clustering was six and a one-way ANOVA showed a significant heterogeneity variation among different MZs. In terms of cost-efficiency analysis, MZ5 and MZ4 showed the lowest and highest total operational cost. Therefore, the simplicity and cost-effectiveness of numerical methods, such as geostatistical tools, PCA, and fuzzy c-means clustering, could be an initial step to utilize fertilizers and reduce the cultural cost in precision agriculture management. However, this has to be performed using an inexpensive technique and input data because the profitability in these types of farms are very low, so owners need a cheap, effective, and reliable method to know which zones have similar production potential. Consequently, remote sensing and topographic variables could be suggested for future studies to reduce the cost of soil MZs more and more at the farm scale when high-resolution (e.g., 1 or 5 m) terrain attributes are available.