Agronomic and Economic Potential of Vegetation Indices for Rice N Recommendations under Organic and Mineral Fertilization in Mediterranean Regions

Rice (Oryza sativa L.) farmers in Mediterranean regions usually apply organic or mineral fertilizers before seeding that are supplemented with mineral nitrogen (N) later in the season. In general, the midseason N is applied without consideration of the actual crop N status, which may lead to over-fertilization and associated environmental problems. Thus, the purpose of this study was to design and evaluate a N recommendation approach using aerial images for Mediterranean paddy rice systems. A two-year rice field experiment was established in northeastern Spain, with different rates of pig slurry (PS) and mineral N fertilizer. Multispectral aerial images were taken at the rice booting stage, and several vegetation indices (VIs) were calculated. The VIs showed strong relationships with yield and the relations significantly differed between the PS and mineral fertilization treatments. The strongest relations with yield were obtained with gMCARINIR, proposed in this study, (R2 = 0.67), GNDVI (R2 = 0.64) and MCARINIR (R2 = 0.64), indicating the importance of including the green band information. The N recommendation approach generated using the VIs information showed a high success (87.5%) in the preliminary evaluation. The economic and environmental analysis showed that this approach provides a useful tool when compared to the usual farmer practices.


Introduction
Rice (Oryza sativa, L.) is the staple food for nearly half of the world's seven billion people [1].Nitrogen (N) is an essential element that is required to obtain high rice yields.Although N application increases rice productivity, poor N use efficiency is characteristic of irrigated rice systems, due to rapid losses of applied N; hence, adjusting N rates and the time of application to crop N requirements is crucial for optimizing N use efficiency and avoiding environmental problems, such as emission of greenhouse gases or surface and groundwater nitrate pollution [2].
In Europe, rice is mostly cultivated in Mediterranean countries, and 17% of the total area is cultivated in Spain [3].The rice extension in northeastern Spain (Aragon and Catalonia regions) represents 24% [4] of the rice lands in the country.Besides, half of the national pig production is concentrated in this area [5]; thus, it is necessary to integrate pig slurry (PS) in the N fertilization schedule for rice fields.The most common practice of local farmers is to apply PS or chemical fertilizers before seeding (~70% of crop N requirements) and to complement it with mineral N-topdressing at the end of the tillering stage or during stem elongation.However, N is typically applied during mid-season without consideration of the crop N status at that time.
Although some studies have shown that a single optimum N application before flooding could allow maximum yield [6][7][8], estimation of the potential yield at the beginning of the season is difficult, because it is strongly influenced by yearly variation in weather conditions, mainly temperature and wind [9]; thus, this N management presents a high risk for over-fertilization.However, if N rates are reduced before seeding and later complemented with topdressing N rates that are adjusted to the crop N status at that time, over-fertilization and related environmental problems can be reduced or avoided and N use efficiency may be increased, as reported by Nishikawa et al. [10] and Xie et al. [11].
One of the most common tools for real-time N management in cereals is a chlorophyll meter (SPAD 502, Minolta Corp., Ramsey, NJ, USA; N-tester, Yara International ASA, Oslo, Norway).This tool has been used successfully to monitor leaf N status and guides N fertilization during the crop season in different cereals [12,13], including rice [14,15].However, this method involves measurements at the leaf level; thus, it is inadequate for large-scale applications.Other approaches include the use of active crop canopy sensors, such as the GreenSeeker (Trimble Navigation Limited, Sunnyvale, CA, USA) or Crop Circle (Holland Scientific Inc., Lincoln, NE, USA), which measure the amount of light reflected by the canopy [16].These sensors reduce the measurement time and are more suitable for field applications; nevertheless, it is still necessary to enter the field and perform the measurements in different areas to obtain a representative value.Aerial or satellite remote sensing images improve these measurements, because they can provide spectral information for whole fields and large areas [17], although they also have limitations [18].Satellites have improved their spatial and temporal resolution over the last years, and cover larger areas at the same time, but are subject to fixed scheduling and strongly depend on cloud cover.Unmanned aerial vehicles (UAVs) offer very high spatial resolution, flexibility on scheduling, and its acquisition are independent on cloud cover conditions, however they are limited for the use in extensive areas and during wind conditions and the cost of a great effort for mosaicking and geocoding.The aircraft platform sits in between satellites and UAVs with more flexibility than satellite and better facility to cover the whole scene than UAVs [18].
Different authors have proposed different approaches for determining the rice N requirements in-season by using the relationships between canopy spectral information (vegetation indices; VIs) and different crop parameters [15,[19][20][21].The key to the success of these approaches is the prediction of yield responsiveness to additional N fertilization.Therefore, the establishment of strong relationships between spectral information and yield or other crop parameters is necessary for generating sound topdressing N recommendations.
The Normalized Difference Vegetation Index (NDVI), which integrates information in the red (R) and near-infrared (NIR) bands, is one of the most widely applied vegetation indices (VIs).The NDVI has been related to the leaf area index (LAI), biomass or yield in rice [22][23][24] and other crops [25][26][27], and it is the most widely used VI for N recommendation approaches [15,20,21,28].However, for high chlorophyll (Chl) concentration or large vegetation coverage values, NDVI losses sensitivity and saturates [16,29].Different authors have tried to handle this saturation phenomenon.Gitelson et al. [30] proposed the Green NDVI (GNDVI), that considers the green (G) instead the red (R) band in its formulation, and found that GNDVI was much more sensitive to Chl concentration, in a wide range of Chl variations, than the original NDVI.Lately, this index has been applied to rice [16].Other authors have proposed three-band VIs to solve the saturation issue [16,31].
Other indices have been formulated to evaluate responses to variation in chlorophyll and N content.One example is the Modified Chlorophyll Absorption in Reflectance Index (MCARI) that was designed to be responsive to chlorophyll variation [32]; later this index was modified by different authors by integrating the NIR wavelength to increase the sensitivity to LAI and aboveground biomass changes [31,33] and was applied to rice [16,33].
Most studies focused on the use of canopy reflectance for estimating crop parameters or adjusting mid-season N in rice paddy fields have been conducted in Asia, and different agricultural practices are used in this region compared to the Mediterranean areas of Europe.Under Mediterranean conditions, Gilabert and Meliá [24] established yield prediction based on VIs obtained from satellite images in Valencia (Spain), and Casanova et al. [22] estimated LAI and biomass from VIs obtained from radiometer measurements in the Ebro Delta (Spain).Recently, efforts have been underway to integrate multispectral information into farm advisory systems.The Earth obseRvation Model based RicE information Service (ERMES) project is being developed in Spain, Italy, and Greece with the primary aim of rice yield prediction [34,35].This project also includes the support of rice growers for fertilization or pest control and management, and two applications (PocketLAI and PocketN) have been developed for estimating in-season the crop status (LAI, leaf and plant N content) from digital photographs acquired with commercial smartphones [36,37].However, although the perspectives based on spectral information to improve N management and increase N use efficiency in rice are promising, more studies focusing on the development of approaches to estimate the N topdressing needs in season are needed, with a special emphasis on the effects of organic fertilizers, such as PS on rice spectral information.Knowledge of the spectral response sensitivity to differences between mineral and organic fertilization is essential to develop sound approaches for N recommendations.Moreover, the use of aerial images, that cover large areas, should be studied, so that recommendations can be made to different farmers or for large farm areas.
Therefore, the main purpose of this study was to design and evaluate a N recommendation approach using multispectral images in a Mediterranean paddy rice system under organic and mineral fertilization.To achieve this goal, three sub-objectives were proposed: 1.
To establish and compare the relationships between different VIs derived from aerial multispectral information and yield at the rice booting stage, and evaluate possible differences in these relationships due to organic and mineral fertilization.2.
To design and evaluate the agronomic performance of a N topdressing recommendation approach based on the information obtained in sub-objective 1.

3.
To compare economically and environmentally different scenarios for N adjustment based on the recommendation approach defined in sub-objective 2.
The work has been developed as a pilot study on the japonica rice cultivar, Guadiamar, which accounts for 70% of the rice surface in the Aragon region in Spain.

Experimental Design and Agricultural Practices
The study was conducted in a flooded rice field (41 •  The experimental design was a split plot with four repetitions (Figure 1).The main plots were assigned to three basal fertilization strategies consisting of two rates of pig slurry (PS) equivalent to 120 kg NH 4 -N•ha −1 (PS120) and 170 kg NH 4 -N•ha −1 (PS170) (Table 1, Pig slurry treatments, PS) and a mineral treatment (ammonium sulfate) at different rates (Table 1, Mineral treatments, M).Secondary plots included different topdressing N rates applied as ammonium sulfate.The experimental plots were 6 m wide by 12 m long for the PS treatments and 6 m wide by 5 m long for the M treatments (Figure 1).The total N rates ranged between 0 and 320 kg N•ha −1 and were applied as basal fertilization or as a combination of basal and topdressing fertilization for a total of 22 treatments (Table 1).
Pig slurry was band-spread on 15 May 2012 and 9 May 2013, and PS rates were established according to the ammonium N concentration of the PS, which was measured in situ by a Quantofix ® N-volumeter [38] and by conductimetry [39].On the same days, the basal mineral N fertilizer (ammonium sulfate) was applied to the plots of the M treatments at the corresponding rates together with P (100 kg P 2 O 5 •ha −1 ) and K (100 kg K 2 O•ha −1 ) to avoid P or K deficiency.
The japonica rice cultivar, Guadiamar, was broadcast-seeded on 16 May 2012 and 15 May 2013 at a seed rate of 180 kg•ha −1 .In 2012, the field was immediately flooded after seeding; however, in 2013, flooding preceded rice seeding.Topdressing N was applied at the end of the tillering stage on 4 July 2012 and 29 July 2013 as ammonium sulfate.The field remained flooded until approximately one month before harvest, except for occasional drainages for the application of herbicides, pesticides and fungicides according to habitual practices in the area.
Rice was harvested from 15-17 October 2012 and 25 October 2013.Grain moisture (PM-600 grain moisture tester, Keller, Japan) was measured to adjust the yield to a moisture content of 140 g•kg −1 .
Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 21 (ammonium sulfate) was applied to the plots of the M treatments at the corresponding rates together with P (100 kg P2O5•ha −1 ) and K (100 kg K2O•ha −1 ) to avoid P or K deficiency.The japonica rice cultivar, Guadiamar, was broadcast-seeded on 16 May 2012 and 15 May 2013 at a seed rate of 180 kg•ha −1 .In 2012, the field was immediately flooded after seeding; however, in 2013, flooding preceded rice seeding.Topdressing N was applied at the end of the tillering stage on 4 July 2012 and 29 July 2013 as ammonium sulfate.The field remained flooded until approximately one month before harvest, except for occasional drainages for the application of herbicides, pesticides and fungicides according to habitual practices in the area.
Rice was harvested from 15-17 October 2012 and 25 October 2013.Grain moisture (PM-600 grain moisture tester, Keller, Japan) was measured to adjust the yield to a moisture content of 140 g•kg −1 .

Spectral Information
Spectral information was collected using a spectral camera from a manned plane (RS Servicios de Teledetección SL, Lleida, Spain).The spectral camera collected data at four wavelengths.The center wavelengths were as follows: Band 1 (blue-B): 450 nm; band 2 (green-G): 550 nm; band 3 (red-R): 675 nm; and band 4 (near infrared-NIR): 780 nm; the bandwidth was 20 nm (for the four bands).The images were collected at solar noon, to minimize shadows effects in the spectral response, on 30 July 2012 and 13 August 2013 at the rice booting stage, a few days before heading.The images were provided pre-processed by SpecTerra Services (WA, Australia) with the differential illumination effect corrected and were ortho-rectified and mosaicked.The spatial resolution was 0.1 m and the radiometric resolution was 16 bits.The spectral information was provided in digital values, thus adimensional VIs were chosen for the study.
Different VIs were evaluated for their relationship to rice yield (Table 2).The indices RVI, GRVI, NDVI and GNDVI were included because they were consistently shown to be related to agronomic parameters, such as yield, and have been used in the development of approaches to recommend N topdressing in different crops.Furthermore, we included three indices, MCARI1, MCARI NIR and gMCARI NIR , which were derived from the Modified Chlorophyll Absorption in Reflectance Index, that was originally developed to be responsive to chlorophyll variation [32].
Haboudane et al. [31] modified the MCARI and proposed the MCARI1 with suppression of the R 700 /R 670 ratio to lower the sensitivity to chlorophyll effects and replacement of the red-edge wavelength (R 700 ) by the near-infrared wavelength (R 800 ) to increase the sensitivity to green LAI variation.Cao et al. [33] modified the MCARI to work with the Crop Circle ACS-470 active sensor.The MCARI modified by Cao et al. [33], also denoted MCARI1 retained the formula of the MCARI, but the R 700 and R 670 were replaced by R NIR (R 760 ) and R red edge (RE, R 730 ), respectively.This index showed consistent relations with rice aboveground biomass (R 2 = 0.79) and plant N uptake (R 2 = 0.83) across growth stages [33].In our study, we have adapted this index to the available bands, NIR 780 nm and red 675 nm, and we have denoted the index as MCARI NIR (Table 2).
The index gMCARI NIR , proposed in this work, retains the MCARI structure, includes the near-infrared wavelength to increase the sensitivity to biomass changes as proposed by Haboudane et al. [31] and incorporates the green reflectance peak information (G-R) to account for the sensitivity to chlorophyll concentrations (Table 2).This index is the product of Green Difference Vegetation Index (NIR-G; [43]) and the Ratio Vegetation Index (NIR/R; [40]).
The VIs were calculated for each pixel and the average values for each plot were extracted.The extraction was performed using a mask for each plot that excluded a width of 1 m from the borders to avoid edge effects.In 2013, due to adverse meteorological conditions, rice seed germination was hindered and five plots were removed from the analyses because of bad emergence.

Relationship between Yield and Vegetation Indices
Many studies have reported improvements in yield estimation of irrigated rice after relativization of VIs with respect to overfertilized plots [21] to obtain sufficiency indices.Although some authors have obtained strong relationships between vegetation indices and yield in different years without converting the yield to relative value [20,21,24,44], in our study, it was not possible to obtain an equation to predict absolute yield values, due to high interannual yield variability.The variability in yield between years (7.8 Mg•ha −1 in 2012 and 5.9 Mg•ha −1 in 2013 in this study) in the study area is high for two reasons.Firstly, the area is at the low limit of temperature for adequate rice cultivation.Secondly, there is strong variability in meteorological conditions between years.Hence, relative yield (R_yield) and the relative VI (R_VI) were calculated for each year (Equations ( 1) and ( 2)). R_V Linear-plateau equations ([7,45]; Equation (3)) were adjusted to model the response of yield to N rates in the PS and the M treatments for each of the two years in order to get Yield max and VI max .
where Y (Mg•ha −1 ) is the yield; N is the applied nitrogen rate (kg N•ha −1 ) (this rate is the sum of N applied before seeding and topdressing); a (intercept) is the yield at 0 kg N•ha −1 ; b is the increase in yield per unit increase in N; and C is the critical (or optimum) N rate, i.e., the minimum N rate above which the maximum yield is obtained.Yield max is the maximum yield or plateau of Equation (3), i.e., the yield for N ≥ than the critical N rate (C).The maximum VI (VI max ) for each year was calculated as the average VI value of treatments with N rates equal to, or higher than, the critical N rate (N ≥ C).The relationships between relativized yield and relativized vegetation indices were established by regression analysis for the years 2012 and 2013 and the pooled data (using 100% of the data).The models tested were linear, multiplicative, logarithmic and exponential.The coefficient of determination and the root mean square error (RMSE) of the regressions were used to evaluate the performance of the seven VIs.
The differences between years and between PS and M treatments were evaluated by a test of equality for regression lines across groups [46].

Design
The N topdressing recommendation approach was based on the calculation of a N sufficiency index (NSI) as indicator of the crop N status, i.e., as quantifiers of N nutrition index.This approach was defined by Denuit et al. [12] using the HydroAgri N-Tester sensor and corrects the readings of the sensor with the readings of the N over-fertilized plot to obtain sufficiency indices (in this case relativized VIs, i.e., NSI is R_VI, Equation ( 4)).Thus, a NSI value below 0 means the plot could present a nutritional deficit, and NSI values of 1 or above 1 means the plot presents 'a priori' a good nutritional status.This type of approach has been successfully tested in rice by Chen et al. [19] or Xue and Yang [20].
The approach is based on the relationship between the variables Delta NSI and Delta N (Equation ( 5)).Delta NSI (Equation ( 6)) is the difference between NSI and 1.Thus, conceptually if Delta NSI is 0 or positive, the plot should have a good nutritional status, and if it is negative, the plot might present a N deficit.Delta N (Equation ( 7)) is the difference between the N applied in the treatment (N T ) and the critical N rate (C) obtained in the yield response to N (Equation ( 3)).
This approach was selected because it does not require absolute yield predictions; instead, N deficit or additional N required for raising NSI from a current level to a target level is calculated based on the relationship of NSI with N rate [20].
If Delta N estimate from spectral information is positive, topdressing N fertilization is not necessary; and on the other hand, if Delta N estimate is negative, the plot must be fertilized.In this case, Delta N gives the N deficit, i.e., the amount of N that needs to be applied to obtain maximum yield.
Three of the four replicates (75% of the data) of the experiment were used for establishing the model (Equation ( 5)).The fourth replicate was used for performance evaluation.The replicate used in the evaluation process was not the same for all treatments and was determined by random draw in each treatment.
Among the seven VIs evaluated, the two indices that presented the strongest relation to yield were chosen for the establishment of the N recommendation model.These two indices were GNDVI and gMCARI NIR , (see Section 3.1.2).

Validation Process
In the validation process, the values of R_GNDVI and R_gMCARI NIR for the validation plots (25% data) in each year (2012 and 2013) were obtained dividing the VIs values obtained in each plot by the average values of the over-fertilized plots of the corresponding year (2012 and 2013).Plots from PS120M150 and PS170M150 were considered overfertilized for PS treatments and plots from M120M90 and M120M120 for M treatments (Table 2).Then, the values of Delta N were obtained using Equation (5).The performance of the model was evaluated by the percent of success.The sign of Delta N and the actual R_yield were compared in the validation plots and the plots were assigned to one of the following three options.

•
Success: Delta N was negative (i.e., the plot would have needed additional N fertilization) and the actual R_yield was below 1; or Delta N was positive (i.e., the plot would not have needed additional N fertilization) and the actual R_yield was equal or higher than 1.

•
Failure by excess: Delta N was negative, but the actual R_yield was equal or above 1 (the approach would have recommended additional N fertilization, but the plot had reached the optimum yield).

•
Failure by defect: Delta N was positive, but the actual R_yield was below 1, (the approach would not have recommended additional N fertilization, but the plot had not reached the optimum yield).
Success and error percentages were calculated considering four strategies: • The use of GNDVI

•
The use of gMCARI NIR • The combination GNDVI & gMCARI NIR : The plot will only be fertilized if both VIs recommend additional N fertilization.

•
The combination GNDVI or gMCARI NIR : If one of the VIs recommends N fertilization, the plot will be fertilized even if the other index does not recommend N fertilization.
The fertilization treatments varied from 0 to 240 kg N•ha −1 and hence, some of the plots presented a high N deficiency.These plots with high N deficiency do not represent a real field situation and give an easy success in the validation; therefore, plots with R_yield below 0.7 were eliminated of the validation.Therefore, the set of 24 plots used for validation represents a real scenario.

Economic and Environmental Analysis
The net benefit due to N topdressing application was calculated according to Equation (8).
where Y Ntop (t•ha −1 ) is the yield hypothetically reached (according to response equation) if the plot had been fertilized; Y real (t•ha −1 ) is the actual yield in the plot, Price grain is the rice grain price ($•t −1 ); Price N is the N (ammonium sulfate) price ($•kg −1 ); N rate (kg•ha −1 ) is the N rate established for each scenario; and Cost Napplication ($•ha −1 ) is the price for N application (fuel + man labor).Prices of the local market (year 2015) [47] were considered for the calculation.The net benefit did not include the costs of the image acquisition and further processing.The net benefit was calculated for each of the 24 plots used in the validation process.Different scenarios were economically compared to a reference scenario.

•
Reference: All plots are fertilized with a fixed predefined N rate (N fix ) without using any recommendation approach.This practice is currently used by farmers in the study area and will be considered as the reference scenario.
• Scenario 1: Plots are fertilized according to Delta N estimates.Topdressing N is applied to the plots when Delta N is negative.The N rate is given by Delta N, (i.e., if Delta N = −30, N rate will be 30 kg N•ha −1 ) (Equation ( 9)).
• Scenario 2: Is a variation of scenario 1, but a minimum topdressing N rate (N m ) is established (Equation ( 10)).This option was considered since machinery is not prepared to apply fertilizers at low rates and farmers do not usually go inside the field to apply small N rates.
• Scenario 3: Is a variation of scenario 2. Plots are fertilized according to Delta N approach establishing a fixed predefined N rate (N fix ), i.e., if Delta N estimate is positive, the plots will not be fertilized, if the Delta N estimate is negative, the plots will be fertilized with N fix (Equation ( 11)).
To consider the environmental impact of overfertilization, N excess was also evaluated.For the plots that would be fertilized according to the different scenarios, N excess was calculated as the difference between the N rate applied and the N rate that would be necessary according to yield response to N equations (Equation (3)).
For all scenarios, different predefined N rates (N fix ) and minimum N rates (N m ) were evaluated in the range 0-200 kg N•ha −1 and the net benefit and the N excess were represented graphically.

Relationships between Yield and Vegetation Indices
3.1.1.Influence of the Year and TYPE of fertilizer Significant relationships were observed between R_yield and the seven R_VIs and those relationships did not significantly differ between the two years (Table 3).Thus, information from years with different crop development and yield potential could be joined if both yield and the indices are relativized to the maximum values of each year.However, the relationship between R_yield and R_VIs significantly differed between the PS and M treatments (Table 3; Figure 2a for NDVI).For example, for the NDVI, for a fixed value of the R_NDVI, the PS treatments had a higher expected yield than the M treatments (Figure 2a), i.e., the PS treatments reached that maximum yield with a lower value of the R_NDVI than in the M treatments.
Remote Sens. 2018, 10, x FOR PEER REVIEW 2 of 21 Thus, information from years with different crop development and yield potential could be joined if both yield and the indices are relativized to the maximum values of each year.However, the relationship between R_yield and R_VIs significantly differed between the PS and M treatments (Table 3; Figure 2a for NDVI).For example, for the NDVI, for a fixed value of the R_NDVI, the PS treatments had a higher expected yield than the M treatments (Figure 2a), i.e., the PS treatments reached that maximum yield with a lower value of the R_NDVI than in the M treatments.These results suggest that the VIs need to be relativized considering the PS and M treatments separately.Therefore, the R_VIs were recalculated using maximum VI values separately for the PS and M treatments within each year.After this process, the relationship between R_yield and R_VIs did not significantly differ between the PS and M treatments (Figure 2b).

Performance and Comparison of the Indices
The R_yield and the seven R_VIs for each year individually and for the pooled data over the two years (2012&2013) were significantly related (Table 4).Linear models fitted well the relationships between R_yield and the R_RVI, R_GRVI, R_NDVI, R_GNDVI and R_MCARI1 (Figure 3a-d) (Table 4).For the R_MCARINIR and R_gMCARINIR (Figure 3e,f), the increase in R_yield per unit increased in the indices decreased as the value of the indices increased, indicating that the relationship was not linear.For these two indices, the best fit was obtained with the multiplicative model (R_yield = a•R_VI b ).The residuals from all regressions were independent and normally distributed.These results suggest that the VIs need to be relativized considering the PS and M treatments separately.Therefore, the R_VIs were recalculated using maximum VI values separately for the PS and M treatments within each year.After this process, the relationship between R_yield and R_VIs did not significantly differ between the PS and M treatments (Figure 2b).

Performance and Comparison of the Indices
The R_yield and the seven R_VIs for each year individually and for the pooled data over the two years (2012&2013) were significantly related (Table 4).Linear models fitted well the relationships between R_yield and the R_RVI, R_GRVI, R_NDVI, R_GNDVI and R_MCARI1 (Figure 3a-d) (Table 4).For the R_MCARI NIR and R_gMCARI NIR (Figure 3e,f), the increase in R_yield per unit increased in the indices decreased as the value of the indices increased, indicating that the relationship was not linear.For these two indices, the best fit was obtained with the multiplicative model (R_yield = a•R_VI b ).The residuals from all regressions were independent and normally distributed.The coefficients of determination ranged between 0.40 and 0.77 (Table 4); in general, the R_NDVI and R_GNDVI improved the relationships in comparison to the R_RVI and R_GRVI (all of them including two bands in their definition).R_gMCARINIR, which includes information on three bands, slightly improved the relationship in comparison to the two-band VIs.
In 2012, R_GNDVI and R_gMCARINIR explained 77% of the R_yield variability and performed better than the other indices.For 2013, R_gMCARINIR was the index that explained the highest percentage (61%) of the R_yield variability.
When data from both years were pooled, R_gMCARINIR, R_MCARINIR and R_GNDVI showed the highest coefficients of determination (R 2 = 0.67, 0.64 and 0.64 respectively) and the lowest RMSE (0.139, 0.145, 0.144) and all of them included the green band information in their definition.The coefficients of determination ranged between 0.40 and 0.77 (Table 4); in the R_NDVI and R_GNDVI improved the relationships in comparison to the R_RVI and R_GRVI (all of them including two bands in their definition).R_gMCARI NIR , which includes information on three bands, slightly improved the relationship in comparison to the two-band VIs.
In 2012, R_GNDVI and R_gMCARI NIR explained 77% of the R_yield variability and performed better than the other indices.For 2013, R_gMCARI NIR was the index that explained the highest percentage (61%) of the R_yield variability.
When data from both years were pooled, R_gMCARI NIR , R_MCARI NIR and R_GNDVI showed the highest coefficients of determination (R 2 = 0.67, 0.64 and 0.64 respectively) and the lowest RMSE (0.139, 0.145, 0.144) and all of them included the green band information in their definition.

Design of the N Topdressing Recommendation Approach
The indices that presented the strongest relation to yield were GNDVI, MCARI NIR and gMCARI NIR .The index gMCARI NIR was selected for the design of the model because it showed the best relation to yield.GNDVI was selected versus MCARI NIR (both with the same strength in their relation to yield) with the aim to test two different functions (linear and multiplicative) and two different indices (two and three bands VIs).
The equations adjusted (Figure 4) to estimate Delta N from Delta NSI (GNDVI) and Delta NSI (gMCARI NIR ) were Equations ( 12) and ( 13 The indices that presented the strongest relation to yield were GNDVI, MCARINIR and gMCARINIR.The index gMCARINIR was selected for the design of the model because it showed the best relation to yield.GNDVI was selected versus MCARINIR (both with the same strength in their relation to yield) with the aim to test two different functions (linear and multiplicative) and two different indices (two and three bands VIs).

Assessment of the N Topdressing Recommendation Approach
The N recommendation approach designed had a high success rate (higher than 83%) and performed better with gMCARINIR, with an 87.5% of success, than with the GNDVI with an 83.3% success (Table 5).The percentage of failure by excess was lower for gMCARINIR than for GNDVI; however, the percentage of failure by defect was higher for gMCARINIR than for GNDVI.The

Assessment of the N Topdressing Recommendation Approach
The N recommendation approach designed had a high success rate (higher than 83%) and performed better with gMCARI NIR , with an 87.5% of success, than with the GNDVI with an 83.3% success (Table 5).The percentage of failure by excess was lower for gMCARI NIR than for GNDVI; however, the percentage of failure by defect was higher for gMCARI NIR than for GNDVI.The combination of both VIs did not improve the ability for N recommendation prediction (Table 5).The GNDVI & gMCARI NIR strategy showed the same percentages of success and failure than gMCARI NIR , and the GNDVI or gMCARI NIR possibility showed the same results than GNDVI.
Since the combination of VIs did not improve the ability of N recommendation in comparison to using only one VI, gMCARI NIR (the VI with the highest percentage of success) was considered the best option and it was evaluated in the economic and environmental analysis.In the framework of the four scenarios analyzed in this work, the net benefit due to topdressing N fertilization ranged between −50 $•ha −1 and 132 $•ha −1 (the additional benefit if these plots had been fertilized with additional N) (Figure 5).
Figure 5 shows the net benefit (a) and N excess (b) according to the fixed predefined N rate (N fix ) for the reference scenario and scenario 3 or to the minimum N rate (N m ) in the case of scenario 2.
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 21 combination of both VIs did not improve the ability for N recommendation prediction (Table 5).The GNDVI & gMCARINIR strategy showed the same percentages of success and failure than gMCARINIR, and the GNDVI or gMCARINIR possibility showed the same results than GNDVI.Since the combination of VIs did not improve the ability of N recommendation in comparison to using only one VI, gMCARINIR (the VI with the highest percentage of success) was considered the best option and it was evaluated in the economic and environmental analysis.In the framework of the four scenarios analyzed in this work, the net benefit due to topdressing N fertilization ranged between −50 $•ha −1 and 132 $•ha −1 (the additional benefit if these plots had been fertilized with additional N) (Figure 5).
Figure 5 shows the net benefit (a) and N excess (b) according to the fixed predefined N rate (Nfix) for the reference scenario and scenario 3 or to the minimum N rate (Nm) in the case of scenario 2. The reference scenario (current farmers' practice), where all plots are fertilized, obtains the lowest net benefit, with a maximum of 84 $•ha −1 for a rate of 75 kg N•ha −1 (Figure 5a), and the highest N excess of the four scenarios reaching 45 kg N excess•ha −1 for the N rate of maximum net benefit (Figure 5b).
The use of the N recommendation approach (Scenario 1) shows a net benefit of 97 $•ha −1 and a N excess of 2.3 kg N•ha −1 .The net benefit and N excess associated with this strategy have been represented as a horizontal green line in Figure 5 for comparative purposes with the rest of scenarios.
The establishment of a threshold minimum N rate (N m ) in Scenario 2, shows a higher net benefit than Scenario 1, reaching 132 $•ha −1 for a threshold N rate of 90 kg N•ha −1 (Figure 5a).However, the increase of the net benefit (35 $•ha −1 ) is counterbalanced with an increase in N excess (20 kg•N•ha −1 for the highest benefit versus 2.3 kg N•ha −1 in Scenario 1) (Figure 5b).
In the case of fixing a predefined N rate (Scenario 3), the maximum net benefit (132 $•ha −1 ) is equal to that in Scenario 2 for the rate of 90 kg•ha −1 (Figure 5a).In this scenario, the net benefit dramatically decreases as the N fix decreases; therefore, the election of the optimum N fix rate is a key factor, i.e., if N fix is too low, the net benefit could be dramatically reduced.For the optimum N fix rate of this study (90 kg•ha −1 ), the N excess was 20 kg•N•ha −1 , the same as in Scenario 2 (Figure 5b).

Influence of the Type of Fertilizer
The relationship between R_yield and R_VIs showed differential effects for the PS and M fertilization treatments for the seven VIs analyzed (Figure 2a).Although studies focusig on the differences in spectral information between different types of fertilization in rice fields are lacking, there are some studies in other crops that support these results.Zhao et al. [48] measured the canopy apparent photosynthesis (CAP), the photosynthetic rate of flag leaves, LAI and yield in a winter wheat experiment in plots fertilized with cow manure, urea and a mixed application.The results showed that during the early growth period, a single application of urea promoted better crop development and resulted in higher values of CAP or LAI; however, during the late growth stage, a single application of cow manure and the mixed application delayed the leaf senescence process compared to a single application of urea.The results suggested that mixed application of organic and inorganic fertilizers delayed leaf senescence and maintained better canopy structure and higher photosynthetic capability at the late grain filling stage, which resulted in a higher grain yield.In our study, we found the same type of behavior as Zhao et al. [48], i.e., for the same value of the R_VIs at the booting stage, the PS treatments reached a higher yield compared to the M plots.This better development during the late stage is usually related to higher availability of micronutrients provided by the organic fertilizers [49][50][51].
The observed differences between the type of fertilization is an important result for the development and use of N recommendation approaches.This result implies that it is necessary to establish over-fertilized plots for each type of fertilization.

Performance and Comparison of the Indices
The R_VIs showed good relationships with R_yield (Table 4).Other studies have also found good relationships between these indices and rice yield [15,16,24,29,52].The coefficients of determination are in the range reported by the abovementioned studies at the same growth stage (booting) [15,16].
One common problem with the two-band VIs is that they become saturated under high values of biomass [16,29] and different authors have proposed three-band VIs to handle this saturation phenomenon.In our study, the R_MCARI1 (three-band VI) did not improve the R_yield prediction in comparison with the traditional two-band indices and the relationship even worsened; this result contrasts the study reported by Haboudane et al. [31] in which the MCARI1 was less sensitive to the saturation phenomenon, although these findings were for different crops than rice (i.e., corn, wheat and soybean).However, the R_MCARI NIR and R_gMCARI NIR had a strong relationship with R_yield and performed equal to or better than the traditional R_NDVI and R_GNDVI (Table 4).The index proposed in this study (R_gMCARI NIR ) was the best for R_yield prediction and improved the relationship obtained with the R_MCARI NIR .
Although the relationships between R_yield and the R_MCARI NIR or R_gMCARI NIR (three-band VIs) were stronger than for the two-band VIs, there was only a 1-7% increase in the variability explained by these indices in comparison to two-band VIs, such as the R_NDVI and R_GNDVI (Table 4).These results contrast with those from the study reported by Cao et al. [16] in which the yield potential variability explained by the three-band indices was between 21 and 26% higher than the variability explained by the NDVI or RVI at the rice booting stage.
These differences in yield prediction improvement with three-band indices may be related to potential yield.The highest yield in this study was approximately 8000 kg•ha −1 ; however, yields reached 10,000 kg•ha −1 in the study conducted by Cao et al. [16] and 14,000 kg•ha −1 for the study conducted by Harrell et al. [29] in which the saturation phenomenon was also observed.Therefore, the saturation phenomenon can be attributed to high biomass conditions in high-yield systems.In these systems, two-band VIs do not perform well, especially at later stages with high biomass conditions; however, three-band VIs can overcome this saturation problem.
Thus, in low-yielding rice systems similar to our system, two-band indices, such as the NDVI and GNDVI, are not subject to the saturation phenomenon, due to low crop biomass.This hypothesis is in agreement with the results of Xue et al. [21], who found strong relationships between the NDVI and rice yield potential across the growing season from tillering to grain filling, with yield levels lower than 8000 kg•ha −1 (similar to this study); thus, the saturation phenomenon, due to high biomass conditions was not observed.
Thus, in our study, the best indices to predict R_yield were the R_GNDVI, R_MCARI NIR and R_gMCARI NIR , both incorporate the green band.This result supports the findings of Cao et al. [16], which demonstrated that the best VIs for estimation of the response index at harvest (R_yield) were green band-based VIs.

Assessment of the N Topdressing Recommendation Approach
The approach evaluated in this study showed success percentages higher than 80% and reaching 87.5% for the best options (Table 5).Thus, these options confer an advantage over the application of N topdressing without any advice.Nevertheless, it is important to point out that this validation process means only a preliminary evaluation of the approach, since it is based on the comparison of the approach's recommendation and the actual yield harvested in the same experiment where the models (equations) were obtained.The approach should be tested further in real field situations, fertilizing the fields following the recommendations of the models and evaluating the yield obtained [15,20,21] compared to other fields with a standard fixed N rate.
In contrast to other studies where only one index was used in the development of these approaches, in our study, the combination of two indices at the same time was tested (the plot is fertilized when both VIs recommend N topdressing or when one of the VIs recommends N fertilization even if the other index does not recommend N fertilization).Although these options did not increase the percentage of success under the studied conditions, it should be considered in further studies, because in other situations it could increase the percentage of success.
The results show that the Reference Scenario, in which all plots are fertilized without spectral information consideration, is the least recommended scenario of the four evaluated scenarios (Figure 5).This scenario shows the lowest net benefit and the highest N excess, suggesting that this option is neither economically nor environmentally viable when compared to the other scenarios.Therefore, the use of a decision tool for in season topdressing N recommendation is clearly advantageous.
The N recommendation approach evaluated in this study seems to be a good strategy.Three possibilities were evaluated: Fertilizing with the N rate according to the approach, Delta N, (Scenario 1), establishing a minimum N rate, N m , (Scenario 2) or establishing a fixed predefined N rate, N fix , (Scenario 3).When a minimum N rate is established (Scenario 2), the maximum net benefit increases in comparison to Scenario 1, the maximum difference is 35 $•ha −1 , for a minimum N rate of 90 kg N•ha −1 , but the N excess increases in 18 kg N•ha −1 (Figure 5).Combining net benefit and environmental impact, we would recommend a N rate of 60 kg N•ha −1 with an increase in net benefit of 27 $•ha −1 and an increase in N excess of 6.1 kg N•ha −1 in comparison to Scenario 1.On the other hand, when a fixed predefined N rate is established (Scenario 3), the maximum net benefit obtained is the same as that obtained in Scenario 2. However, in Scenario 3, a careful selection of the N fix rate is crucial, since the net benefit can dramatically lower with small modifications of the N rate (Figure 5); this N rate will depend on the year and it is difficult to estimate.Hence, Scenario 1 and Scenario 2 seem to be better than Scenario 3. Scenario 2 was considered because farmers do not usually enter the field to apply small N rates or even machinery is not able to apply tiny amounts of fertilizer, thus a minimum N rate should be considered.This option allows increasing the net benefit under a large interval of N m rates (from 0 to 150 kg N•ha −1 ) (Figure 5), in which the net benefit is the same or higher than in Scenario 1.Therefore, in Scenario 2 the net benefit is not at risk; however, the N excess increases as N m increases.The cost of the risk for contamination, the N excess, has not been considered in the economic analysis and should be incorporated in future works.As a first approach in Scenario 2, N m rates between 50-60 kg N•ha −1 are believed to be optimal considering both economic return and N excess.
Thus, the economic analysis suggests that the best option is the use of the N recommendation approach fertilizing with the recommended N rate (Scenario 1) or establishing a minimum N rate (Scenario 2) in order to increase the net benefit, although the increase of N excess should be economically evaluated.
These findings should be validated in a different set of field experiments, where some of the plots were fertilized according to the approach's N recommendation and other plots according to the local practices in order to evaluate the real benefits.
The inclusion of the cost of the images' capture and further processing is crucial to evaluate whether these tools are economically feasible, although multispectral information arises as a useful tool for increasing net benefit and decreasing N excess in the agro-systems of Northeast Spain.

Conclusions
The seven studied VIs showed good relationships with yield.However, differences between the PS and M treatments were detected, i.e., fields fertilized with pig slurry had a different spectral response than those fertilized with synthetic fertilizers.This is an important result for the development and use of N recommendation approaches since it is necessary to establish over-fertilized plots for each type of fertilization.
The best relationships between yield and the VIs were obtained with indices including the green band.Additionally, three-band VIs performed better than the traditional two-band VIs, but the improvement was minor compared to other studies with higher yields when the saturation of two-band indices was observed.The best relationships were obtained with GNDVI, MCARI NIR and gMCARI NIR ; the latter was proposed in this study, all three included the green band, and the MCARI NIR and gMCARI NIR included 3 bands.
The N recommendation approach evaluated was useful for N recommendation.The best option was the use of the index gMCARI NIR achieving an 87.5% of success and the combination of VIs did not improve the ability for N recommendation prediction.The economic analysis showed that the use of a N recommendation approach clearly increases the net benefit and lowers the N excess in comparison to fertilization without any recommendation.The best option in this study was the use of the recommended N rate (Delta N) or establishing a minimum N rate (optimum minimum N rate between 50 and 60 kg N•ha −1 ).
Therefore, the use of aerial remote sensing is a promising tool for developing strategies for advising rice farmers.However, more research is needed, including a validation of this approach to field level, the inclusion of the cost of the recommendation system in the evaluation of the net benefit and the response of yield to spectral information in earlier crop development stages to adjust the N topdressing as much as possible to the usual practices in the area.
45 31.87N, 0 • 2 18.16 W), with three fertilization strategies during two consecutive years (2012 and 2013) [7].The climate of the region is semiarid continental Mediterranean with high temperatures during the summer and low precipitation (15.0 • C annual average temperature and 349 mm annual precipitation; average period 1980-2010).

Figure 1 .
Figure 1.Experimental plots layout (yellow lines) superimposed on the images (false color RGB: Near-Infrared band/ Green band/ Red band) taken with a multispectral camera at the booting stage, 30 July 2012 (a) and 13 August 2013 (b).The shaded areas represent additional treatments excluded from the analysis in this paper.

Figure 1 .
Figure 1.Experimental plots layout (yellow lines) superimposed on the images (false color RGB: Near-Infrared band/ Green band/ Red band) taken with a multispectral camera at the booting stage, 30 July 2012 (a) and 13 August 2013 (b).The shaded areas represent additional treatments excluded from the analysis in this paper.

Figure 2 .
Figure 2. Relationship between R_yield and R_NDVI (relativized to maximum values of each year) in the PS and M treatments (a) and relationship between R_yield (relativized to maximum value of each year) and R_NDVI (relativized to maximum value of each year and within each year individually for the PS and M treatments) for PS and M treatments (b); 2012 and 2013 pooled data.

Figure 2 .
Figure 2. Relationship between R_yield and R_NDVI (relativized to maximum values of each year) in the PS and M treatments (a) and relationship between R_yield (relativized to maximum value of each year) and R_NDVI (relativized to maximum value of each year and within each year individually for the PS and M treatments) for PS and M treatments (b); 2012 and 2013 pooled data.

Figure 5 .Figure 5 .
Figure 5. Net benefit ($•ha −1 ) (a) and N excess (kg N•ha −1 ) (b) according to the fixed predefined N rate (Nfix, kg N•ha −1 ) or minimum N rate (Nm, kg N•ha −1 ).Each point is the average of the 24 plots used in the economic analysis.

Table 1 .
Amounts of N applied (kg N ha −1 ) before seeding (BS) and topdressing (TS) in the different treatments.For the PS treatments, amount indicates target N rates.

Table 3 .
Best models for the relationships between R_yield (y) (relativized to maximum value of each year) and R_VIs (x) (relativized to maximum value of each year) for years 2012 and 2013, for the PS and M treatments.The last two rows give the significance of the comparison between years and the type of fertilizer.

Table 4 .
Coefficients of determination (R 2 ) of the relationships between R_yield and R_VIs for years 2012 and 2013 and the pooled data, and RMSE for the pooled data.

Table 4 .
Coefficients of determination (R 2 ) of the relationships between R_yield and R_VIs for years 2012 and 2013 and the pooled data, and RMSE for the pooled data.

Table 5 .
Percentage of success and failure by excess and defect using the indices GNDVI, gMCARI NIR and the combinations GNDVI & gMCARI NIR and GNDVI or gMCARI NIR (25% data, excluded plots with R_yield below 0.7, total number of plots used for validation: 24).

Table 5 .
Percentage of success and failure by excess and defect using the indices GNDVI, gMCARINIR and the combinations GNDVI & gMCARINIR and GNDVI or gMCARINIR (25% data, excluded plots with R_yield below 0.7, total number of plots used for validation: 24).