Economic and Environmental Assessment of Variable Rate Nitrogen Application in Potato by Fusion of Online Visible and Near Infrared (Vis-NIR) and Remote Sensing Data

Abstract

Addressing within-field spatial variability for nitrogen (N) management to avoid over and under-use of nitrogen is crucial for optimizing crop productivity and ensuring environmental sustainability. In this study, we investigated the economic, environmental, and agronomic benefits of variable rate nitrogen application in potato (Solanum tuberosum L.). An online visible and near-infrared (vis-NIR) spectroscopy sensor was utilized to predict soil moisture content (MC), pH, total organic carbon (TOC), extractable phosphorus (P), potassium (K), magnesium (Mg), and cation exchange capacity (CEC) using a partial least squares regression (PLSR) models. The crop’s normalized difference vegetation index (NDVI) from Sentinel-2 satellite images was incorporated into online measured soil data to derive fertility management zones (MZs) maps after homogenous raster and clustering analyses. The MZs maps were categorized into high fertile (VR-H), medium–high fertile (VR-MH), medium–low fertile (VR-ML), and low fertile (VR-L) zones. A parallel strip experiment compared variable rate nitrogen (VR-N) with uniform rate (UR) treatments, adjusting nitrogen levels based on fertility zones as 50% less for VR-H, 25% less for VR-MH, 25% more for VR-ML, and 50% more for VR-L zones compared to the UR treatment. The results showed that the VR-H zone received a 50% reduction in N fertilizer input and demonstrated a significantly higher crop yield compared to the UR treatment. This implies a potential reduction in negative environmental impact by lowering fertilizer costs while maintaining robust crop yields. In total, the VR-N treatment received an additional 1.2 Kg/ha of nitrogen input, resulting in a crop yield increase of 1.89 tons/ha. The relative gross margin for the VR-N treatment compared to the UR treatment is 374.83 EUR/ha, indicating substantial profitability for the farmer. To further optimize environmental benefits and profitability, additional research is needed to explore site-specific applications of all farm resources through precision agricultural technologies.

1. Introduction

The essential role of nitrogen (N) in ensuring a sufficient food supply for 8 billion people is undeniable. However, the overuse of N fertilizers leads to significant release of N into the environment, contaminating air and water. Globally, 40% of total nitrogen applied to agricultural lands is taken up by crops, yielding a nitrogen use efficiency (NUE) of 40% [1]. On the other hand, about half of N is released to the environment through the leaching and runoff of N compounds into waterbodies and greenhouse gas emissions [1,2]. This not only harms the environment but also poses a risk to human health. In addition, ongoing escalation in the N fertilizer use in recent years has not been matched by proportional improvements in crop yield for essential crops. This suggests that the continuous increase in nitrogen pollution is linked to the overuse of nitrogen fertilizers. Consequently, agronomists are actively seeking a strategy that achieves a ‘win–win’ outcome by improving crop production while simultaneously mitigating nitrogen pollution.

A positive plant response to nitrogen fertilizer application depends on the N threshold level, beyond which additional nitrogen input does not increase the crop yield. On the contrary, excessive N input can increase the fertilizer input cost and decrease yield due to crop lodging problems [3] in the case of growing cereal crops. Despite the threshold N limit, farmers frequently apply N fertilization more than plant biological needs with the aim of increasing crop yield potential [4]. This over-application not only increases unnecessary costs and environmental pollution but also deteriorates the soil quality by increasing soil acidification [5,6,7]. Having said that, it is suggested that advanced and sustainable nutrient management practices can help to enhance crop production and mitigate N pollution by minimizing the N loss [8].

In the context of financial constraints faced by farmers and growing environmental concerns, a range of strategies are being implemented to address the challenges associated with the overuse of N fertilizers. These strategies include monitoring of manure production, mandating adherence to nitrogen fertilizer recommendations, and performing post-harvest analysis of mineral N in soil [9,10]. However, convincing farmers to adopt the most efficient N fertilization practices proves challenging because of the crucial significance of N in affecting crop yield. Therefore, there is an urgent need to develop an innovative and effective approach to N fertilizer use for sustainable crop production.

In traditional agricultural practices, spatial variability within a field is ignored, and fields are considered as a single unit whose soil characteristics, topography, and environmental factors are spatially homogeneous [11]. Consequently, applying a uniform rate of N fertilization on the entire field area may result in over and under-application of N, potentially leading to significant environmental pollution and poor crop yield [12]. Tackling the spatial variability within the field for variable rate fertilization is a key component of precision agriculture and emerges as a sustainable nutrient management approach. In variable rate application, N fertilizers are applied based on spatial variability of soil and crop characteristics within a field [12]. Site-specific application of fertilizer has been demonstrated to increase nutrient use efficiency and net returns while also addressing the issue of excessive use of fertilizers [13]. Variable rate N (VR-N) application affects the soil system by optimizing soil nutrient distribution based on spatial variability in soil fertility and crop demand. This precision nutrient management can avoid the overapplication risk, which can lead to nutrient imbalance, soil acidification, and disruption of soil microbial communities [14,15]. Adjusting VR-N application based on soil nutrient level and crop need contributes to maintaining soil health and fertility. Successful implementation of VR-N fertilization requires high-resolution data on soil and crop characteristics to enable mapping the spatial variability at sufficient resolution. This can be achieved through the use of online proximal soil sensing and remote crop sensing technologies, coupled with chemometrics and machine learning [16]. These measurement techniques should be rapid, cost-effective, and convenient, and they require minimum labor effort. The most reported proximal technologies that fulfill these requirements are online soil sensors, which play a crucial role in the effectiveness of variable rate fertilizer application. A good example is the visible and near-infrared (vis-NIR) spectroscopy-based soil sensor [17,18]. However, the potential of these online soil sensors for VR-N fertilization was rarely explored and only for cereal crops, e.g., wheat, barley, and oil seed rape [3,19]. However, to the best of our knowledge, no work on adopting online soil sensing for VR-N fertilization in potato was reported in the literature.

Potato (Solanum tuberosum L.) is an important food and high cash crop whose yield quantity and quality are significantly affected by N fertilization. Therefore, improvement of NUE is of high priority to enhance potato tuber yield and quality. It can be hypothesized in this work that online vis-NIR soil sensing-based VR-N fertilization in potatoes increases potato yield at a reduced N application rate compared to the uniform rate N fertilization. Therefore, this work evaluates the economic, environmental, and agronomic benefits of variable rate nitrogen application in potatoes, based on management zones, derived from the fusion of multiple soil properties measured with an online vis-NIR sensor [18] with remote sensing measured crop normalized difference vegetation index (NDVI).

2. Materials and Methods

2.1. Experimental Field

The experiment was conducted in one field (Blondel-3) of an area of 6.62 ha on a commercial farm in the northwest of Belgium near the French border (51°01′47.3″ N 2°33′31.6″ E) (Figure 1). This region has a mean annual temperature of 6 to 10  °C and a mean annual precipitation ranging from 700 to 850 mm. The soil texture type of the experimental field is a sandy loam determined by the texture-by-feel method by Soil Service of Belgium (Heverlee, Belgium). The topography is flat and has almost no undulation. The experiment was carried out in the 2023 cropping season, before which the field was cultivated with barley. The overall methodology followed in this study is shown in Figure 2.

2.2. Online Soil Scanning and Sampling

Online measurement of soil spectra was carried out in August 2022 by using the online vis-NIR multi-sensor platform [18] (Figure 3). The platform consists of a subsoiler fitted to a metal frame attachable to the tractor’s three-point linkage. The subsoiler creates 15–20 cm deep trenches, whose bottom is smoothened by a downward force acting on the subsoiler chisel. An optical probe hosted by a mild steel lens holder is attached to the backside of the subsoiler heel. Soil reflectance spectra are measured while the tractor is driven at a speed of 3.5 km/h. Online spectra with a spectral range of 305–1700 nm and resolution of 1 nm are recorded by a fiber-type spectrophotometer (CompactSpec Tec5 Technology, Oberursel, Germany). A 100% ceramic disk with a 50 mm diameter was used as a white reference. Recalibration of this reference was performed at regular intervals of 30 min. The position of soil spectra was recorded by a Differential Global Positioning System (DGPS) (Trimble AG., Trimble Navigation Ltd., Sunnyvale, CA, USA). Each soil spectrum with georeferenced data was recorded at a rate of 1 Hz in a computer using MultiSpec Pro-II software from Tec5 Technology, Germany. The sensing platform was driven along parallel transects spaced 12 m apart. During online soil scanning, 13 soil samples were collected randomly from the bottom of the 12 trenches. On average, the data logger captured more than 800 online soil spectra per hectare. These samples were merged with selected samples from our spectral library to produce calibration models for several soil properties, as detailed below.

2.3. Laboratory Scanning and Analysis

Each soil sample was mixed well and cleaned manually to remove visible plant roots and stones, then ground to pass through a 2 mm sieve. Employing the quartering method proposed by Mukhopadhyay and Maiti [20], each soil sample was reduced to 150 g, which was used for the chemical analysis in the laboratory. The soil moisture content (MC) of each soil sample was determined by using the oven-drying method at 105 °C for 24 h. Subsequently, soil chemical analyses were performed at the soil laboratory of Soil Service of Belgium (HEverlee, Belgium) by using standard protocol (

Table 1. Summary statistics of laboratory-measured soil properties for the 13 soil samples collected from the target field and for the 97 samples from the soil library.

Samples	Soil Properties	Min	Max	Median	Mean	SD	Range	Skewness	Kurtosis
Blondel-3	MC	4.05	6.09	4.76	4.76	0.55	2.03	0.99	3.47
	pH	7.40	7.70	7.60	7.61	0.100	0.30	−0.96	−0.29
	TOC	0.74	1.56	0.97	1.00	0.20	0.82	1.39	1.63
	P	32.00	47.00	41.00	41.08	3.93	15	−0.64	−0.12
	K	14.00	26.00	17.00	18.23	3.49	12.00	0.85	−0.45
	Mg	43.00	52.00	47.00	47.00	2.52	9.00	0.17	−0.97
	Ca	2400.00	3500.00	3200.00	3162.00	287.34	1100.00	−1.36	1.22
	CEC	125.40	180.20	164.80	163.00	14.33	54.78	−1.31	1.10
Online spiked from spectral library	MC	4.94	24.81	14.79	14.64	6.14	19.87	−0.14	−1.35
	pH	6.10	8.00	7.70	7.55	0.37	1.90	−1.68	3.51
	TOC	0.89	2.50	1.74	1.74	0.38	1.61	−0.23	−0.50
	P	10.00	69.00	24.00	27.38	11.68	59.00	1.45	2.52
	K	8.00	30.00	14.00	16.72	7.34	22.00	0.49	−1.23
	Mg	33.00	74.00	44.00	46.85	9.96	41.00	1.21	0.77
	Ca	1670.00	6610.00	3330.00	3665.50	1398.00	4940.00	0.81	−0.50
	CEC	91.59	335.92	171.91	189.33	68.80	244.32	0.84	−0.46

Abbreviations: Moisture content, MC; total organic carbon, TOC; phosphorus, P; potassium, K; magnesium, Mg; calcium, Ca; cation exchange capacity, CEC and standard deviation, SD.

). Soil pH was measured in a 1:2.5 soil–water suspension. Soil total organic carbon (TOC) was determined by dry combustion according to the Dumas principle, followed by the removal of inorganic carbon with HCl treatment. The extractable soil potassium (K), phosphorus (P), magnesium (Mg), and calcium (Ca) were determined in an ammonium lactate extract and analyzed by using inductively coupled plasma atomic emission spectroscopy [21]. Cation exchange capacity (CEC) was calculated from base cations (K⁺, Mg²⁺, Ca²⁺, Na⁺) following the method suggested by Gupta and Gupta [22].

2.4. Online Vis-NIR Spectral Prediction Models and Maps

Vis-NIR calibration models for soil MC, pH, TOC, P, K, Mg, Ca, and CEC were developed using partial least squares regression (PLSR) after spectra preprocessing using R studio (RStudio Inc., Boston, MA, USA) with free available libraries [23]. The 13 soil samples collected from the experimental field were insufficient to develop the PLSR calibration models. Therefore, 97 soil spectra from our spectral library, collected from different fields on the same farm, were used for the development of calibration models. Adding samples from the experiment yielded a total of 110 soil samples that were used in modeling.

Spectral modeling was started with the preprocessing of spectra by using several algorithms, and the best-performing combination of algorithms was used. Before preprocessing, edge trimming was executed to eliminate the noisy boundaries within the range of 305–429 nm and 1666–1700 nm, resulting in a final spectral range of 430–1665 nm that was used for subsequent analysis. Different spectra ranges were selected for different soil properties to ensure the best possible prediction accuracy. The correction of spectral jump occurring at the junction between the visible and NIR detectors at 1045 nm was performed following the approach proposed by Mouazen et al. [24]. The parameters applied for the different preprocessing steps are shown in

Table 2. Parameters adopted for different algorithms used for the preprocessing of visible and near-infrared spectra for the different soil properties.

Order →	Moving Average	Normalization	SG Derivative			GS Derivative			Smoothing
Soil properties	w	Type	w	p	m	m	w	s	w	p	m
MC, pH, P	3	0 to 1	7	2	1	1	5	3	5	2	0
TOC	3	0 to 1	-	-	-	1	5	3	5	2	0
K	5	0 to 1	3	1	0	1	5	3	5	2	0
Mg	3	0 to 1	3	2	1	1	5	3	5	2	0
Ca, CEC	15	0 to 1	-	-	-	1	7	3	5	2	0

Abbreviations: Savitzky and Golay, SG [26]; gap segment, GS; size of the moving window, w; degree of polynomial fitting, p; t order of derivatives, m; gap size, s; moisture content, MC; total organic carbon, TOC; phosphorus, P; potassium, K; magnesium, Mg; calcium, Ca; cation exchange capacity, CEC. The right arrow (→) indicates the sequential preprocessing steps, moving from the left with the moving average towards the right with smoothing.

for moving average (MV), normalization, Savtizky and Golay (SG) derivative, gap segment (GS), and SG smoothing. The processed spectra were randomly divided into a 70% calibration set and a 30% prediction set. The PLSR analysis with cross-validation was conducted on the calibration dataset, and output models were validated using the prediction dataset. The prediction accuracy was evaluated by means of the coefficient of determination (R²), the root mean square error (RMSE), the ratio of prediction to deviation (RPD), and the ratio of performance to interquartile distance (RPIQ). Spectral analysis and PLSR modeling were performed using the “pls” package, as provided by Mevik and Wehrens [25], and accessible on R CRAN. Models providing the best accuracy using the prediction set were used for the prediction of all vis-NIR spectra collected during online measurement.

After the successful prediction of soil properties using online soil spectra, the prediction dataset containing DGPS coordinates was used to develop maps of individual soil properties. High-resolution maps necessitate fitting a semi-variogram model for each soil property, followed by interpolation using ordinary kriging (OK) available in ArcGIS (ESRI ArcGIs v10.7, Redlands, CA, USA).

2.5. NDVI Data Collection

To account for crop growth characteristics, NDVI was calculated for barley grown in the previous cropping season. It was not possible to use the NDVI of potatoes using data from the current cropping season, as N fertilizer has to be applied during potato seeding. We used the Band 4 and Band 8 of multi-band images collected on 26 March 2022 from the Sentinel-2 open data hub to calculate the NDVI index [27] using the following equation:

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\frac{(B8-B4)}{(B8+B4)}$$

(1)

2.6. Management Zones Delineation and Application Map

Crop NDVI data and online predicted soil properties were used to delineate management zones (MZ) for nitrogen application maps. Before MZs delineation, the different layers of online soil parameters and crop NDVI with different resolutions have to be resampled into the same resolution. This was performed by a raster analysis to develop a common grid of 5 m by 5 m resolution. This raster analysis was followed by K-mean cluster analysis, which was performed on the normalization data to divide the field into a limited number of clusters by using R version 4.4.3 software (RStudio Inc., Boston, MA, USA). Each zone has similar soil and crop characteristics. Since no historical yield data were available for this field, it was necessary to discuss with farmers the ranking of MZ in relation to different soil fertility levels for each management zone based on soil and crop NDVI maps. The output of this discussion allowed ranking the four MZs into high soil fertility (VR-H), medium–high fertility (VR-MH), medium–low fertility (VR-ML), and low fertility zone (VR-L).

In conventional UR farming, farmers estimate nitrogen fertilization rate based on a composite soil sample, analyzed by laboratory chemical methods to include nitrate and texture. For this study, the N fertilization rate for potatoes was 100 L/ha standard liquid fertilizer containing 39% of N, which was used in this study for the UR N fertilization treatment. A strip map for nitrogen application was developed based on fertility management zones using the fishnet function in an open-source ArcGIS software for comparison. Treatments were designed as variable rate (VR) and uniform rate (UR) parallel to strips. Each strip was 20 m wide, parallel to the tramline. Nitrogen fertilizer rate was assigned according to fertility MZs for VR application treatment. The UR treatment received the recommended N rate by the farmer, which was applied uniformly over the entire treatment. For the VR-N treatment, we followed the Robin Hood principle [3,28], which recommends more N for poor zones (feeding the poor principle) and vice versa for the fertile zones. Therefore, the VR-H zone (plot) received 50% less N fertilizer than UR, VR-MH received 25% less N, VR-ML received 25% higher N fertilizer, and VR-L received a 50% higher N fertilizer rate than the UR treatment (Figure 4).

2.7. Crop Management and Yield Measurement

This experiment followed standard crop management practices routinely adopted by the farmer for potato production. Prior to seeding, compost and chicken manure were spread across the entire field at the rate of 10 t/ha and 15 t/ha, respectively. Standard liquid N fertilizer (containing 39% N) was applied at a variable rate according to MZ recommendation after crop germination on 15 May 2023. The applied potato seed rate was 2200 kg/ha, while different pesticides were applied during the cropping season. Overall, the field received 3 irrigations by means of a gun-type sprinkler irrigation system, applying about 30 mm of water each time.

Crop yield was assessed by harvesting potato tuber manually, following the protocol suggested by Munnaf et el. [29]. Sampling points for harvesting were allocated in each polygon (treatment plot), from each of which 11 plants were harvested manually. Yields were then calculated per plot based on the area covered by 11 plants, and yield data were extrapolated to strip and treatment (3 strips). Potatoes were cleaned manually, and any rotten potatoes were removed prior to weighing.

2.8. Cost–Benefit and Environmental Analysis

Nitrogen fertilizer cost was calculated in euros (EUR) from the fertilizer’s current price and percentage of N in the fertilizer. Liquid N fertilizer with 39% N was applied. The cost of N fertilizer was 400 EUR/t. The area for each treatment was calculated from the field map. The current market price for potato yield obtained from the farmer is 200 EUR/t. The gross margin was calculated from the fertilizer price and the current market price of potato tuber for consumption potato. The cost of field scanning was not included at this stage because soil scanning cost evaluated at EUR 25 per year should be divided among several VR applications of farming input resources (e.g., manure, seeding, irrigation, tillage). Simulated profit was calculated for the whole field, assuming one treatment is applied in the entire field. The marginal value for N fertilization (MV_N) in agriculture refers to the additional benefit of applying an extra unit of N fertilizer. This concept involves assessing the impact of a small change in nitrogen application on the output of interest, typically crop yield. MV_N was calculated by the following equation:

$${\mathrm{M}\mathrm{V}}_{\mathrm{N}}=\frac{\Delta \mathrm{N}}{\Delta \mathrm{Y}}$$

(2)

where: ΔN is a change in N input, and ΔY is a change in yield.

The yield response index (YRI) measures the relative yield response to a variable rate N application (VR-N), compared to a uniform rate N (UR-N). YRI was calculated by following equation.

$$\mathrm{Y}\mathrm{R}\mathrm{I}=\frac{\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \mathrm{V}\mathrm{R}-\mathrm{N} \mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}}{\mathrm{Y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \mathrm{U}\mathrm{R}-\mathrm{N} \mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}}$$

(3)

Environmental analysis for N use was performed by calculating NUE, which quantifies how effectively N was utilized by the crop in response to crop yield. It was calculated by dividing crop yield by N fertilizer applied in the respective treatment [30]. Net N impact was determined by the change in N input in VR-N treatment with respect to the UR treatment.

2.9. Statistical Analysis

Statistical analyses were performed on laboratory-measured soil properties to identify the outlier and assess the overall quality and consistency of soil data by using R studio (RStudio Inc., USA). Significant differences among the various treatments on crop yield were assessed using analysis of variance (ANOVA), followed by Tukey’s post hoc HSD test at a significance level of p = 0.05. The relative influence of predictors (soil properties) on crop yield was calculated using the randomforest package and importance function in R studio. This analysis provides a measure of variable importance based on the mean decrease in accuracy when each predictor variable is permuted. The random forest model extracts variable importance scores using the importance function and creates a variable importance plot using the varImpPlot function.

3. Results

3.1. Accuracy of Prediction Models of Online Vis-NIR Sensor

The accuracy of PLSR models demonstrated fair to very good performance across both the calibration set and online prediction of soil properties (

Table 3. The results of partial least squares regression (PLSR) models for the prediction of soil properties, shown for cross-validation (calibration set) and validation (prediction set).

Soil Properties	Calibration (70%)				Independent Validation (30%)
	R2	RMSECV	RPD	RPIQ	R2	RMSEP	RPD	RPIQ
MC	0.94	1.32	4.18	8.62	0.87	1.82	2.86	2.69
pH	0.84	0.12	2.56	3.27	0.72	0.21	1.93	1.16
TOC	0.70	0.23	1.86	2.69	0.69	0.24	1.85	2.89
P	0.87	3.74	2.85	3.80	0.80	6.38	2.31	2.58
K	0.70	4.68	1.85	2.24	0.71	5.06	1.92	2.81
Mg	0.96	1.81	5.36	5.51	0.73	3.57	2.00	2.58
CEC	0.91	12.15	3.53	4.21	0.64	36.85	1.73	2.79

Abbreviations: Moisture content, MC; total organic carbon, TOC; phosphorus, P; potassium, K; magnesium, Mg; calcium; Ca; cation exchange capacity, CEC; coefficient of determination, R²; root mean square error of cross-validation, RMSECV; root mean square error of prediction, RMSEP; ratio of prediction to deviation, RPD; ratio of performance to inter-quartile distance, RPIQ.

). Significant variations were observed in model accuracy for different soil properties. Accuracy in the calibration was relatively higher than that of the online prediction for most of the soil properties. For instance, in the calibration set, the highest accuracy was observed for soil Mg with an R² value of 0.94, an RPD value of 5.36, and an RPIQ value of 5.51. In the online prediction set, Mg accuracy exhibited an R² value of 0.73, an RPD value of 2.00, and an RPIQ value of 2.58. Soil MC, pH, TOC, P, K, and CEC accuracies in calibration sets exhibited R² values of 0.94, 0.84, 0.70, 0.87, 0.70, and 0.91 and RPD values of 4.18, 2.56, 1.86, 2.85, 1.85, and 3.53, respectively. Moreover, in online prediction sets, MC, pH, TOC, P, K, Mg, and CEC accuracies exhibited R² values of 0.87, 0.72, 0.69, 0.80, 0.71, and 0.64 and RPD values of 2.86, 1.93, 1.85, 2.31, 1.92, and 1.73, respectively. These findings highlight the model’s efficacy in predicting all soil properties, with particularly notable performance in the calibration set. Overall, the online prediction results evaluated as RPD [19] indicated excellent performance (RPD > 2.5) for MC; very good (RPD = 2.5–2.0) for P and Mg; and fair (RPD = 1.9–1.4) for K and TOC, pH, and CEC.

3.2. Spatial Variation in Soil Fertility

Summary statistics of measured soil properties are shown in Table 1. Online soil-predicted maps and NDVI maps showed notable spatial variability, necessitating variable rate management of farming input resources using precision agricultural technologies (Figure 5). For instance, soil pH, P, and CEC showed similar patterns of fertility variation, indicating a strong correlation among these properties, which is evidenced in NDVI and yield map. While K and Mg showed a similar pattern, soil TOC demonstrated a unique spatial distributed pattern. These spatial variations emphasize the importance of employing advanced proximal and remote sensing technology to enhance soil management practices.

Four MZ classes were determined by K-means clustering and categorized into different fertility levels, namely, H, MH, ML, and L (Figure 6). The quality of the MZ seems convincing, as it mirrors the spatial variability pattern of a few soil properties and NDVI (Figure 6). The classification was based on soil fertility status and the experience of the farmer in this field, so it was not a straightforward decision to make. According to fertility classification, the middle part of the field has MH and H fertility levels, the southern part of the field was of L fertility (low P, NDVI in particular), and the northern part was classified as ML fertile zone. Nitrogen application rates were estimated based on these fertility MZs.

3.3. Crop Yield Response to Variable Nitrogen Application

Crop yield response was different among the different VR-N application rates (Figure 7) within each MZ. In contrast to the UR-N treatment, the VR-H treatment, which received 50% less N than the UR treatment, exhibited an 8.1% increase in potato yield. Additionally, the VR-ML and L treatments, which received 25% and 50% higher N than the UR treatment, demonstrated increases in potato yield by 14.3% and 2%, respectively. The only zone that had a smaller yield than the control UR was the VR-MH. The impact of predictors on crop yield was assessed using a random forest model. The analysis revealed TOC, P, and CEC as the predominant factors influencing crop yield, with relative influences of 18.02%, 17.33%, and 16.81%, respectively. Meanwhile, Mg, K, and pH exhibited relative influences of 16.53%, 16.05%, and 15.22%, respectively (Figure 8).

3.4. Cost–Benefit and Environmental Assessment of Variable Rate Nitrogen

Results of the cost–benefit analysis demonstrated that VR-N application significantly affected yield throughout the different MZs in the field (

Table 4. Calculation of cost–benefit analysis for the variable rate nitrogen (N) fertilization compared to the uniform rate N fertilization.

Treatment	Area (ha)	N-Fertilizer Application Rate (N kg/ha)	Fertilizer Cost (EUR/ha)	Yield (ton/ha)	Revenue (EUR/ha)	Gross Margin (EUR/ha)	Relative Gross Margin (EUR/ha)	Simulated Field Relative Gross Margin (EUR)	MVN (t/kg)	YRI
UR	2.00	39.00	39.78	58.66	11,712.91	11,693.13
VR-H	0.79	19.50		63.44						1.08
VR-MH	1.39	29.25		55.44						0.95
VR-ML	1.30	48.75		67.06						1.14
VR-L	1.14	58.50		57.32						0.98
Total VR	4.63	40.24	41.08	60.55	12,109.05	12,067.96	374.83	2481.37	0.68	1.03

Abbreviations: Uniform rate, UR; variable rate high fertile zone, VR-H; medium–high fertile, VR-MH; medium–low fertile, VR-ML; low fertile zone, VR-L; marginal value of nitrogen, MV_N; yield response index, YRI.

). The yield of VR-H and VR-ML MZs exceeded that of the UR treatment, whereas the yield for the VR-MH and VR-L MZs was smaller. Overall, the VR-N treatment increased the average yield (60.55 t ha⁻¹), compared to the UR-N treatment (58.66 t ha⁻¹), which resulted in a corresponding increase in overall gross margin. Among different variable rate treatments, the highest yield (67.06 t ha⁻¹) obtained by the VR-ML MZ treatment followed by VR-H (63.44 t ha⁻¹) has compensated for the reduced yield in the VR-MH, and VR-L MZs, resulting in an overall higher relative gross margin of the VR-N treatment of 374.83 (EUR/ha), compared to that produced by the UR treatment. The simulated relative gross margin of the VR-N treatment for the field was EUR 2481.37. The highest yield response index (YRI—the ratio of yield of VR-N/yield of UR treatments) was in VR-ML MZ (1.14), followed by VR-H MZ (1.08), VR-L (0.98) and VR-MH (0.95), respectively. Overall, YRI for the VR-N treatment was 1.03, indicating the positive yield response to the VR treatment application compared to the UR application. Marginal more N of 0.68 t/kg was consumed in the VR-N treatment, which has resulted in overall potato yield, indicating improved N use efficiency per unit of N fertilizer applied.

Environmental assessment of N use in VR application was performed by determining N use efficiency (NUE) and net N impact (ΔN) (

Table 5. Environmental assessment of nitrogen (N) applied by the variable rate N (VR-N) fertilization, compared to the uniform rate N fertilization in potato (S. tuberosum L.).

Treatments	Area (ha)	N-Fertilizer Application Rate (N kg/ha)	Yield (t/ha)	NUE (%)	Net N Impact ΔN (kg/ha)
UR	2.00	39.00	58.66	58.70
VR-H	0.79	19.50	63.44	126.87	19.50
VR-MH	1.39	29.25	55.44	73.92	9.75
VR-ML	1.30	48.75	67.06	53.64	−9.75
VR-L	1.14	58.50	57.32	38.21	−19.5
Total VR	4.63	40.24	60.55	60.81	−1.28

Abbreviations: Uniform rate, UR; variable rate high fertile zone, VR-H; medium–high fertile, VR-MH; medium–low fertile, VR-ML; low fertile zone, VR-L; nitrogen use efficiency, NUE.

). The results showed that NUE was observed highest in VR-H (126.87%), followed by VR-MH (73.92%), UR (58.70%), VR-ML (53.64%), and VR-L 38.2%, respectively. On average, the overall NUE for the VR-N treatment was relatively higher (60.81%) than that of the UR treatment (58.70%), confirming the more efficient use of N by the VR-N treatment compared to the UR treatment. Net N impact (ΔN), indicating the difference in N input in the VR-N treatment with respect to the UR N input, demonstrated a negative value of −1.28 (Table 5). This result explains that the VR-N treatment consumed a larger N input compared to the UR treatment. The values of ΔN were negative in poor fertility MZs (VR-ML and VR-L), indicating a higher N input compared to UR treatment, whereas positive values can be observed in the VR-H and VR-MH zones.

4. Discussion

4.1. Models Accuracy

Despite encountering low moisture levels (4.05–6.09%) that resulted in a higher noise generated by vibration during online soil scanning, the vis-NIR models demonstrated a satisfactory performance ranging from fair to excellent prediction of various soil properties (Table 3). The model’s performance in the cross-validation was slightly higher than that in prediction for most of the soil properties, findings which are comparable with previous studies’ results [31,32]. The model’s performance closely resembled that documented in previous studies [32,33]. In this study, the accuracy obtained for the majority of vis-NIR models met the criteria established by Maleki et al. [34], requiring an R² > 0.70 to be considered suitable for the implementation of site-specific applications. Furthermore, current models showed higher prediction accuracies than those models utilized in other variable rate applications employing the same online vis-NIR sensing system. For instance, Zhang et al. [35] and Guerrero et al. [3] successfully implemented VR manure and N applications, respectively, using the same online vis-NIR models, reporting R² values ranging from 0.51 to 0.81. Similarly, Mouazen and Kuang [36] reported a successful site-specific P management using the same online vis-NIR sensor, using (in advance developed) prediction model for P with R² = 0.60, compared to the P model of this work with an R² = 0.8 (Table 3). Thus, it can be inferred that current vis-NIR models exhibit sufficient accuracy to determine MZs for VR-N application.

4.2. Within Field Spatial Variability

Considerable spatial variability was identified in the online-predicted soil properties maps (Figure 7). Zones in the field characterized by low pH exhibited similarities with reduced P content, affirming the impact of pH on P. This correlation may be attributed to lower Ca levels in this specific field area. Additionally, the areas with low pH demonstrated elevated Mg and K content, indicating a negative correlation between pH and the levels of Mg and K. The influence of soil pH is crucial for P availability, as P accessibility can be limited in both acidic and alkaline soils [37]. In the current field, pH is in the alkaline range, and this should be corrected for a better P uptake by the crop. The distinct spatial pattern observed in NDVI closely resembled the distribution of P, pH, and CEC. Consequently, special consideration is directed in this work and perhaps in future work toward P, CEC, and NDVI maps in making informed crop management decisions. Furthermore, the NDVI maps confirmed the significance of P, CEC, and pH as pivotal factors limiting yield, as they exhibited similar spatial patterns.

In previous studies, a lack of similarity between crop yield patterns and NDVI was observed [19,38], highlighting the potential unreliability of relying solely on NDVI for ranking fertility classes in soil management decisions. Since other studies reported a strong correlation between crop yield and NDVI [39], we opted to use NDVI as an alternative to yield (as historical yield data were not available), along with other soil properties, in classifying MZ into different fertility classes (Figure 5). The determination of N application rates relied on the fertility level of each zone and the farmer’s experience with the current field. Similar practices were employed in previous studies to optimize manure and P management, resulting in successful variable rate applications [35,36].

4.3. Variable Rate Nitrogen Impact on Crop Yield

Applying the right amount of N fertilizer at the right time in the right place is critical for the effective management of N to meet the N needs of plants. Variable rate N application significantly affected crop yield (Figure 5). Nitrogen is the most limiting nutrient in crop production, and in most plants, it is higher in concentration than all other mineral nutrients [40]. Nitrogen is applied in greater amounts than other nutrients, which indicates its critical role in potato production [41,42]. According to Moller et al. [43], N availability poses a significant limitation to yield in potato cultivation. To optimize the potential of a crop in a given location, it is crucial to apply the appropriate and optimal amount of N, ensuring vitality and maximizing crop performance. In the VR-ML fertility zone, using a 25% higher N rate led to a higher crop yield compared to the standard rate (UR). Similarly, in the VR-H zone, a 50% lower N application than the standard rate still resulted in a relatively higher crop yield. These results align with findings by Munnaf et al. [19], who reported that a high-fertility zone exhibited increased crop yield despite lower N input. This could be attributed to the possibility that there was a high natural availability of N through mineralization. Even with a lower external N application, the soil’s inherent fertility was sufficient to support a relatively higher crop yield. But, to confirm if this was the case for this output, it would have been necessary to measure the N mineralization rate and available N throughout the cropping sensor. However, this was not performed in the current work, as no sensor technology is available to collect the necessary data on nitrate and mineralization rates. Conversely, in the VR-ML zone, N was a limiting factor, hence, increasing N application by 25% improved crop yield compared to the standard rate. However, in the VR-L zone, even applying a 50% higher N input than the UR treatment did not improve crop yield (Table 5). Indeed, the extra N rate used in the VR-L zone was not utilized by the crop. This can result in a risk of increasing residual N in the soil without a corresponding increase in crop yield [44]. Therefore, future work should indeed consider monitoring nitrate and measurement of mineralization rates for the calculation of the recommendation of VR-N application rather than relying on the % increase or decrease in N by the Robin Hood method adopted in this work. Increasing N levels in the plant tissue can positively affect crop yield by improving the absorption of solar radiation through enhanced tillering, leaf enlargement, and enhanced photosynthesis [45]. However, it is essential to note that excessive N application can have negative consequences, causing problems like lodging in cereal crops [46] and foliar infections [45], which can ultimately reduce overall yield. For potatoes, the latter issue might be the cause of reducing tuber yields when applying N fertilizer beyond the crop response. Moreover, soil properties such as soil TOC, pH, and CEC significantly influence the N availability and crop yield. It is well known that TOC serves as a source of N through microbial decomposition. The breakdown of TOC releases N in various forms, including ammonium (NH₄⁺) and organic nitrogen compounds [47]. CEC crucially affects N availability by enhancing the retention of ammonium in high CEC soils [48]. In the low-fertility zone (VR-L) shown in Figure 6, with the lowest pH, CEC, and P, the limited crop yield response to increased N applied (+50% of that of UR) can be attributed to nutrient imbalances and interactions [49]. Phosphorus deficiency constrains the plant’s ability to utilize nitrogen effectively due to its impact on plant root development and nitrogen metabolism regulation. Moreover, the random forest model confirms that TOC, P, and CEC were major influencing factors on crop yield (Figure 8). Thus, the study underscores the importance of VR-N application in precision agriculture, emphasizing the need to balance inputs across varied soil fertility zones and adopt a holistic soil management approach for sustainable crop production.

4.4. Cost–Benefit and Environmental Analysis

Variable rate applications, based on recommendations developed from high-resolution soil data, could be economically and environmentally beneficial. An online multi-sensor platform for the collection of high-resolution soil data offers small- to medium-sized farmers the opportunity to access scanning services without the requirement of investing in this specialized equipment. They can acquire scanning services at a rate of approximately 25 EUR per hectare per year. This includes online scans, maps illustrating predicted soil properties, and recommendation maps for VR application. Implementation of VR-N application based on MZ delineation through the fusion of high-resolution soil and remote sensing data can increase profitability by increasing crop yield and reducing excessive N use [12,19]. Variable rate N input significantly enhanced crop yield compared to the UR treatment, leading to increased gross margins (Table 4). These findings align with the results reported by Basso et al. [50]. In the VR-H zone, despite receiving 50% less N than the UR treatment, there was a higher crop yield response. Conversely, in the VR-ML zone, receiving 25% higher N input than the UR treatment, the highest yield (67.76 t/ha) was achieved. However, in the VR-L zone, where the crop did not respond to a 50% increase in N compared to the UR rate, the excessive N use in this class exceeded the maximum crop response threshold to N fertilizer [51]. Therefore, it is not surprising to observe almost equal or smaller potato yield compared to the UR treatment, which is most likely attributed to increased crop sensitivity to foliar infection due to biotic stresses. This suggests that a 50% increase in N than the recommended N fertilizer should be avoided. Consequently, a recommendation with a 25% change in N input rate proved effective and would not exceed N over the crop N response limit. Overall, the VR N application increased the gross margin by EUR 373.84 per hectare, with a marginal N increase of 0.68 (t/kg), signifying an increase in crop yield with each additional kilogram of N input. This study demonstrates that VR-N fertilization can offer additional profits to farmers compared to UR-N applications.

Compared to the UR treatment, N use per hectare increased by 3.25% with the VR-N treatment; however, N remains below the environmental threshold level while contributing to increased crop yield and overall profitability (Table 4). Environmental analysis of N use revealed a significant increase in NUE, compared to the UR-N treatment, attributed to a better crop yield response to site-specific N application. It is reported that precision nitrogen management is an environmentally friendly approach that can help mitigate soil nitrogen pollution and greenhouse gas (GHG) emissions while maintaining soil health, preserving fertility, and safeguarding the environment by preventing the overuse of fertilizers to sustain crop yields [15]. Precision N application based on MZs is shown to be both environmentally and economically profitable [3,52]. However, recommending a 25% change in the N input rate instead of a 50% change could potentially yield greater profits and environmental benefits compared to the UR treatment.

5. Conclusions

Implementation of VR-N fertilization in a potato field based on MZ delineation using high-resolution data on soil fertility parameters and crop NDVI measured by means of an online vis-NIR spectroscopy sensor and Sentinel-2 satellite imagery, respectively, has increased crop yield and profitability indicated as the gross margin. The highest fertile zone (VR-H) receiving 50% lower N input and the medium-low fertile (VR-ML) zone receiving 25% more N rate than the UR treatment noticeably increased the crop yield, hence increasing the gross margin of the VR-N treatment over the UR treatment. This finding suggests a possibility of mitigating the negative environmental effects by reducing fertilizer expenses while sustaining the crop yield.

A marginal N value of 0.68 (t/kg) applied in the VR-N treatment has resulted in increased yield and gross margin. This suggests that the crop has used this minor extra N to produce more yield, minimizing the chance for excess N leached to the environment. Moreover, a simulated gross margin of EUR 2481 was observed from the VR-N fertilization for the whole field. These results need to be validated further for potatoes by running similar work in several fields to confirm the potential profitability, with a minor increase in applied N that has no negative effect on the environment as long as the extra N applied is consumed by the crop to produce more yield compared to the traditional UR fertilization. Furthermore, historical yield data should be incorporated into the decision-making process to improve VR-N recommendations. Similarly, data on nitrate and nitrogen mineralization rate per MZ should also be monitored and used in the decision-making to overcome the Robin Hood ± % change in N rate adopted in this work.