Potential of Precise Fertilization through Adoption of Management Zones Strategy to Enhance Wheat Production

Abstract

The variability in soil properties and crop yield can be overcome by adoption of smart farming practices through interpolation and mapping of spatial variability patterns. Geospatial technologies can be utilized to determine the cause of spatial variability in fields for site-specific application of fertilizer. This study was designed to quantify and identify the spatial variation in soil properties and wheat (Triticum aestivum L.) yield and to delineate prescription maps for precise application of fertilizer in a semi-arid subtropical region of Pakistan. To examine the variability in soil properties on the production of the considered crop, this study comprised two different fields and each field was divided into (20 × 20 m) grids. The samples of soil were collected at 15 cm and 30 cm soil depths before the fertilization to analyze the different soil characteristics i.e., nitrogen (N), electrical conductivity (EC), potassium (K), soil organic matter (SOM), phosphorus (P), and pH. The boundaries of selected fields and grid points were established with a real-time kinematics-global positioning system (RTK-GPS). The soil data were acquired with a soil proximal sensor at a depth of 7 cm after fertilization. The statistical analysis coefficient of variation (CV), geostatistical-analysis-nugget-to-sill ratio (N:S), and the interpolated maps (ArcGIS pro 2.3) were used to characterize the least to moderate variability of soil parameters and yield, demanding site-specific management of fertilizer application. Cluster analysis was conducted using Minitab 21, which classified soil and yield characteristics into five categories: “very good”, “very low”, “good”, “poor”, and “medium”, with an external heterogeneity and internal homogeneity both more than 60%. Significant relationships (p < 0.05) between soil and crop properties were used to develop the management zones (MZs) for the precise application of fertilizer in wheat fields. Significant differences (p < 0.05) in soil nutrients were found in the very high and very low productivity zones at both sampling times, which suggest delineating the MZs for precise application of fertilizer according to the need of crop and soil properties. The results revealed that the optimum number of MZs for the wheat fields was five and there was heterogeneity in the soil nutrients in five MZs. The findings of this study also highlight the necessity of predicting the crop and soil factors by using precision technologies to develop the prescription maps, because sampling and analysis of soil are expensive and time-consuming. Based on the demand of the soil and crops, site-specific fertilization can increase economic and environmental efficiency.

1. Introduction

The world population is projected to increase by 35% in 2050, which increases food production demand; therefore, sustainable agricultural production systems must be developed to feed the growing population [1]. Wheat (Triticum aestivum L.) is the major crop grown globally and the most consumed staple crop in Pakistan with an annual production of approximately 26 million tons in 2021 [2,3]. Pakistan’s economy strongly depends on the wheat crop as it contributes 9.6% of the value added in agriculture and 1.9% of the gross domestic product (GDP) [4]. However, farmers are still adherent to use traditional farming practices of fertilizer application without considering the soil and crop variability within the field that can reduce crop performance and yield, especially in the arid areas (dry). Fertilizer consumption has been significantly increased in the last few decades in order to increase productivity and meet production goals [5,6]. However, excessive fertilization cannot only increase the cost of farm inputs, but it can also pollute the environment. Fertilizer application rates are dependent on soil condition, field moisture, crop types, and other factors [7]. Soil texture variability plays an important role in controlling within-field moisture variability and nutrients holding capacity especially in the arid regions [8]. Therefore, site-specific crop and soil management plays a crucial role in nutrients use efficiency (NUE), profitability, and environmental protection [9].

In Pakistan, fertilizers are applied uniformly without consideration of the spatial variability in soil properties and crop parameters that could result in yield variation. Therefore, the soil and crop variability issues can be addressed by precise application of fertilizer [10]. However, proper characterization, quantification, identification, and development of the management zones (MZs) for site-specific applications of fertilizers are the basic steps for addressing spatial variability using precision agriculture technologies [11,12,13].

The physicochemical properties of soil have a great effect on plant growth and the variability of yield [14,15]. Nitrogen (N), phosphorus (P), and potassium (K) are three primary macronutrients [16] essential for plant growth. In an arid cropping system, due to unpredictable rainfall, in terms of both amount and intensity, the unavailability of water, excessive use of fertilizer, and soil erosion can decrease the availability of macronutrients in the soil and increase the environmental risks [17]. Soil in Pakistan is unable to support significant levels of productivity due to low organic matter that is essential for water holding capacity and nutrients, and improves soil structure [18,19]. Recently, soil maps have been developed for all the districts of Pakistan that is a feasible strategy for enhancing crop production and improving soil health [20]. The variable-rate (VR) fertilization could control the main factors, which are responsible for this variation [21].

Digital mapping is an easy and practical method to classify the soil parameters into low to high ranges, which can be visualized easily using different color patterns. These soil properties are used to develop the productivity zones for the precise application of fertilizer. Site-specific nutrient management (SSNM) has been developed on the bases of crop and soil properties requirements. SSNM is the main component of PA as it is used to improve NUE, especially in crops such as rice, maize, and wheat. Soil proximal sensors, GIS, remote sensing, simulation models, decision support systems, and machine learning for variable application of nutrients play an important role in SSNM [22]. This is consistent with SSNM, which promotes the 4R’s of nutrient stewardship (Right product, Right time, Right place, and Right rate), as shown in Figure 1. Studies have revealed that 40% of agrochemicals can be saved by correctly distinguishing variability based on the MZs, and by using these MZ zones, 75% of pesticides can also be saved [23,24].

The fertility of the soil is decreasing due to conventional management and agricultural intensification in the last few decades. Many studies have reported the effects of spatial variation on crop production. However, there is very little work conducted on addressing the issue of spatial variability through management zone delineation especially in wheat crop. It is necessary to provide the right amount of fertilizer according to the need of the crop and it can be applied using prescription maps, which represents the variation within the field. These precision technologies are much needed in developing countries such as Pakistan. However, the spatial variability in soil and crop and its impact on yield need to be studied in arid regions. Therefore, the objective of this research was to characterize, quantify, and identify the spatial variability in soil properties and crop yield and to delineate the MZs for precise application of fertilizer. These developed MZs would help in avoiding the over- and under-application of fertilizer. The crop production will be enhanced in the Potohar region by applying fertilizer according to developed MZs.

2. Materials and Methods

2.1. Study Area

Two wheat (Triticum aestivum L.) fields were selected in (33°11′66″ N, 73°01′11″ E) Koont Research Farm, in PMAS-Arid Agriculture University Rawalpindi, Pakistan. The selected fields with two different varieties of wheat were Field 1—Barani 17 (0.4 ha, 33°6′54″ N, 73°0′50.4″ E) and Field 2—Fakhae e Bhkkar (0.4 ha, 33°6′58″ N, 73°0′59″ E) during the wheat growing season 2021–2022 to assess the variation in crop yield and soil properties (Figure 2). The study area was located in an arid region, which only depends on rainfall as the source for crop water requirement. The wheat has been growing in the selected fields for the past few years. The soils at the selected sites are sandy loam with a low organic matter and the pH ranges between 7.0 and 8.27, which is optimum for maximum crop growth [25]. The field boundaries and grid were demarked by a real-time kinematics global positioning system (RTK-GPS) (SinoG NSS^® T30 IMU, ComNav Technology Ltd., Shanghai, China) at each experimental site to collect soil and crop yield samples (Figure 2).

2.2. Sampling and Analysis

The grid pattern (20 m × 20 m) of sampling points was established at the selected fields to collect the soil and yield samples (n = 13 and n = 21 for Field 1 and Field 2 sites, respectively). The samples were collected twice before fertilization (November 2021) at 15 cm and 30 cm depths, and after fertilization (February 2022) at 7 cm depth. The first sampling before fertilization was performed using a hand auger from each georeferenced point and collected soil samples were placed in bags and labeled. The collected samples were first air-dried and then oven-dried for 24 h at 65 °C and sieved through a 2 mm polyethylene sieve to remove large particles. The collected samples were analyzed for Total Nitrogen (TN) [26], available phosphorus (P) [27], available potassium (K) [28], pH [29], EC [30], and soil organic matter [31]. The Hydrometer 152H [32] was used to measure the texture of the soil. Soil texture and pH do not frequently change, so these parameters were measured only once at the initiation of an experiment. However, for 2nd sampling, the proximal sensor was used for the analysis of soil, and it gives the value of pH, N, P, K, EC, humidity, and temperature with great accuracy, quick response speed, and stable output [10,33]. The MZs were developed for the crop’s need after these analyses were used to evaluate the variation in soil characteristics.

2.3. Yield Estimation

The wheat crop was harvested in the last week of April with the help of georeferenced coordinates at each grid by using a hand sickle of a quadrant of 1 m × 1 m. Grains from the debris were separated, including leaves, spikes, and spikelets from each grid, and weighed at the time of harvesting. Yield was obtained from each grid of crop fields and the total average yield was estimated.

2.4. Statistical and Geostatistical Analysis

The normality test was conducted by using the 5% significant level of the Shapiro–Wilk test and non-normal data were normalized by the z-score method. The minimum, maximum, mean, coefficient of variance (CV), standard deviation (SD), and skewness were computed by statistical software (Minitab, LLC 21.1.0) for descriptive statistics. The result of descriptive statistics provided details about the field’s soil variability but did not provide any information about spatial trends or patterns. Wilding et al. [34] classified the variability based on the CV values of soil properties: if CVs were <15%—least variation; from 15 to 35%—moderate variation; >35%—high variation. A semivariogram was used in ArcGIS pro 2.3 to perform geostatistical analysis for the investigation of spatial patterns. The author [35] assessed the variability pattern by using Equation (1) to calculate the semivariograms.

$$\mathsf{\Upsilon }(\mathrm{h})=\frac{1}{2\mathrm{n}\left(\mathrm{h}\right)} \mathrm{n} (\mathrm{h}) \sum \mathrm{i} = 1 {[\mathrm{V}(\mathrm{X}\mathrm{i} + \mathrm{h}) - \mathrm{V}(\mathrm{X}\mathrm{i})]}^2$$

(1)

where ϒ (h) denotes the semi-variance at a particular distance h; V (Xi) denotes the variable V value at the Xi location; n (h) is the number of pairs of sample points separated by the lag distance h. The semivariograms were used as indicators in this investigation, and the anisotropy of semivariogram parameters was calculated as described by using semivariogram modeling [36,37]. Therefore, anisotropy models were fitted theoretically by using ArcGIS pro software. The coefficient of (r) between soil parameters and crop yield was determined by using a correlation matrix.

2.5. Spatial Variability Mapping

The variability in the collected soil and crop data was analyzed using the Kriging interpolation technique in the ArcGIS environment. Exponential, Gaussian, linear, spherical, and circular models were used for the identification of the optimum Kriging model used for characterizing and quantifying the soil parameters, and crop yield data, the spatial correlation were represented using the root-mean-square error (RMSE) method [38]. The residual sum of squares (RSS) and coefficient of determination (R²) were used to confirm the best-fitted model according to the soil parameters. The RMSE was calculated using Equation (2).

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}\sum {\left\{ z\left(xi\right)-z\left(xi\right)\right\}}^2}$$

(2)

The prediction of the spatial values in soil properties of an un-sampled location was estimated by Equation (3) [39].

Z(X₀) = n∑_i = 1 λ i Z (X_i)

(3)

The spatial dependency in the soil parameter was defined by the nugget-to-sill ratio (N:S). If the (N:S) < 25%, the parameter was considered strongly spatially dependent; if the (N:S) was between 25 and 75%, the parameter was considered moderately spatially dependent; if the (N:S) was >75%, the variable was considered weakly spatially dependent [40]. The Kriging interpolated maps were produced to allow easier comparison and to observe the pattern of variability in the fields. The collected soil and crop data were arranged into natural clusters/groupings with prior knowledge of productivity potential [18].

The cluster analysis results were displayed using dendrograms, where the hierarchical relationship among the parameters is based on similarity levels and depicts the various groups as well as the differences between them and their nearest neighbors. For the comparison of the mean average of distinct productivity zones, the analysis of variance (ANOVA) statistical method was used to compare the mean with the least significant difference (LSD) at the 5% significance level. The results were used to organize yield data and spatial soil properties into naturally occurring clusters, to minimize the variation within each cluster and maximize the variance between clusters and ultimately developing MZs in ArcGIS pro-2.3.

3. Results and Discussion

3.1. Descriptive Statistics of Soil Properties and Wheat Yield

Results of the Shapiro–Wilk normality test of collected soil properties data for both sampling datasets (before and after fertilization) showed that all the collected parameters were normally distributed (p > 0.05) except EC, SOM, and total N. These parameters can be normalized by z-score, which is the deviation in data points that would lie above or below the mean [41]. The overall variation within the field assessed using the coefficient of variation (CV) [34] and the statistical analysis showed the least to a higher variation in yield and soil properties. The pH of the soil had a low variation (CV < 15%) for both sampling rounds. The SOM, N, and P showed high variations (CV > 35%) at the 30 cm sampling depth, as shown in

Table 1. Descriptive Statistics of Soil Properties and Crop Yield for Field 1.

Parameters	Min	Max	Mean	S.D	C.V (%)	Skewness	Kurtosis
First Sampling, November 2021 (0–15 cm) Before Fertilization
pH	6.9	7.7	7.3	0.3	4.2	0.5	−1.4
EC (dS m−1)	0.1	0.4	0.3	0.1	35.6	−0.1	−0.4
SOM (%)	0.4	0.7	0.5	0.1	18.0	0.5	−0.74
N (%)	0.04	0.07	0.06	0.007	12.7	0.21	−1.3
P (mg Kg−1)	3.7	9.0	6.3	1.9	29.5	0.0	−1.6
K (mg Kg−1)	22.2	27.1	25.0	1.9	7.5	−0.3	−1.5
M.C	16.3	30.2	22.9	4.7	20.5	0.02	−1.33
First Sampling, November 2021 (0–30 cm)
pH	6.4	7.6	7.1	0.4	6.1	−0.1	−1.5
EC (dS m−1)	0.3	0.6	0.5	0.1	20.1	0.0	−1.1
SOM (%)	0.0	0.2	0.1	0.1	53.2	0.0	−1.6
N(%)	0.002	0.042	0.050	0.021	42.0	−1.01	0.27
P (mg Kg−1)	2.5	7.8	5.2	1.8	35.4	0.0	−1.5
K (mg Kg−1)	7.2	12.1	10.1	1.9	18.7	−0.5	−1.5
Second Sampling, February 2022 (0–7 cm) After Fertilization
N(mg Kg−1)	12.1	18.2	14.6	1.9	12.8	−0.4	−0.9
P (mg Kg−1)	19.9	25.6	22.9	2.1	9.0	0.0	−1.6
K (mg Kg−1)	49.9	55.6	52.6	1.8	3.5	0.3	−1.5
pH	6.8	7.5	7.0	0.2	2.5	1.6	2.5
EC (dS m−1)	0.1	0.3	0.2	0.1	45.0	0.3	−1.0
M.C	17.7	37.1	24.1	5.6	23.3	1.14	0.6
Crop Yield sampling, April 2022
Yield (Kg ha−1)	925	996	971	23.57	2.4	−0.95	−0.01

M.C: Volumetric Moisture Content, EC: Electrical Conductivity, SOM: Soil Organic Matter, N: Total Nitrogen, N: Nitrogen, P: Phosphorus, K: Potassium.

, and N showed high variability for Field 2 for 15 cm and 30 cm depths. Mean values of M.C and available forms of N, P, and K values were higher in Field 1 as compared to Field 2 that might be due to the variation in soil properties and crop’s physiological growth parameters amid different varieties of wheat. The total N mean value was higher at 15 cm than 30 cm for both fields. The results of this study are supported by other studies conducted in wheat and wild blueberry cropping systems, where the spatial variation also found moderate to high CVs for these soil properties except pH [35,42], which may be due to the logarithmic nature of pH.

Soil EC at 30 cm depth (before fertilization) and 7 cm depth (after fertilization) showed high variations while, on the other hand, M.C showed moderate variation before and after fertilization. SOM at 15 and 30 cm depths showed moderate variation (

Table 2. Summary of Statistics of Soil Properties and Crop Yield for Field 2.

Parameters	Min	Max	Mean	S.D	C.V (%)	Skewness	Kurtosis
First Sampling, November 2021 (0–15 cm) Before Fertilization
pH	7.6	8.3	8.0	0.3	3.2	−0.2	−1.6
EC (dS m−1)	0.2	0.2	0.2	0.0	12.0	0.3	−1.2
SOM (%)	0.5	1.0	0.7	0.1	20.2	−0.3	−0.7
N (%)	0.002	0.07	0.04	0.02	36.4	−0.6	0.21
P (mg Kg−1)	0.7	7.8	4.4	1.9	42.1	−0.1	−0.7
K (mg Kg−1)	64.5	114.4	84.2	13.2	15.7	0.1	−0.5
M.C	15.3	31.2	22.2	4.16	18.7	0.24	−0.5
First Sampling, November 2021 (0–30 cm)
pH	5.7	7.8	6.9	0.5	7.4	−0.7	0.2
EC (dS m−1)	0.2	1.2	0.5	0.3	57.9	2.2	4.1
SOM (%)	0.1	0.6	0.4	0.1	39.8	−0.3	−0.7
N (%)	0.003	0.07	0.04	0.02	46.3	−0.5	−0.4
P (mg Kg−1)	0.3	5.2	2.4	1.3	55.6	0.8	−0.2
K (mg Kg−1)	49.5	99.4	69.5	13.0	18.8	0.1	−0.4
Second Sampling, February 2022 (0–7 cm) After Fertilization
N(mg Kg−1)	6.9	20.9	13.0	3.1	23.9	0.6	1.0
P (mg Kg−1)	9.9	27.9	17.8	4.0	22.7	0.6	0.8
K (mg Kg−1)	23.9	65.9	41.6	9.8	23.6	0.8	0.7
pH	6.5	8.0	7.0	0.5	6.6	1.0	−0.2
EC (dS m−1)	0.0	0.3	0.1	0.1	42.1	0.4	2.1
M.C	18.2	27.4	22.4	3.4	15.0	0.23	−1.5
Crop Yield sampling, April 2022
Yield (Kg ha−1)	702	928	793	53	14.9	0.55	0.51

M.C: Volumetric Moisture Content, EC: Electrical Conductivity, SOM: Soil Organic Matter, N: Total Nitrogen, N: Nitrogen, P: Phosphorus, K: Potassium.

). The EC showed high variation at 15 cm and 7 cm depths before and after fertilization (Table 2). The total N and SOM showed high variations at 30 cm (Table 1). The wheat yield in Table 1 and Table 2 had a CV between 2.4% and 14.2% that showed a weak variation. The average yield was poor in both fields due to the presence of more weeds and dry areas, different varieties of wheat, or less nutrient availability for plant growth and development. The main problem in the arid region is a lack of water availability, and water stress might have reduced the crop yield [43]. This study investigated the spatial variations in soil properties and yields. The M.C showed a similar pattern of change at both stages indicated by their CVs. All the parameters showed positive skewness for all depths except EC and K at 15 cm, pH and K at 30 cm, and N from 7 cm, while P results showed the same level of data skewness at both sampling rounds (Table 1). The EC was highly skewed at 30 cm for Field 2 and pH was highly skewed at 7 cm as depicted in Table 2. The SOM at the sampling depth of 0–30 cm showed a high variation (CV = 53.2%) in Field 1 and Field 2, and P was highly variable (CV = 55.6%) at the 0–30 cm depth. The temporal variation also affected the results, but this research focused on spatial variation in the arid region.

3.2. Geostatistical Variation in Soil Properties and Grain Yield

Geostatistical analysis was used to represent the spatial autocorrelation between crop yield and soil properties. The results of geostatistical analysis revealed that the low value of RSS and high value of R² were used to consider the best models as the literature suggested that the best-fit models for soil properties are exponential, circular, spherical, and Gaussian. The theoretical models for different soil properties and grain yield with the best-fit models were estimated as shown in Figure 3 and Figure 4, while the results are presented in

Table 3. Semivariogram Parameters of Soil Properties and Grain Yield for wheat Field 1.

Parameters	Nugget Co	Partial Sill C	Sill Co + C	Range (m)	Nugget Sill Ratio	R2	N/S	Model
First Sampling, November 2021 (0–15 cm)
pH	0.09	0.06	0.15	39.28	0.62	0.37	62.07	Exponential
EC (dS m−1)	0.01	0.01	0.02	60.8	0.35	0.03	34.73	Exponential
SOM (%)	0.00	0.00	0.00	12.56	0.54	0.06	54.42	Exponential
N (%)	0.67	1.58	2.25	16.65	0.30	1.42	29.74	Circular
P (mg Kg−1)	1.46	0.46	1.92	5.59	0.76	1.56	76.09	Spherical
K (mg Kg−1)	43.32	80.92	124.24	16.7	0.35	7.34	34.87	Linear
M.C	11.60	12.50	24.10	15.54	0.48	5.22	48.13	Exponential
First Sampling, November 2021 (0–30 cm)
pH	0.03	0.08	0.11	11.9	0.27	0.28	27.11	Exponential
EC (dS m−1)	0.01	0.01	0.01	42.6	0.39	0.05	38.73	Circular
SOM (%)	0.00	0.00	0.00	46.56	0.04	0.05	3.91	Exponential
N (%)	2.81	0.54	3.35	11.1	0.84	1.83	83.89	Gaussian
P (mg Kg−1)	2.01	0.07	2.08	15.5	0.97	1.55	96.63	Exponential
K (mg Kg−1)	50.95	59.19	110.15	16.7	0.46	7.91	46.26	Linear
Second Sampling, February 2022 (0–7 cm)
N(mg Kg−1)	0.36	0.50	0.86	15.86	0.42	1.25	41.91	Exponential
P (mg Kg−1)	5.54	0.08	5.62	13.32	0.99	2.49	98.58	Spherical
K (mg Kg−1)	8.34	30.21	38.55	8.6	0.22	4.75	21.63	Spherical
pH	0.11	0.04	0.14	11.3	0.75	0.38	75.47	Exponential
M.C	18.69	15.40	34.09	17.4	0.55	6.50	54.82	Gaussian
EC (dS m−1)	0.01	0.00	0.01	28.5	0.95	0.10	95.26	Spherical
Crop Yield sampling, April 2022
Yield (Kg ha−1)	242.00	589.86	831.86	18.7	0.29	0.79	29.09	Spherical

M.C: Volumetric Moisture Content, EC: Electrical Conductivity, SOM: Soil Organic Matter, N: Total Nitrogen, N: Nitrogen, P: Phosphorus, K: Potassium.

and

Table 4. Semivariogram Parameters of Soil Properties and Grain Yield for wheat Field 2.

Parameters	Nugget Co	Partial Sill C	Sill Co + C	Range (m)	Nugget Sill Ratio	R2	N/S	Model
First Sampling, November 2021 (0–15 cm)
pH	0.31	0.05	0.36	11.5	0.86	0.25	86.1	Exponential
EC (dS m−1)	0.00	0.00	0.0006	48.6	0.06	0.50	6.05	Spherical
SOM (%)	0.02	0.01	0.02	14.5	0.16	0.16	15.68	Exponential
N (%)	1.26	1.83	3.09	19.4	0.41	1.70	40.71	Gaussian
P (mg Kg−1)	3.56	0.87	4.43	13.17	0.80	1.98	80.37	Spherical
K (mg Kg−1)	24.00	142.93	166.93	16.14	0.14	1.10	14.38	Linear
M.C	9.66	11.6	21.3	9.71	0.45	2.28	45.35	Spherical
First Sampling, November 2021 (0–30 cm)
pH	0.2	0.09	0.27	18.05	0.67	0.53	67.11	Linear
EC (dS m−1)	0.0	0.09	0.10	41.2	0.04	0.24	4.20	Exponential
SOM (%)	0.0	0.02	0.03	33.5	0.04	0.16	3.99	Exponential
N (%)	2.1	1.44	3.55	18.05	0.59	2.13	59.37	Exponential
P (mg/Kg)	1.3	0.69	1.97	17.1	0.65	1.39	64.79	Spherical
K (mg/Kg)	74.9	157.62	232.57	16.6	0.32	1.10	32.23	Circular
Second Sampling, February 2022 (0–7 cm)
N(mg/Kg)	4.95	3.81	8.76	18.0	0.56	1.30	56.47	Spherical
P (mg/Kg)	10.45	4.75	15.20	17.52	0.69	1.31	68.73	Exponential
K (mg/Kg)	67.55	28.48	96.03	11.06	0.70	1.45	70.34	Gaussian
pH	0.21	0.04	0.25	16.6	0.84	0.48	84.14	Linear
M.C	28.69	7.00	35.69	22.4	0.80	1.26	80.39	Exponential
EC (dS m−1)	0.00	0.00	0.00	58.05	0.83	0.65	83.34	Exponential
Crop Yield sampling, April 2022
Yield (Kg ha−1)	2007	1069	3076.0	11.35	0.65	4.9	65.25	Spherical

M.C: Volumetric Moisture Content, EC: Electrical Conductivity, SOM: Soil Organic Matter, N: Total Nitrogen, N: Nitrogen, P: Phosphorus, K: Potassium.

. The variation in soil properties and yield might occur due to extrinsic and intrinsic sources. The intrinsic variation occurred naturally due to parent material but extrinsically occurred due to different management practices used within the field [44,45].

The geostatistical analysis showed that semivariogram models were exponential, spherical, and Gaussian models for Field 1, and for Field 2, exponential, spherical, linear, and Gaussian were found the best-fit models for their soil properties and yield. The properties of soil were inconsistent because the range of influence was greater than the distance of the grid interval. The results indicate that there were moderate to large variations within the field except pH, EC, and SOM for Field 1 (Table 3), and EC, SOM, and M.C for Field 2 (Table 4). The (N:S) ratio indicated the spatial dependency for the (0–15 cm) depth of soil samples, and that all parameters showed moderate spatial dependence (N:S > 25%) for Field 1 except P, as shown in Table 3. Total N from the first sampling (November 2021) and N from the second sampling (February 2022) exhibited medium spatial dependence, and this variation could be controlled by applying the fertilizer via a variable-rate application. The second sampling showed a moderate spatial dependence within the field as indicated by N: S except for SOM, which indicated a high variability as the N:S ratio was less than 25% and the range of influence was 46.56 m, as indicated in Table 3 for Field 1. The soil properties including EC, SOM, and K had less than 25% N:S, indicating high spatial dependence for Field 2 (Table 4). Total N and M.C both were moderately variable, and pH and P showed a weak dependence indicated by the N: S ratio for Field 2 (Table 4). The same pattern of variation in soil and yield was observed for the 30 cm and 7 cm depths. The yield showed moderate spatial dependence because N:S was greater than 25%, and the yield could be controlled by managing the spatial variation in soil properties.

The variation in soil properties occurred due to applications of fertilizer by the broadcast method. The distance (h) at which the semivariance reached its highest value (sill), which equated to the sample variance, was the semivariograms “range.” The high spatially dependent factors were controlled by intrinsic soil properties and the moderate and low spatial dependences were controlled by extrinsic soil properties. The yield represented high spatial dependency (N:S 29.09%), as depicted in Table 3 for Field 1, and moderate spatial dependency (N:S 65.25%), as shown in Table 4 for Field 2, of the selected area of the experimental site.

The soil variability was the major factor that affected the regional yield variation within the fields and the sampling technique and interpolated maps determined their spatial dependency. The sampling interval should be less than the range of the semivariogram [46] but sampling intervals in this study were more than the semivariogram range. Field 2 showed a large variation and had a low yielding value as compared to Field 1, indicated by their CVs.

3.3. Correlation of Soil Properties and Grain Yield

The significant relation between the grain yield and soil properties of the selected fields is presented in

Table 5. Correlation Matrix among the Soil Properties (0–15 cm) and Grain Yield for Field 1.

	pH	EC	N	P	K	TN	M.C
EC	0.26 NS
N	0.09 NS	−0.18 NS
P	−0.07 NS	0.48 **	−0.62 ***
K	−0.20 NS	0.41 **	−0.43 **	0.85 ***
TN	0.49 **	0.35*	−0.21 NS	0.32 *	0.26 NS
M.C	−0.24 NS	−0.54 ***	0.16 NS	−0.33 *	−0.02 NS	−0.22 NS
Yield	−0.11 NS	0.39 *	0.78 ***	0.81 ***	0.68 ***	0.41 **	0.31 *

* Significance of correlation, p = 0.05, ** Significance of correlation, p = 0.01, *** Significance of correlation, p = 0.001.

and

Table 6. Correlation Matrix among the Soil Properties (0–15 cm) and Grain Yield for Field 2.

	pH	EC	N	P	K	TN	M.C
EC	0.10 NS
N	−0.04 NS	0.85 ***
P	−0.03 NS	0.82 ***	0.90 ***
K	0.00 NS	0.80 ***	0.86 **	0.97 ***
TN	−0.26 NS	0.10 NS	0.67 ***	0.06 NS	0.07 NS
M.C	−0.12 NS	0.13 NS	−0.08 NS	−0.06 NS	−0.01 NS	0.30 *
Yield	−0.06 NS	0.47 **	0.68 ***	0.74 ***	0.58 ***	0.49 **	0.30 *

* Significance of correlation, p = 0.05, ** Significance of correlation, p = 0.01, *** Significance of correlation, p = 0.001.

. The soil properties were significantly correlated with yields except for the soil pH of both fields. The correlation results indicated a linear trend between EC with P, K, and N that indicates that the EC values were influenced by the soil properties. The EC negatively correlated with M.C (r = −0.54) for Field 1, whereas the relationship of EC with M.C was non-significantly correlated for Field 2. There was a positive correlation between EC and yield, which indicates that higher values of EC had high yielding areas and vice versa. Mann et al. [47] investigated that the spatial variation in soil physiochemical properties affects the citrus yield and the finding of this research article agrees of this study.

The yield was significantly related to EC (r = 0.39 to 0.47), N (r = 0.68 to 0.78), P (r = 0.74 to 0.81), K (r = 0.58 to 0.68), Total N (r = 0.41 to 0.49), and M.C (r ~0.30–0.31) for both fields. The relationship indicated that the yield of crop depends on the abovementioned soil parameters. The pH had a non-significant relation with other soil properties for Field 2 and pH had a positive correlation with N in Field 1. However, the correlations of soil parameters revealed that N, P, and K could be used to evaluate the soil’s fertility status, perceive how they affect yield, and improve both productive and unproductive areas within the field. Based on the variations in soil parameters that could be quantified, identified anddelineated the prescription maps for variable-rate fertilizer application, which will improve farm profitability and lower environmental contamination. The correlation matrix results indicated that soil properties are one of the main components that limits the yield, while many other factors also affect the yield; for example, climate change and damage from disease and insects are two clear instances that are not discussed in this research. Crop yield can also be negatively impacted by weed competition and bare spots in wheat fields, seasonal variability, winter kill, and water scarcity.

The interpolated maps were generated using Kriging interpolation techniques in Arc GIS pro 2.3 software as the literature suggested that the Kriging interpolation is more accurate than inverse distance weighting (IDW) due to the nature of soil sampling variations within the field and its application for the generation of prescription maps [48]. The interpolated maps showed a variation with significantly different values across the field. The N, P, and K maps for Field 1 (Figure 5a–c) showed almost similar variation patterns with high values in the northwest region, lower values in the southern region, and medium values in the center of the field. The M.C, pH, EC, SOM, Total N, P, and K maps indicated the large spatial variability within Field 1 (Figure 5). The semivariogram results of these soil properties showed the large spatial dependency on other soil properties (Table 3). The matrix of correlation among these soil properties also explained the similar pattern of variations within the field (Table 5). The EC and total N showed high values in the southeast region. The interpolated maps of pH revealed less variation as compared to the other soil properties. The low values of pH were observed in the western part of the field, higher values were observed in the southwest region, and medium values were observed in the center of the field. The pH map showed that most of the wheat crop had a pH ranging from 7.3 to 7.36 (Figure 5d). A similar pattern of variation was observed for 30 cm and 7 cm depths. The statistics, geostatistical analysis, correlation matrix (Figure 5, Table 1 and Table 3), and the interpolated maps of soil attributes and yield all demonstrated that there is spatial variability in the selected fields.

The soil properties and grain yield Krigged maps as shown in Figure 6 suggested substantial variability within Field 2. The spatial patterns of variation showed higher values in the north-central and north regions and moderate values in the southeastern region of the field, except N, which showed a high value on the south side. The low values were in the center of the field and some were in the south region side. A map of the soil pH shown in Figure 5d showed the substantial variability across the field, and grain yield showed its lower values on the north side. Most Field 2 sites had a pH ranging from 7.09 to 7.41.

The Krigged maps of the yield indicated substantial variation across the field (Figure 5h and Figure 6h). The low-yielding areas were on the north side, as shown in Figure 5h, and high-yielding areas were on the south side of Field 1. The high-yielding areas were in the southeast, and medium-yielding areas were in the center of Field 2. The fields inspection revealed that the low yielding was due to weeds, grass, and shortage of M.C. According to CVs, the combined findings of traditional statistics suggested that there was a moderate to high degree of variation in soil characteristics and crop yield. Except for soil pH, the geostatistical analysis also showed moderate to high spatial dependence in soil characteristics and crop production, which varied significantly within the fields. The findings of the descriptive and geostatistical analyses were consistent with the maps of soil characteristics and grain yield, indicating significant variance within the field. The significant relationships between the variables M.C., EC, SOM, Total N, N, P, and K, and grain production agreed with the relationships illustrated on the maps.

3.4. Cluster Analysis and Delineation of MZs for Wheat Crop

The interpolated maps showed large variations in soil properties and crop yield in the selected fields of wheat that requires site-specific nutrient management. The cluster analysis (Figure 7) was used to develop MZs (Figure 8) that reflect areas of very good, good, medium, low, and very low productivity. The results of clustered analysis grouped the soil properties and yield into the same class on the basis of similarity level. The dendrograms were created by performing observation cluster analysis on the field’s natural productive potentials. The soil properties and crop yield data were categorized into five classes based on their degree of similarity (Figure 6). Crop yield data were used to determine productivity levels for management zones, which included very good yield (>990 kg ha⁻¹), good yield (983 to 990 kg ha⁻¹), medium yield (971 to 982 kg ha⁻¹), low yield (951 to 970 kg ha⁻¹), and very low yield (<970 kg ha⁻¹). Crop yield had substantial relationships with soil parameters, which was supported by these developed zones. Most of the wheat crop yield fell in each productivity potential evenly, according to the grouping of the developed management zones (Figure 7), with a few data points falling in medium-productivity zones indicating average yielding areas.

The clustered observations in each group had internal homogeneity and external heterogeneity at a similarity level of more than 60%. The analysis was used to describe the difference between the developed MZs, which defined the status of soil properties in the selected fields. Soil properties and yield were higher in the very good zone, moderate levels in the medium zones, and low in the poor zones.

The wheat production might be low due to the insufficient application of fertilizer. The grasses and weeds in low-lying locations may worsen the water quality and also encouraged the weed development. These weeds limit the nutrients availability to the wheat crop, which will ultimately diminish output and raise production costs [49]. The yield, soil properties, and topographic features data were used to delineate MZs, and these developed MZs depend on spatial information that is stable or predictable over time and is related to crop yield [50].

The cluster analysis results were imported into ArcGIS software for the delineation of MZs (Figure 7). The means comparison in prescribed zones (Figure 8) for both fields of wheat based on soil properties and crop yield revealed that the highest yield and nutrient value was in the good zone and the lowest levels were in the low-productivity zones (

Table 7. Comparison of mean crop yield and soil properties for management zone of Field 1.

Crop Yield (kg ha−1) Management Zone
Soil Properties	Zone 1	Zone 2	Zone 3	Zone 4	Zone 5
	Yield	Yield	Yield	Yield	Yield
	<950	951–970	971–982	983–990	>991
	(Very Poor)	(Poor)	(Medium)	(Good)	(Very Good)
Wheat Yield (kg ha−1)	936.0 e	956.8 d	970.4 c	982.0 b	989.9 a
θv	21.8 c	23.8 ab	23.2 c	23.0 b	24.6 ab
pH	7.4 a	7.5 a	7.6 a	7.6 a	7.5 a
EC (dS m−1)	0.16 c	0.19 c	0.20 bc	0.24 b	0.3 a
SOM (%)	0.15 b	0.17 b	0.20 ab	0.21 ab	0.25 a
TN (%)	0.053 cd	0.055 c	0.056 bc	0.057 a	0.055 ab
N (mg kg−1)	15.00 cd	15.12 c	15.42 b	15.38 b	15.67 a
P (mg kg−1)	20.77 d	22.01 cd	22.99 b	23.73 a	23.78 a
K (mg kg−1)	51.89 a	52.17 d	52.53 c	52.78 b	52.96 a

The Means are represented by different letters according to the significance level of 0.05.

and

Table 8. Comparison of mean crop yield and soil properties for management zone of the Field 2.

Crop Yield (kg ha−1) Management Zone
Soil Properties	Zone 1	Zone 2	Zone 3	Zone 4	Zone 5
	Yield	Yield	Yield	Yield	Yield
	<950	951–970	971–982	983–990	>991
	(Very Poor)	(Poor)	(Medium)	(Good)	(Very Good)
Wheat Yield (kg ha−1)	758.92 d	784.61 c	801.55 b	809.26 b	833.75 a
θv	24.15 a	24.14 a	23.89 ab	23.96 ab	23.73 c
pH	6.97 b	6.99 b	7.06 a	7.00 ab	7.01 ab
EC (dS m−1)	0.14 b	0.15 a	0.14 b	0.15 a	0.15 a
SOM (%)	0.72 c	0.74 b	0.73 b	0.75 a	0.74 b
TN (%)	0.04 b	0.04 b	0.05 a	0.05 a	0.05 a
N (mg kg−1)	12.69 c	13.31 bc	12.01 d	12.95 b	13.83 a
P (mg kg−1)	17.27 cd	18.20 b	16.41 d	17.66 c	19.01 a
K (mg kg−1)	39.96 c	42.52 b	37.73 c	41.03 b	44.57 a

The Means are represented by different letters according to the significance level of 0.05.

). The rate of application will be modified as per developed management zones, and in this way, an excess amount of fertilizer will be saved. The cluster analysis, mean comparison, and interpolated maps of soil and yield revealed that there is a need of development of these prescription maps. These findings also suggested that these developed MZs could be used for the variable rate of fertilizer application.

The statistical and geostatistical results revealed the existence of variability in the selected fields and, with the help of interpolated maps, the pattern was identified. A moderate to high variation was shown with the help of all statistical results. The ANOVA and mean comparison were used to develop the zones. The MZs could be used to implement site-specific nutrient management with the help of a variable spreader, thus helping improve crop yield, reduce the cost of farms, and save the environment.

4. Conclusions

The focus of this research was to characterize and quantify the variations among the soil properties and crop yield in the selected fields of wheat that can be the basis of delineating the MZs for VR fertilizer application in arid regions. The statistical results indicated that (CV > 15%) medium to high variability exists in selected fields of wheat, which required site-specific fertilization application to improve productivity. Moreover, the analysis of geostatistical indicated (N:S < 25%) that large spatial dependency and the soil sampling grid selection should be based on the sampling interval suggested by the semivariogram range of the easily measured soil properties. Results of geostatistical showed that the 15 to 20 m grid size would provide reliable predictions suggested by the range of influence. The interpolated maps of soil properties and wheat grain yield were used to understand the fertility status in the monitoring fields. The correlation matrix results showed a significant relation (p < 0.05) between M.C, EC, SOM, N, P, K, and crop yield. The characterization and quantification of variations in soil properties would help us in the future for planning soil sampling in the fields that have soil and crop variability. Hence, the identification of factors that are responsible for a lower yield in certain areas may increase our understanding of soil properties relationships with crop yield. Based on the results, the selected fields were grouped into five different classes on the basis of internal homogeneity and external heterogeneity >60%. The development of MZs provided the pattern to group and manage the spatial variability in the selected fields of wheat, and the results showed that soil N, P, K, and EC were influential factors to develop MZs in wheat yield. The soil properties can be easily measured by a proximal sensor such as N, P, and K for VR fertilization application, and this would help to avoid the over and under-application of fertilizers to improve the productivity of crop and reduce environmental risks in wheat fields. In future studies, it would be good to link these variations among different soil and crop properties with climate change and develop the MZs according to the requirement. A web-based application can be developed to automatically delineate the MZs in the arid regions for site-specific nutrient management.