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

Abstract

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 gMCARI_NIR, proposed in this study, (R² = 0.67), GNDVI (R² = 0.64) and MCARI_NIR (R² = 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.

1. 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:

• 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.
• To design and evaluate the agronomic performance of a N topdressing recommendation approach based on the information obtained in sub-objective 1.
• 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.

2. Materials and Methods

2.1. Experimental Design and Agricultural Practices

The study was conducted in a flooded rice field (41°45′31.87′′ N, 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).

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₄-N·ha⁻¹ (PS120) and 170 kg NH₄-N·ha⁻¹ (PS170) (

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

Pig Slurry Treatments (PS)			Mineral Treatments (M)
	BS	TP		BS	TP
	kg NH4-N·ha−1	kg N·ha−1		kg N·ha−1
PS120M0	120	-	M120M0	120	-
PS120M30	120	30	M120M30	120	30
PS120M60	120	60	M120M60	120	60
PS120M90	120	90	M120M90	120	90
PS120M120	120	120	M120M120	120	120
PS120M150	120	150
PS170M0	170	-	Control (M0)	-	-
PS170M30	170	30	M30	30	-
PS170M60	170	60	M60	60	-
PS170M90	170	90	M90	90	-
PS170M120	170	120	M120 = M120M0	120	-
PS170M150	170	150	M150	150	-

, 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⁻¹ 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₂O₅·ha⁻¹) and K (100 kg K₂O·ha⁻¹) 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⁻¹. 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⁻¹.

2.2. 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. Vegetation indices (VIs) evaluated in this study (Green-G: 550 nm, Red-R: 675 nm and Near Infrared-NIR: 780 nm).

Indices (VIs)	Formula	Reference
RVI Ratio Vegetation Index	NIR/R	[40]
GRVI Green Ratio Vegetation Index	NIR/G	[41]
NDVI Normalized Difference Vegetation Index	(NIR-R)/(NIR+R)	[42]
GNDVI Green Normalized Difference Vegetation Index	(NIR-G)/(NIR+G)	[30]
MCARI1 Modified Chlorophyll Absorption in Reflectance Index1	1.2[2.5(NIR-R)-1.3(NIR-G)]	[31]
MCARINIR Modified Chlorophyll Absorption in Reflectance IndexNIR	[(NIR-R)-0.2(NIR-G)](NIR/R)	Adapted from Cao et al. [33]
gMCARINIR Green peak Modified Chlorophyll Absorption in Reflectance IndexNIR	[(NIR-R)-(G-R)](NIR/R) = (NIR-G)(NIR/R)	Proposed in this study

). 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₇₀₀/R₆₇₀ ratio to lower the sensitivity to chlorophyll effects and replacement of the red-edge wavelength (R₇₀₀) by the near-infrared wavelength (R₈₀₀) 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₇₀₀ and R₆₇₀ were replaced by R_NIR (R₇₆₀) and R_{red edge} (RE, R₇₃₀), respectively. This index showed consistent relations with rice aboveground biomass (R² = 0.79) and plant N uptake (R² = 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]).

These vegetation indices have been successfully applied to monitor field crops at different spatial resolutions; spatial resolution of centimeters using field devices, as radiometers or Crop Circle sensors [16,21,22,23,33], spatial resolution of meters with aerial imagery [18,31], or spatial resolution of several meters with satellite images [17,18,24,25].

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.

2.3. 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⁻¹ in 2012 and 5.9 Mg·ha⁻¹ 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\_yiel{d}_{\left(plot\right)}=\frac{Yiel{d}_{\left(Plot\right)}}{Yiel{d}_{max}}$$

(1)

$$R\_V{I}_{\left(plot\right)}=\frac{VIvalu{e}_{\left(Plot\right)}}{V{I}_{max}}$$

(2)

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.

$$\begin{array}{c}IfN<C;Y=a+b·N\\ IfN\ge C;Y=a+b·C\end{array}$$

(3)

where Y (Mg·ha⁻¹) is the yield; N is the applied nitrogen rate (kg N·ha⁻¹) (this rate is the sum of N applied before seeding and topdressing); a (intercept) is the yield at 0 kg N·ha⁻¹; 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 ($\mathrm{N}\ge \mathrm{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].

2.4. N Topdressing Recommendation Approach

2.4.1. 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].

$$NSI=R\_V{I}_{\left(Plot\right)}=\frac{VIvalu{e}_{\left(Plot\right)}}{V{I}_{max}}$$

(4)

$$Delta N=f\left(Delta NSI\right)$$

(5)

$$DeltaNSI=NSI-1$$

(6)

$$DeltaN=N_T-C$$

(7)

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).

2.4.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⁻¹ 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.

2.4.3. Economic and Environmental Analysis

The net benefit due to N topdressing application was calculated according to Equation (8).

$$Netbenefit=\left(Y_{Ntop}-Y_{real}\right)·Pric{e}_{grain}-Pric{e}_N·N_{rate}-Cos{t}_{Napplicattion}$$

(8)

where Y_Ntop (t·ha⁻¹) is the yield hypothetically reached (according to response equation) if the plot had been fertilized; Y_real (t·ha⁻¹) is the actual yield in the plot, Price_grain is the rice grain price ($·t⁻¹); Price_N is the N (ammonium sulfate) price ($·kg⁻¹); N_rate (kg·ha⁻¹) is the N rate established for each scenario; and Cost_Napplication ($·ha⁻¹) 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⁻¹) (Equation (9)).

  $$\begin{array}{cc}IfDeltaN\ge 0;& Nrate=0\\ IfDeltaN<0;& Nrate=\left|DeltaN\right|\end{array}$$

  (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.

  $$\begin{array}{cc}IfDeltaN\ge 0;& Nrate=0\\ IfDeltaN<0\mathrm{and}\left|DeltaN\right|\ge N_m;& Nrate=\left|DeltaN\right|\\ IfDeltaN<0\mathrm{and}\left|DeltaN\right|<N_m;& Nrate=N_m\end{array}$$

  (10)
• 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)).

  $$\begin{array}{cc}IfDeltaN\ge 0;& Nrate=0\\ IfDeltaN<0;& Nrate=N_{fix}\end{array}$$

  (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⁻¹ and the net benefit and the N excess were represented graphically.

3. Results

3.1. 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. 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.

	R_RVI	R_GRVI	R_NDVI	R_GNDVI	R_MCARI1	R_MCARINIR	R_gMCARINIR
2012	y = 0.22 + 0.74x	y = 0.04 + 0.93x	y = −0.18 + 1.13x	y = 0.29 + 0.68x	y = −0.22 + 1.18x	y = 0.98x0.43	y = 1.00x0.33
2013	y = 0.22 + 0.75x	y = 0.14 + 0.84x	y = −0.39 + 1.36x	y = 0.22 + 0.76x	y = −0.17 + 1.12x	y = 0.99x0.48	y = 1.01x0.39
PS	y = 0.25 + 0.75x	y = 0.20 + 0.80x	y = −0.26 + 1.26x	y = 0.31 + 0.70x	y = −0.12 + 1.12x	y = 1.03x0.44	y = 1.04x0.35
M	y = 0.19 + 0.72x	y = −0.01 + 0.94x	y = −0.26 + 1.17x	y = 0.22 + 0.71x	y = −0.33 + 1.21x	y = 0.93x0.46	y = 0.95x0.37
Year	n.s.	n.s.	n.s.	n.s.	n.s.	n.s.	n.s
Fertilizer	***	***	***	**	***	***	***

n.s = not significant; ** p < 0.01; *** p < 0.001.

).

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).

3.1.2. 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. Coefficients of determination (R²) of the relationships between R_yield and R_VIs for years 2012 and 2013 and the pooled data, and RMSE for the pooled data.

	Model Type †	2012 n = 88	2013 n = 83	Pooled 2012+2013 n = 171
		R2	R2	R2	RMSE
R_RVI	L	0.70 ***	0.56 ***	0.62 ***	0.149
R_GRVI	L	0.74 ***	0.53 ***	0.61 ***	0.151
R_NDVI	L	0.74 ***	0.56 ***	0.63 ***	0.148
R_GNDVI	L	0.77 ***	0.56 ***	0.64 ***	0.144
R_MCARI1	L	0.69 ***	0.40 ***	0.52 ***	0.168
R_MCARINIR	M	0.74 ***	0.58 ***	0.64 ***	0.145
R_gMCARINIR	M	0.77 ***	0.61 ***	0.67 ***	0.139

*** p < 0.001. ^† L = linear, M = multiplicative.

). 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_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² = 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.

3.2. N Topdressing Recommendation Approach

3.2.1. 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):

$$DeltaN=2.843+211.491·DeltaNNI\left(GNDVI\right)$$

(12)

$$Delta N+200=207.291·{\left(Delta NNI\left(gMCAR{I}_{NIR}\right)+1\right)}^{0.517}$$

(13)

3.2.2. 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. 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).

Strategy	SUCCESS	EXCESS	DEFECT
R_GNDVI	83.3	12.5	4.2
R_gMCARINIR	87.5	4.2	8.3
R_GNDVI&R_gMCARINIR	87.5	4.2	8.3
R_GNDVI or R_gMCARINIR	83.3	12.5	4.2

). 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⁻¹ and 132 $·ha⁻¹ (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.

The reference scenario (current farmers’ practice), where all plots are fertilized, obtains the lowest net benefit, with a maximum of 84 $·ha⁻¹ for a rate of 75 kg N·ha⁻¹ (Figure 5a), and the highest N excess of the four scenarios reaching 45 kg N excess·ha⁻¹ 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⁻¹ and a N excess of 2.3 kg N·ha⁻¹. 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⁻¹ for a threshold N rate of 90 kg N·ha⁻¹ (Figure 5a). However, the increase of the net benefit (35 $·ha⁻¹) is counterbalanced with an increase in N excess (20 kg·N·ha⁻¹ for the highest benefit versus 2.3 kg N·ha⁻¹ in Scenario 1) (Figure 5b).

In the case of fixing a predefined N rate (Scenario 3), the maximum net benefit (132 $·ha⁻¹) is equal to that in Scenario 2 for the rate of 90 kg·ha⁻¹ (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⁻¹), the N excess was 20 kg·N·ha⁻¹, the same as in Scenario 2 (Figure 5b).

4. Discussion

4.1. Relationships between Yield and Vegetation Indices

4.1.1. 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.

4.1.2. 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⁻¹; however, yields reached 10,000 kg·ha⁻¹ in the study conducted by Cao et al. [16] and 14,000 kg·ha⁻¹ 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⁻¹ (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.

4.2. 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⁻¹, for a minimum N rate of 90 kg N·ha⁻¹, but the N excess increases in 18 kg N·ha⁻¹ (Figure 5). Combining net benefit and environmental impact, we would recommend a N rate of 60 kg N·ha⁻¹ with an increase in net benefit of 27 $·ha⁻¹ and an increase in N excess of 6.1 kg N·ha⁻¹ 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⁻¹) (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⁻¹ 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.

5. 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⁻¹).

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.