Influence of Soil Properties on Maize and Wheat Nitrogen Status Assessment from Sentinel-2 Data

Abstract

Soil properties variability is a factor that greatly influences cereals crops production and interacts with a proper assessment of crop nutritional status, which is fundamental to support site-specific management able to guarantee a sustainable crop production. Several management strategies of precision agriculture are now available to adjust the nitrogen (N) input to the actual crop needs. Many of the methods have been developed for proximal sensors, but increasing attention is being given to satellite-based N management systems, many of which rely on the assessment of the N status of crops. In this study, the reliability of the crop nutritional status assessment through the estimation of the nitrogen nutrition index (NNI) from Sentinel-2 (S2) satellite images was examined, focusing of the impact of soil properties variability for crop nitrogen deficiency monitoring. Vegetation indices (VIs) and biophysical variables (BVs), such as the green area index (GAI_S2), leaf chlorophyll content (Cab_S2), and canopy chlorophyll content (CCC_S2), derived from S2 imagery, were used to investigate plant N status and NNI retrieval, in the perspective of its use for guiding site-specific N fertilization. Field experiments were conducted on maize and on durum wheat, manipulating 4 groups of plots, according to soil characteristics identified by a soil map and quantified by soil samples analysis, with different N treatments. Field data collection highlighted different responses of the crops to N rate and soil type in terms of NNI, biomass (W), and nitrogen concentration (Na%). For both crops, plots in one soil class (FOR1) evidenced considerably lower values of BVs and stress conditions with respect to others soil classes even for high N rates. Soil samples analyses showed for FOR1 soil class statistically significant differences for pH, compared to the other soil classes, indicating that this property could be a limiting factor for nutrient absorption, hence crop growth, regardless of the amount of N distributed to the crop. The correlation analysis between measured crop related BVs and satellite-based products (VIs and S2_BVs) shows that it is possible to: (i) directly derive NNI from CCC_S2 (R² = 0.76) and either normalized difference red edge index (NDRE) for maize (R² = 0.79) or transformed chlorophyll absorption ratio index (TCARI) for durum wheat (R² = 0.61); (ii) indirectly estimate NNI as the ratio of plant nitrogen uptake (PNUa) and critical plant nitrogen uptake (PNUc) derived using CCC_S2 (R² = 0.77) and GAI_S2 (R² = 0.68), respectively. Results of this study confirm that NNI is a good indicator to monitor plants N status, but also highlights the importance of linking this information to soil properties to support N site-specific fertilization in the precision agriculture framework. These findings contribute to rational agro-practices devoted to avoid N fertilization excesses and consequent environmental losses, bringing out the real limiting factors for optimal crop growth.

1. Introduction

Crop nutritional stress is usually assessed considering N deficiency on the base of plant analysis methods with no direct measurements or in consideration of soil properties information. If not considered, soil variability can prevent a correct diagnosis of the nutritional status of crops. The influence of soil properties variations in relation to crop nitrogen status assessment performed with proximal or remote sensing technique is a topic that has not been much studied in literature. As reported by [1], several studies have demonstrated separately the potential advantages of crop spatial variability analysis with soil-based and plant-based methods in order to drive variable N fertilization, but only few studies assessed the interaction between the two, hence the possibility to combine this information in a decision process [1,2].

Cereals cultivation requires sustainability in terms of production and environmental impacts. Nitrogen is one of the most important nutrients for plants. Soil N availability varies spatially and temporally within a field due plant, soil–atmosphere interaction almost guided by the variations of soil properties and land morphology. Crop N demand changes during the season depending on growing conditions related to the presence of limiting factors. This requirement is met by soil N supply capacity and N fertilizations. N mineral fertilizers are globally the most used substances for fertilization, with a world supply of ammonia (NH₃) of ~180 million tons [3]. Unfortunately, the large quantities used do not correspond to high levels of nitrogen use efficiency (NUE) [4,5,6,7]. Sustainable agriculture requires optimizing N management during the crop season, in order to achieve profitability and a healthy environment [8,9,10]. N supply from soil and fertilizations and N accumulation in crops are very dynamic processes that should be considered through an integrated approach [11,12] and a spatial relationship analysis [13]. The uncertainty in both plant N demand and soil N supply due to seasonality and growth potential often leads to non-optimal nitrogen management strategies [14].

One of the most widely used crop N status assessment methods relies on the estimation of the nitrogen nutrition index (NNI), defined as the ratio of the actual N concentration in the dry above-ground biomass (Na%) to a critical N concentration (Nc%), which is in relation to the specific above-ground biomass (W) expressed in t ha⁻¹ [14]. NNI can also be derived considering crop N content (kg N ha⁻¹), instead of concentration (%). In this case, NNI is computed using the ratio between actual plant nitrogen uptake (PNUa) and critical plant nitrogen uptake (PNUc), both expressed in kg N ha⁻¹ [15,16]. This is because the amount of N in the plant, expressed as a concentration or quantity, and biomass accumulation, have a relationship based on the theory of N dilution during plant growth [17]. The critical nitrogen dilution curve derived from this theory [18] determines the minimum N concentration (Nc) or N content (PNUc) [16] in the crop, for a specific moment, which allows for a normal growth with no biomass accumulation reduction. Other authors identified N dilution curves as a function of development stages [19,20] or Leaf Area Index (LAI) [16,21,22,23] instead of using W. The use of LAI to obtain NNI estimation is indeed a promising approach, because it can be performed with non-destructive measures of W using in field or remote indirect methods [24,25]. Besides the field estimation of NNI, this indicator can also be calculated by empirical relations with remote sensing (RS) vegetation indices (VIs) or biophysical variables (BV) using proximal or satellite data. RS data can be used to derive both Na% and Nc% (via their relationships with W) or PNUa and PNUc. Many authors suggest that the latter option is preferable, because the relationship between RS data and total quantities of nitrogen present in the canopy is more robust [14,26].

From the literature, two approaches have been used to estimate NNI from satellite data: (1) direct methods (DM) and (2) indirect methods (IM). DM rely on the possibility to estimate NNI by “direct” relationships with Vis. According to [27], we can use a (i) “mechanistic approach” (DM_1) that requires first to estimate Na% and W to calculate NNI or (ii) a “semi–empirical approach” (DM_2) based on the direct definition of a parametric regression between VIs and NNI values that are usually estimated with in situ data. Many authors have successfully used DM of NNI from remote sensing using multispectral [23,28,29,30,31] or hyperspectral data [32,33] for optimizing the timing and the rate of N fertilizer applications. The study of [34] tested a DM_1 approach in an operational workflow for producing NNI maps based on the combined use of high-resolution satellite images and ground-based estimates of crop parameter using smart apps.

IM are based on a different paradigm: RS data are firstly used to retrieve biophysical variables from RS, such as leaf area index (LAI) and canopy chlorophyll content (CCC). Then, from these estimates, a relation with PNUa and W is adopted to derive PNUc according to a crop-specific dilution curve. Once PNUa and PNUc are estimated, NNI is calculated. IM approaches have the advantage to exploit RS data to estimate optically related biophysical variables (BV), hence to be a physically-based solution. BV estimation can be performed using (i) parametric regression methods using Vis [15,35,36,37], (ii) non-parametric machine learning algorithm (e.g., neural network, SVR, PLSR etc.) or by inverting radiative transfer models (RTM) adopting look up table approaches [38] or hybrid methods [39,40,41]. For a detailed review of the different available approaches, see [42]. Recently, the European Space Agency (ESA) Sentinel Application Platform (SNAP) toolbox, reached an operational stage in providing biophysical variables (LAI, fAPAR. fcover, CCC) through a biophysical processor tool based on neural networks [43]. This greatly facilitates the operational adoption of RS products in agricultural management decision support systems.

Under this framework, the objectives of this study were to investigate the robustness and representativeness of nutritional status assessment through the computation of NNI from Sentinel-2 for maize and durum wheat in relation to different nitrogen fertilization levels and soil properties. A field experiment for two consecutive crop seasons was set up to interpret the potential and limits of NNI in identify nitrogen deficiency.

More specifically, Sentinel-2 acquisitions were exploited to analyze spatial and temporal variability of crop status from remote sensing with the following objectives: (i) compare the performances of direct and indirect NNI estimation approaches; (ii) assess the added value of the operational S2-BV product with respect to ad hoc calculated VIs.

2. Materials and Methods

2.1. Study Area and Experimental Design

The study area is located in northern Italy, in the eastern part of the Ferrara province, in the Po valley. This flat area, between the city of Ferrara and the Adriatic Sea, derives from extensive land reclamation works concluded at the beginning of the 20th century. This operation generated new arable lands from zones that were previously permanently flooded marshland. This explains a great deal of heterogeneity in soil properties, affecting yields and crop growth.

A field experiment was carried out at the Bonifiche Ferraresi farm, located in Jolanda di Savoia (FE) (lat. 44°52′26″N, long. 11°56′24″E, elev. −3 m), in the 2018 and 2019 seasons, respectively, on maize and durum wheat. Previously, in 2015, the farm conducted a survey for the determination of soil properties using the automatic resistivity profiler (ARP) geophysical sensor [44] as well as soil sampling. From these data, a soil map of the whole farm was derived by the surveying company (SOING, Livorno, Italy). This information was used to identify areas with contrasting soil properties potentially related to different productivity.

For the 2018 and 2019 experiments, a 16 ha field (Figure 1) was chosen and sown with maize, cv SY HYDRO (Syngenta FAO 700) in 2018 and with durum wheat ODISSEO (Syngenta) in 2019 (

Table 1. Agronomic details of field experiments for maize in 2018 and wheat in 2019.

	Variety	Density	Plowing Date	Planting Date	Harvest Date	Irrigation
Maize	SY	8.8 seeds m2	2-Apr 2018	17-Apr 2018	25-Aug 2018	Sub-irrigation
Wheat	ODISSEO	200 kg ha−1	20-Oct 2018	12-Nov 2018	15-Jun 2019	No irrigation

). Four different soil classes, coded as CDS2, JOL3, JOL1, FOR1 were present within the field as detected by a pedological survey conducted in 2015 (Figure 1b). A total of 20 plots were established in the field, in which 5 N fertilizer rates were applied with 4 replicates. Each replicate was located inside one of the identified soil classes (Figure 1b) with the goal of maximizing the intra-field nutritional variability. Plot size, 35 m width by 50 m length, was established based on the width of a variable rate (VR) fertilizer spreader, i.e., Sulky model x36 (Sulky-Burel, Chateaubourg, France), and for being large enough to contain at least 10 pixels of a 10 m resolution S2 image.

The N treatments applied to maize and wheat respectively in the 2018 and 2019 season are reported in

Table 2. Experimental design details.

Year	Crop	N Treat. Code	N Rate (kg N ha−1)	BBCH a	N Application Date	BBCH b	Field Campaign Date	BBCH c	S2 Acquisition Date
2018	Maize	N0	10
		N1	90	00 d	17-Apr	12	09-May	13	19-May
		N2	180	13	21-May	15	29-May	16	03-Jun
		N3	300	16	01-Jun	30	18-Jun	50	23-Jun
		N4	350	60	26-Jun	65	03-Jul	65	08-Jul
		Standard	240
2019	Durum Wheat	N0	10
		N1	60	00 d	20-Oct
		N2	120	22	04-Mar	25	19-Mar	25	20-Mar
		N3	180	237	06-Apr	31	19-Apr
		N4	240	51	02-May
		Standard	160

^a Crop growth stage BBCH [45] at the date of the N fertilization. ^b Crop growth stage BBCH at the date of the field campaign. ^c Crop growth stage BBCH at the date of the S2 acquisition. ^d 00 refers to the N application at the sowing date (basal dressing).

. For both crops, two N treatments (N1 and N2) were below the standard N prescription rate of the farm, in order to create N limiting conditions, whereas two treatments (N3 and N4) were above, to create theoretically non-limiting N conditions. Plots with the N0 treatment were considered as non-fertilized controls, but actually received 10 kg N ha⁻¹ because a uniform basal dressing was applied for both crops at sowing with 125 kg ha⁻¹ of 8-18 fertilizer (i.e., 10 kg N ha⁻¹) as a standard practice at the farm. A uniform basal dressing with 200 kg ha-1 Pheoscor was also applied before sowing for both crops, whereas no K fertilizer was applied because of good soil content.

During the season, top dressing fertilizations with urea using prescription maps were carried out according to phenological key stages (Table 2), with the exception of N0 plots. Application maps (as-applied maps) recorded by the VR spreader system, allowed to check the amount of N provided and the accuracy of the spreading in terms of geolocation. In the 2018 season, the last N application at flowering was carried out with a uniform distribution of N, due to the absence of VR machinery for liquid fertilizer distribution. The field was irrigated by means of a drainage sub-irrigation system (Table 1). Other practices for pest, diseases, and weeds control were adequately carried out according to usual farm management practice for both crops. An effort was made to obtain non-limiting conditions for all factors, except for nitrogen fertilization.

2.2. Field Data Collection

In April 2018, an intensive soil sampling campaign was carried out, stratified in the more contrasting areas identified by the resistivity map of the experimental field, collecting a total of 36 soil samples, from a 0–30 cm depth. The number of samples was chosen in order to have 9 points for each of the soil classes identified from the resistivity map, as a compromise between accuracy and resources available (Figure 1). Soil texture, soil organic matter (SOM), total nitrogen, total organic carbon (TOC), pH, C/N were measured in the laboratory using official standard analytical protocols. Along the season, field measurements were carried out to acquire ground data of crop BVs (Table 2).

An ESU (elementary sample unit) [46] of 20 by 20 m size, was defined at the center of each plot and a Trimble GNSS receiver was used to conduct the measurements at the center of each ESU. A Canon digital camera with a fish-eye lens was used to acquire 16 digital hemispherical photographs (DHP), randomly spread in each ESU, at each date and processed with the Can_Eye V6.494 software [47] to estimate crop green area index (GAI). GAI includes all the green parts of the crop and not just the leaves, as for the leaf area index (LAI). The effective_GAI (GAI_eff), computed by the Can_Eye software, was used and not the true_GAI, considering the fact that the former is more theoretically related to the LAI estimated from satellite images [39,43,48]. Leaf chlorophyll content (Cab – mg cm² leaf) readings from 60 plants in each ESU were acquired with the Dualex (Force-A, Orsay, France) at each measurement date. GAI_eff and Cab were used to compute the chlorophyll content (CCC) by multiplying these two measured variables. Above ground biomass (W) was measured by harvesting 8 plants for each ESU for maize and 1 m² of wheat crop, which were then dried in the oven at 80 °C until constant weight. The whole plant samples were then finely ground and used for N content (Na%) determination using the Kjeldhal method. The plant phenological stages were identified and registered using the BBCH scale (Biologische Bundesanstalt, Bundessortenamt and CHemical industry) [45].

2.3. Satellite Data

Sentinel-2 SAFE (Standard Archive Format for Europe) level-2A (Bottom-of-Atmosphere BOA) cloud free images, both from S2-A and S2-B satellites, were downloaded and processed using the sen2r R toolbox [49] to produce L2A products and obtain the VIs listed in

Table 3. List of vegetation indices used in the present study.

Index	Name	Reference
NDVI	Normalized Difference Vegetation Index	[52]
NDRE	Normalized Difference Red Edge Index	[53]
MSAVI	Modified Soil-Adjusted Vegetation Index	[54]
GNDVI	Green Normalized Difference Vegetation Index	[55]
MCARI	Modified Chlorophyll Absorption in Reflectance Index	[56]
MCARI2	Modified Chlorophyll Absorption in Reflectance Index 2	[57]
MTVI2	Modified Triangular Vegetation Index 2	[57]
OSAVI	Optimized Soil Adjusted Vegetation Index	[58]
TCARI	Transformed Chlorophyll Absorption Ratio Index	[59]
TCARI/ OSAVI	Transformed Chlorophyll Absorption Ratio Index/Optimized Soil Adjusted Vegetation Index	[59]
TCI	Triangular Chlorophyll Index	[60]

. The ESA SNAP biophysical processor toolbox [43] was used to retrieve the satellite-derived biophysical variables (S2_ BVs) LAI_S2 and CCC_S2. LAI_S2 estimated from remote sensing includes all the green parts of the plant and therefore can be called GAI [48]. This variable is similar to the effective LAI output of CAN_EYE [39,43,48]. CCC_S2 is a quantitative variable (g of chlorophylls a + b per m² of soil) calculated at canopy level, strongly related to leaf nitrogen content [26,39] and it is directly estimated from the SNAP toolbox (LAIxCab). S2Toolbox level2 biophysical variables algorithms implemented in SNAP [43] allow to estimate GAI and CCC and not Cab. Since Cab estimation is fairly complicated due to the gradient of Cab content in the plant and the complexity of canopy structure that influence the radiometric response, the CCC variable is provided, to avoid ambiguities between the GAI and Cab estimation during the inversion process [43]. The spectral domain where Cab can be estimated with high accuracy is in the visible and near infrared and also GAI is influencing the reflectance in this domain, especially for dense canopies [31]. For this reason, CCC is used as an indicator of N stresses [25,39,50], because it includes information on the leaf nitrogen content and the photosynthetic capability [51].

In the 2018 season, four S2 images acquired after N applications and close to field measurement dates were processed, whereas only one useful S2 image close to the field campaign dates was available in 2019 (Table 2), due to cloud occurrence and lodging problems that affected wheat crop before flowering. A vector layer with plots boundaries was used to extract VIs/BVs values from images considering an internal buffer area. Pixels values were averaged for each plot and compared with ground-measured variables.

2.4. NNI Determination

NNI calculation is based on the concept of the dilution curve [18] derived by the theory of decline of N concentration (Nc%) in plants in relation to their growth [17]. Nc is defined as the minimum N concentration that is necessary to achieve potential non limiting production [61]. The critical N dilution curves of both maize and durum wheat crops were calculated using Equation (1), adopting the crops specific coefficients (a, b) proposed for maize and wheat by [16] (a = 3.4, b = 0.37) and [18] (a = 5.35, b = 0.442), respectively:

$$\mathrm{Nc}\%={\mathrm{aW}}^{-\mathrm{b}},$$

(1)

where W is above-ground biomass (t ha⁻¹), a represents the plant N concentration when W is 1 t ha⁻¹ and b is a dimensionless coefficient. Nc % is constantly considered 3.4% for W < 1 t ha⁻¹ for maize and 4.4% for W < 1.55 t ha⁻¹ for wheat according to [16,18]. The NNI was calculated using Equation (2) as the ratio between the actual plant N concentration (Na%) and the critical plant N concentration (Nc%):

$$\mathrm{NNI}=\frac{\mathrm{Na}}{\mathrm{Nc}},$$

(2)

and was used to define conditions of N deficiency (NNI < 0.9), N surplus (NNI > 1.1), and optimal status (0.9 < NNI < 1.1). NNI can also be determined as the ratio between the actual plant N uptake (PNUa, kg ha⁻¹) and the critical plant N uptake (PNUc, kg ha⁻¹) exploiting the theory of the critical nitrogen absorption curve [15,16]. PNUa and PNUc were calculated by multiplying W by Na% and Nc%, respectively. Excess or deficit of nitrogen uptake can be determined by the relationship between the nitrogen absorbed by the crop (PNUa kg ha⁻¹) and its biomass [15]. Then, a quantitative indicator ΔPNU (kg ha⁻¹) [26] can be determined using Equation (3) and directly used to calculate the N rates (increase or reduction of planned dose) to be considered for site-specific fertilization purpose:

$$\mathsf{\Delta }\mathrm{PNU}=\mathrm{PNUa}-\mathrm{PNUc}.$$

(3)

2.5. NNI Estimation from Remote Sensing

In this study, two approaches have been tested for NNI calculation from RS data. For method 1, (DM) linear regressions between selected Vis and field measured Na% and NNI were derived, while for method 2, (IM) S2_BVs were first calculated from SNAP biophysical toolbox and then linear regressions to measured PNUa and W was exploited in order to determine NNI. The studies of [15,26] proposed three IM to estimate PNUa and PNUc from remote sensing data for then calculating ΔPNU or NNI, evaluating the third scheme as the best in terms of prediction. This scheme consists of calculating PNUa directly from CCC using experimental relationships and then obtaining PNUc by estimating W from GAI. The NNI values obtained from the two methods were then compared with field measured NNI. In order to further test the information content of NNI as a proxy for nitrogen deficiency diagnosis, an NNI map for all the 2018 maize cultivated fields was produced with S2 data and compared to field level reported yield. The regression model between VIs and BVs calculated from S2 images and crop variables collected during field experiments were analyzed using statistical parameters using R 3.5 version software [62]. The performance of the different regression models was estimated by analyzing the coefficient of determination (R²), mean absolute error (MAE), root mean square error (RMSE) (

Table 4. Statistics used to compare the model’s performance in regression analysis; p is value derived from remote sensing (Vis and S2_BVs) and o is the observed value (BVs). Ō is the mean value of o based on N observations.

Statistic	Equation
Coefficient of Determination	${\mathrm{R}}^2=1-\frac{{{\displaystyle \sum }}^{}{\left({\mathrm{p}}_{\mathrm{i}}-{\mathrm{o}}_{\mathrm{i}}\right)}^2}{{{\displaystyle \sum }}^{}{\left({\mathrm{p}}_{\mathrm{i}}-{\mathrm{ō}}_{\mathrm{i}}\right)}^2}$
Mean Absolute Error	$\mathrm{MAE}=\frac{1}{\mathrm{N}}{\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}}}\left|{\mathrm{p}}_{\mathrm{i}}-{\mathrm{o}}_{\mathrm{i}}\right|$
Root Mean Square Error	$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{N}}{\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}}}{\left({\mathrm{p}}_{\mathrm{i}}-{\mathrm{o}}_{\mathrm{i}}\right)}^2}$

).

The results of soil laboratory analyses were assessed using one way analysis of variance (ANOVA) and Tukey HSD post hoc test (p < 0.05) in statistical software R 3.5 version [62] to evaluate differences in soil properties among soil classes.

3. Results

3.1. Soil Variability in the Experimental Field

In general, laboratory analyses revealed that the field had a clay soil (63%) with generally good fertility, as indicated by high total organic carbon (TOC) and soil organic matter (SOM) (6%), moderate total N content (0.34%), and good C/N (10.25%). The comparison among soil classes showed a statistically significant difference in pH and in total N (%) for the FOR1 class with respect to the others (

Table 5. Results of the soil sampling carried out in the experimental field.

Soil Class	N Total (%)	TOC (%)	C/N (%)	SOM (%)	Clay (%)	Silt (%)	Sand (%)	pH (-)
CDS2	0.33 a	3.44	10.39	5.93	63.72	24.01	12.28	7.56 a
JOL3	0.33 a	3.30	10.07	5.7	63.84	22.78	13.39	7.24 a
JOL1	0.32 a	3.40	10.68	5.87	62.28	25.66	12.06	7.63 a
FOR1	0.40 b	4.00	9.86	6.89	62.24	26.85	10.91	5.48 b
Mean	0.34 *	3.54	10.25	6.1	63.02	24.82	12.16	6.98 *

* Significant differences among means (p < 0.05). The letters a, b, represent statistically significant difference among groups (p < 0.05) in post hoc test; * significant differences among means (p < 0.05).

). The differences in pH, sub-acidic for FOR1 (pH = 5.45) compared to the sub-alkaline of the other classes (pH~7.5), respect to N total content, are considered quite important because of: (i) the general high amount of N total in the soil; (ii) the potential impact of pH on elements availability in the soil for plant nutrition.

3.2. Crop Response to Fertilization and Soil Influence

3.2.1. Plant Growth and Yield

During the 2018 season, maize was harvested at the phenological phase of the milky stage to produce silage of maize mash (grains and maize cob). The harvest was weighed at the farm Weighbridge showing a yield of 14.2 t ha⁻¹. In 2018, the farm harvested 278 ha of silage maize in total, with an overall average yield of 18.22 t ha⁻¹, therefore the experimental field produced less than the farm average (considering FAO 600–700 classes) and less than the variety yield potential (HYDRO FAO 700). The 1st of June 2018 during phase BBCH 15 the whole farm and the experimental field were affected by a hailstorm that caused serious damage to plant foliage with different intensity for each field.

For the maize experiment, crop response to variable N fertilization showed different patterns depending on N rate and soil classes. In general, it is expected that maize biomass would increase with increasing N rates [1,18,31]. In our case, biomass response was strongly influenced by soil properties as shown in Figure 2 and the biomass trends are not always correlated with increasing N rates, indicating that other limiting factors might have affected crop N uptake.

On the first measurement date (09 May, BBCH 12), the effect of N fertilization on maize biomass was still not visible (Figure 2a). On the other dates, whereas for most soil classes, a rate of 90 or 180 kg N ha⁻¹ allowed to reach the maximum biomass, the plots with soil type FOR1 showed a remarkably lower response to N (Figure 2b–d). This effect was probably due to other nutritional limitations, caused by the low pH of this soil class (Table 5). For example, it can be seen from Figure 2, that the JOL1 soil class is characterized by increasing biomass values according to N application on both phases, BBCH 15 (Figure 2b) and BBCH 65 (Figure 2d), whereas the FOR1 soil class remained stable with a lower biomass accumulation. In general, a limited or negative effect on biomass production was obtained for the highest N rates in this experiment. In particular, the range of measured biomass values at flowering (

Table 6. Maize above ground biomass (t ha⁻¹) measured during the last field campaign of 2018 (BBCH 65) and analyzed by the ANOVA method (different letters are significantly different at p < 0.05).

Soil Class	Max	Min	Mean	Tukey’s Post-Hoc Test
JOL1	10.2	5.0	8.3	a
CDS2	10.0	3.4	7.9	ab
JOL3	8.0	3.7	6.1	ab
FOR1	6.3	1.6	4.2	b

), which are expected to be related to the final yield, highlights how the biomass values reached in JOL1 almost doubled the biomass level achieved in the FOR1 soil class and how the biomass in the different soil classes are statistically different.

Different N treatments showed a substantial effect on Na% changes across different growth stages, except for the first measurement date (Figure 2). During subsequent growth phases, similar patterns of response of Na% to N rates were observed among the soil classes. Unlike biomass, Na% showed increasing concentration values for increasing N rates in all soil classes. The FOR1 class achieved high concentrations of nitrogen in highly fertilized plots as opposed to what was seen for biomass (Figure 2f–h).

Concerning the durum wheat 2019 experiment, marked effects of different soil classes of the field, on crop response to N rates, were also observed. In Figure 3, the results of only 2 field measurement dates are presented, due to lodging problems caused by wind occurring later in the season. At the BBCH 31 stage, plots in JOL1 registered the highest biomass values (Figure 3b) in all N treatments (reaching 3.24 t ha⁻¹ at the highest N4 rate), while plots in FOR1 yielded the lowest biomass, less than half of the total biomass of JOL1 at the highest N rate (1.1 t ha⁻¹). The interpretation of the Na% trend in wheat is more complex with respect to maize, since no clear response to N rate is visible (Figure 3), especially for the FOR1 soil class, which showed low Na% values for high levels of N fertilization. These results confirm differences of crop response to different soil classes to the N fertilization in both years. In some zones of the field (JOL1, CDS2), the crop manages to reach substantial W and Na% values, even if it is less fertilized, whereas in others, is not able to produce high biomass even if it is highly fertilized (FOR1). An apparent non responsiveness of crop biomass and N content to N fertilization rates in a specific zone of the field, was observed for both seasons and crops, and can be associated to static variables of soil properties, such as in this case, the pH.

3.2.2. Nutritional Status–NNI Values

Figure 4 provides a representation of crop nutritional status for each measured plot in relation to N supply, data are plotted in the W (x-axis) and Na space (y-axis) and compared to the critical N from the dilution curves for maize [16] and durum wheat [18].

The calculated NNI values ranged between 0.42 and 1.12 and 0.43 and 1.15 for maize and wheat respectively, revealing the condition of plot with N deficiency (NNI < 0.9). NNI values under varied N applications did not show a clear gradient as expected and observed in other studies [30].

In Figure 4, the grayscale of the symbols, from light grey to black, is proportional to the N application. It would be expected that brighter grey points (N0 and N1 treatments) are more distant from the critical curve with respect to the darker ones (N3 and N4 treatments). This is visible only for the points corresponding to the low doses (N0-N1-N2), but it is not equally respected for the high doses (N3 and N4), evidencing how other factors might have limited nitrogen uptake and plant growth even in non-limiting N conditions. The NNI values for N0 and N1 generally decreased during all BBCH phases, whereas for N3 and N4 treatment, NNI remained substantially stable (data not shown). NNI values for N3 and N4 treatments in the FOR1 soil class did not reach the value of 1 for both the considered seasons and crops. This finding indicates a sub optimal N condition (Figure 4), whereas CDS2 and JOL2 were very often higher than 0.9. Substantial differences were observed under different soil classes, where classes FOR1 and JOL1 markedly differed in NNI values, especially in N4 treatment. FOR1 had a constant limiting factor that influenced N uptake regardless of N applications.

The boxplot of Figure 5 shows how NNI varied during the growth seasons for maize (2018) and wheat (2019). The general non-optimal nutritional condition of maize and wheat crop in the experiment, especially late in the growth season, is highlighted by many Na% values lower than the critical concentration represented by the dilution curve. Only few plots reached optimal nutritional condition with NNI values greater than 0.9. The effects of soil classes properties on NNI were statistically non-significant (p < 0.05), but FOR1 class, constantly had values unequivocally lower than 0.9 during the season independently of the fertilization rate. Instead, all the plots in classes CDS2 and JOL1 had higher values and no stress for maize (wheat) at BBCH 12 and 30 (25) was detected. The different behavior emerged in the different areas of the field regardless of the nitrogen fertilization class, suggested that the capability of the plants to absorb nutrients and therefore nitrogen depends on different properties of the soil and not only on the quantity of nitrogen that is given with fertilizers. In this case, the sub-acid pH of the FOR1 class limited more the others classes the absorption of nitrogen available in soil and provided by fertilization.

3.3. Remote Sensing Estimation of NNI

3.3.1. S2 BVs Accuracy

Before considering NNI estimation with indirect methods, S2_BVs retrieval performance was analyzed (Figure 6). From our experimental data, GAI_S2 for maize are well correlated with field estimates (R² = 0.87 and RMSE = 0.58). Lower performances are obtained for wheat (R² = 0.41 and RMSE = 0.81), however Sentinel-2 retrievals generally produced higher values than the field measurements, i.e., over-estimation (Figure 6). The obtained accuracy metrics are comparable to those reported by [39], and also in their study, the trend of GAI_S2 values showed some bias (over-estimation).

For what concerns CCC_S2 for maize, the estimations are well correlated with field measurements (R² = 0.81 and RMSE = 0.67) even though there is an overestimation at the beginning and an underestimation towards the end of the season. Worse performances were obtained for wheat, with R² = 0.38 and RMSE = 0.72. It is important to note that in 2019 for wheat, a shorter phenological period was monitored due to cloud contamination in EO data that did not allow to analyze images when field data were acquired. These results are in general in agreement with [39,63].

3.3.2. Crop Nutritional Status Estimation

Direct Method

The analyses for maize crop 2018 data (

Table 7. Coefficients of determination for the linear regressions between field measured or calculated biophysical variables (BVs) and satellite (S2) derived Vis and S2_BVs for maize crop 2018. Only statistically significant values are reported.

Variable	W	GAI DHP	GAI DHP	Cab	CCC1 a	CCC2 b	Na	NNI	PNUa
Unit	t ha−1	Effective	True	µg cm2	g m2	g m2	%		Kg ha−1
GAI_S2	0.68	0.87	0.91	0.25	0.77	0.81	0.56	0.74	0.75
CCC_S2 (g m2) c	0.65	0.84	0.88	0.28	0.77	0.81	0.46	0.76	0.77
Cab_S2 (g m2) d	0.54	0.77	0.83	0.30	0.69	0.76	0.47	0.82	0.70
GNDVI	0.56	0.78	0.86	0.26	0.67	0.74	0.58	0.77	0.66
MCARI	0.34	0.48	0.56		0.32	0.37	0.71	0.33	0.29
MCARI2	0.53	0.76	0.84	0.23	0.64	0.71	0.59	0.74	0.63
MSAVI	0.55	0.77	0.85	0.23	0.65	0.72	0.63	0.73	0.63
MTVI2	0.53	0.76	0.84	0.23	0.64	0.71	0.59	0.74	0.63
NDRE	0.51	0.75	0.83	0.26	0.64	0.71	0.56	0.79	0.64
NDVI	0.48	0.72	0.80	0.21	0.59	0.66	0.56	0.76	0.60
OSAVI	0.53	0.75	0.84	0.24	0.63	0.70	0.61	0.74	0.62
SAVI	0.56	0.78	0.86	0.24	0.66	0.73	0.61	0.73	0.65
TCARI	0.43	0.64	0.72	0.28	0.54	0.62	0.54	0.74	0.54
TCARI/OSAVI							0.47

^a CCC1 expressed as GAI_DHP_effxCab. ^b CCC2 expressed as GAI_DHP_truexCab. ^c CCC_S2 expressed as LAIxCab (SNAP BV product). ^d Cab_S2 expressed as GAI_S2/CCC_S2.

) show that there were highly significant relationships between NNI and NDRE (R² = 0.79; RMSE = 0.26; MAE = 0.2) and GNDVI (R² = 0.77; RMSE = 0.27; MAE = 0.22) (Figure 7). The highest R² and the smallest RMSE confirm the usefulness of VIs that exploit the red edge spectral region (e.g., NDRE) to provide prediction of leaf N content as mentioned in many studies [36,64,65,66]. We also found a good relationship between Na% and MCARI (R² = 0.71; RMSE = 2.48; MAE = 2.29), as previously shown for the same crop by [32], but in general, VIs achieved worse results as compared to S2_BV for retrieving BV as W, CCC, and PNUa. For this reason, these results suggest to use VIs to directly derive NNI instead of calculating W and Na% and then indirectly deriving NNI values.

An analysis of the relationships between S2 BV and nutritional status indicators was also performed (Table 7). No index reached the performance of the GAI_S2 for biomass estimation (R² = 0.68), confirming the correctness in the use of this variable to derive W. Chlorophyll related variables (CCC_S2 and Cab_S2) resulted well correlated with NNI, comparably to VIs analysis, but not so well with Na%. A better relationship was found between NNI and leaf chlorophyll Cab_S2 (R² = 0.82) where Cab is estimated as the ratio of the two variables output (S2_BVs) from the SNAP Sentinel-2 biophysical processor, i.e., CCC_S2 and GAI_S2. Being Cab_S2 a derived secondary variable, the use of CCC_S2 (R² = 0.77) in NNI estimation could be preferred to that of the leaf chlorophyll content. In general, good results were found for BV retrieval (W, CCC, PNUa) which suggests the possibility to both use direct or indirect methods to derive NNI.

For wheat crop 2019 data (

Table 8. Coefficients of determination for the relationship between field measured or calculated BVs and satellite (S2) derived Vis and S2_BVs for wheat crop 2019. Only statistically significant values are reported.

Variable	W	GAI DHP	GAI DHP	Cab	CCC1 a	CCC2 b	Na	NNI	PNUa
Unit	t ha−1	Effective	True	µg cm2	g m2	g m2	%		Kg ha−1
GAI_S2	0.36	0.42	0.39		0.37	0.35			0.33
CCC_S2 (g m2) c	0.35	0.37	0.37		0.37	0.38			0.27
Cab_S2 (g m2) d	0.28	0.21	0.27		0.30	0.38
GNDVI	0.45	0.45	0.44		0.46	0.46			0.39
MCARI		0.32	0.20		0.18	0.10	0.48	0.47	0.31
MCARI2	0.36	0.46	0.36		0.36	0.29	0.35	0.35	0.42
MSAVI	0.36	0.46	0.36		0.36	0.29	0.35	0.35	0.42
MTVI2	0.47	0.46	0.45		0.48	0.48			0.39
NDRE	0.42	0.50	0.43		0.41	0.36	0.32	0.31	0.45
NDVI	0.39	0.48	0.39		0.38	0.31	0.36	0.36	0.44
OSAVI	0.36	0.46	0.36		0.35	0.28	0.38	0.38	0.43
SAVI				0.34			0.61	0.61
TCARI				0.42			0.40	0.40
TCARI/OSAVI		0.29		0.23			0.54	0.53	0.29

^a CCC1 expressed as GAI_DHP_effxCab. ^b CCC2 expressed as GAI_DHP_truexCab. ^c CCC_S2 expressed as LAIxCab (SNAP BV product). ^d Cab_S2 expressed as GAI_S2/CCC_S2.

), only one S2 image was available during the key phenological phases and close to the field measurements, due to cloud occurrence and the lodging problems that affected the crop before flowering. The partial results for wheat were affected by dependence of the relationship with crop phenology and environmental conditions [39,67] and cannot be a valid exportable predictive model. Despite this, TCARI obtained the best performance in Na% (R² = 0.61, RMSE = 4.28) and NNI (R² = 0.61 RMSE = 0.87) retrieval for direct estimation.

Indirect Method

Indirect methods require the capacity to directly estimate different biophysical variables related to crop status. The advantage of using PNUa besides NNI indicator is because it is possible to quantitatively estimate crop N deficit/surplus [15,26] and then use this amount to make N prescription maps. This analysis was performed only on maize data because S2 data were available across the whole 2018 crop season. Figure 8 shows the relationship between GAI_S2 and W for maize, described by a logarithmic function, as expected, with good accuracy R² = 0.84 and RMSE = 2.60. The variable GAI_S2 can be used to estimate through empirical relationships W. From these estimates, PNUc can be derived using the crop specific critical canopy nitrogen content curve [16,26]. To further calculate ∆PNU, it is crucial to accurately estimate PNUa. The correlation between CCC_S2 and PNUa shows a R² of 0.77 that is in line with what achieved by [39] and [31], but with an higher RMSE (RMSE = 61.12; MAE = 41.26).

NNI estimates derived from the S2 image of 8th of July with the two methods (direct and indirect) showed similar results if compared with NNI field data. For a single image with direct method and NDRE index, we achieved a R² of 0.62 (RMSE = 0.21), while when using NNI obtained by the ratio between PNUa and PNUc, we obtained a R² of 0.63 (RMSE = 0.27).

3.3.3. NNI Map and Yield Relation

A NNI map for all the 2018 maize fields of the farm was produced and compared to the reported yield at field level. The S2 image of the 8th of July, corresponding to the end of maize vegetative phase, was used to generate an NNI map using the NDRE relation (Figure 8). The NNI map represents, at the beginning of flowering, when the fertilization period is concluded, a proxy of the potential end-of-season yield if other limiting factors are not occurring (e.g., water stress, pest, disease, or lodging).

The NNI map derived using NDRE showed the variability and the low NNI values of the crop in the experimental field (Figure 9), especially for FOR1 plots (top left of Figure 9). Anova test confirmed that field level NNI values were well correlated with the final yield and maize group (FAO class 600–700) according to nutritional conditions at flowering. The most productive fields had the higher NNI values and vice versa and the NNI map generally predicted well the final yield (

Table 9. Yield and NNI mean values (directly derived using S2 image of 08 July) of each maize field of the farm analyzed using the ANOVA method (different letters are significantly different at p < 0.05).

Maize Field ID	NNI Mean Value	Tukey’s Post-Hoc Test	Yield t ha−1
9	0.90	a	22.1
175	0.90	a	18.5
10	0.89	a	24.2
11	0.88	ab	24.1
258	0.86	bc	19.5
8	0.83	cd	17.8
308	0.83	d	18.0
312	0.82	d	13.0
368	0.78	e	19.5
Exp. Field	0.77	e	14.2
331	0.69	f	19.5
467	0.68	f	16.0
92	0.58	g	10.5
91	0.5	h	16.0

). Correlation analysis between NNI and yield revealed a significant positive correlation (Pearson’s r = 0.6).

4. Discussion

It is recognized that seasonal N status assessment can be performed using the NNI indicator [68], however, not much investigation is present in literature about the interaction between detected crop N deficiency and actual source of crop grow limiting factors. NNI can be used to highlight which zones of a field present crop anomalous conditions: NNI > 1 “luxury consumption” and NNI < 1 “nitrogen deficiency”. Usually this information is used to tactically modulate nitrogen fertilization with VRT machinery, providing the appropriate amount of fertilizer. However, in some cases, plant N deficiency is not related to nitrogen availability, but to other factor that must be identified and taken into account, to support the most appropriate fertilization strategy. The present study highlighted that NNI is an appropriate tool for the diagnosis of the nutritional status of the crops, however, soil properties can have an impact on crop N uptake, showing an effect independent from N rates supply. The capacity of NNI to assess crop nutritional status in relation to soil properties was confirmed for two consecutive crop seasons. In our case, measured NNI showed that the experimental plots located in an area of the field with acidic soil (FOR1) had lower NNI values than all the other plots in the other soils’ types both for maize and wheat. High N rates did not increase Na% and consequently crops biomass in the problematic soil area of the field (FOR1), as expected and as demonstrated by other studies [16,18,29]. In addition to this, the high rate plots in the FOR1 class reached only half of the biomass accumulation (4 t ha⁻¹) as compared to the best performance of JOL1 and CDS2 soil classes (8 t ha⁻¹) and the two classes are statistically different considering this variable. It is well known that most plant nutrients are optimally available to plants within the soil pH range between 6.5 and 7.5 and soil classes of the field experiment belonged to this range except FOR 1. Nutrients absorption depends on elements availability in the root zone and the growth rate potential, even in the N excessive condition, is regulated by the soil N supply [30]. N fertilization can temporarily reduce soil pH in the surface soil layer after the supply, especially at high amounts as demonstrated by [69]. In presence of acidic soils, this effect may be even higher, causing even greater problems in nutrients absorption, with a combination of yield limiting factors. The study of [11] underlined that nitrogen assimilation needs to be considered in the context of inter-regulation of multiple crop physiological processes (i.e., C N assimilation and allocation) and soil N availability.

Besides this consideration, overall, the estimated NNI values during the season have revealed areas with a general sub-optimal crop growth in the experimental field, which resulted in an average yield (14.18 t ha⁻¹) much lower than the average of all the maize fields of the farm (18.22 t ha⁻¹), confirming the goodness of NNI to indicate crop yield potential [14]. NNI calculated from remote sensing during the season resulted in a dynamic diagnostic tool to assess spatial explicit crop N status and limitation in N uptake [14]. NNI calculated from remote sensing during the season resulted in a dynamic diagnostic tool to assess spatial explicit crop N status and limitation in N uptake [14]. Chlorophyll content estimation for crop stress status monitoring has been considered in many studies and for site-specific N management [66,70], but the relationship between chlorophyll content and Na% remains highly empirical [26,71] because it is usually influenced by variety, phenology, and study area, limiting the generalization capability of the method using only chlorophyll content information [39,50]. RS DM are proposed as a rapid and robust way for nutritional status assessment to be used as an alternative to direct field measurements [29,31]. Information on actual nitrogen content or concentration, or direct indication of nitrogen deficiency (NNI) from remote data are of great help in crop monitoring.

From the 2018 maize experiment, we found that the direct estimation of NNI using VIs or S2_BVs derived by Sentinel-2 is feasible as demonstrated by general good correlation metrics comparable to other studies. The best performance was obtained using Cab_S2 secondary variable (R² = 0.82) slightly better than using CCC_S2 (R² = 0.77) or NDRE (R² = 0.79), this performance is similar to what was reported by [65]. Our results indicated also that Na% can be directly derived using MCARI with an agreement of R² = 0.71, confirming the conclusions of [32]. However, biomass estimation directly from VIs resulted not reliable. So, the direct “mechanistic approach” (DM_2) was not performing as expected. These results confirmed the possibility to directly derive NNI (DM_1) with VIs and/or S2_BVs as a suitable solution for operational crop monitoring to support site-specific nitrogen fertilization. The NNI derived maps showed comparable results in terms of correlation with DM (R² = 0.62, RMSE = 0.21) and IM (R² = 0.63, RMSE = 0.27) approaches.

Indirect method (IM) derive NNI from the separate estimations of W and PNUa. The advantage of using such approach is that estimated PNUa can be used to compute ∆PNU, this is a promising quantitative approach to define nitrogen absorption deficit of a crop [14,15,26]. PNUA (kg ha⁻¹) is a quantitative variable preferred over a variable expressed in percentage as Na because it avoids an intermediate step to convert N contents from a mass basis to an area basis to a mass basis [72] as reported by [14]. Experimental results demonstrated that the biophysical variable CCC estimated from S2 is well related with PNUa (R² = 0.78), confirming the results of previous studies [31,39]. S2 estimated GAI can be used to estimate W (R² = 0.84), since this variable is necessary to calculate PNUc according to crop specific dilution curves. The quantitative deficit or surplus values given by the indicator ∆PNU is thus suitable for adjusting N supply in site-specific fertilization management.

5. Conclusions

The present study highlighted that soil types have an effect on the diagnosis of the nutritional status of the crops in response to N rates supply. Sentinel-2 data can be exploited to estimate nutritional deficit or surplus computing NNI with direct and indirect methods. In addition to this indicator using indirect methods, it is possible to quantitatively defined the required N amount to supply with VRT machinery by calculating ∆PNU. However, this information is not sufficient and it is necessary to take in account the potential causes of crop nutritional status such as soil properties variability. In the experiment presented here, a limiting factor for optimal crop growth was identified in the soil properties (pH) that reduce the capacity of the crop to uptake N, also in conditions of high N availability in the soil and supplied by fertilization. This influences the interpretation of N status, because even though a nitrogen deficit is detected (i.e., NNI < 0.9) a supplementary provision of fertilizer (∆PNU) would not be efficient. The experiment confirmed that different soils types within a cultivated field can have an effect on growth independently of N rates provided and consequently impact on dry matter production and final yield. Indeed, in the location of a field where yield potential is low, added N fertilizer profitability can be reduced as demonstrated, for example by [73]. The results confirm that, since crop N absorption is a very dynamic and complex process, site-specific management in a precision farming paradigm requires a combination of spatial explicit information: static maps describing variability of soil properties and seasonal information related to actual crop nutritional status achievable from RS data. The information from soil characterization is potentially complementary to crop N status monitoring, and could improve NUE and site-specific N fertilization strategies of cereals cropping systems. Further experiments and tests are needed in different study areas and exploiting data acquired on different seasons to deeply investigate the effect of different soil properties and N supply on crop N deficiency diagnosis and to explore the possibility of using soil and crop derived information in combination.