Performances of Vegetation Indices on Paddy Rice at Elevated Air Temperature, Heat Stress, and Herbicide Damage

: Spectral reﬂectance-based vegetation indices have sensitive characteristics to crop growth and health conditions. The performance of each vegetation index to a certain condition is di ﬀ erent and needs to be interpreted, correspondingly. This study aimed to assess the most suitable vegetation index to identify the crop response against elevated air temperatures, heat stress, and herbicide damage. The spectral reﬂectance, yield components, and growth parameters such as plant height, leaf area index (LAI), and above-ground dry matter of paddy rice, which was cultivated in a temperature gradient ﬁeld chamber to simulate global warming conditions, were observed from 2016 to 2018. The relationships between the vegetation indices and the crop parameters were assessed considering stress conditions. The normalized di ﬀ erence vegetation index (NDVI) represented the changes in plant height (R-square = 0.93) and the LAI (R-square = 0.901) before the heading stage. Furthermore, the NDVI and the cumulative growing degree days had a Sigmoid curve and an R-square value of 0.937 under the normal growth case, but it decreased signiﬁcantly in the herbicide damage case. This characteristic was useful for detecting the damaged crop growth condition. Additionally, to estimate the grain yield of paddy rice, the medium resolution imaging spectrometer (MERIS) terrestrial chlorophyll index was better: R-square = 0.912; root mean square error = 95.69 g / m 2 . Photochemical reﬂectance index was sensitive to physiological stress caused by the heatwave, and it decreased in response to extremely high air temperatures. These results will contribute towards determining vegetation indices under stress conditions and how to e ﬀ ectively utilize them.


Introduction
Remote sensing is a useful technique to directly monitor crop conditions non-destructively [1,2]. Crop spectral reflectance response varies depending on the distribution of crop biomass, crop conditions (e.g., disease, nutrition), and soil background [3][4][5][6][7]. This response constitutes the vegetation indices developed to represent the structural, biochemical, and physiological characteristics of vegetation [8]. Moreover, certain vegetation indices have been used to estimate the crop yield using an empirical relationship [9][10][11]. However, this method, with only one input of a vegetation index, has some limitations in various spatial and temporal applications [12,13]. Current crop yield models not only use meteorological variables but also the vegetation indices instead of the input parameters related to crop condition [14]. Thus, these crop models, unlike other crop models that only use meteorological variables and the fixed crop parameter values, can accurately represent the influence of the cultivation environment, such as agrometeorological disaster, disease stress, and nutriture

Meteorological Data
Four instrument shelters hanging in the middle of the TGFC measured the air temperature ( • C) in each TGFC. A temperature probe (107, Campbell Scientific Inc., Logan, UT, USA), which is a type of thermistor with a tolerance of approximately ±0.2 • C, was used to measure the air temperature in the TGFC. The four temperature probes and instrument shelters were installed, and the position of the four shelters denoted as AT, AT + 1, AT + 2, and AT + 3 • C. TGFC systems are automatically an effort to maintain these conditions using the heater and three exhaust fans. In this study, the positions were divided into 16 equal intervals in each TGFC, and they were labeled from position 1 (P1) up to position 16 (P16). The four positions located in instrument shelters were position 1 (P1), position 6 (P6), position 11 (P11), and position 16 (P16). All meteorological data were recorded by the CR10X data logger (Campbell Scientific Inc., Logan, UT, USA) every 5 min and averaged based on 4 s interval scan data. Solar radiation information such as solar insolation and photosynthetically active radiation (PAR) was measured by a pyranometer (LI-200R, LI-COR Inc., Lincoln, Nebraska, USA) and a quantum sensor (LI-190R, LI-COR Inc., Lincoln, Nebraska, USA). These radiation sensors were installed inside and outside of TGFC to consider TGFC transmittance. Figure 3 is an example of a time series for a daily mean air temperature and solar radiation for the three-year period in one TGFC. A solid blue line indicates the AT. Green dash-dot line, orange dot line, and red dash line show AT + 1, AT + 2, and AT + 3 • C, respectively. The difference in air temperature between two adjacent stations was not always 1 • C because these facilities were automatically controlled based on the AT and AT + 3 • C positions. Thus, analyses were conducted using cumulative air temperature in each position. The growing degree days (GDD) were calculated as follows: GDD = daily mean air temperature -base air temperature. Here, the daily mean air temperature was calculated based on data recorded in the data logger every 5 min. The base air temperature of the paddy rice is 10 • C [42]. GDD was accumulated for GDDs > 0 to evaluate the effect of air temperature under the same solar radiation condition in each position.
A heatwave occurred in this study area during the booting and heading stages in 2016, but this meteorological situation did not occur in 2017. According to the Korea Meteorological Administration (KMA), the worst heatwave in meteorological history happened in 2018. The heatwave lasted for one month from mid-July to mid-August in 2018. The highest solar radiation and air temperature were recorded during this period ( Figure 3c). Therefore, the response of crop due to high air temperature was analyzed during this heatwave period. Although the optimum daily mean air temperature of paddy rice for the growth and development changes depending on the cultivars and the growth stage [43], it is generally 25~32 • C [44].

Ground Measurements of Growth Parameters
Growth parameters such as plant height, LAI, above-ground dry matter (AGDM), and grain yield were measured to evaluate their responses to elevated air temperature, heat stress, and herbicide damage. Plant height was measured at all chambers every week. For three years, the number of chambers used for plant height was 8, and measurement positions were 22. Finally, the 227 observations of plant height were used in this study. LAI was measured by LAI-2000 (LI-190R, LI-COR Inc., Lincoln, Nebraska, USA) at all chambers every week. The number of chambers used for observing LAI observation was nine for the three year period. The number of positions and data on LAI was 26 and 279, respectively. The AGDM was measured to accurately comprehend the growth of paddy rice. In the first year (2016), AGDM was only measured at mature stage to analyze yield components of rice. The number of AGDM positions and data, during the three years, was 38 and 113, respectively. To measure AGDM, the root of paddy rice was removed, and the AGDM was dried at 70 • C in the dryer for four days to maintain a constant humidity. The yield components such as grain yield (g/m 2 ), AGDM (g), harvest index, number of panicles per area, number of spikelets per panicle, ripening ratio (%), and 1000 grain weight (g) were measured in 9 chambers and 38 positions during the harvesting period. Grain yield and ripening ratio among the yield components were analyzed with vegetation indices. The number of ground measurements for growth parameters in 2017 and 2018 was adjusted based on the 2016 results. Information on ground measurements is summarized in Table 1.  (3,15,20) 38 (3,15,20) LAI: Leaf area index; AGDM: Above-ground dry matter.

Optical Measurements for Vegetation Indices
Optical measurements were carried out to monitor the conditions of paddy rice. In this study, two types of devices were used to calculate vegetation indices. First, a field-spectrometer was used to measure spectral reflectance-based vegetation indices. Analytical spectral devices (ASD) Fieldspec 3 spectrometers (hereafter Fieldspec; Malvern Panalytical Ltd., Malvern, UK) were used in 2016. The wavelength range on Fieldspec spectrometer is 350~2500 nm. Spectral resolution and sampling interval in wavelength for calculating vegetation indices are 3 at 700 nm and 1.4 at 350~1050 nm. In 2017 and 2018, Avantes AvaSpec-UL2048L spectrometers (hereafter Avantes; Avantes BV, Apeldoorn, Netherlands) were used. The wavelength range of the Avantes spectrometer is 300 to 1100 nm, and this device observes visible and NIR wavelengths. The sampling interval is about 0.576 at 700 nm. Field of view (FOV) of the spectrometers is 23 • . In this study, two spectrometers were used to measure solar irradiance and radiance at solar noon time that has the minimum solar zenith angle during the day. The radiance was the value reflected from paddy rice, and the observation radius and areas of the radiance were about 0.41 m and 0.52 m 2 at a 1-m height from the target. Solar irradiance was the value reflected from a white reference panel (RS50 White Reflectance Standard, StellarNet Inc., Florida, USA; >97% reflectance between 300 and 1700 nm). Measurements were carried out inside the TGFC to determine the influence of the film as the micro-climate conditions inside the TGFC were different from the conditions outside [45]. After conducting an inter-calibration between two spectrometers, Remote Sens. 2020, 12, 2654 8 of 25 the reflectance at each wavelength was calculated from a ratio between radiance and solar irradiance. Vegetation indices were then computed.
Three vegetation indices were selected among various others based on three biophysical types, i.e., structural, biochemical, and physiological types ( Figure 4, Table 2). The NDVI represents the structure of vegetation; it usually uses red and NIR wavelengths. It has been thoroughly studied compared to other structural vegetation indices such as SAVI and OSAVI and measured using a satellite and drone. The EVI was excluded to avoid the influence of the weak transmittance of blue wavelength on the ETFE film of the TGFC. The MTCI was chosen as a structural-biochemical index. It is produced using the Sentinel-2 satellite and the data obtained can be utilized in precision agriculture [46,47]. Last, the PRI was selected as a biochemical-physiological vegetation index. We hypothesized that the PRI can detect the crop stress caused by heatwave because it is highly related to the photosynthetic performance [48,49]. Three vegetation indices were selected among various others based on three biophysical types, i.e., structural, biochemical, and physiological types ( Figure 4, Table 2). The NDVI represents the structure of vegetation; it usually uses red and NIR wavelengths. It has been thoroughly studied compared to other structural vegetation indices such as SAVI and OSAVI and measured using a satellite and drone. The EVI was excluded to avoid the influence of the weak transmittance of blue wavelength on the ETFE film of the TGFC. The MTCI was chosen as a structural-biochemical index. It is produced using the Sentinel-2 satellite and the data obtained can be utilized in precision agriculture [46,47]. Last, the PRI was selected as a biochemical-physiological vegetation index. We hypothesized that the PRI can detect the crop stress caused by heatwave because it is highly related to the photosynthetic performance [48,49].   Second, a photodiode-based spectral reflectance sensor (hereafter SRS; METER group Inc., Pullman, WA, USA) was used to continually monitor paddy rice at AT and AT+3 • C positions in 2017 and 2018. The NDVI and PRI were measured using SRS, and they were processed as daily mean values. In this process, data measured from 11 to 14 o'clock were used to utilize the seasonal signal of PRI that reflects the diurnal and seasonal changes [50]. The processed PRI was used to evaluate the response of paddy rice under the heatwave. SRS was installed in the TGFC ( Figure 5), and the distance between the Remote Sens. 2020, 12, 2654 9 of 25 sensor and plants was maintained consistently. SRS data used in this study was saved in the EM50 data logger (METER group Inc., Pullman, WA, USA) every one minute.
Remote Sens. 2020, 10, x FOR PEER REVIEW 9 of 26 Second, a photodiode-based spectral reflectance sensor (hereafter SRS; METER group Inc., Pullman, WA, USA) was used to continually monitor paddy rice at AT and AT+3 °C positions in 2017 and 2018. The NDVI and PRI were measured using SRS, and they were processed as daily mean values. In this process, data measured from 11 to 14 o'clock were used to utilize the seasonal signal of PRI that reflects the diurnal and seasonal changes [50]. The processed PRI was used to evaluate the response of paddy rice under the heatwave. SRS was installed in the TGFC ( Figure 5), and the distance between the sensor and plants was maintained consistently. SRS data used in this study was saved in the EM50 data logger (METER group Inc., Pullman, WA, USA) every one minute.

Time Series of NDVI, MTCI, and PRI Under Elevated Air Temperature
Time series of the vegetation indices, which reflect the elevated air temperature and measured using spectrometers, were evaluated during the whole cultivation period. The NDVI, MTCI, and PRI were selected as representative structural, structural-biochemical, and biochemical-physiological vegetation indices, respectively. These vegetation indices included the fraction of vegetation and each characteristic of vegetation indices. In 2018, NDVI increased after the transplanting stage until the heading stage (Figure 6a), and the gradient of NDVI was sharp at a high air temperature before the heading stage. The NDVIs at the eight positions were saturated in the heading stage, although they were exposed to different air temperatures. This saturation effect was known to occur in the vegetation index using a red wavelength such as NDVI [51,52]. After the heading stage, the NDVI decreased at the P1 (AT) position, but it had consistent value at P15 (AT+2.37 °C). The time series of MTCI sensitively expressed the growth of paddy rice ( Figure 6b). The value of MTCI at P15 (AT+2.37 °C) was higher than that at P1 (AT). Unlike NDVI, the MTCIs exposed to various air temperatures were largely different among each temperature gradient condition in the heading stage because the red-edge band, which was used to calculate MTCI, tended to reduce saturation [53]. After the heading stage, the MTCI decreased at all positions. The PRI showed variability in the first optical measurement due to background effects, such as soil and water (Figure 6c), and it rapidly increased after that. The PRI decreased after the heading stage at the AT position due to senescence of the leaf [37]. However, decreasing trends in PRI occurred under high air temperatures before the heading stage, unlike NDVI and MTCI. According to Cao et al. [53], PRI was the most sensitive to heat stress among 17 published vegetation indices; this sensitivity is due to the relation between PRI and photosynthesis activity. In addition, given that seasonal change of PRI is associated with the ratio of

Time Series of NDVI, MTCI, and PRI Under Elevated Air Temperature
Time series of the vegetation indices, which reflect the elevated air temperature and measured using spectrometers, were evaluated during the whole cultivation period. The NDVI, MTCI, and PRI were selected as representative structural, structural-biochemical, and biochemical-physiological vegetation indices, respectively. These vegetation indices included the fraction of vegetation and each characteristic of vegetation indices. In 2018, NDVI increased after the transplanting stage until the heading stage (Figure 6a), and the gradient of NDVI was sharp at a high air temperature before the heading stage. The NDVIs at the eight positions were saturated in the heading stage, although they were exposed to different air temperatures. This saturation effect was known to occur in the vegetation index using a red wavelength such as NDVI [51,52]. After the heading stage, the NDVI decreased at the P1 (AT) position, but it had consistent value at P15 (AT + 2.37 • C). The time series of MTCI sensitively expressed the growth of paddy rice ( Figure 6b). The value of MTCI at P15 (AT + 2.37 • C) was higher than that at P1 (AT). Unlike NDVI, the MTCIs exposed to various air temperatures were largely different among each temperature gradient condition in the heading stage because the red-edge band, which was used to calculate MTCI, tended to reduce saturation [53]. After the heading stage, the MTCI decreased at all positions. The PRI showed variability in the first optical measurement due to background effects, such as soil and water (Figure 6c), and it rapidly increased after that. The PRI decreased after the heading stage at the AT position due to senescence of the leaf [37]. However, decreasing trends in PRI occurred under high air temperatures before the heading stage, unlike NDVI and MTCI. According to Cao et al. [53], PRI was the most sensitive to heat stress among 17 published vegetation indices; this sensitivity is due to the relation between PRI and photosynthesis activity. In addition, given that seasonal change of PRI is associated with the ratio of chlorophyll and carotenoid [49,54], PRI can accurately detect the stress. Chou et al. [55] showed that the PRI had a relatively low value under the water stress and it was related to the change in the chlorophyll-to-carotenoid ratio. However, the chlorophyll content at the leaf scale did not significantly change. Given that the MTCI, which is related to chlorophyll, did not decrease under the extremely high air temperature in this study (Figure 6b), the PRI decrease is assumed to be due to the carotenoid change rather than chlorophyll. Carotenoid levels can increase to protect plants from heat stress [56]. The PRI had a positive correlation with the chlorophyll/carotenoid ratio, which decreased under heat shock [55,57].
Remote Sens. 2020, 10, x FOR PEER REVIEW 10 of 26 chlorophyll and carotenoid [49,54], PRI can accurately detect the stress. Chou et al. [55] showed that the PRI had a relatively low value under the water stress and it was related to the change in the chlorophyll-to-carotenoid ratio. However, the chlorophyll content at the leaf scale did not significantly change. Given that the MTCI, which is related to chlorophyll, did not decrease under the extremely high air temperature in this study (Figure 6b), the PRI decrease is assumed to be due to the carotenoid change rather than chlorophyll. Carotenoid levels can increase to protect plants from heat stress [56]. The PRI had a positive correlation with the chlorophyll/carotenoid ratio, which decreased under heat shock [55,57]. In this study, vegetation indices in the four positions, P1 (AT), P5 (AT+0.91), P11 (AT+1.71), and P13 (AT+2.00 °C), were analyzed during the heatwave period using detailed air temperature data (Figure 6d,e,f). The four positions were selected to consider the change of spectral reflectance. The spectral reflectance in visible wavelength decreased during the period from July 17 2018 to August 12 2018 in all positions. However, the spectral reflectance in the NIR wavelength was different. The NIR increased in P1 and P5 from July 17 2018 to August 12 2018 due to the growth of paddy rice. However, NIR wavelength had a consistent value for about one month in P11; this value decreased in P13 when exposed to the extremely high air temperatures. As a result, NDVI in P1 and P5 increased, and the value in P11 remained constant during the heatwave period. Although the NIR decreased in P13, NDVI was similar to NDVIs in the other positions during the heatwave period due to a low reflectance of the red band. In addition, MTCI showed a similar trend to the NDVI. Conversely, PRI exposed to high air temperatures decreased faster than that in AT position (Figure 6f). The first PRI decrease during the heatwave period (from July 10 2018 to August 17 2018) occurred in P13, with a daily mean and maximum air temperature of 33.35 and 35.69 °C, respectively. In P1 and P5, the PRI decreased after the heading stage, and the daily mean air temperature was below 32.0 °C. Positions were divided into 16 at equal intervals in each temperature gradient field chamber (TGFC) and they were named as position 1 (P1)~position 16 (P16). P1 indicates the time series of vegetation indices exposed to ambient temperature (AT), and P15 indicates the time series of vegetation indices exposed to high air temperature (about AT+2.37 • C). P5 and P11 were averagely higher temperatures as +0.91 and +1.71 • C than AT.
In this study, vegetation indices in the four positions, P1 (AT), P5 (AT + 0.91), P11 (AT + 1.71), and P13 (AT + 2.00 • C), were analyzed during the heatwave period using detailed air temperature data (Figure 6d,e,f). The four positions were selected to consider the change of spectral reflectance. The spectral reflectance in visible wavelength decreased during the period from July 17 2018 to August 12 2018 in all positions. However, the spectral reflectance in the NIR wavelength was different. The NIR increased in P1 and P5 from July 17 2018 to August 12 2018 due to the growth of paddy rice. However, NIR wavelength had a consistent value for about one month in P11; this value decreased in P13 when exposed to the extremely high air temperatures. As a result, NDVI in P1 and P5 increased, and the value in P11 remained constant during the heatwave period. Although the NIR decreased in P13, NDVI was similar to NDVIs in the other positions during the heatwave period due to a low reflectance of the red band. In addition, MTCI showed a similar trend to the NDVI. Conversely, PRI exposed to high air temperatures decreased faster than that in AT position (Figure 6f). The first PRI decrease during the heatwave period (from July 10 2018 to August 17 2018) occurred in P13, with a daily mean and maximum air temperature of 33.35 and 35.69 • C, respectively. In P1 and P5, the PRI decreased after the heading stage, and the daily mean air temperature was below 32.0 • C.

Performance of PRI Under the Heat Stress
PRI and NDVI were continually observed using SRS at AT and AT + 3 • C in 2017 and 2018. These two vegetation indices were processed based on 1 day. SRS-based PRI (hereafter PRI SRS ) in 2018 was analyzed because no heatwave occurred in 2017. SRS-based NDVI did not analyze because it did not respond to the heatwave like the results of spectrometers. Figure 7 shows the time series of PRI SRS at two positions in 2018. PRI SRS increased after the transplanting stage as PRI was measured by the spectrometer. PRI SRS at AT + 3 • C was higher than that at AT before the heading stage. The decreasing time of PRI SRS at AT appeared in the heading stage, but at AT + 3 • C, the decrease occurred early on the heatwave before the heading stage. This pattern of PRI SRS was most noticeable in the AT + 3 • C position of chamber A; PRI SRS did not increase in the AT + 3 • C position of chamber B. In other words, the decreasing characteristics of PRI in plot spatial scales showed the same behavior under heat stress as the previous study PRI characteristics at leaf spatial scale [53].

Performance of PRI Under the Heat Stress
PRI and NDVI were continually observed using SRS at AT and AT+3 °C in 2017 and 2018. These two vegetation indices were processed based on 1 day. SRS-based PRI (hereafter PRISRS) in 2018 was analyzed because no heatwave occurred in 2017. SRS-based NDVI did not analyze because it did not respond to the heatwave like the results of spectrometers. Figure 7 shows the time series of PRISRS at two positions in 2018. PRISRS increased after the transplanting stage as PRI was measured by the spectrometer. PRISRS at AT+3 °C was higher than that at AT before the heading stage. The decreasing time of PRISRS at AT appeared in the heading stage, but at AT+3 °C, the decrease occurred early on the heatwave before the heading stage. This pattern of PRISRS was most noticeable in the AT+3 °C position of chamber A; PRISRS did not increase in the AT+3 °C position of chamber B. In other words, the decreasing characteristics of PRI in plot spatial scales showed the same behavior under heat stress as the previous study PRI characteristics at leaf spatial scale [53]. PRI normally decreases after the heading stage due to the leaf senescence [37,58]. Thus, this study assumes that the difference (hereafter, delta date; unit: Days) between the decline date in PRI and the heading stage date was related to heat stress. The delta date includes the error because the optical measurement using the spectrometer was averagely conducted once per week. Figure 8 shows the relationship between cumulative hours when the air temperature is above 35 °C and delta date. Here, 35 °C is the threshold air temperature of paddy rice for heat stress [59]. The data collected from spectrometer-based PRI and PRISRS from 2016 to 2018 were used in this relationship with cumulative hours. The correlation value in this relationship was -0.75 (p-value < 0.0001). The decline in PRI started early when the heat stress was maintained for a long time. This means that PRI responded to heat stress. PRI normally decreases after the heading stage due to the leaf senescence [37,58]. Thus, this study assumes that the difference (hereafter, delta date; unit: Days) between the decline date in PRI and the heading stage date was related to heat stress. The delta date includes the error because the optical measurement using the spectrometer was averagely conducted once per week. Figure 8 shows the relationship between cumulative hours when the air temperature is above 35 • C and delta date. Here, 35 • C is the threshold air temperature of paddy rice for heat stress [59]. The data collected from spectrometer-based PRI and PRI SRS from 2016 to 2018 were used in this relationship with cumulative hours. The correlation value in this relationship was −0.75 (p-value < 0.0001). The decline in PRI started early when the heat stress was maintained for a long time. This means that PRI responded to heat stress.

Relationship between Cumulative Growing Degree Days and Plant Height
The paddy rice in 2017 suffered due to the herbicide treatment; the abnormal growth and development on paddy rice were evaluated using cumulative GDD. Cumulative GDD is known to relate to LAI, dry matter, and grain yield, and it has been used to simulate crop growth in crop models [13,14,60]. Plant height (number of samples = 167) remains constant after the heading stage, so the relationship between cumulative GDD and plant height has a sigmoid trend (Figure 9). The coefficient of determination (hereafter R-square) was 0.958 (p-value < 0.0001), and the root mean square error (RMSE) was 5.82 cm under the normal growth states. These results were similar to the ones from previous studies, where the R-square between cumulative GDD and plant height was 0.88-0.91 for Sesamum indicum [61] and 0.99 for sweet corn [62]; these relationships follow a Sigmoid growth curve like the one obtained in this study. In general, the growth and development of paddy rice were fast under high air temperatures [63], and the relationship between cumulative GDD and plant height was stable.
The LAI has been estimated using GDD that is one of the variables in the rice growth and productivity model. The R-square (RMSE) was 0.97 (0.47) [64]. Williams and Lindquist [62] reported an R-square (RMSE) of 0.92-0.98 (0.0224-0.2138) for two years. In this study, R-square between cumulative GDD and LAI before the heading stage was 0.861 (p-value < 0.0001). These results mean that cumulative GDD has a high correlation with growth parameters of plant height and LAI under normal growth conditions. However, the relationship between cumulative GDD and LAI under the abnormal growth conditions had different trends with that under the normal growth condition. The growth and development of paddy rice were late due to herbicide influence. Thus, the relationship between cumulative GDD and growth parameters can effectively express the low and late growth and development of paddy rice due to herbicide damage.
This concept can be applied to other abnormal growth conditions as well as the effect of herbicide. For example, the damages of growth and development caused by disease and harmful insects can be detected using cumulative GDD and growth parameters. Conversely, rapid growth and development of crops due to the bountiful fertilizer can be detected because of the relatively high values of growth parameters compared to the normal relationship between cumulative GDD and

Relationship between Cumulative Growing Degree Days and Plant Height
The paddy rice in 2017 suffered due to the herbicide treatment; the abnormal growth and development on paddy rice were evaluated using cumulative GDD. Cumulative GDD is known to relate to LAI, dry matter, and grain yield, and it has been used to simulate crop growth in crop models [13,14,60]. Plant height (number of samples = 167) remains constant after the heading stage, so the relationship between cumulative GDD and plant height has a sigmoid trend (Figure 9). The coefficient of determination (hereafter R-square) was 0.958 (p-value < 0.0001), and the root mean square error (RMSE) was 5.82 cm under the normal growth states. These results were similar to the ones from previous studies, where the R-square between cumulative GDD and plant height was 0.88-0.91 for Sesamum indicum [61] and 0.99 for sweet corn [62]; these relationships follow a Sigmoid growth curve like the one obtained in this study. In general, the growth and development of paddy rice were fast under high air temperatures [63], and the relationship between cumulative GDD and plant height was stable.
The LAI has been estimated using GDD that is one of the variables in the rice growth and productivity model. The R-square (RMSE) was 0.97 (0.47) [64]. Williams and Lindquist [62] reported an R-square (RMSE) of 0.92-0.98 (0.0224-0.2138) for two years. In this study, R-square between cumulative GDD and LAI before the heading stage was 0.861 (p-value < 0.0001). These results mean that cumulative GDD has a high correlation with growth parameters of plant height and LAI under normal growth conditions. However, the relationship between cumulative GDD and LAI under the abnormal growth conditions had different trends with that under the normal growth condition. The growth and development of paddy rice were late due to herbicide influence. Thus, the relationship between cumulative GDD and growth parameters can effectively express the low and late growth and development of paddy rice due to herbicide damage. growth parameter. Although the relationship between cumulative GDD and growth parameters was useful to detect plant stress on the crop growth and development, this method for detecting crop stress conditions needs micro-meteorological data, as well as crop growth information, and its application is limited to small areas. Thus, the relationships between growth parameters and vegetation indices under normal and abnormal growth states were evaluated in order to replace growth parameters with vegetation indices before the heading stage.

Relationship Between Plant Height and Vegetation Indices
Plant height increases after the transplanting stage, and then reaches its maximum at the heading stage. Conversely, vegetation indices increase until the heading stage, and then they start to decrease. For this reason, the relationship between plant height and vegetation indices was analyzed using data before the heading stage (number of samples = 121). The highest value of the R-square among the three vegetation indices was 0.938 in NDVI (Figure 10a), and the lowest value as 0.832 in PRI ( Figure  10c). The relationship between plant height and MTCI was linear (Figure 10b). All vegetation indices accurately explained the plant height until the heading stage. Structural vegetation index such as NDVI was saturated in high values of plant height, but MTCI could sensitively express high plant height values. In the literature, the plant height of durum wheat had R-square values of 0.518~0.828 and root mean squared error values of 0.083~0.105 in the NDVI. In the OSAVI the same values were 0.519~0.817 and 0.102~0.127, respectively [65]. This result was lower than ours because it included the data measured after heading. In the ripening phase, the plant height has a constant value, but the value of vegetation index decreases. Thus, the estimation of the plant height of paddy rice using vegetation indices before the heading stage was recommended.
The influence of data related to the abnormal growth was evaluated in the relationship between plant height and vegetation indices. When data on abnormal growth (number of samples = 38) were included with the normal growth data, the R-square decreased approximately 0.003~0.015 in three vegetation indices (Table 3). Unlike PRI, the R-square values in NDVI and MTCI decreased by just This concept can be applied to other abnormal growth conditions as well as the effect of herbicide. For example, the damages of growth and development caused by disease and harmful insects can be detected using cumulative GDD and growth parameters. Conversely, rapid growth and development of crops due to the bountiful fertilizer can be detected because of the relatively high values of growth parameters compared to the normal relationship between cumulative GDD and growth parameter. Although the relationship between cumulative GDD and growth parameters was useful to detect plant stress on the crop growth and development, this method for detecting crop stress conditions needs micro-meteorological data, as well as crop growth information, and its application is limited to small areas. Thus, the relationships between growth parameters and vegetation indices under normal and abnormal growth states were evaluated in order to replace growth parameters with vegetation indices before the heading stage.

Relationship Between Plant Height and Vegetation Indices
Plant height increases after the transplanting stage, and then reaches its maximum at the heading stage. Conversely, vegetation indices increase until the heading stage, and then they start to decrease. For this reason, the relationship between plant height and vegetation indices was analyzed using data before the heading stage (number of samples = 121). The highest value of the R-square among the three vegetation indices was 0.938 in NDVI (Figure 10a), and the lowest value as 0.832 in PRI (Figure 10c). The relationship between plant height and MTCI was linear (Figure 10b

Relationship between Above-Ground Dry Matter and Vegetation Indices
AGDM indicates biomass information of paddy rice, and it increased continually until the harvesting stage. The characteristics of stem, leaf, and panicle in the AGDM were separated into before and after the heading stage. AGDM on stem and leaf increased before the heading stage. After the heading stage, AGDM of leaf decreased, and the panicle AGDM increased. In addition, the stem's weight decreased due to assimilating translocation from stem to panicle. The ratio of panicle weight and all AGDM was called the harvest index. Given that the AGDM and harvest index are directly related to yield, the estimation of AGDM is critical. According to previous studies, SR and NDVI have an exponential relationship with the paddy rice AGDM, with an R-square of 0.68 and 0.64 (p- The influence of data related to the abnormal growth was evaluated in the relationship between plant height and vegetation indices. When data on abnormal growth (number of samples = 38) were included with the normal growth data, the R-square decreased approximately 0.003~0.015 in three vegetation indices (Table 3). Unlike PRI, the R-square values in NDVI and MTCI decreased by just 0.008 and 0.003, and they were similar regardless of stress, which means that vegetation indices, especially related to the structure of the crop, include the abnormal growth response of paddy rice. Table 3. The relationship between vegetation indices and cumulative growing degree days (GDD), plant height, leaf area index (LAI) or above-ground dry matter (AGDM) before the heading stage. Difference indicates the difference of R-square between "normal growth samples" and "normal growth + abnormal growth samples". Li, Lo, S, and E are abbreviations for linear, logarithm, Sigmoid, and exponential function (p-value < 0.0001).

Relationship between Leaf Area Index and Vegetation Indices
LAI that indicates 3-dimension crop growth information had a high correlation with plant height that indicates 2-dimensional crop growth information. The correlation coefficient for this relationship between plant height and LAI was 0.946 (p-value < 0.0001), and the equation on linear trend was "plant height (cm) = 19.802 × LAI (m 2 /m 2 ) + 28.566" (number of samples = 158). This paddy rice result was similar to the garlic result with an R-square of 0.96 [66]. Unlike the plant height, LAI increased after transplanting and then decreased after the heading stage [67].
The relationship between LAI and vegetation indices was evaluated (number of samples = 155), excluding the abnormal growth. In this study, NDVI and PRI maintained an exponential relationship with LAI (Figure 10d,f). MTCI, which include biochemical characteristics, had a linear relationship with LAI (Figure 10e). The NDVI relationship showed the highest R-square (0.903; p-value < 0.0001), similar to its relationship with plant height. The relationship between LAI and MTCI was a linear relationship like with plant height, and it was effective for estimation of high LAI values [68]. Structural vegetation indices had a non-linear relationship with LAI, but they had a higher R-square than the ones in previous studies. The obtained R-square for the LAI were (1) 0.87, 0.91, and 0.87 for the SR, NDVI, and PRI correlations [69]; (2) 0.98, 0.99, and 0.99 for the SR, NDVI, and OSAVI correlations [67]; (3) 0.94 and 0.92 for the NDVI and OSAVI correlations [70]. Most previous studies have used NDVI to estimate LAI; NDVI was also found to be the best variable to predict the LAI in this study. Given that NDVI was saturated when the canopy was closed (LAI > about 2 m 2 /m 2 ), vegetation index using the red-edge band (e.g., MTCI) can be used to diagnose the crop growth, because it does not saturate under high LAI values [68,71].
Data on abnormal growth (number of samples = 23) did not affect the relationship between LAI and vegetation indices that include structural characteristics such as the NDVI and the MTCI. The R-square in the NDVI and the MTCI decreased from 0.002 to 0.004. However, the values of LAI and vegetation indices under the abnormal growth were lower than those under normal growth because of the slow growth and development of paddy rice. Conversely, the R-square increased by 0.017 for PRI. These results mean that the relationship between LAI and vegetation indices remained similar independent of the stress conditions; this was the same as the plant height correlations behavior. Vegetation indices indicated which of the growth conditions were affected by stress and environmental factors. Accordingly, structural variables such as plant height and LAI can be replaced with vegetation indices.

Relationship between Above-Ground Dry Matter and Vegetation Indices
AGDM indicates biomass information of paddy rice, and it increased continually until the harvesting stage. The characteristics of stem, leaf, and panicle in the AGDM were separated into before and after the heading stage. AGDM on stem and leaf increased before the heading stage. After the heading stage, AGDM of leaf decreased, and the panicle AGDM increased. In addition, the stem's weight decreased due to assimilating translocation from stem to panicle. The ratio of panicle weight and all AGDM was called the harvest index. Given that the AGDM and harvest index are directly related to yield, the estimation of AGDM is critical. According to previous studies, SR and NDVI have an exponential relationship with the paddy rice AGDM, with an R-square of 0.68 and 0.64 (p-value < 0.0001) for SR and NDVI [72]. The relationship (p-value < 0.0001) between dry biomass and NDVI or OSAVI in barley shows an exponential regression with an R-square of 0.61 (NDVI) and 0.68 (OSAVI) [73]. Zheng et al. [74] reported the statistical results on the relationship between AGDM and structural vegetation indices such as the NDVI, green NDVI, and OSAVI before the heading stage. R-square values were 0.62 (NDVI), 0.66 (green NDVI), and 0.62 (OSAVI).
In this study, the number of samples on AGDM (number of samples = 22) was lower than that of other growth parameters because AGDM values can only be obtained using a destructive method. In addition, the number of samples matched with vegetation indices on the same day is lower than the number of AGDM measurements. In this study, NDVI was the best vegetation index to estimate AGDM before the heading stage (Figure 10g). The R-squares (p-value < 0.0001) for data on normal growth and data including abnormal growth samples were 0.914 and 0.910 in NDVI. MTCI had the second-best result with AGDM ( Figure 10h); the R-square was 0.863 (p-value < 0.0001), and it decreased to 0.816 (p-value < 0.0001) when the abnormal growth data (number of samples = 15) were included. PRI showed an exponential and linear relationship with AGDM ( Figure 10i). The statistic results of these vegetation indices were lower than the others. In addition, when abnormal growth data were added, degrees of change were relatively large, like MTCI. The relationship between AGDM and vegetation indices before the heading stage reflected the fraction effect of the crop before canopy closure [72].

Relationship Between Cumulative Growing Degree Days and Vegetation Indices
The relationships between the growth parameters of plant height, LAI, and AGDM and the vegetation indices were evaluated. The relationship between the two factors was divided into normal growth and abnormal growth conditions in vegetative and reproductive period. According to previous results, the growth parameters can be estimated using the structural vegetation index and MTCI, regardless of stress before the heading stage. Thus, the growth parameters were replaced with the vegetation index to detect abnormal growth of the paddy rice.
NDVI, MTCI, and PRI had a Sigmoid, logarithm, and Sigmoid trend with cumulative GDD (Figure 10j,k,l). The R-square (p-value < 0.0001) is 0.937 in the relationship between cumulative GDD and NDVI under the normal growth conditions (number of samples = 151), and it is the highest among the three indices. However, statistical results decreased when the abnormal growth data (number of samples = 41) was included. The trend in the abnormal growth data was different from that in the normal growth data. NDVI and PRI were effective in detecting the abnormal growth of paddy rice. The PRI statistics results were lower than those in NDVI because PRI includes physiological and biochemical characteristics.
Kimball et al. [75] reported the curve between cumulative GDD and NDVI during the cultivation period on wheat that was cultivated under the different air temperature conditions. The curves of cumulative GDD-NDVI were stable regardless of the different air temperature, but they were different when the crop was damaged by frost. The curve of cumulative GDD-NDVI in the frost-damaged samples was different with that in the heated samples. As a result, the grain yield of wheat damaged from frost was close to zero. It meant that the curve of cumulative GDD and vegetation indices could detect the abnormal status.

Estimation of Grain Yield Using Vegetation Indices
The grain yield of paddy rice has been estimated using vegetation indices. Vegetation indices are directly utilized as input data for the crop model, or they are used to estimate specific variables to use as input data for crop models, such as LAI [13]. In this study, vegetation indices were evaluated using the maximum vegetation index value during the cultivation period to estimate grain yield [9][10][11], because this study focuses on finding a suitable vegetation index instead of in the accurate estimation of grain yield. In addition, the ripening ratio of paddy rice was considered. Finally, the relationship between the maximum vegetation index value during the cultivation period multiplied by the ripening ratio (hereafter maximum vegetation index × ripening ratio) and grain yield of paddy rice was evaluated using three years of data ( Figure 11). The samples on yield components were divided into common condition, herbicide damage, and heat stress samples to confirm the effect of stress. Additionally, the sPRI was used instead of PRI in order to remove the negative value, and it was calculated as (PRI+1)/2.
The R-square between the vegetation index × ripening ratio and the grain yield of paddy rice on all data ranged from 0.861 to 0.912 (p-value < 0.0001) ( Table 4). MTCI had the best statistical results. Additionally, grain yield could be estimated after excluding the heat stress samples. The R-square values decreased in the three vegetation indices because the influence of ripening ratio on the estimation of grain yield was larger than that of the vegetation indices. In particular, the maximum NDVI and sPRI could not be used to estimate the grain yield. The maximum MTCI × ripening ratio had R-square values of 0.481 (p-value = 0.0002). The RMSE in MTCI was 96.04 g/m 2 when heat stress samples were excluded. divided into common condition, herbicide damage, and heat stress samples to confirm the effect of stress. Additionally, the sPRI was used instead of PRI in order to remove the negative value, and it was calculated as (PRI+1)/2. Figure 11. Relationship between maximum vegetation index value during cultivation period × ripening ratio and grain yield.
The R-square between the vegetation index × ripening ratio and the grain yield of paddy rice on all data ranged from 0.861 to 0.912 (p-value < 0.0001) ( Table 4). MTCI had the best statistical results. Actually, grain yield dependence on the degree of stress was significant for p-value < 0.1 (Figure 12a). The grain yield in samples with herbicide damage was lower than that in samples in common conditions, and the yield in heat stress samples was notably low. The significant difference was assessed in NDVI, MTCI, and sPRI considering the ripening ratio. The maximum NDVI and sPRI × ripening ratio showed no significant difference between the normal growth and abnormal growth data (Figure 12b,d). Heat stress-induced damage of grain yield was reflected by the NDVI and sPRI.
Conversely, the maximum MTCI × ripening ratio had a significant difference depending on the degree of stress for p-value < 0.1.
Remote Sens. 2020, 10, x FOR PEER REVIEW  21 of 26 can be estimated using mean air temperature (maximum air temperature) for 40 days after the heading stage (during the flowering period) [81,82].

Conclusions
The phenology of paddy rice was monitored using structural, biochemical, and physiological vegetation indices under the elevated air temperature conditions for three years. The extremely high temperature-induced heat stress occurred in 2016 and 2018, and the herbicide damage occurred in 2017. The most suitable vegetation indices were selected based on their performance under the future meteorological environment, physiological stress, and herbicide damage. The conclusions in the study were summarized as follows: (1) The cumulative GDD had a stable relationship with the plant height and LAI under normal growth conditions; (2) The relationship between vegetation indices and growth parameters such as plant height, LAI, and AGDM was constant regardless of abnormal growth data until the heading stage. This means that vegetation indices can replace growth parameters; These results corresponded with previous studies. Kanke et al. [51] showed that red-edge-based vegetation indices were better than red-based vegetation indices to predict grain yield. Zhang and Liu [76] discussed that MTCI measured at the heading stage was the best indicator to predict grain yield (R-square = 0.652). In addition, MTCI is the best vegetation index to predict leaf nitrogen concentration in rice [77]. Although red-edge-based vegetation indices were useful, the NDVI was the selected representative vegetation index to estimate the crop grain yield of the crop in numerous studies [78][79][80]. NDVI at the specific growth stage of crops was utilized to predict the grain yield instead of maximum NDVI. Grain yield of rice was estimated using NDVI at the booting and heading stages (R-square = 0.76) [78], and it had a significant correlation with NDVI measured at 91 days after seeding (R-square = 0.605) [79]. In a corn crop, the relationship between grain yield and NDVI at the V8 stage was high (R-square = 0.77) [80]. However, to use only vegetation index values is not suitable for estimating the grain yield. Although biomass increases, grain yield can decrease due to external factors such as disease infection, heat stress, and lodging [51]. In this study, the difference of grain yield due to herbicide damage was detected by the vegetation indices; the decreased grain yield caused by physiological stress was not detected by the vegetation indices. If the decreasing of the ripening ratio, which happened due to heat stress, was not considered, none of the vegetation indices could have estimated the grain yield of paddy rice. MTCI in AT+3 • C was higher than that in AT, but actual grain yield was lower due to the low ripening ratio. Therefore, an additional study to estimate the ripening ratio of paddy rice using vegetation index is demanded, or meteorological variables should be utilized to the ripening ratio or fertility percentage. The ripening ratio (fertility percentage) can be estimated using mean air temperature (maximum air temperature) for 40 days after the heading stage (during the flowering period) [81,82].

Conclusions
The phenology of paddy rice was monitored using structural, biochemical, and physiological vegetation indices under the elevated air temperature conditions for three years. The extremely high temperature-induced heat stress occurred in 2016 and 2018, and the herbicide damage occurred in 2017. The most suitable vegetation indices were selected based on their performance under the future meteorological environment, physiological stress, and herbicide damage. The conclusions in the study were summarized as follows: (1) The cumulative GDD had a stable relationship with the plant height and LAI under normal growth conditions; (2) The relationship between vegetation indices and growth parameters such as plant height, LAI, and AGDM was constant regardless of abnormal growth data until the heading stage. This means that vegetation indices can replace growth parameters; (3) The structural vegetation index, i.e., NDVI, can detect the slow growth and development of the crop caused by stress using the relationship with cumulative GDD. This concept can explain the abnormal condition of the crop, even if measurement data for several decades is not available; (4) PRI is the most useful vegetation index under physiological stress caused by heat stress. Unlike other vegetation indices, PRI decreased under the extremely high air temperatures even before the heading stage. When PRI decreases before the heading stage, this phenomenon might be the effect of the physiological stress and can then take the following actions: (1) Put silicate and potassium fertilizer; (2) let heatwave-induced the heated water flow to outside, and supply new water in order to lower the water temperature; (5) MTCI, which uses a red-edge band, is one of the most useful indices to predict the grain yield.
MTCI showed higher sensitivity to the growth and development of paddy rice than structural vegetation indices because it did not saturate at the heading stage. In addition, the MTCI did not decrease under the extremely high air temperatures, indirectly indicating that the chlorophyll levels of paddy rice did not decrease before the heading stage.
In other words, structural vegetation indices are useful to distinguish the normal and abnormal growth and development states (Category 1). Although structural damage does not exist, the crop stress under the extremely high air temperature is detected through biochemical-physiological vegetation index (Category 2). Lastly, the structural-biochemical vegetation index is more effective to grasp the crop growth and development than the vegetation index which has the only structural characteristic. These results suggest that it is important to select a reasonable vegetation index in accordance with the development intention of the index and the purpose of its application.
In this study, the relationship between growth parameters and vegetation indices was assessed using many observation data under various environments of air temperature. These results can be helpful to remote sensing researchers in the agricultural domain for efficient monitoring of crop growth as well as indicating crop stress. However, although this study investigates not only the vegetative phase but also the grain yield, it focused on the performance of vegetation indices before the heading stage under the various environmental conditions. In further research, the temporal response of paddy rice after the heading stage will be analyzed considering different states of rice panicle, particularly, the ripening ratio.
Author Contributions: All authors contributed this paper. J.-H.R., the main author, designed the research, measured and analyzed the remote sensing data, and wrote the manuscript; H.J., the co-author, reviewed the paper and contributed to the discussion; J.C., the corresponding author, designed the research and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.