Comparative Analysis of Energy and GHG Emissions Using Fixed and Variable Fertilization Rates

: Reducing the overuse of mineral fertilizers in crop production is a key factor related to healthier soil, healthier food, and more economical, efﬁcient, and cleaner agricultural production. The aim of this study was to investigate and to compare the effects of ﬁxed and variable rate fertilization with fertilizer consumption, energy consumption, and environmental impact. A 4-year experimental ﬁeld study using crop rotation of four plants (spring barley, winter oilseed rape, winter wheat, and faba bean) was performed. Fertilization with phosphorus and potassium at a variable rate were performed and applied based on the soil properties analyzed before the start of the research and completion of the fertilization maps. Nitrogen fertilization was performed by additional fertilization using a proximal N-sensor, which gives the accurate nitrogen uptake in plants in real-time. This was followed by a comparative evaluation analysis of variables and conventional ﬁxed-rate fertilization methods to assess fertilizer consumption, energy consumption, environmental impact, and economic efﬁciency. The results of the study showed that an application of a variable fertilization rate can reduce the total amount of nitrogen, phosphorus, and potassium fertilizers by 24.9%, energy consumption by 3463.1 kg ha − 1 , and emissions of greenhouse gases (GHG) by 341.5 kg CO 2eq ha − 1 compared to ﬁxed fertilization rate. The method of fertilization with a variable application rate reduced the costs of and at the same time, the costs of the total plant production by € 168.0 ha − 1 , on average.


Introduction
Plant nutrition is one of the most important factors on which the yield and the quality of agricultural products depend [1]. Therefore, organic and mineral fertilizers are used to enrich the soil with nutrients and to achieve the highest soil fertility. Over the last few decades, fertilizer use has been growing exponentially around the world. In 2018, global consumption of chemical and mineral fertilizers amounted to 188 million tons, of which nitrogen (N) averaged 58%, phosphorus (P 2 O 5 ) 22%, and potassium (K 2 O) 21%. In a period of 18 years (2000-2018), fertilizers' consumption increased by about 40% [2]. In Europe, fertilizers per crop area increased by 10% during this period and reached 77 kg ha −1 , N-49, P 2 O 5 -13, and K 2 O-15 kg ha −1 , respectively [2]. Intensive use of mineral fertilizers and application of conventional fertilization methods have negative effects on the soil, the environment, and human health. For example, ammonia emission from over-fertilized lands can damage organisms living in lakes and reservoirs, and the amount of nitrate concentration in drinking water can increase as well. Excessive amounts of fertilizers are related to higher GHG emissions [1,3]. Traditionally, fertilizers are applied to all arable land, regardless of soil variation in the field area [4]. Variable rate fertilization (VRF) is essential for the implementation of precision farming and for the efficient use of mineral fertilizers and nutrient management in plants adapted to the conditions of individual field sites [5]. It should be taken into consideration that the application of VRF usually used in conventional farming practice and variable rate, applied after analyzing the soil and assessing the changes in fertilizer consumption, energy consumption, GHG emissions, and economic costs related to fertilizers. The soil was described as sandy to silty loam according to the soil granulometric classification [27]. Experimental studies were performed comparing 2 different fertilization methods. The fixed fertilization rate (FRF) was determined using the recommended fertilizer rates based on the assumed yield for each crop. The variable fertilization rate (VRF) was determined based on soil characteristics in site-specific areas, crop assumed yield, and nutrient requirements for each crop.  Field boundaries were measured on 16 September 2015, with an off-road vehicle Toyota Hilux (Toyota Motor Corporation, Toyota, Japan). The total field size was determined to be 11.37 ha by using the GPS coordinate positioning equipment Trimble EZ-Guide 250 with GPS antenna (Trimble Navigation Ltd., Alpharetta, USA). John Deere 6630 tractor (John Deere GmbH & Co., Mannheim, Germany), spreader Rauch Axis-H (Rauch Landmaschinenfabrik GmbH, Sinzheim, Germany), Rauch spreader terminal CCI 100 (Müller GmbH, Rheinfelden, Germany) were used to spread the fertilizer according to variable rate fertilization maps. The crop was harvested with a John Deere S685 harvester (John Deere GmbH & Co.KG, Zweibrücken, Germany) with a working width of 9.14 m.

Measurements of Soil Electrical Conductivity, Field Subdivision, and Soil Sampling
Prior to the start of the experiment, soil apparent electrical conductivity (ECa) was measured using an EM-38 MK-2 (Geonics Ltd., Mississauga, Canada) at a depth of 0-1.5 m in the vertical position to determine differences in the field soil's properties and to establish an accurate soil sampling plan. EM-38 readings were applied while driving a Toyota off-road vehicle towing EM38-MK2 mounted on plastic sledges at a speed of 10-15 km h −1 in parallel technological tracks, separated by 30 m, with a total of 1046 readings (92 readings ha −1 ) for the field of 11.37 ha. These readings were stored on a computer Panasonic Toughbook CF 19 (Panasonic Corporation, Osaka, Japan).
ECa readings were converted to a CSV file format using Convert EM38-MK2 software, and subsequently, an ECa map ( Figure 3) was created using the Open Source Geographic Information System (QGIS) program. A total of 4 areas (each approximately 3 ha) with different ranges of ECa were delineated ( Figure 1). The delineated plots in a digital shape (SHP) file was sent back to the field computer, and soil samples were taken. Field boundaries were measured on 16 September 2015, with an off-road vehicle Toyota Hilux (Toyota Motor Corporation, Toyota, Japan). The total field size was determined to be 11.37 ha by using the GPS coordinate positioning equipment Trimble EZ-Guide 250 with GPS antenna (Trimble Navigation Ltd., Alpharetta, GA, USA). John Deere 6630 tractor (John Deere GmbH & Co., Mannheim, Germany), spreader Rauch Axis-H (Rauch Landmaschinenfabrik GmbH, Sinzheim, Germany), Rauch spreader terminal CCI 100 (Müller GmbH, Rheinfelden, Germany) were used to spread the fertilizer according to variable rate fertilization maps. The crop was harvested with a John Deere S685 harvester (John Deere GmbH & Co.KG, Zweibrücken, Germany) with a working width of 9.14 m.

Measurements of Soil Electrical Conductivity, Field Subdivision, and Soil Sampling
Prior to the start of the experiment, soil apparent electrical conductivity (ECa) was measured using an EM-38 MK-2 (Geonics Ltd., Mississauga, Canada) at a depth of 0-1.5 m in the vertical position to determine differences in the field soil's properties and to establish an accurate soil sampling plan. EM-38 readings were applied while driving a Toyota off-road vehicle towing EM38-MK2 mounted on plastic sledges at a speed of 10-15 km h −1 in parallel technological tracks, separated by 30 m, with a total of 1046 readings (92 readings ha −1 ) for the field of 11.37 ha. These readings were stored on a computer Panasonic Toughbook CF 19 (Panasonic Corporation, Osaka, Japan).
ECa readings were converted to a CSV file format using Convert EM38-MK2 software, and subsequently, an ECa map ( Figure 3) was created using the Open Source Geographic Information System (QGIS) program. A total of 4 areas (each approximately 3 ha) with different ranges of ECa were delineated ( Figure 1). The delineated plots in a digital shape (SHP) file was sent back to the field computer, and soil samples were taken.  Soil (at a depth of 20-30 cm) was sampled on 17 September 2015 and 2 September 2019 using an automatic soil sampling equipment ( Figure 4) developed by Agricon (GmbH, Jahna, Germany) and manufactured by Adigo AS, Langhus, Norway). The subsamples were collected while driving non-stop (50% more efficient than sampling at a stop [28]) in the trajectory of the letter Z and combined in one sample in each of the 4 homogeneous areas delineated previously. Four samples were taken from 4 areas formed according to soil ECa differences. One sample consisted of 15-20 subsamples (Figure 1). These samples were collected into boxes containing 300-500 g of soil. After all the boxes were filled, soil samples from each box were placed into separate plastic bags, and labeled with a bar code including information, such as the coordinates, date, and time of sampling. Samples were sent to an accredited laboratory Agrolab GmbH (Leinefelde-Worbis, Germany), where their granulometric Soil (at a depth of 20-30 cm) was sampled on 17 September 2015 and 2 September 2019 using an automatic soil sampling equipment ( Figure 4) developed by Agricon (GmbH, Jahna, Germany) and manufactured by Adigo AS, Langhus, Norway). The subsamples were collected while driving non-stop (50% more efficient than sampling at a stop [28]) in the trajectory of the letter Z and combined in one sample in each of the 4 homogeneous areas delineated previously.  Soil (at a depth of 20-30 cm) was sampled on 17 September 2015 and 2 September 2019 using an automatic soil sampling equipment ( Figure 4) developed by Agricon (GmbH, Jahna, Germany) and manufactured by Adigo AS, Langhus, Norway). The subsamples were collected while driving non-stop (50% more efficient than sampling at a stop [28]) in the trajectory of the letter Z and combined in one sample in each of the 4 homogeneous areas delineated previously. Four samples were taken from 4 areas formed according to soil ECa differences. One sample consisted of 15-20 subsamples (Figure 1). These samples were collected into boxes containing 300-500 g of soil. After all the boxes were filled, soil samples from each box were placed into separate plastic bags, and labeled with a bar code including information, such as the coordinates, date, and time of sampling. Samples were sent to an accredited laboratory Agrolab GmbH (Leinefelde-Worbis, Germany), where their granulometric Four samples were taken from 4 areas formed according to soil ECa differences. One sample consisted of 15-20 subsamples (Figure 1). These samples were collected into boxes containing 300-500 g of soil. After all the boxes were filled, soil samples from each box were placed into separate plastic bags, and labeled with a bar code including information, such as the coordinates, date, and time of sampling. Samples were sent to an accredited laboratory Agrolab GmbH (Leinefelde-Worbis, Germany), where their granulometric structure, magnesium (Mg) and pH (CaCl 2 method), potassium (K) and phosphorous (P) (CAL method) were determined.

Determination of Fixed and Variable Fertilization Rates for Phosphorus and Potassium
In order to apply the fertilizer at a variable rate, fertilization plans were drawn up ( Figure 5). Results of the soil analysis were imported into an agronomic software Agriport. Then, variable fertilizer application rates for a 4-year crop rotation were calculated using the data from soil analysis, the uptake of nutrients P and K for each crop and its assumed yield (t ha −1 ), respectively: Spring barley (7 t ha −1 ) P 2 O 5 -33 kg ha −1 , K 2 O-25.7 kg ha −1 ; winter rape (5 t ha −1 ) P 2 O 5 -33 kg ha −1 , K 2 O-25.7 kg ha −1 ; winter wheat (8 t ha −1 ) P 2 O 5 -33 kg ha −1 , K 2 O-25.7 kg ha −1 ; faba beans (4 t ha −1 ) P 2 O 5 -33 kg ha −1 , K 2 O-25.7 kg ha −1 . The maps for P and K fertilization for the computer of Rauch spreader were created in a SHP file also using Agriport online program. In order to generate a variable rate fertilization (VRF) map for each element, soil analysis was loaded into the program. Then, crop rotation and assumed yield of each crop were entered into the program, which automatically generated VRF maps for each fertilizer. Fertilizers P and K were spread at a variable rate according to the fertilization map, distributing 25% of the planned 4-year fertilization rate in each year. Magnesium was enough to satisfy crop demand in the whole rotation and was not applied. structure, magnesium (Mg) and pH (CaCl2 method), potassium (K) and phosphorous (P) (CAL method) were determined.

Determination of Fixed and Variable Fertilization Rates for Phosphorus and Potassium
In order to apply the fertilizer at a variable rate, fertilization plans were drawn up ( Figure 5). Results of the soil analysis were imported into an agronomic software Agriport. Then, variable fertilizer application rates for a 4-year crop rotation were calculated using the data from soil analysis, the uptake of nutrients P and K for each crop and its assumed yield (t ha −1 ), respectively: Spring barley (7 t ha −1 ) P2O5-33 kg ha −1 , K2O-25.7 kg ha −1 ; winter rape (5 t ha −1 ) P2O5-33 kg ha −1 , K2O-25.7 kg ha −1 ; winter wheat (8 t ha −1 ) P2O5-33 kg ha −1 , K2O-25.7 kg ha −1 ; faba beans (4 t ha −1 ) P2O5-33 kg ha −1 , K2O-25.7 kg ha −1 . The maps for P and K fertilization for the computer of Rauch spreader were created in a SHP file also using Agriport online program. In order to generate a variable rate fertilization (VRF) map for each element, soil analysis was loaded into the program. Then, crop rotation and assumed yield of each crop were entered into the program, which automatically generated VRF maps for each fertilizer. Fertilizers P and K were spread at a variable rate according to the fertilization map, distributing 25% of the planned 4-year fertilization rate in each year. Magnesium was enough to satisfy crop demand in the whole rotation and was not applied. Fixed rates of fertilization (FRF) for nutrients P and K were calculated using the data from the agronomic software Agriport only for nutrient uptake for each crop and assumed yield, respectively: For spring barley (7 t ha −1 ) P2O5-56 kg ha −1 , K2O-42 kg ha −1 ; for winter rape (5 t ha −1 ) P2O5-89.3 kg ha −1 , K2O-49.8 kg ha −1 ; for winter wheat (8 t ha −1 ) P2O5-64 kg ha −1 , K2O-48 kg ha −1 ; for faba beans (4 t ha −1 ) P2O5-47.6 kg ha −1 , K2O-56 kg ha −1 .

Yara N-Sensor ALS Calibration
Before starting to fertilize winter wheat in the field, a manual proximal sensor Yara N-tester (manufactured by Konica Minolta INC, Japan with Yara software) was used to calibrate the proximal Yara N-Sensor ALS (manufactured by Tec5 GmbH, Germany and developed by Yara International ASA, Oslo, Norway). Calibration was needed to set the Fixed rates of fertilization (FRF) for nutrients P and K were calculated using the data from the agronomic software Agriport only for nutrient uptake for each crop and assumed yield, respectively: For spring barley (7 t ha −1 ) P 2 O 5 -56 kg ha −1 , K 2 O-42 kg ha −1 ; for winter rape (5 t ha −1 ) P 2 O 5 -89.3 kg ha −1 , K 2 O-49.8 kg ha −1 ; for winter wheat (8 t ha −1 ) P 2 O 5 -64 kg ha −1 , K 2 O-48 kg ha −1 ; for faba beans (4 t ha −1 ) P 2 O 5 -47.6 kg ha −1 , K 2 O-56 kg ha −1 .

Yara N-Sensor ALS Calibration
Before starting to fertilize winter wheat in the field, a manual proximal sensor Yara N-tester (manufactured by Konica Minolta INC, Japan with Yara software) was used to calibrate the proximal Yara N-Sensor ALS (manufactured by Tec5 GmbH, Germany and developed by Yara International ASA, Oslo, Norway). Calibration was needed to set the N application rate for the Skagen variety of winter wheat at certain growth stages (BBCH 30; 37) in a repeatable part of the field. Yara N-Sensor ALS can collect data only to determine N uptake in the crop and cannot recommend N application rate. An Ntester helps to determine N optimum application rate for different wheat varieties. The calibration location for the device was selected after a visual assessment of the winter wheat field after entering the 4th technological track from the field edge. Measurements were performed by randomly selecting the middle part of the fully developed last leaves of 30 different plants in the 20 m technological track section ( Figure 6). The area measured with the N-tester was scanned with Yara N-Sensor ALS device to complete the calibration process. Then, the N application rate for the winter wheat was entered into the terminal using the N application program in the tractor and this software adjusted the N rate in site-specific areas in real-time while driving in the field. N application rate for the Skagen variety of winter wheat at certain growth stages (BBCH 30; 37) in a repeatable part of the field. Yara N-Sensor ALS can collect data only to determine N uptake in the crop and cannot recommend N application rate. An N-tester helps to determine N optimum application rate for different wheat varieties. The calibration location for the device was selected after a visual assessment of the winter wheat field after entering the 4th technological track from the field edge. Measurements were performed by randomly selecting the middle part of the fully developed last leaves of 30 different plants in the 20 m technological track section ( Figure 6). The area measured with the N-tester was scanned with Yara N-Sensor ALS device to complete the calibration process. Then, the N application rate for the winter wheat was entered into the terminal using the N application program in the tractor and this software adjusted the N rate in sitespecific areas in real-time while driving in the field.

Fixed and Variable Rate Application for Nitrogen Fertilization
The FRF method was not performed in the field as an actual comparison. FRF for N nutrient was calculated based on conventional farming experience for spring barley N-50 kg ha −1 , winter oilseed rape N-188 kg ha −1 , winter wheat N-210 kg ha −1 , and faba beans N-70.2 kg ha −1 , respectively. For the first N application for winter wheat N-40 kg ha −1 (BBCH 23), and for winter oilseed rape N-60 kg ha −1 (BBCH 18), fixed rates were calculated based on conventional farming experience also when using VRF method. For the first N application, it was not possible to use the Yara N-Sensor ALS since it could not collect data from the poor biomass of the crops. For spring barley and faba beans, fixed N nutrient rates were calculated according to conventional farming experience and applied only once before seeding using a fixed rate.
Yara N-Sensor ALS agronomic software calculated the actual N-uptake of the crop. Using N application software optimum rates for VRF in site-specific areas were derived from the N-uptake data in the crop. Calculated rates were sent to the spreader controller, which adjusted the shutters of the implement. At that time, the N-sensor was installed on the roof of the tractor, scanning the crop in 4 m stripes on both sides when driving at 30 m width tramlines.
Information about the fertilizers applied and the amount of N absorbed by the plants in individual parts of the field was transmitted online from the computer inside the tractor to the Agriport program into the office computer. Fertilizer rates for each part of the field and an average could be seen in this program. In our case, N fertilizers YaraBela AXAN

Fixed and Variable Rate Application for Nitrogen Fertilization
The FRF method was not performed in the field as an actual comparison. FRF for N nutrient was calculated based on conventional farming experience for spring barley N-50 kg ha −1 , winter oilseed rape N-188 kg ha −1 , winter wheat N-210 kg ha −1 , and faba beans N-70.2 kg ha −1 , respectively. For the first N application for winter wheat N-40 kg ha −1 (BBCH 23), and for winter oilseed rape N-60 kg ha −1 (BBCH 18), fixed rates were calculated based on conventional farming experience also when using VRF method. For the first N application, it was not possible to use the Yara N-Sensor ALS since it could not collect data from the poor biomass of the crops. For spring barley and faba beans, fixed N nutrient rates were calculated according to conventional farming experience and applied only once before seeding using a fixed rate.
Yara N-Sensor ALS agronomic software calculated the actual N-uptake of the crop. Using N application software optimum rates for VRF in site-specific areas were derived from the N-uptake data in the crop. Calculated rates were sent to the spreader controller, which adjusted the shutters of the implement. At that time, the N-sensor was installed on the roof of the tractor, scanning the crop in 4 m stripes on both sides when driving at 30 m width tramlines.
Information about the fertilizers applied and the amount of N absorbed by the plants in individual parts of the field was transmitted online from the computer inside the tractor to the Agriport program into the office computer. Fertilizer rates for each part of the field and an average could be seen in this program. In our case, N fertilizers YaraBela

Assessment of Energy Use and GHG Emissions
Fertilizer energy consumption (MJ ha −1 ) using variable rate fertilization (VRF) and fixed rate fertilization (FRF) was calculated by multiplying the amount of fertilizer used by the energy equivalent the N fertilizer (60.6 MJ kg −1 ) [29]. Energy consumption equivalents for P 2 O 5 and K 2 O are presented in Table 1. In order to compare fertilization methods, the difference in energy consumption between VRF and FRF was calculated. Another very important benchmark that has been identified was the reduction of fertilizers energy consumption per ton of crop. Fixed and variable fertilizer rates were calculated for an assumed average yield using yield historical data in this field from the last 4 years for each crop: 7.0 t ha −1 for spring barley; 5.0 t ha −1 for winter oilseed rape; 8.0 t ha −1 for winter wheat; 4.0 t ha −1 for faba beans. Therefore, such yields were used in further calculations. In this field, the last historical yields were received using a fixed application rate based only on nutrient uptake with the yield. For our study, it was important in 4 year crop rotation to get the same yield with a lower application rate, which was based not only on nutrient uptake with the yield but also using data about nutrient stocks in the soil. In the farming business, the main goal should be to get more profit, but not the highest yields. This can be achieved using modern precision farming technologies like VRF. The environmental impact assessment identified GHG emissions (kg CO 2eq ha −1 ) from different VRF and FRF fertilization methods. GHG emissions were calculated by multiplying the amount of fertilizer used for each method by the GHG equivalents (kg CO 2eq kg −1 ) given in Table 1. According to Lal et al. [29,30], the CO 2 emission equivalent per kilogram of N fertilizers was 4.96, for P 2 O 5 and K 2 O were 1.35 and 0.58 kg CO 2 kg −1 , respectively. The difference in GHG emissions for each fertilizer type between VRF and FRF methods was also calculated. In order to assess the environmental impact of different fertilization methods, the reduction of GHG emissions per ton of crop yield was calculated based on yield (for spring barley, winter oilseed rape, winter wheat, and faba beans 7.0, 5.0, 8.0, 4.0 t ha −1 , respectively).

Statistical Analysis
To ensure the statistical reliability of the data, the soil properties and yield studies of each crop were performed in 4 replicates of random samples. Each soil sample consisted of randomized 15-20 subsamples. Experimental study data were processed using one-way analysis of variance program ANOVA. The data was evaluated by calculating the smallest significant difference with the level of statistical significance (p < 0.05) using Tukey's HSD method [31,32]. The different letters in Table 2 indicate the significant difference between the values.  Figure 1). * Recommended interval means the optimum nutrient amount available for the most row crops in the soil according to its granulometric structure and soil group. The different letters indicate the significant difference between the values.

Soil Analysis
At the beginning of the research in 2015, the average of each element level in analyzed soil samples was above the optimum level (Table 2). At the end of the studies in 2019, the average K element concentration in the re-analyzed soil samples was reduced significantly (35.6%) to the optimum level, but the rest of the elements stayed above. The main reason for K element reduction was due to the sampling area No. 3 because of higher biomass and yield. Increasing yield leads to higher K element uptake. Phosphorus concentration in all samples in the soil was in or above the optimum level. The highest concentration for elements P and K was found in sample No. 4.
Mg and pH levels in all samples were above the optimum level at the beginning and at the end of the study.
The amount of elements in the soil changed in all the samples throughout the 4-year period: pH increased on average by 1.4%, Mg increase-1.7%, P decreased-12.5%, K decreased-35.6%, respectively ( Figure 7). Results received from the laboratory were imported into agronomic software Agriport. All nutrient distribution maps were generated automatically using this program. Agriport has algorithms, which helps to interpolate the spaces between the different sampling areas with different soil sample results. This feature helps to get a better visual view of the nutrient distribution in the field.
Analogous studies by Kulczycki and Grocholski [26] at the start of experimental field research found very high levels of potassium in the soil, and the element potassium exceeded the recommended levels of 82% of the total area. After six years of potassium fertilization at a variable rate, marked changes were observed-38% of the total field area was already attributed to the optimal element content in the soil, and very high potassium levels were found in only 18% of the field area. Meanwhile, the variable rate of phosphorus fertilization increased the field area with a low P element level by 9% and the medium element level by 6%, but the field area with a high P element level decreased by 7%, and the very high P element level decreased by 8%, respectively. Agronomy 2021, 11, x FOR PEER REVIEW 10 of 19 Figure 7. Distribution of soil elements (pH, P, K, Mg) according to the classes and the area. The class C or the optimal number of elements is yellowish and light green color; classes D and E are dark green and blue, indicating an excess of soil elements; letters A and B are red and brown, indicating the deficiency of elements. "+" next to a letter means higher and "−"-lower nutrient amount in the same element distribution class.

Additional N-Variable Rate Fertilization of Winter Wheat
This study presents additional fertilization of winter wheat at growth stages BBCH 30 and BBCH 37. In spring 2018, N uptake ( Figure 8) and N fertilization maps ( Figure 9) were generated with additional N fertilization at a variable rate in real-time with the Yara N-Sensor ALS. N uptake was highest in 15.2% of the total area (56-67 kg ha −1 ), and N application rate was lowest in 16.6% of the total area (66-76 kg ha −1 ) (Figure 8a). In the areas where crop N uptake was below the biomass threshold (17 kg ha −1 at BBCH 30), the fertilization rate was reduced from 66 to 0 kg ha −1 (1.6% of the total area).
During variable fertilization at BBCH 37 (Figure 8b), the average N uptake of winter wheat was 108 kg ha −1 and the average N fertilizer rate was 75 kg ha −1 . In the last additional fertilization, the field became more homogeneous. For example, when the N uptake went from 105 to 117 kg ha −1 , the field uniformity reached 84.3%. The uniformity of the field (71% of the total area) also affected the fertilizer application rate, which ranged from 68 to 78 kg ha −1 . After 16 days, N fertilization change in N uptake at BBCH 37 was on average 60 kg ha −1 . Figure 7. Distribution of soil elements (pH, P, K, Mg) according to the classes and the area. The class C or the optimal number of elements is yellowish and light green color; classes D and E are dark green and blue, indicating an excess of soil elements; letters A and B are red and brown, indicating the deficiency of elements. "+" next to a letter means higher and "−"-lower nutrient amount in the same element distribution class.

Additional N-Variable Rate Fertilization of Winter Wheat
This study presents additional fertilization of winter wheat at growth stages BBCH 30 and BBCH 37. In spring 2018, N uptake ( Figure 8) and N fertilization maps ( Figure 9) were generated with additional N fertilization at a variable rate in real-time with the Yara N-Sensor ALS. N uptake was highest in 15.2% of the total area (56-67 kg ha −1 ), and N application rate was lowest in 16.6% of the total area (66-76 kg ha −1 ) (Figure 8a). In the areas where crop N uptake was below the biomass threshold (17 kg ha −1 at BBCH 30), the fertilization rate was reduced from 66 to 0 kg ha −1 (1.6% of the total area).
During variable fertilization at BBCH 37 (Figure 8b), the average N uptake of winter wheat was 108 kg ha −1 and the average N fertilizer rate was 75 kg ha −1 . In the last additional fertilization, the field became more homogeneous. For example, when the N uptake went from 105 to 117 kg ha −1 , the field uniformity reached 84.3%. The uniformity of the field (71% of the total area) also affected the fertilizer application rate, which ranged from 68 to 78 kg ha −1 . After 16 days, N fertilization change in N uptake at BBCH 37 was on average 60 kg ha −1 . Agronomy 2021, 11, x FOR PEER REVIEW 11 of 19

Crop Yield
For calculation of the demand of mineral fertilizer, energy, and pollution indicators, the average yields were determined using historical yield data in this field from the last 4 years for each crop in the crop rotation. The average yields for spring barley, winter oilseed rape, winter wheat, and faba beans were 7.0, 5.0, 8.0, and 4.0 t ha −1 , respectively.

Crop Yield
For calculation of the demand of mineral fertilizer, energy, and pollution indicators, the average yields were determined using historical yield data in this field from the last 4 years for each crop in the crop rotation. The average yields for spring barley, winter oilseed rape, winter wheat, and faba beans were 7.0, 5.0, 8.0, and 4.0 t ha −1 , respectively.

Crop Yield
For calculation of the demand of mineral fertilizer, energy, and pollution indicators, the average yields were determined using historical yield data in this field from the last 4 years for each crop in the crop rotation. The average yields for spring barley, winter oilseed rape, winter wheat, and faba beans were 7.0, 5.0, 8.0, and 4.0 t ha −1 , respectively. These yields were taken as assumed and used in subsequent calculations as yields obtained using the FRF method. The average yields for the whole field using the VRF method for spring barley, winter oilseed rape, winter wheat, and faba beans were 7.3, 5.2, 8.1, and 3.9 t ha −1 , respectively.
During the threshing of winter wheat, a precise yield map was created ( Figure 10) and yield data were recorded: An average of 8.1 t ha −1 , a minimum of 5.4 t ha −1 , a maximum of 9.2 t ha −1 . After comparing the yield map with the N uptake map, it can be seen that the highest yield was in the places with the lowest N application rate, because in this site, specific areas contained more available N nutrients than other areas. During the application of variable N fertilization in the places where the maximum rate of 220 kg ha −1 was applied, the lowest yield was 5.4 t ha −1 . The reason could be that these site-specific areas had lower N availability in the soil. The N uptake map made at growth stage BBCH 37 showed more even nitrogen distribution in the crop matching to the yield map with smaller yield variability in the field. In 64% of the total area, the yield uniformity ranged from 6.8 to 8.4 t ha −1 . These yields were taken as assumed and used in subsequent calculations as yields obtained using the FRF method. The average yields for the whole field using the VRF method for spring barley, winter oilseed rape, winter wheat, and faba beans were 7.3, 5.2, 8.1, and 3.9 t ha −1 , respectively. During the threshing of winter wheat, a precise yield map was created ( Figure 10) and yield data were recorded: An average of 8.1 t ha −1 , a minimum of 5.4 t ha −1 , a maximum of 9.2 t ha −1 . After comparing the yield map with the N uptake map, it can be seen that the highest yield was in the places with the lowest N application rate, because in this site, specific areas contained more available N nutrients than other areas. During the application of variable N fertilization in the places where the maximum rate of 220 kg ha −1 was applied, the lowest yield was 5.4 t ha −1 . The reason could be that these site-specific areas had lower N availability in the soil. The N uptake map made at growth stage BBCH 37 showed more even nitrogen distribution in the crop matching to the yield map with smaller yield variability in the field. In 64% of the total area, the yield uniformity ranged from 6.8 to 8.4 t ha −1 . The threshing speed of the combine was uniform and accounted for 68% of the total threshing area in the range from 3.5 to 4.1 km h −1 . The average threshing rate of winter wheat was 3.7 km h −1 with average grain moisture of 14.4% and a yield of 8.1 t ha −1 . Feiffer et al. [33] found that a more even crop of winter wheat resulted in a higher threshing rate of 9-33%.
Grain moisture was also uniform and accounted for 75% of the total area in the range of 13.8-14.9% moisture. More evenly dried grain can be threshed earlier, and more even threshing speed results in higher harvester's performance.
The highest yield of winter wheat 9.2 t ha −1 was obtained by performing two additional fertilizations at growth stages BBCH 30 and BBCH 37 in the places where the lowest amount of N fertilizer was applied at a variable rate of 139 kg ha −1 . The threshing speed of the combine was uniform and accounted for 68% of the total threshing area in the range from 3.5 to 4.1 km h −1 . The average threshing rate of winter wheat was 3.7 km h −1 with average grain moisture of 14.4% and a yield of 8.1 t ha −1 . Feiffer et al. [33] found that a more even crop of winter wheat resulted in a higher threshing rate of 9-33%.
Grain moisture was also uniform and accounted for 75% of the total area in the range of 13.8-14.9% moisture. More evenly dried grain can be threshed earlier, and more even threshing speed results in higher harvester's performance.
The highest yield of winter wheat 9.2 t ha −1 was obtained by performing two additional fertilizations at growth stages BBCH 30 and BBCH 37 in the places where the lowest amount of N fertilizer was applied at a variable rate of 139 kg ha −1 .

Fertilizer Use and Savings
The results of comparison of VRF based on field soil analysis and crop nutrient demand with FRF, which takes into account only crop nutrient removal with the assumed average yield (Table 3). The results showed that less fertilizers were consumed when using the VRF method in 4-years compared to FRF: N 24 kg ha −1 (4.6%), P 2 O 5 124 kg ha −1 (48.6%), K 2 O 93.1 kg ha −1 (47.5%), respectively. The amounts of P element in the soil before the start of the experimental studies were higher (on an average three times in the whole field) than the recommended optimal amount. Therefore, in the phosphorus fertilization maps for variable fertilization, the rate of phosphorus fertilizers was reduced. A fixed rate of phosphorus fertilization was determined at the traditional recommended rate for each crop as shown in Table 3. These reasons resulted in more phosphorus savings than other fertilizers (N and K 2 O). Other researches stated that applying VRF of phosphorus fertilizer during the study period decreased the variability of the amount of this element in the soil [26]. The second-highest savings were achieved using potassium fertilizers, and average K 2 O level stayed at the optimum level at the end of the study. However, it needs to be clarified that the smallest N savings were achieved because variable N application software was created not for all crops. In our study, N application software was used only for winter wheat and oilseed rape. This could be the main reason why so little savings was found in N-fertilizer during the 4-year crop rotation. Other research have recorded 4 to 37% of nitrogen fertilizer savings in wheat in Turkey [4]. Yara International Research Centre (Hanninghof, Germany) applied tested winter wheat fertilization with N fertilizer at a variable rate with Yara N-Sensor in 2001-2005. Their test resulted in 12% lower fertilizer consumption compared to uniform fertilization [34]. Koch et al. [35] found that variable rates of maize fertilization can save from 6 to 46% of the fertilizer compared to a uniform rate. Other researchers compared urea fertilization between FRF (105 kg ha −1 ) and VRF (64 kg ha −1 ) in maize production. According to Greenseeker's recommendation, they obtained an average yield of 3.9 t ha −1 in both fertilization methods [36]. Notes: VRF-variable rate fertilization; FRF-fixed-rate fertilization.
Our research found that the total amount of saved NPK fertilizers was 241.6 kg ha −1 . In 4-years, from the total field were saved: 1418.8 kg of P 2 O 5 , 1058.3 kg of K 2 O, 272.9 kg of N. Depending on the crop type, NPK fertilizer savings were uneven. The highest NPK fertilizer savings were in winter oilseed rape-1061.8 kg (28.6%), spring barley-446.7 kg (40.1%), winter wheat-731.0 kg (20.0%), and faba beans-510.4 kg (43.3%). The research showed that while maintaining the optimal average nutrient content in the field, less fertilizer was used to grow different crops.

Energy Use Reduction with VRF
Compared VRF with FRF in spring barley (2016), energy consumption related to P 2 O 5 and K 2 O fertilizers was lower at 41.0% and 38.9%, respectively (Table 4). N variable fertilization was not applied for spring barley. Notes: VRF-variable rate fertilization; FRF-fixed-rate fertilization; total field area-11.37 ha.
In 2017 for winter oilseed rape applying VRF, energy consumption related to fertilizer P 2 O 5 was lower 63.0%, K 2 O-48.4% and N-6.9%, compared to FRF.
In 2018 for winter wheat using VRF, energy consumption related to P 2 O 5 fertilizers was lower 48.4%, K 2 O-46.5% and N-5.2%, compared to FRF practice.
In 2019 for faba beans applying VRF, energy consumption related to P 2 O 5 fertilizers was lower 30.6%, K 2 O-54.1%, compared to FRF method. N variable fertilization was not applied for faba beans.
Comparing VRF with FRF, the lowest energy savings were obtained by fertilizing winter wheat with N fertilizers, the highest-faba beans with K 2 O.
Summarizing the crop rotation in 4-years, when the VRF fertilization method was applied, energy consumption for phosphorus fertilizers was 48.6% lower and for potassium fertilizers-47.5% lower compared to the FRF fertilization method. Energy consumption related to nitrogen fertilizers was 4.6% lower compared to FRF. With variable fertilization, the total energy consumption for all fertilizers was 9.7% lower compared to FRF. The results showed that energy savings by reducing the application of P and K were generally lower compared to N fertilizer with respect to the energy use per ha due to fertilizer application. Summarizing the study's data on the energy consumption associated with NPK fertilizers over 4-years in the entire experimental field, it was found that the energy consumption of VRF was reduced by 9.7% Field −1 compared to the FRF method.

Changes in Fertilization Reduce Environmental Pollution
By applying VRF according to a pre-arranged fertilization plan and using fertilizer application maps and efficient distribution equipment, we can reduce both fertilizer costs as well as GHG emissions. Comparison of a 4-year VRF with FRF method allowed concluding that lower GHG emissions were obtained due to the use of fertilizers: 119.0 kg CO 2eq ha −1 of N, 168.5 kg CO 2eq ha −1 -P 2 O 5 , 54.0 kg CO 2eq ha −1 -K 2 O, respectively (Table 5). Similar research in Latvia showed that a variable N fertilizer application in winter wheat reduced GHG emissions by 46.8 kg CO 2eq ha −1 or 6% compared to a fixed rate N fertilizer application technology [37]. Notes: VRF-variable rate fertilization; FRF-fixed-rate fertilization; total field area-11.37 ha.
Our study found an overall reduction in GHG emissions of 341.5 kg CO 2eq ha −1 . Most greenhouse gas emissions reductions using the VRF method were achieved with P 2 O 5 (168.5 kg CO 2eq ha −1 ). Fertilization with NPK fertilizers in the whole field at a variable rate over 4-years reduced the GHG emissions by 3882.6 kg CO 2eq (11.27%). This study shows that GHG emissions in the experimental field decreased over 4-years, depending on the type of fertilizer: 48.6% for P 2 O 5 , 47.5%-K 2 O and 4.6%-N, respectively. Balafoutis et al. [16] reported that the use of VRF technology can reduce the GHG emission potential of 5% attributed to the use of N mineral fertilizers without yield effects. Italy scientist's findings suggest that VRF is able to decrease CO 2 emission in comparison with FRF [38].
According to the crop rotation, different GHG emissions were obtained by growing different crops and applying VRF, compared to FRF method in the whole study field: 460.2 kg CO 2eq Field −1 for spring barley, 1755.9 kg CO 2eq Field −1 -winter rape, 1242.9 kg CO 2eq Field −1 -winter wheat, 423.6 kg CO 2eq Field −1 -faba beans, respectively. The study showed that the greatest savings in GHG emissions were achieved in winter oilseed rape.

Economic Benefits
An evaluation of the economic benefits of changing a fixed fertilization rate to a variable rate in order to save the cost for fertilizers for the entire research field over 4 years was performed. It was revealed that when the price of N was 0.9€ kg −1 and the amount of fertilizer used was 272.9 kg, the economic benefit was achieved 21.6€ ha −1 . When the price of the fertilizer P 2 O 5 was 0.8€ kg −1 and K 2 O 0.5€ kg −1 and their consumption was 1418.8 and 1058.3 kg, then the savings were 99.8 and 46.5€ ha −1 , respectively. The total 4-year economic benefit for the whole field due to the savings on NPK fertilizers was 168.0€ ha −1 . Other researchers [34], who compared fertilization at a fixed rate with the one at a variable rate, calculated an economic benefit of 19$ (~16.4€) per hectare from savings on nitrogen fertilizers. Diacono et al. [3] analyzed the cost-effectiveness of variable rate N fertilizer management methods in wheat crops and concluded that real-time N-based management had the higher profitability, ranging from 5 to 60$ ha −1 (~4.3 to~51€ ha −1 ) compared to FRF. Dobermann et al. [39] found that fertilization of rice with NPK fertilizers using VRF can increase the yields by 12% compared to FRF.
A comparative analysis of the yields showed that VRF had an advantage over FRF due to higher crop yields. In our studies, no actual yield results were obtained using fixed rate fertilization. Yields for all four plants were assumed on the basis of average values from previous years. Using the VRF system, the yields of different crops obtained were very similar to the assumed yields used for the FRF method. Other authors also found that the average wheat grain yield was 4.24 t ha −1 and the mean increase in wheat grain yield applying VRF was 0.8% when compared with the FRF treatment averaged over 10 sites and two years [40]. In the study of scientists from Italy, there were no significant differences of the barley grain yield between FRF and VRF technologies observed [38]. The method used to apply variable N fertilization has been proved to be effective, leading to a similar barley yield to conventional practice while using less fertilizer (75 kg N ha −1 ). The VRF method has been proved to be a valid alternative to conventional fertilization, considering economic evaluation and leading a saving of 266€ ha −1 [38].
All crops yielded economic benefits, except for faba beans. A 4-year benefit due to the increased yield was 108.0€ ha −1 . Comparing VRF and FRF methods of spring barley, winter oilseed rape, winter wheat and faba beans, the economic benefits from the yield were 45, 70, 17 and −24€ ha −1 , respectively. Economic benefit from yield was calculated using the grain price for each crop: 150€ ha −1 for spring barley, 350€ ha −1 -winter oilseed rape, 170€ ha −1 -winter wheat and 240€ ha −1 -faba beans, respectively.
The greatest savings of energy input and GHG emissions were obtained with N and P 2 O 5 fertilizers at a variable rate in winter oilseed rape. Comparing different fertilizers, the best cost savings results were achieved using P 2 O 5 fertilizers for winter oilseed rape, winter wheat, and spring barley 45, 24.8, and 18.4, respectively. For all fertilizers, energy inputs, GHG emissions, and cost savings using VRF comparing to FRF are presented in Figure 11.  Figure 11. Fertilizers, energy inputs, GHG emissions and cost savings using VRF comparing to FRF.

Conclusions
A 4-year study of soil elements showed that even after 4-years of using VRF, excess phosphorus of 48.5 mg per 100 g of soil and potassium with 27.6 mg per 100 g of soil were still present in certain field areas. Repeated soil element analysis showed that application of VRF, taking into account the existing nutrient content in the soil and demand of elements for the crop, can save fertilizers and reduce environmental impact without overfertilizing P2O5 and N in the areas where there was already an excess of these elements. In the experimental field, the application of VRF maps allowed avoiding over-fertilization in separate field locations, thus saving fertilizers and making the obtained yield closer to the planned one. Application of this method allowed reducing the costs of fertilizers: 4.6% for N, 48.6%-P2O5 and 47.5%-K2O. The total consumption of NPK fertilizers was lower by 24.9% compared to the FRF. Thus, variable rate fertilization results in higher fertilizer use efficiency.
Lower consumption of mineral fertilizers resulted in lower energy consumption throughout the 4-year crop rotation. With VRF, the total energy input for all fertilizers was 3463.1 MJ ha −1 or 9.74% lower than with FRF. The largest reduction in energy consumption was found for nitrogen and phosphorus fertilizers, 1454.4 and 1385.1 MJ ha −1 , Figure 11. Fertilizers, energy inputs, GHG emissions and cost savings using VRF comparing to FRF.
Summarizing the research results obtained in this study, it can be stated that variable fertilization is an important technological process that reduces the consumption and the costs of NPK fertilizers, saves energy costs, and contributes to the preservation of a cleaner environment. In order to achieve even stronger positive energy, as well as economic and environmental impact, future research needs to investigate the technological processes of not only fertilization but also of sowing and spraying of plant protection products at a variable rate.

Conclusions
A 4-year study of soil elements showed that even after 4-years of using VRF, excess phosphorus of 48.5 mg per 100 g of soil and potassium with 27.6 mg per 100 g of soil were still present in certain field areas. Repeated soil element analysis showed that application of VRF, taking into account the existing nutrient content in the soil and demand of elements for the crop, can save fertilizers and reduce environmental impact without over-fertilizing P 2 O 5 and N in the areas where there was already an excess of these elements. In the experimental field, the application of VRF maps allowed avoiding over-fertilization in separate field locations, thus saving fertilizers and making the obtained yield closer to the planned one. Application of this method allowed reducing the costs of fertilizers: 4.6% for N, 48.6%-P 2 O 5 and 47.5%-K 2 O. The total consumption of NPK fertilizers was lower by 24.9% compared to the FRF. Thus, variable rate fertilization results in higher fertilizer use efficiency.
Lower consumption of mineral fertilizers resulted in lower energy consumption throughout the 4-year crop rotation. With VRF, the total energy input for all fertilizers was 3463.1 MJ ha −1 or 9.74% lower than with FRF. The largest reduction in energy consumption was found for nitrogen and phosphorus fertilizers, 1454.4 and 1385.1 MJ ha −1 , respectively. Lower spreads of fertilizers all over the field resulted in lower GHG emissions to the environment during the study period. Phosphorus and potassium VRF application was estimated to reduce the GHG emissions by almost twice compared to FRF. For the total amount of NPK fertilizers, 11.3% lower GHG emissions were obtained. At the end of the study, GHG emission reductions of 341.5 kg CO 2eq ha −1 or 3882.6 kg CO 2eq Field −1 were determined.
VRF, according to the application plan and by using VRF maps, allowed a reduction in NPK fertilizers by 241.9 kg ha −1 or 2750.0 kg Field −1 . A significant decrease in the need of mineral fertilizers reduced production costs and contributed to cleaner agricultural production. The total economic benefit coming from the saved funds for fertilizers due to an application of VRF was 168.0€ ha −1 or 1909.7€ Field −1 .
Based on the results of our study, we can say that the application of the VRF method allows farmers to benefit from lower costs for fertilizers and, at the same time, helps them to meet the EU's requirements for environmental pollution in agriculture. The energy, economic, and environmental benefits of variable rate fertilization are a great recommendation to continue this type of research in the future, including more details on precision farming such as variable rate sowing, spraying of growth regulators, and plant protection products.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.