Delineation of Soil Fertility Management Zones for Site-specific Nutrient Management in the Maize Belt Region of Nigeria

Site-specific nutrient management can reduce soil degradation and crop production risks related to undesirable timing, amount, and type of fertilizer application. This study was conducted to understand the spatial variability of soil properties and delineate spatially homogenous nutrient management zones (MZs) in the maize belt region of Nigeria. Soil samples (n = 3387) were collected across the area using multistage and random sampling techniques, and samples were analyzed for pH, soil organic carbon (SOC), macronutrients (N, P, K, S, Ca and Mg), micronutrients (S, B, Zn, Mn and Fe) content, and effective cation exchange capacity (ECEC). Spatial distribution and variability of these parameters were assessed using geostatistics and ordinary kriging, while principal component analysis (PCA) and multivariate K-means cluster analysis were used to delineate nutrient management zones. Results show that spatial variation of macronutrients (total N, available P, and K) was largely influenced by intrinsic factors, while that of S, Ca, ECEC, and most micronutrients was influenced by both intrinsic and extrinsic factors with moderate to high spatial variability. Four distinct management zones, namely, MZ1, MZ2, MZ3, and MZ4, were identified and delineated in the area. MZ1 and MZ4 have the highest contents of most soil fertility indicators. MZ4 has a higher content of available P, Zn, and pH than MZ1. MZ2 and MZ3, which constitute the larger part of the area, have smaller contents of the soil fertility indicators. The delineated MZs offer a more feasible option for developing and implementing site-specific nutrient management in the maize belt region of Nigeria.


Introduction
Sustainable management of soil is of prime importance for increasing and sustaining crop yields globally and more importantly in sub-Saharan Africa (SSA), where food production is not at pace with the high food demand that resulted from the rapidly growing population [1,2]. Substantial area with high production potential in the SSA region has been degraded due to continuous cropping with poor nutrient management [3]. In the same vein, increasing demographic pressure has also resulted in the cultivation of marginal lands that are prone to fertility decline and other environmental degradations [3].
Use of the extended traditional fallow periods to restore soil fertility and enhance/sustain crop productivity is no longer feasible in the region. Consequently, the use of inorganic fertilizer to increase crops yield gained momentum among farmers in the region [4]. However, use of fertilizer in the region is being promoted through blanket recommendations in which the same rate is prescribed for heterogeneous soil and agro-ecological zones [5]. This practice leads to nutrient imbalance [6] and low fertilizer use efficiency [7], which affects the sustainability and productivity of farmlands [8].
Tropical soils (including those in SSA) are characterized by high spatial variability on both macroand/or micro-scales due to the combined effects of intrinsic (e.g., bio-physical and chemical processes) and extrinsic factors (e.g., crop management, fertilizer and tillage, among others) operating at different intensities and on different spatiotemporal scales [9,10]. Spatial variability in soils and the corresponding variability of yield response of various crops to nutrient application have already been observed across SSA [11][12][13][14]. Therefore, uniform fertilizer recommendation and application might result in over application in zones, fields, or soils with high nutrient status and under application in those with low nutrient status. Several studies in SSA [11][12][13]15] have recommended that nutrient management and fertilizer recommendation strategies should be tailored towards field-, site-, or soil-specific conditions to achieve balanced and effective fertilizer use and close nutrients related yield gaps. Matching the right fertilizer with the right recommendation is therefore regarded critical in optimizing and sustaining crop yield and economic returns.
The soil management (or soil fertility management) zoning approach is a popular and feasible approach to developing and implementing site-specific nutrient management strategies to manage spatial variability on macro-and micro-scales [10]. In this approach, sites (fields) are divided into variable zones that are relatively homogenous and can be used as zones or clusters for the development of site-specific nutrient management and recommendations. The nutrient management or soil fertility management zone (MZ) is defined as a technique utilizing the spatial management tools of precision agriculture, which represent a cost-effective approach to improving crop productivity, nutrient management, and reducing detrimental environmental impact [16]. Spatial delineation of MZs with their specific soil management recommendations, has proved capable of providing sustainable management decision solutions for natural resource management in many crop productions zones [17].
Various qualitative approaches, such as detailed soil survey maps [18,19], yield maps [20], and interpretation of remotely sensed images [21] have been used for delineation of MZs. However, these approaches were criticized as being either expensive, less precise, laborious, or localized. Currently, a quantitative approach using geostatistics to delineate the MZs is recognized for its ability to provide precise spatial information about soil properties, which aids decision-making on appropriate use of farming inputs and soil management decisions [17,22]. With the use of geo-statistics, soil property values at sampled locations are used to predict those of unsampled locations based on the spatial correlation behavior of that property [23,24]. When multiple soil factors are involved, the development of MZs using geostatistical analysis also requires the use of a multivariate statistical tools. Principal component analysis (PCA) and clustering are the multivariate methods commonly used to delineate MZs in most natural resource studies [25]. As a dimension reduction technique, the PCA identifies an orthogonal linear recombination of correlated variables that summarizes the principal sources of variability in the data [26]. Some researchers used either PCA [26] or multivariate cluster analysis [27,28] for delineating MZs. However, Tripathi et al. [10], Mohamed et al. [17], and others combined both techniques to delineate MZs by performing the multivariate clustering on the principal components (PCs) scores of soil fertility variables and demonstrated that this technique improves the accuracy and representativeness of MZ delineation.
In Nigeria, the maize belt region is the most intensified farming area [29]. The region is characterized as the most suitable zone for arable crop production (especially maize) in the country due to its intense solar radiation, adequate and well-defined rainy season, and low incidences of pests and diseases [30]. In contrast, the small content of organic matter and the sandy nature of the soils in the region increases vulnerability to degradation and nutrient leaching [31]. Currently, maize and other crops produced in Sustainability 2020, 12, 9010 3 of 18 the region are being largely managed by uniform "blanket" fertilizer recommendation(s) developed in the 1970s [30]. With our discussion above in mind, this results in nutrient imbalances in the soil through over-or under-application of fertilizer. To our knowledge, there is a dearth of information about the soil spatial variability and corresponding MZs covering the entire maize belt region of Nigeria. Moreover, site-specific nutrient recommendation decision support tools such as the Quantitative Evaluation of the Fertility of Tropical Soils (QUEFTS) [6] and the Nutrient Expert (NE) tool [5] have been developed for maize in the Nigerian maize belt. The challenge remains to develop and integrate the geospatial soil MZs so that the models can be efficiently used among peasant small-scale farmers who, in most cases, cannot afford soil testing of their fields. As such, this study aimed to develop a soil fertility information guide for informed policy making and input provision, especially for fertilizer and extension services in this region, by (1) assessing the spatial variability of soil properties in the maize belt cropping region of Nigeria using geostatistical analysis, and (2) delineating nutrient management zones using the combination principal component and multivariate cluster analyses.

Description of Study Area
The study covered 60 of 133 administrative Local Government Areas (LGAs) across 8 administrative States that fall within the maize belt cropping region of Nigeria ( Figure 1). The area encompasses five major agro-ecologies as follows: (i) derived savanna (DS), which is a transitional area found in the southern part of the area bordering the forest vegetation zone and the savannas; (ii) the southern Guinea savanna (SGS); (iii) northern Guinea savanna (NGS); (iv) Sudan savanna (SS); and (v) mid-altitude (MA). Altitude increases from less than 100 m in the DS to 1400 m in the MA. Dominant soils in the area are Entisols/Inceptisols, Alfisols, Ultisols, and Vertisols. Oxisols are also present but on a much smaller scale and restricted only to the extreme southern parts of the area [30]. Rainfall is unimodal and decreases in amount from 1800 mm in the derived savanna to 700 mm in the Sudan savanna, in the duration of seven to four months in the same direction as rainfall amount [30]. Mean annual minimum and maximum temperatures are 16 and 33 • C, respectively, throughout the entire region. and diseases [30]. In contrast, the small content of organic matter and the sandy nature of the soils in the region increases vulnerability to degradation and nutrient leaching [31]. Currently, maize and other crops produced in the region are being largely managed by uniform "blanket" fertilizer recommendation(s) developed in the 1970s [30]. With our discussion above in mind, this results in nutrient imbalances in the soil through over-or under-application of fertilizer. To our knowledge, there is a dearth of information about the soil spatial variability and corresponding MZs covering the entire maize belt region of Nigeria. Moreover, site-specific nutrient recommendation decision support tools such as the Quantitative Evaluation of the Fertility of Tropical Soils (QUEFTS) [6] and the Nutrient Expert (NE) tool [5] have been developed for maize in the Nigerian maize belt. The challenge remains to develop and integrate the geospatial soil MZs so that the models can be efficiently used among peasant small-scale farmers who, in most cases, cannot afford soil testing of their fields. As such, this study aimed to develop a soil fertility information guide for informed policy making and input provision, especially for fertilizer and extension services in this region, by (1) assessing the spatial variability of soil properties in the maize belt cropping region of Nigeria using geostatistical analysis, and (2) delineating nutrient management zones using the combination principal component and multivariate cluster analyses.

Description of Study Area
The study covered 60 of 133 administrative Local Government Areas (LGAs) across 8 administrative States that fall within the maize belt cropping region of Nigeria ( Figure 1). The area encompasses five major agro-ecologies as follows: (i) derived savanna (DS), which is a transitional area found in the southern part of the area bordering the forest vegetation zone and the savannas; (ii) the southern Guinea savanna (SGS); (iii) northern Guinea savanna (NGS); (iv) Sudan savanna (SS); and (v) mid-altitude (MA). Altitude increases from less than 100 m in the DS to 1400 m in the MA. Dominant soils in the area are Entisols/Inceptisols, Alfisols, Ultisols, and Vertisols. Oxisols are also present but on a much smaller scale and restricted only to the extreme southern parts of the area [30]. Rainfall is unimodal and decreases in amount from 1800 mm in the derived savanna to 700 mm in the Sudan savanna, in the duration of seven to four months in the same direction as rainfall amount [30]. Mean annual minimum and maximum temperatures are 16 and 33 °C, respectively, throughout the entire region.  Crop production in the region is significantly at subsistence level and is dominated by smallholders (75%) with an average land holding between 0.8 and 1.5 ha [32]. Major crops grown in the area are cereals, which include maize, sorghum, millet, and rice cultivated as sole, intercropped, or in Sustainability 2020, 12, 9010 4 of 18 rotation mostly with legumes such as cowpea, soybean, and groundnuts [33,34]. Production is rainfed, with little commercial irrigation for production of high value crops like vegetables, rice, and wheat.

Soil Sampling and Laboratory Analyses
A stratified multistage sampling procedure on a 1 km 2 resolution Discrete Global Grid (DGG) of the study area was used for soil sampling site selection. Sixty (60) of the 133 LGAs in the study area were selected based on the criterion that they contain more than 3 unique cropland sentinels. In each of the 60 selected LGAs, a 10 × 10 km sentinel was further selected at random. In each of the 60 selected sentinels (1 in each LGA), 10 grid cells of 1 × 1 km were then randomly selected. Finally, 5 (100 × 100 m) grid cells were randomly selected from the geographical centroid of the 10, 1 × 1 km grid cells in each of the 60 sentinels (resulting in 3000 geo-referenced points) and used as sites for soil sampling. Soil samples were collected in 2016 by sampling from 4 points at each identified site at a depth of 0-20 cm using a soil auger. Other soil data used in this study were obtained from published (95 data points described by Shehu et al. [13]) and unpublished data sets (292 points) collected between 2016 and 2018 by the project "Taking Maize Agronomy to Scale in Africa" (TAMASA) in Nigeria. Predictions based on AfSIS Bruker-Alpha KBr Mid InfraRed (MIR) spectral models were used for individual Mehlich-3 extractable nutrient contents as well as pH, effective cation exchange capacity (ECEC), SOC, and total N (N tot ) for the 3000 soil samples. The procedure for soil sampling and laboratory analyses of both the published and unpublished datasets obtained from the TAMASA project remains the same as described by Shehu et al. [13].

Descriptive Statistics
The data for each of the measured soil parameters were first subjected to descriptive statistical analysis (which include mean, standard deviation, minimum, maximum, lower and upper quartiles, skewness, kurtosis, and coefficient of variation (CV)). Correlation analysis was also performed to understand the relationships among the soil properties using JMP Pro version 14 statistical software.

Geostatistical Analysis
Geostatistical analysis of the soil parameters was performed in Geostatistical Analyst extension module of ArcGIS 10.4.1 software. Ordinary kriging procedure was used to estimate interpolated values of each soil parameter at non-measured points based on Equation (1) [35]: where γ(h) represents semivariance at lag distance (h), N (h) is the number of measured pairs separated by (h), and Z(xi) and Z (xi + h) are sampled measured values at i th sampling point separated by vector (h). Before fitting the semi-variogram model, the skewed values of soil properties were log-transformed to assume a near normal distribution. To select the best fitting model, different variogram models, namely, spherical, exponential, stable, circular, and Gaussian, were fitted to each soil parameter. Each model was evaluated for interpolation bias and accuracy based on three cross validation indices: root mean square error (RMSE), index of agreement (d-index), and coefficient of determination (R 2 ) (Equations (2)-(4)).
Sustainability 2020, 12, 9010 where M and P are respectively the measured and predicted values at ith point, M is the mean of the measured values, P is the mean of the predicted values, and n is the number of observations. The RMSE calculates the degree of model error of estimation; thus, a lower RMSE value indicates better model performance [36]. The d-index, as recommended by Willmott [37], and R 2 are non-dimensional values in the range of 0-1. Values of zero (0) represents no agreement between measured and predicted values, while 1 represents perfect agreement between measured and predicted values [38]. The d-index represents the ratio of mean square error to potential error. The coefficient of determination (R 2 ) estimates the combined dispersion against the single dispersion of the measured and predicted series [39]. The best fitting model for each parameter is that with both lowest RSME, highest d, and highest R 2 values and is thus selected because it has the highest prediction accuracy. Therefore, the model with higher prediction accuracy for individual soil property was selected for the spatial analysis.
C 0 is the nugget effect; C 1 is the partial sill of the semivariogram.

Principal Component Analysis
Because of the large volume of data points, multiple numbers, and co-linearity of the variables involved in this study, there might be the tendency for less stable variables to result in less reliable eigenvector values. This may consequently compromise expressions of the results. Based on that, principal component analysis (PCA) was used to summarize major variations contained in the whole dataset to reduce the number of uncorrelated dimensions, called the principal components (PCs). Use of PCA to discard correlated or redundant variables was described in McCabe [41]. Given that each variable is standardized to obtain a variance of 1 in PCA, only PCs with eigenvalues ≥1 were selected, so that the variability explained by each selected PC is greater than that attributed to individual parameters. The scores of the selected PCs were used to delineate the management zones.

Multivariate K-means Clustering
The scores of the selected PCs were used to develop MZs based on K-means cluster analysis performed using JMP Pro version 14. As an unsupervised machine learning technique, a cluster algorithm identifies sub-groups of variables in a dataset that have similar patterns and groups them as clusters. The major aim of clustering is therefore to create homogeneity among members of the same cluster while maximizing heterogeneity among clusters. In K-means clustering, the optimum number of clusters (K) that defines the dataset is selected when the algorithm stops detecting new trends in the data. At this point, the cubic clustering criterion (CCC) or the stopping criterion is high [42]. Finally, analysis of variance (ANOVA) was performed to compare the soil contents of each parameter in the different clusters (MZs). Mean values with significant differences (p ≤ 0.05) among the MZs were compared and separated using Tukey's HSD (Honestly Significant Difference) test.

Exploratory Statistics of the Soil Properties
The descriptive or exploratory statistics summarizing the status and spatial variability of the soil properties are presented in Table 1. Except the soil pH, significant variability in the spatial distribution of the soil properties within the study area (as indicated by the high coefficient of variability) have been observed. The average soil pH (5.99) in the study area fell within the moderately acidic category, and this level is considered suitable for most crops [43]. The average content of SOC in the study area (i.e., 0.79%) is rated low (according to the Esu [44] soil fertility classification of Nigerian Savanna soils), despite the SOC content varying from 0.11 to 3.85%. The low average level of SOC might be linked to the low application of organic manure as most farmers in the area seldom leave crop residues on their fields during harvest.
Average total N {N tot } (i.e., 0.65 g/kg) and available P (i.e., 4.29 mg/kg) in the area were classified as low according to the National Special Programme for Food Security (NSPFS) [45] soil fertility classification. Low N and P contents in some parts of the study area have been similarly reported in previous studies by Ekeleme et al. [46], Kamara et al. [47], and Shehu et al. [48]. The coefficient of variation (CV) of exchangeable cations (Ca, Mg, Na and K) in the area ranged between 52 and 67%. The average contents of exchangeable cations fall within the moderate fertility class (except K, which falls within high fertility status) according to Esu [44] fertility ratings of Nigerian savanna soils. This finding corroborates those of Shehu et al. [6] and Jimoh et al. [49] who also reported moderate to high levels of exchangeable cations, which can be related to the soil development process of the area. Another study in this area by Møberg and Esu [50] reported similar findings and attributed the results to the presence of large deposits of K-bearing Feldspar minerals in most soils in the area.
Distribution of micronutrients Fe, Mn, Zn, Cu, and B exhibits moderate to high spatial variability in the area based on CV > 30% according to Wilding [51] classification. Average values for Fe (129 mg/kg) and Mn (68.34 mg/kg) are very high, while that of Zn (1.18 mg/kg) is moderate according to Esu [44] soil fertility classification. These findings are contrary to those of Shehu et al. [6], where all these nutrients were observed to be high in the area. Kparmwang and Malgwi [52] similarly found moderate levels of these nutrients in the same area. The content of B in the study area is small (0.085 mg/kg) and varies over a range of 0.01 to 0.77 mg/kg. Copper also shows high variability (ranging between 0.05 and 30.57 mg/kg). The moderate to high spatial variability observed in the distribution of most of the studied soil nutrients re-justified the adoption of site-specific nutrient management to ensure balanced nutrient supply, improve nutrient limited yield, and maintain good soil health [53].

Relationships between the Soil Properties
The result of correlation analysis among the soil properties is presented in Table 2. Significant positive correlation was observed between SOC and other soil properties (with the exception of soil pH), though the strength of the correlation strongly varied among the soil properties. The soil properties with strong positive correlation (i.e., correlation coefficient >0.50) with SOC are N tot , S, and B. This implies the importance of soil organic matter in soil nutrient supply and bioavailability, especially in areas with high sand content. Moreover, the ECEC is largely dependent on the soil's organic matter (correlation coefficient 0.48). Studies by Zingore et al. [54] and Vanlauwe et al. [55] have recommended the rotation of cereal crops with legumes, crop residue preservation, and the application of well managed manure among others as measures to increase the SOC in SSA. A significant but weak negative correlation was observed between soil pH, SOC, and N tot as well. Likewise, a weak positive correlation was found between the soil pH and Ca, Mg, and available P (AP). This might indicate the potential for rainfall (which commonly comes in high intensity in the study area) to substantially leach out Ca and Mg, thus decreasing the soil pH. Meanwhile, the moderately acidic nature of the soils (as explained above) and increase in soil pH can enhance the AP, thus justifying the positive relationship between soil pH and AP. It has been documented that only at a pH above 7.5 can an increase in soil pH reduce AP due to the formation of calcium phosphates.  SOC: soil organic carbon; N tot : total nitrogen; AP: available phosphorus; ECEC: effective cation exchange capacity. ** Correlation is significant at 1% probability level. * Correlation is significant at 5% level of probability.

Geospatial Analysis of the Soil Properties
As shown in Table 1, the data distribution of all the soil properties was skewed, and this indicates a lack of normality in the distribution. Based on that, the data was log-transformed before performing geostatistical analysis. Before choice of model, different models (stable, exponential, Gaussian, spherical, and circular) were evaluated for prediction accuracy based on cross validation using RMSE, d-index, and R 2 as earlier presented above in the materials and methods section. Summary of the geostatistical characteristics of the soil properties is presented in Table 3, while cross validation results are shown in Figure 2. For most of the properties, the stable model function had the highest prediction accuracy and was therefore selected. Other model functions selected were exponential for S, Cu, Fe, and Na; circular for P; and spherical for Zn. However, studies of Jiang et al. [56], Tripathi et al. [10], and Mohamed et al. [17] found that the spherical function is the best fitting for most environmental variables. The nugget-to-sill ratio measures the spatial dependence of variables [40]. A ratio of ≤0.25 means that a large portion of variance was introduced due to intrinsic/localized soil factors such as parent material and other characteristics inherent to the soil. Ratio ≥0.75 indicates that variance was largely a result of extrinsic factors, such as cropping practice, fertilizer application, and farming systems, which vary among farmers [10], while the ratio between 0.25 and 0.75 indicates that variance was a result of the interaction of intrinsic and extrinsic factors. As presented in Table 3, overall, the nugget-to-sill ratio ranged from 0.000 (for available P) to 0.719 (for Ca), while overall range values varied from the smallest (1100 m) for K to the largest (9515 m) for Ca. This indicates that the spatial distribution of most studied soil properties possesses high spatial dependence across the area of study. The range value of semivariogram determines the distance after which spatial dependence ceases to exist. Therefore, the smaller range and nugget-to-sill ratio values of K mean that K variations in the soil are strongly due to inherent soil and environmental factors. Moreover, the high range and nugget-to-sill ratio values of Ca hint to artificially or anthropogenically induced variability [40], probably due to farming practices and soil management history.
A map of the study area showing the spatial distribution of the soil variables is presented in Figure 3. As observed through the significant correlation between the variables (Table 2) earlier presented above, the distribution pattern in the study area among the soil properties likewise followed a similar trend (i.e., properties with positive correlation showed a similar spatial distribution pattern in the study area, and vice versa). For instance, the highest correlation between SOC and N tot resulted in a very similar distribution pattern. Equally, variables that have a negative correlation with pH showed high values where pH was low, and vice-versa. This is consistent with results reported by Mohamed et al. [17], which showed a negative correlation between SOC, N, P, and K with pH. Except for Cu and B, most of the soil properties values were lower in the northern part of the study area, which constitutes parts of northern Guinea and Sudan savannas. Levels of SOC, N tot , S, K, Ca, ECEC, and Mg were higher around the central part of the study area, which is characterized by mid-altitude/mountainous vegetation. The higher content of these nutrients (i.e., SOC, N tot , S, K, Ca, ECEC, and Mg) around the mountainous vegetation area could be related to its higher altitude (i.e., >1400 m above sea level). The higher altitude area is associated with high rainfall and low temperature that slows the rate of soil biological activities and thus creates rich nutrient contents in the soil. Lower nutrient levels in the extreme south-western part of the area could be explained by their proximity to the river, low altitude (<100 m above sea level), and higher rainfall (>1500 mm), which encourage heavy nutrient leaching beyond the crop rooting zone. A map of the study area showing the spatial distribution of the soil variables is presented in Figure 3. As observed through the significant correlation between the variables (Table 2) earlier presented above, the distribution pattern in the study area among the soil properties likewise followed a similar trend (i.e., properties with positive correlation showed a similar spatial distribution pattern in the study area, and vice versa). For instance, the highest correlation between SOC and Ntot resulted in a very similar distribution pattern. Equally, variables that have a negative correlation with pH showed high values where pH was low, and vice-versa. This is consistent with results reported by Mohamed et al. [17], which showed a negative correlation between SOC, N, P, and K with pH. Except for Cu and B, most of the soil properties values were lower in the northern part of the study area, which constitutes parts of northern Guinea and Sudan savannas. Levels of SOC, Ntot, S, K, Ca, ECEC, and Mg were higher around the central part of the study area, which is characterized by mid-altitude/mountainous vegetation. The higher content of these nutrients (i.e., SOC, Ntot, S, K, Ca, ECEC, and Mg) around the mountainous vegetation area could be related to its higher altitude (i.e., >1400 m above sea level). The higher altitude area is associated with high rainfall and low temperature that slows the rate of soil biological activities and thus creates rich nutrient contents in the soil. Lower nutrient levels in the extreme south-western part of the area could be explained by their proximity to the river, low altitude (<100 m above sea level), and higher rainfall (>1500 mm), which encourage heavy nutrient leaching beyond the crop rooting zone.

Principal Component Analysis (PCA)
The results of correlation analysis in Table 2 indicated significant correlation among the soil variables. For this reason, the PCA was performed to summarize the variability in the data into principal components (PCs). The number of PCs produced (Table 4) is always equal to the number of variables involved in the analysis. Fifteen PCs were produced, out of which the first 4 that together explained 62.14% of the variation were retained for further analysis based on having an eigen value ≥1. As explained by Sharma [57], a PC with eigen value ≥1 explains more variance than an individual attribute. The maps of the PCs are shown as Figure 3. The result of component is presented in Tables 4 and 5. Principal component 1 explained 34.01% of the total variance in the data (Table 4), and is associated with SOC and related soil parameters, such as available N, K, Ca, S, ECEC, Mg, and B ( Table 5). The additional 11.31% of the variance was explained by PC 2, which seems to be associated with pH and related parameters (such as Ca). PC 3 explained another additional 9.26% and was loaded by available P but with negative correspondence to K and positive correspondence to Cu. PC 4 was most loaded by available P but with positive correspondence to K and negative correspondence to Cu, explaining the additional 7.55%. The kriging interpolation map of the 4 retained PCs is shown in Figure 4.

Multivariate Clustering and Delineation of MZs
K-means cluster analysis was done on the scores of the four selected PCs. The optimum number of clusters was selected based on attainment of high cubic clustering criterion (CCC). High CCC in this analysis was attained by four clusters. The clusters show clearly segregated zones when mapped apart from MZ4, with distinct characteristics ( Figure 5). MZ1 is the third largest and is located predominantly around the central part of the study area. MZ2 has the largest land proportion of all the MZs and seems to cover the central parts, whereas MZ3, the second largest, is found at the outer sections of the study area. MZ4 comprises of only a very small area and seems to be in the same form as spots within MZ3. The ANOVA results (Table 6) show that MZ1 is characterized by a relatively high content of most soil variables especially SOC, N tot , S, Ca, Mg, K, and ECEC (except AP). MZ4 is characterized by relatively high pH, AP, and Zn contents. MZ1 and MZ4 both have relatively high micronutrient concentrations (Zn, B, and Mn) compared with the other two MZs (i.e., MZ2 and MZ3). Both these zones (i.e., MZ1 and MZ4) represent the more fertile parts of the study area, with MZ4 probably being the more fertile because of the higher availability of P but representing a very small portion of the study. MZ2 and MZ3 do not differ significantly. Both zones have low SOC and ECEC, and are consequently relatively low in K, Ca, and Mg. MZ2 has a slightly lower pH and clearly higher Cu (1.20 mg/kg) and Fe content (164.6 mg/kg) in the soil compared with MZ3. The overall soil organic matter content (SOC) is low in all MZs, and thus farmers in these areas are recommended to increase soil organic matter content through cereal-legume rotation, application of manure, and preserving crop residues on the field after harvest. For MZ2 and especially MZ3, the SOC is extremely low, portraying high degradation and low productivity of the soils in these two zones. This can be related partly to inherently low nutrient contents in the soils and continuous cropping of high nutrient demanding crops without adequate application of organic or inorganic fertilizers [58].
In the same trend, AP content is generally small across all the MZs. Consequently, high P rates are required for P demanding crops like maize [6,13]. However, available compound fertilizer for most cereals is in the form of either NPK 15:15:15 or 20:10:10, i.e., containing the same ratios of P and K. Therefore, using these fertilizer blends to target high P and low K application rates, as in the case of MZ2 which has a lesser K requirement, is not advisable. Therefore, compound fertilizer with low K and high P contents is recommended for this zone. The low P content in the area can be related to intrinsic factors, such as the result of soil formation from low weatherable parent materials [48], or extrinsic factors, such as intensive crop cultivation [29], limited use of P fertilizers [59], and complete removal of crop residue that is commonly practiced by farmers in the area [32]. Sulphur content in MZ3 (the second largest zone) is low, and to enhance crop productivity and ensure nutrient balance, application of S might be required in this zone. Additionally, micronutrients (Zn, B and Mn) are low, especially in MZ2 and MZ3. Response of maize to S and other micronutrients has been reported in some parts of SSA, including the study area [56][57][58][59][60][61][62]. This might also point to the need for application of these micronutrients in MZ3 and MZ4 to sustain yield and prevent nutrient imbalances. However, we first recommend further research in these zones to understand the response and factors affecting the bioavailability of these three micronutrients (i.e., Zn, B, and Mn). 020, 12, x FOR PEER REVIEW Figure 5. Soil management zones for the study area.

Conclusions
Moderate to high spatial variability was observed in most of the studied soil properties in the study area. Spatial variation of pH and macronutrients (N, P, and K) was largely influenced by indigenous or inherent bio-physical factors, while variation of most micronutrients was significantly influence by both indigenous and artificially human-induced soil management factors. This emphasizes the need for the development of site-or area-specific soil management advises. Combined use of geostatistical, principal component, and multivariate K-means cluster analyses was successful in delineating soil nutrient management zones with a unique nutrient requirement in this study. Four MZs were delineated and differed significantly in spatial extent and distribution, as well as respective nutrient levels. Overall, soil organic carbon, total nitrogen, and available phosphorus are the most deficient across all the MZs, and therefore their management is of critical importance for sustainable crop production in the area. Therefore, the identified geospatial management zones, if integrated into the already calibrated and validated site-specific nutrient management (SSNM) and decision support tools such as nutrient expert (NE) and QUEFTS (Quantitative Evaluation of Fertility of Tropical Soils) in the study area when field soil information is not available, can offer a more feasible means of implementing SSNM. This is because most of the farmers in the study area are small holders with inadequate access to laboratory facilities to periodically conduct soil tests on their fields in a quest to use developed SSNM decision support tools. Moreover, the differences in nutrient contents of the MZs indicates the need for informed decisions by policy makers and other stakeholders regarding nutrient management.