Use of Multispectral Airborne Images to Improve In-Season Nitrogen Management, Predict Grain Yield and Estimate Economic Return of Maize in Irrigated High Yielding Environments

: Vegetation indices (VIs) derived from active or passive sensors have been used for maize growth monitoring and real-time nitrogen (N) management at ﬁeld scale. In the present multilocation two-year study, multispectral VIs (green- and red-based), chlorophyll meter (SPAD) and plant height (PltH) measured at V12–VT stage of maize development, were used to distinguish among the N status of maize, to predict grain yield and economic return in high yielding environments. Moreover, linear plateau-models were performed with VIs, SPAD and PltH measurements to determine the amount of N needed to achieve maximum maize grain yields and economic return. The available N in the topsoil (0–30 cm) was measured, and its relationship with VIs, maize yield and maize N requirements was analyzed. Green-based VIs were the most accurate indices to predict grain yield and to estimate the grain yield optimum N rate (GYON r ) (216.8 kg N ha − 1 ), but underestimated the grain yield optimum N available (GYON a ) (248.6 kg N ha − 1 ). Red-based VIs slightly overestimated the GYON r and GYON a , while SPAD highly underestimated both of them. The determination of the available N did not improve the accuracy of the VIs to determine the grain yield. The green chlorophyll index (GCI) distinguished maize that would yield less than 84% of the maximum yield, showing a high potential to detect and correct maize N deﬁciencies at V12 stage. The economic optimum nitrogen rate (EON r ) and economic optimum nitrogen available (EON a ) were determined below the GYON r and the GYON a , demonstrating that maximum grain yield strategies in maize are not normally the most proﬁtable for farmers. Further research is needed to ﬁne-tune the response of maize to N applications when deﬁciencies are detected at V12 stage, but airborne imagery could be useful for practical farming implementation in irrigated high yielding environments.


Introduction
Productivity and resource-use efficiency are desirable agronomic, economic and environmental goals [1,2]. In most crops, nitrogen (N) is the most limiting nutrient for crop production in many of the world's agricultural areas [3], and its efficient use is important for economic and environmental sustainability [4]. Maize (Zea mays L.) is a high N-demanding crop, where insufficient application of N can result in important economic losses, whereas an excessive fertilization implies wasting resources and an increase of the environmental pollution [5]. Nitrogen fertilization of maize can account for up to 30% of the total production cost [6]. Hence, it is important to optimize its application to maximize the economic return.
High-yielding environments, such as those of the Ebro Valley (northeastern Spain), where maize could achieve more than 16 Mg ha −1 of grain, usually require high N rates (N rates ) to cover crop necessities [7]. The N rates are the amount of N fertilizer applied to corn during the growing season. In this area, maize N extractions are about 250-300 kg ha −1 [8][9][10]. Farmers apply high N rates to cover maize N requirements and therefore, the N fertilization strategy has a great influence in the economic return. Multiple assessment tools for N have been used for N management. However, temporal and spatial variation of crop N requirements [11], together with within-field soil N supply variability, can influence the assessment of N-field availability [12]. Nitrogen fertilization of crops, such as maize, has been mostly underpinned by the simplification of fieldwork management and by avoiding under-fertilization risk. Conducting farms to the common practice of over-fertilization by the application of high N rates at early stages of the crop.
The worldwide N recovery in crops is usually less than 50% [3], with the impact that this supposes on N resource efficiency and the pollution of agroecosystems [13]. Knowing that approximately one-third of the total maize N uptake occurs after pollination (under favorable soil moisture conditions) [14], better N practices should be developed to increase N use efficiency.
Determination of within-field soil spatial patterns is a useful tool for N management but requires the collection and analysis of a large number of samples. Soil [15], plant or chlorophyll concentration [16] sampling is costly, time-consuming and barely captures the spatial variability often present within fields. Plant analysis and chlorophyll meter readings have been used for confirmation of N responsive sites, whereas soil tests for available N (N available ) have predicted successfully the crop N needs in arid regions, but not in humid regions [17]. The N available is the sum of the N applied (N rate ) plus the NO 3 − -N present in the soil at sidedress time, which gives an estimate of the N that is going to be accessible for the crop during the growing season. Therefore, there is a need to develop a faster, more accurate and possibly more economical method to gather crop information and to estimate and adjust N requirements [18]. Due to the link between net photosynthesis and steady-state fluorescence at airborne spectral level, field-spectral imaging is a very useful technology for crop growth monitoring and real-time management at field scale [19][20][21]. Vegetation indices (VIs) derived from active or passive sensors have been used to distinguish temporal patterns in crop development [22][23][24][25][26]. For example, VIs have been used for detecting N stress in maize at early development stages (V4-V7 stage, 4-7 leaves with visible leaf collar), though with relatively low accuracy [27][28][29]. This was attributed to the strong background soil reflectance. Hence, maize ground cover is vital if the pixel resolution does not allow removal of soil pixels [30]. However, image acquisition after V8 stage seems to be consistently useful to determine maize N status and to predict yield [23,[31][32][33][34][35], increasing progressively its accuracy up to VT stage (tasseling). In areas where N can be applied via irrigation systems, the image acquisition at V12 stage is generally preferred to other alternatives. At V12, maize has already taken up about 40% of the total N [36] and crop response to N fertilizer is still high if N deficiencies are detected. Image acquisition after VT stage could be affected by color disturbance of the tassels, reducing the correlation between VIs and N status or yield [37]. Moreover, after VT, the maize response to extra N application is further limited because the crop has less time to absorb the applied N (compared with V8 or V12 stage) [38].
Despite the usefulness of multispectral images to aid N management, the adoption of remote sensing technologies by farmers is still limited [39]. Different technologies have been used for this purpose. Satellite images normally have lower spatial and temporal resolution and, in some areas, they can be disturbed by cloud cover and/or sprinkler irrigation during the image acquisition [40]. Unmanned aerial vehicles (UAV) have tremendous potential for high-resolution requirements for detailed site-specific weed control treatments in early post-emergence [41] or N recommendation [35]. However, such high-resolution is probably not necessary for determining maize N status at V12 (when the canopy covers the ground entirely), and airborne images could be useful for practical farming implementation because of their high potential to cover large cropping areas.
While numerous studies have been conducted using remote sensing techniques to determine N status and predict maize grain yield [19,23,[31][32][33][34], limited research have been made comparing the effectiveness of these technologies in irrigated high-yielding environments.
The objectives of the present study were to use VIs derived from multispectral aerial images: (1) to distinguish among N status of maize at V12 stage, (2) to predict grain yield and economic return, and (3) to determine the amount of N needed to achieve maximum maize grain yields and economic return in high yielding environments.  The experimental fields comprised an area of 108 × 60 m 2 and 110 × 130 m 2 with plot sizes of 180 m 2 and 150 m 2 in AL and GI, respectively. However, in the present study only 20 of the 36 and 49 experimental plots were analyzed for AL and GI, respectively. These 20 plots were selected because they had the same N treatments. It should be noted that plots had been continually used for the same purpose throughout the previous 13 and 6 years in GI and AL, respectively. Same N fertilization treatments were maintained over the aforementioned period, so an emphasizing effect of the different N treatments could be expected.

Study Area
The study area is characterized by a semiarid climate with low precipitation (192 mm) and high average temperature (19.1 • C) during the maize growing period [7]. The summer of 2015 (June-August period) was hotter than the same period of 2014 (1.0, 1.5 and 2.1 • C for minimum, mean and maximum temperatures, respectively). Each experimental field received 750-800 mm of water (lacking nitrate) per growing season, in order to provide the maize with optimum yield conditions. The soils described in the experimental fields were Typic Calcixerept (AL) and Petrocalcic Calcixerept (GI) [42], in both cases well-drained. Soil quality indicators and physicochemical parameters were analyzed using standard methods [43]: soil texture, pH, electrical conductivity (EC), bulk density, available P (Olsen P) and extractable K (NH 4 Ac) ( Table 1). In both locations (AL, GI) and years (2014, 2015) conventional tillage was applied, including disc ploughing and cultivation to a depth of 30 cm to incorporate maize stover and to prepare the soil for the next sowing. In early April, the hybrid PR33Y72 (FAO cycle 600) was planted in both fields at a rate of 90,000 plants ha −1 with 71 cm between rows. Two herbicides treatments were applied: one at pre-emergence to control the majority of weeds (S-Metolachlor 40% and Terbuthylazine 18.75%, at 3 L ha −1 ) and the other at post-emergence to specifically control Abutilon theophrasti M. and Sorghum halepense L. (Dimethylamine salt of dicamba 48.2% at 1 L ha −1 and Nicosulfuron 6% at 0.75 L ha −1 ).
Maize was harvested the last week of October with an experimental small-plot combine. Grain yield was measured in the two central rows (1.5 × 10 m) of each plot and subsequently, the it was adjusted to 14% moisture content by the use of a 250 g sample in each plot (GAC II, Dickey-John, Auburn, IL, USA).
The economic return of each plot was calculated as the difference between the income produced by the selling of the grain yield and the cost of the N fertilizer applied. The N:maize price ratio is defined as the price per kilogram of N divided by the price per kilogram of maize (price ratio = price of fertilizer N, € kg −1 N/price of maize, € kg −1 maize) [44]. In the present study, the N:maize price ratio was 5.3:1, considering a N price (N fertilizer plus application cost) of 0.90 € kg -1 and a maize grain price of 0.17 € kg -1 .

Fertilizer Treatments
Five different inorganic fertilizer N rates : 0, 100, 200, 300 and 400 kg N ha −1 (N0, N100, N200, N300 and N400, respectively), with four replications under a randomized block design were considered. The N fertilization treatments were randomized at the beginning of the experiments (2002 and 2010, in GI and AL, respectively) and applied in the same plots the following seasons. The N fertilizer (34.5% N, as ammonium nitrate) was manually applied and split into two equal sidedress applications (50% at V3-V4 and 50% at V6-V7 stage) for all plots except for N400 treatment, where a third sidedress was applied at V10 stage. Thus, the N distribution in the N400 treatment was 37.5% at V3-V4, 37.5% at V6-V7 and 25% at V10 stage. The N400 treatment had different N management to increase N use efficiency (NUE) and to reduce the risk of pollution by nitrate leaching associated with high N rates at early stages. Phosphorus and potassium were also applied every two years in both locations before planting, at rates of 150 kg P 2 O 5 ha −1 and 250 kg K 2 O ha −1 , to avoid deficiencies of those elements.

Remote and Proximal Sensing: Data Acquisition and Analysis
A digital multi-spectral camera (DMSC) with a DMSC-2k System sensor (Specterra-Services, Leederville, WA, Australia), mounted on an airplane (operated by RS Servicios de Teledetección, Lleida, Spain), was used to acquire multispectral aerial images of the experimental fields at V12 stage. The surveys were conducted under optimum flight conditions at 880 m (AGL) above both fields on 30th June 2014 and 25th June 2015. The time of flight was less than 1 h between 12:30 and 1:30 h (GTM) on sunny days, without clouds. The DMSC-2k System sensor consisted of a 4-interline transfer, a 2048 × 2048 pixel charge-coupled device (CCD) with a Nikon F mount, and 24-28 mm fixed focal length lenses. The spatial and radiometric resolutions were 0.25 m and 14 bit (recorded as 16 bit), respectively. The spectral resolution of the camera consisted in four independent and replaceable narrow bandwidth spectral filters. Twenty nm range width bands centred at 450 nm (blue), 550 nm (green), 675 nm (red) and 780 nm (near infrared) were captured. The spectral bands were pre-processed by the provider to compensate for mis-registration due to lens distortion (less than 0.2 pixels) and for scene brightness due to the bi-directional reflectance distribution function (BRDF).
Seven different VIs reported in the literature (NDVI, GNDVI, GCI, SAVI, GSAVI, WDRVI and EVI) were computed (Table 2). Then, to summarize the response of the different N treatments to each VI, fifty points within each individual experimental plot (excluding 1 m buffers from the borders), were sampled and the VI values extracted. Soil distortion did not affect the VIs because soil pixels were avoided in the analysis.  [45] (2) GNDVI (Green NDVI) NIR−Green NIR+Green [46] (3) GCI (Green Chlorophyll Index) NIR Green − 1 [47,48] (4) ‡ SAVI (Soil Adjusted Vegetation Index)

NIR−Red
(NIR+Red+L) ·(1 + L) [49] (5) ‡ GSAVI (Green SAVI) [50] . † Green, Red, Blue and NIR are the reflectance of the green, red, blue and near infrared light, respectively. L is a correction factor (calculated by the formula presented at the bottom of the table), where s is the slope of the soil line. α is a weighting coefficient (α = 0.2). In the EVI calculation, G is a gain factor, C 1 , C 2 are the coefficients of the aerosol resistance term and L 2 functions as the soil-adjustment factor (L 2 = 1; C 1 = 6; C 2 = 7.5 and G = 2.5).
In addition to the VIs, the heights of five plants were manually measured in each experimental plot at VT stage (10 days after image acquisition). Plant height (PltH) was measured at the last leaf to avoid differences caused by tassel size and shape, which vary with both plant population and variety [52]. Non-destructive chlorophyll readings were taken at VT using a small, lightweight, portable, hand-held meter (SPAD-502, indirect chlorophyll meter; MinoltaCorp, Ramsey, NJ, USA). In agriculture, the SPAD meter is widely accepted to improve N management by predicting N status and determining fertilization requirements [53]. SPAD values indirectly evaluate the leaf chlorophyll content based on the amount of light transmitted by the leaves at two different wavelengths: red (650 nm) and near infrared (940 nm). In the present study, SPAD values were determined by sampling three different parts of the ear leaf (base, middle and top) in five central row plants of each experimental plot (15 measurements per plot). Then, SPAD values were averaged per each plot.
Soil NO 3 − -N was determined in each plot before planting at a depth of 0-30 cm by a composite of five individual cores. Soil NO 3 − -N was extracted using deionized water and measured using test strips with a Nitrachek ® device calibrated according to the standard procedure [54]. The N available for maize during the growing season was determined by the sum of the N rate applied in the experimental plots and the NO 3 − -N determined at sidedress: N available = N soil at sidedress + N rate .
In order to compare the results of the different N rates and N available among fields and years, indices ratios were calculated in each field and year by dividing the plot average value (VIs, SPAD, PltH, grain yield and economic return) by the maximum value of each field and year. Normally, the highest values were obtained in the considered overfertilized treatment (400 kg N ha −1 ). Furthermore, the calculated ratios of VIs, SPAD and PltH were fitted to linear-plateau models with the N rates , N available , and grain yield to determine the accuracy of these variables to determine maize N status and to predict grain yield.

Statistical Analysis
The experiment was analyzed as a split-plot in time with completely randomized blocks and four replications. The location (GI or AL) was the main plot, and the N rates (N0, N100, N200, N300 and N400) were subplots. The ratios calculated from the extracted values in the sampling points were subjected to analysis of variance using the Mixed Model of the Statistical Analysis System (JMP Pro 12, SAS Institute, Cary, NC, USA), considering N rates , the year and the location as fixed factors and the replication as random effect. In addition, the grain yield and economic return means were separated by the LSMeans Tukey HSD test (p < 0.05), with levels not connected by the same letter considered significantly different.
Linear-plateau regression analyses were also carried out: (1) between the ratios of the VIs, SPAD, PltH and grain yield with the N rates or N available to estimate the grain yield optimum N rate (GYON r ) and the grain yield optimum N available (GYON a ) [55], (2) between the ratios of the economic return with the N rates or N available to determine the economic optimum nitrogen rate (EON r ) and the economic optimum nitrogen available (EON a ), and (3) between the ratios of the VIs, SPAD and PltH with the ratios of grain yield to determine their usefulness to predict grain yield at the V12 stage. Analysing together locations and years, the GYON r , the GYON a , the EON r and the EON a were determined where the plateau was reached, indicating the non-responsiveness amount of N (saturation point) to extra N application. Table 3 presents the mean values and the ANOVA of the VIs at V12 stage (n = 50 for each plot), PltH and SPAD at VT stage (n = 15 for each plot), grain yield and economic return (n = 4) for the different N rates applied in the experimental plots. The rates are sorted according to the total amount of applied N.

Vegetation Indices, SPAD and Plant Height Responses to N Rates and Available N
As expected, the control treatment (0 kg N ha −1 ) presented the lowest VIs, SPAD and PltH values in both locations (p-value < 0.05). The differences were more accentuated in GI than in AL. Probably, the longer period of time without any N application (13 years) and the lower soil organic matter (OM) content, reduced the soil capacity to provide N to the crop during the growing season. In GI, the VIs and PltH values of the N0 treatment were higher in 2015 than in 2014. The higher summer temperatures in 2015 could have contributed to the mineralization of the OM [56,57] and consequently, to N availability.
The N100 treatment presented significantly higher VIs, SPAD and PltH values than the control treatment (N0), but lower than the N rates above 200 kg N ha −1 . In fact, there were very little and non-significant differences for most VIs, SPAD and PltH at N rates between 200-400 kg N ha −1 .
Previous studies in the area determined the need of 290 kg N available ha −1 (0-90 cm) to achieve maximum maize grain yields [8], with extractions up to 386 kg N ha −1 [58]. Then, the N400, and probably the N300 treatment, could provide enough N to cover maize's N requirements. The N200 seemed a non-sufficient amount of N to cover maize N extractions. However, the N200 rate, together with the initial soil N content at planting and the N mineralization through the growing season, seemed to provide enough N for maize development. At least until V12-VT, when maize has already absorbed around 40-50% of the total N [14].
All VIs, SPAD and PltH were able to differentiate among different N rates of the experiment (p-value < 0.01), showing their usefulness for differentiating N maize status at V12 stage [37,59,60]. The interaction between location and N rates (L × N) was significant, probably due to the higher values of the indices observed in AL with respect to GI in the N0 and N100 treatment. These differences between locations could be due to field conditions that contribute to soil fertility [61], as well as to different soil OM contents in AL than in GI, which would have produced higher N mineralization in Al than GI ( Table 1). As mentioned above, the difference in the quantity of N mineralization in the two fields was evident and can be observed with the grain yields obtained in the control treatments (N0). The yields of AL (9.0 Mg ha −1 ) nearly doubled those of GI (4.8 Mg ha −1 ). The higher yields in 2015 provoked a significant interaction between year and N rates (Y × N).
The VIs, PltH and SPAD were highly correlated with the N rates applied ( Figure 2) and with the N available ( Figure 3). All indices were fitted to a linear-plateau model when compared with the applied N rates or N available . The GYON r and GYON a were determined when the plateau was reached (saturation point), indicating where the VIs, SPAD and PltH did not increase with additional N inputs.
Green-based VIs obtained better correlations than red-based VIs with the N rates or N available (Figures 2 and 3). Indeed, green-based VIs have been normally considered more useful for assessing leaf chlorophyll variability when the leaf area index is moderately high [62], as it is in maize. Sripada et al. [18] showed the usefulness of green-based VIs at predicting the GYON r (R 2 = 0.67) with aerial CIR (colour infrared) imagery. Bausch and Khosla [32] used the GNDVI to determine the GYON r (R 2 = 0.91) with QuickBird panchromatic imagery (satellite). In our study, the highest correlations at V12 between N rates and VIs were found for GCI (R 2 = 0.80), GSAVI and GNDVI (R 2 = 0.73). Otherwise, at VT stage the SPAD correlated up to R 2 = 0.89. Red-based VIs, as well as PltH, showed low correlation with the N rates . Indeed, the R 2 coefficient was 0.44 for PltH, 0.51 for NDVI and 0.45 for EVI, SAVI, WDRVI ( Figure 2).
The response of the VIs, SPAD and PltH to the N available was quite similar to that observed for the N rates . However, a consistent increase in the R 2 coefficient when fitting linear-plateau models was observed due to the extra information (initial soil NO 3 − -N) added to the model (Figure 3). The use of the N available explained better some variations of the indices values for the same N rates applied. Thus, in-season nitrate analysis can be used for the fine-tuning of remote sensing prediction of N sufficiency rates [63]. The GYON r and the GYON a undoubtedly gives interesting information for N management since they allow to differentiate between N-deficient and N-sufficient plots [23]. All VIs and PltH achieved saturation between applied N rates of 202.9 and 235.0 kg N ha −1 . However, the SPAD was not able to differentiate above 142.1 kg N ha −1 . This findings are similar to the reported by Liu et al. [64] in a 5 N rates study at lower yielding conditions. In their study, SPAD was able to distinguish between N rates up to 135 kg N ha −1 . In high-yielding environments, SPAD has shown saturation at 173 kg N ha −1 , which is consider not enough to achieve maximum yields [65]. Therefore, the usefulness of SPAD as a predictor of maize N status in high-yielding environments is probably limited.  Generally, the GYON r determined by the green-based VIs (≈210 kg N ha −1 ) was lower than that determined by the red-based VIs (≈230 kg N ha −1 ). The GYON a was similar to the GYON r for the green-based VIs, however, the GYON a predicted with red-based VIs increased to 250-270 kg N available ha −1 . SPAD increased the GYON a by 27% compared with its determined GYON r , but even so, the GYON a was highly underestimated. The relationship between N available and the studied indices could also help to explain the better N status seen in 2015 compared to 2014 in GI, when higher indices values in the N0 treatment were correlated with higher N available (Figure 3).
In our conditions, the GYON a determined by the VIs for maize production was between 200 and 270 kg N ha −1 . These amounts of available N were similar to those reported in other studies in the Ebro Valley [7,8,66,67]. Thus, the tested VIs enabled the spatial characterization of N status and were able to assess crop performance even under low-N stress [68].

Grain Yield and Economic Return Responses to N Rates and Available N
Grain yield and economic return varied significantly according to N rates (p-value < 0.01), location (p-value < 0.01) and year (p-value < 0.05) ( Table 4). Moreover, interactions between location and N rates (L × N), year and location (Y × L), and year and N rates (Y × N) were detected and can be explained by the same reasons described above.
Averaging both years, there were not significant yield differences among the highest N rates (200-400 kg N ha −1 ), however, AL showed high average yields than GI (Table 4).
Consistently with high grain yields, higher economic returns were found at high N rates . The maximum economic return was obtained in an N200 treatment in AL (3068 € ha −1 ) and in an N300 treatment in GI (2578 € ha −1 ) but there were no statistical differences among the treatments of 200, 300 and 400 kg N ha −1 . The economic return was mostly determined by grain yield at the studied N:Maize price ratio of (5.3:1), so a slight reduction in grain yield was drastically penalized by the economic return. It was evident that, at higher N:Maize price ratios (worse relation of prices for farmers), the efficiency of the N will affect more the economic return of the crop. Sripada et al. [44] tested an interval of price ratios from 4:1 to 14:1 (N:Maize), so the 5.3:1 price ratio used in our study was close to the optimum for farmers and highly affected by the grain yield.
The economic return is important not only for N management in maize production system, but also for the reduction in N contamination of groundwater [11]. As EON r for maize is usually consistent with good environmental stewardship, it could therefore be used as a tool to determine maize N requirements [44], and could be considered essential for responsible N management of maize crops.
In the present study, the grain yield and economic return were highly correlated with the N rates (R 2 = 0.82 and 0.77, respectively) and N available (R 2 = 0.89 and 0.87, respectively) in linear-plateau models ( Figure 4). As observed for the VIs, knowledge of the N available in the topsoil layer slightly increased the accuracy of the model. The GYON r and GYON a were respectively determined at 216.8 and 248.6 kg N ha −1 , similar to the N rates reported in previous studies in the area [8,67,69,70]. However, the maize N extractions were probably not completely covered by the determined GYON r and GYON a in the linear-plateau model. Figure 4. Response of (a) yield to nitrogen rates (N rates ), (b) yield to available nitrogen (N available ), (c) economic return to N rates and (d) economic return to N available as tested and determined in the experiments. Saturation N dose was determined when there was no response to higher N fertilization. The results of Almacelles (AL) and Gimenells (GI) are presented as ratios dividing the yield and economic return values by the maximum observed in each year. 'Both' indicates that the two fields were analyzed together.
In similar environments, Maresma et al. [35] and Yagüe and Quílez [70] reported respectively a GYON r of 239.8 and 300 kg N ha −1 , where the N extraction of maize is around 250-300 kg N ha −1 [8][9][10]. These N rates can be considered low compared to the extractions, but if N is provided when needed, maize could uptake the N and translates it into yield. Moreover, Biau et al. [58] quantified the maize N uptake up to 386 kg N ha −1 for grain yields of 16 Mg ha −1 , in similar conditions. Probably, the difference between applied N and N uptake was supplied by soil OM mineralization during the growing season, which has been quantified in around 100 kg N ha −1 [8] for our conditions. The EON r was lower than GYON r , showing the non-economic justification of the maximum yield strategies. The maize yield per kg of applied N is reduced when increasing the N rates and, consequently, the last yield increment normally entails a decrease in crop profitability. This is because the income obtained per kg of grain is lower than the cost of the N fertilizer input. Di Paolo and Rinaldi [71] reported a reduction of 25 kg of maize per kg of N fertilizer when applying high amounts of N fertilizer (300 kg N ha −1 ) in irrigated Mediterranean conditions (14 Mg maize grain ha −1 ). At a N:Maize price ratio of 5.3:1, and considering both site location and year, the EON r was 176.6 kg N ha −1 while the EON a was 209.4 kg N ha −1 . Although EON r could vary depending on the prices of maize and N [44], Schlegel et al. [72] concluded that in same fields, the EON r was relatively insensitive to changes in these prices. In their study, the EON r was similar for low-, medium-, and high-yielding years, concluding that application of an extra N as insurance reduced crop profitability. Therefore, as EON r and EON a could be considered stable in same fields, the EON r (176.6 kg N ha −1 ) and the EON a (216.8 kg N ha −1 ) reported in this study could be useful for the guidance of the N management strategies in maize fields located in high yielding environments.

Correlation between Vegetation Indices, SPAD and Plant Height with Grain Yield
The calculated ratios of VIs at V12, SPAD and PltH at VT stage showed moderate to high correlations with the grain yield ratios (R 2 between 0.45 and 0.90) ( Figure 5). Several studies have proven the usefulness of red-based VIs, as NDVI, for predicting grain yield in different conditions [18,68,[73][74][75][76][77]. However, in the present study the green-based VIs showed higher correlation with grain yield than the red-based VIs. In a two-year study and under similar conditions, Bausch et al. [31] demonstrated the usefulness of green-based indices at V12 and later growth stages for predicting grain yield (R 2 = 0.81) and its variability within a field for in-season N management. Isla et al. [33] also reported better grain yield predictions with GNDVI and SPAD (R 2 = 0.93) than with red-based VIs (R 2 = 0.65).
The higher indices values determined at V12 stage resulted in higher grain yields. However, there was variability among the highest indices values that was not clearly translated into grain response ( Figure 5). Linear-plateau models were fitted to identify the saturation of the indices for predicting grain yield. This saturation provided information about the percentage of maximum grain yield up to which the index can differentiate. Therefore, the indices that presented saturation at higher grain yields were determined to be more useful due to their capacity to distinguish between nearly optimal maize N statuses. In this case, green-based VIs performed better than red-based VIs in determining maize vigor. This fact could be explained by problems of saturation associated with red-based indices for some types of vegetation during their later growth stages [33]. The highest saturation of VIs for grain yield prediction corresponded to SAVI (at 86% of maximum yield). The usefulness for predicting leaf chlorophyll content and N status in maize in different growth stages of VIs derived from SAVI has been reported [76,78]. However, its low correlation with grain yield (R 2 = 0.48) demonstrated that its usefulness for grain yield prediction is still uncertain.The best VI to predict grain yield was GCI, which was able to differentiate up to 84% of the maximum yield and showed high correlation with the linear-plateau model (R 2 = 0.87). Hence, small N deficiencies in maize at V12 stage could be detected by GCI and corrected with extra N application.
It is known that intra-field yield spatial variability can be occasioned by different stress problems, resulting in a grain yield variability of 15-38% with respect to non-stress conditions [79,80]. Therefore, under N deficits these differences would probably be higher. Detection of maize N status that would yield less than 84% of maximum yield would improve the average maize grain yields and, consequently, N management. Similarly, the other green-based VIs tested at V12 showed similar trends to GCI. GNDVI and GSAVI saturated at 80% and 79% of the maximum grain yield, respectively.
High variability was observed among the red-based VIs. The NDVI predicted grain yield accurately (R 2 = 0.73), but was able to differentiate only up to 64% of the maximum yield. WDRVI and EVI were respectively saturated at 45% and 58% of the maximum yield, and, together with the PltH (56%), were the least accurate at differentiating among the highest grain yields. These results partially disagree with those presented by Maresma et al. [35], who in a one year-site study found the WDRVI to be the best grain predictor using a high-resolution UAV service. Nevertheless, although in the present study the plateau was not reached for high relative grain yield, there was a clear tendency of increasing WDRVI values when increasing the yield (Figure 5f). This trend might corroborate the potential of this index to distinguish yield at high N rates with high resolution imagery.
SPAD and PltH measured at VT proved their potential to predict grain yield. SPAD differentiated up to 86% of the maximum yield, confirming the usefulness of chlorophyll meters for predicting grain yield [37,53,81]. However, PltH was partially useful for grain yield prediction (differentiating up to 56% of the maximum grain yield). Although Liu and Wiatrak [64] did not find a correlation between PltH and grain yield at V8 stage (R 2 = 0.03), our results do concur with Yin et al. [82], who predicted grain yield with PltH at VT stage (R 2 between 0.52 and 0.86, depending on the year).

Vegetation Indices, SPAD and Plant Height to Predict GYON r , GYON a , EON r and EON a
The studied indices were useful determining the GYON r (216.8 kg N ha −1 - Figure 4). Linear-plateau models (R 2 = 0.73-0.80) showed saturation of green-based VIs at rates of 210 kg N ha −1 , which makes them the best predictors of GYON r at V12 stage ( Figure 2). The red-based VIs (at V12 stage) or PltH (at VT stage) also predicted the GYON r (206-235 kg N ha −1 ) but the lower correlation with the linear-plateau model (R 2 = 0.44-0.51) causes uncertainty. The SPAD widely underestimated (142.1 kg N ha −1 ) the GYON r , agreeing with the results of Maresma et al. [35]. Thus, despite the potential of SPAD to predict yield, it should not be used to estimate the GYON r in high yielding environments.
The GYON a (in the first 30 cm of soil) was 248.6 kg N ha −1 (Figure 4). The red-based VIs predicted accurately the GYON a (252-271 kg N ha −1 ) ( Figure 3) and slightly increased their correlation with the linear-plateau model (R 2 = 0.48-0.58) compared with the GYON r . However, the effectiveness of green-based VIs, WDRVI and PltH at determining the GYON a was negatively impacted by adding the N available to the model. SPAD, as was observed for the GYON r , highly underestimated the GYON a (181.9 kg N ha −1 ).
The EON r (and EON a ) were determined below the GYON r (and GYON a ). Therefore, since the VIs proved their usefulness for determining the GYON r and GYON a , they would over-estimate the EON r (and EON a ). However, at low N:maize price ratios, EON r (and EON a ) will be more accurately predicted by the indices since the economic return is mainly affected by the grain yield.

Conclusions
Vegetation indices derived from multispectral airborne images showed their usefulness to distinguish maize N status and yield potential in high yielding environments. The VIs, SPAD and PltH were useful to predict maize grain yield (R 2 between 0.48-0.90). Similarly, VIs have proven to be an outstanding tool to accurately differentiate between N responsive and N nonresponsive areas in the field at V12 stage, when N deficits can still be corrected by late sidedress applications. Green-based VIs, and especially GCI, were the most accurate at predicting the GYON r , whereas red-based VIs slightly overestimated it and SPAD highly underestimated it. The EON r was overestimated by the VIs, confirming their higher potential to predict grain yield than economic return.
The determination of the N available did not improve the accuracy of the VIs to determine the GYON a . However, in our study N fertilization treatments were maintained during a long-term period and probably the residual N in the soil was stabilized for each N treatment. Commercial fields might have larger differences in N patterns within fields, and N available calculation could be very useful to determine N responsive and N non-responsive areas.
Overall, VIs could contribute to enhancing N management in farming practices for maize cropping. The use of VIs could reduce the risk of contamination by over-fertilization, improve the economic return, and increase the NUE by reducing N inputs while maintaining grain yield. Further research is needed to fine-tune the response of maize to N applications when deficiencies are detected at V12 stage.

Acknowledgments:
The authors would like to thank Jordi Voltas for his help with the statistical analyses, Helga Ochagavia for her help in the graphic presentation of the results, Sergio Serrano for his help in the data processing, IRTA Research Station (Gimenells, Lleida, Spain) and Camilo Solsona for allowing the research to take place, and the GIS&Remote Sensing Laboratory of the University of Lleida for the use of their facilities to carry out image processing and analysis. This work was supported by the Spanish Ministry of Science and Innovation (Project AGL2012-35122) and the University of Lleida (Ángel Maresma held a PhD scholarship).
Author Contributions: Ángel Maresma contributed to the creation of the split-plots, supervised the experiment in the field, acquire field data and he wrote the bulk of the paper. Jaume Lloveras conceived and designed the experiment and contributed to the writing of the paper. José A. Martínez-Casasnovas worked on data acquisition, methodology and the analysis of the remote sensing data as well as writing the paper.

Conflicts of Interest:
The authors declare no conflict of interest. The funding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.