Evaluation of a Phenology-Dependent Response Method for Estimating Leaf Area Index of Rice Across Climate Gradients

Abstract

Accurate estimate of the seasonal leaf area index (LAI) in croplands is required for understanding not only intra- and inter-annual crop development, but also crop management. Lack of consideration in different growth phases in the relationship between LAI and vegetation indices (VI) often results in unsatisfactory estimation in the seasonal course of LAI. In this study, we partitioned the growing season into two phases separated by maximum VI (${\mathrm{VI}}_{\max }$) and applied the general regression model to the data gained from two phases. As an alternative method to capture the influence of seasonal phenological development on the LAI-VI relationship, we developed a consistent development curve method and compared its performance with the general regression approaches. We used the Normalized Difference VI (NDVI) and the Enhanced VI (EVI) from the rice paddy sites in Asia (South Korea and Japan) and Europe (Spain) to examine its applicability across different climate conditions and management cycles. When the general regression method was used, separating the season into two phases resulted in no better estimation than the estimation obtained with the entire season observation due to an abrupt change in seasonal LAI occurring during the transition between the before and after ${\mathrm{VI}}_{\max }$. The consistent development curve method reproduced the seasonal patterns of LAI from both NDVI and EVI across all sites better than the general regression method. Despite less than satisfactory estimation of a local ${\mathrm{LAI}}_{\max }$, the consistent development curve method demonstrates improvement in estimating the seasonal course of LAI. The method can aid in providing accurate seasonal LAI as an input into ecological process-based models.

1. Introduction

Leaf area index (LAI) is one of the key parameters in estimating ecosystem productivity of various process-based models and is strongly related to plant phenology and vegetation dynamics [1,2,3]. LAI influences many biological and physical processes driving the exchange of matter and energy flow [4]. LAI serves as a useful indicator to characterize the condition of vegetation owing to its rapid response to different stress factors and changes in climatic conditions [5]. Therefore, estimation of LAI is an essential step in most of the process-based models for carbon and water fluxes in vegetative ecosystems [6,7,8]. Use of inaccurately estimated LAI as an input variable for process-based models will propagate errors in estimating ${\mathrm{CO}}_2$ and ${\mathrm{H}}_2$O exchange in vegetative ecosystems. Kucharik et al. (1997) [9] found that 45% underestimation of LAI led to 39% underestimation of ${\mathrm{CO}}_2$ assimilation. Coarse LAI led to a skewed pattern of ${\mathrm{CO}}_2$ exchange and crop growth in a crop growth model [10]. A regional eco-hydrology model, which uses LAI as a key canopy parameter, resulted in underestimation of evapotranspiration due to the use of inaccurately estimated LAI [11].

A remote sensing-based vegetation index (VI) has been used to monitor phenological and seasonal changes in vegetation development [5,12]. The normalized difference vegetation index (NDVI) is most widely used in estimating LAI over different spatiotemporal coverage [5,13,14,15]. This depends on the amount of chlorophyll and other pigments exposed to the view of the satellite [12]. The enhanced vegetation index (EVI) is often used for the crop mapping and crop phenology detection [16,17]. EVI has improved sensitivity to vegetation canopy structure and exhibits high correlation with the vegetation cover in areas with dense vegetation such as forested area [16,18].

Studies indicate that LAI and VI have either a linear or exponential correlation in forests, shrublands, and grasslands [19,20,21,22,23]. Higher correlation of the general regression between LAI and VI is observed when the vegetation is not extremely dense. Fan et al. (2008) [24] reported good general linear and exponential relationships (${\mathrm{R}}^2$ > 0.79) between LAI and NDVI over grasslands. Xiao et al. (2002) [25] reported a linear relationship (${\mathrm{R}}^2$ = 0.78) at the rice paddy. During the leaf production period, LAI and NDVI had a linear relationship (${\mathrm{R}}^2$ = 0.80) over a deciduous forest [26].

This correlation, however, is constrained by nonlinearity between LAI and VI as LAI increases, and this nonlinearity considerably varies among vegetation types [5,15,27]. Wang, et al. [15] showed an unconvinced linear relationship between LAI and NDVI during the middle of the growing season in a deciduous forest due to the saturation of NDVI at high LAI. In an alpine meadow, Lu, et al. [3] reported a difficulty in estimating LAI using an exponential relationship. Saturation of NDVI at high LAI results in a weak relationship of crops between LAI and NDVI [28]. Gupa, et al. [29] attempted to establish a general relationship between LAI and NDVI for wheat and onion crops using polynomial regression, but a complex crop growth cycle produced difficulty in generalizing the relationship. Brando, et al. [30] demonstrated an insensitive relationship between LAI and EVI, but showed the meaningful relationship with leaf phenology. Due to the strong dependence of VI on biophysical factors (e.g., vegetation canopy structure and growth states; [31]) that change at different growth stages, the general regression approach, using the entire growing season, produces less satisfactory estimates in seasonal dynamics of LAI [12,25,32].

Cultivated cropland covers 12% (1.5 billion hectare) of the Earth’s land area [33], and is important in its role as a strong modifier of regional carbon balance [34,35,36]. Growth stages of crops are primarily divided into vegetative/growth and reproductive phases. Each phase has different reflectance features, leading to different relationships between LAI and VI in each phase. Capturing the influence of seasonal phenological development on the relationship between LAI and VI is critical to improve applicability of the seasonal LAI of crops. However, few studies have considered the effect of different growth phases on this relationship [37]. In addition, most of the relationship has been derived from single or more localized sites. Estimating the seasonal LAI of crops across different climate conditions remains a challenge due to different timings of growth stages under varying management practices and local crop-planting schedules [38].

Here, we examine the relationships between LAI and NDVI and between LAI and EVI of rice, which is one of the major crops global and is cultivated in a wide range of locations with diverse climatic conditions [39]. Our objectives were (1) to estimate the seasonal course of LAI based on the traditional approach by partitioning the entire growing season into two growth phases; and (2) to develop an alternative method using a consistent development curve to estimate seasonal LAI by identifying the time at which the maximum NDVI (${\mathrm{NDVI}}_{\max }$) and EVI (${\mathrm{EVI}}_{\max }$) occur. To examine the applicability of this method across different climatic conditions and management cycles, we considered rice paddy sites in Asia (South Korea and Japan) and Europe (Spain).

2. Materials and Methods

2.1. Study Sites

The analysis was conducted using available rice paddy data from South Korea (Haean Catchment) and Japan (Mase) in Asia, and Spain (El Saler-Sueca) in Europe. Data used in the analysis cover the year of 2010 for Haean (HK), 2002–2005 for Mase (MSE), and 2007–2008 for El Saler-Sueca (ESES2). In order to further validate the alternative method, we applied it to an independent dataset (2006 for MSE, 2000 for Shizukuishi, and 2003 for Aso in Japan), which were not used to draw the consistent development curve. All of the Asian sites were under the influence of a monsoon climate, with >50% of annual precipitation occurring during the summer monsoon period. Meteorological conditions and site characteristics are summarized in

Table 1. Site information and meteorological conditions at Haean (HK, S. Korea), Mase (MSE, Japan), Shizukuishi (Japan), Aso (Japan), and El Saler-Sueca (ESES2, Spain). ${\mathrm{R}}_{\mathrm{g}}$ is the total global radiation, ${\mathrm{T}}_{\mathrm{a}}$ is the mean air temperature, and P is the sum of rainfall over the crop growth period (from transplanting date to harvest date). ${\mathrm{LAI}}_{\max }$ indicates maximum LAI.

Site	Year	${\mathbf{R}}_{\mathbf{g}}$	${\mathbf{T}}_{\mathbf{a}}$	P	Transplanting Date	Harvest Date	${\mathbf{LAI}}_{\mathbf{max}}$
		(MJ/${\mathbf{m}}^2$)	($°\mathbf{C}$)	(mm)	(DOY)	(DOY)
Haean	2010	1713	20.4	1165	144	290	5.8
Mase	2002	2356	22.0	593	122	262	5.5
	2003	2049	20.3	545	122	262	5.1
	2004	2384	22.7	547	123	254	4.9
	2005	2237	21.8	647	122	256	4.4
	2006	1989	21.5	632	122	141	6.0
El Saler	2007	3224	22.8	437	134	270	5.7
	2008	3263	22.1	121	132	278	6.1
Shizukuishi	2000	2057	21.4	615	143	263	4.5
Aso	2003	2176	21.4	1491	138	268	3.5

.

2.1.1. Haean, South Korea (HK)

Haean Catchment is a typical erosion mountain basin in South Korea located northeast of Chuncheon, Gwangwon Province in Yanggu County (38°17′N, 128°8′E, 450–1200 m above sea level (a.s.l.)). The total area of the catchment is 64 ${\mathrm{km}}^2$, consisting of 58% forested mountain area, 30% agricultural area, and 12% as residential, riparian, field margins, and farm road area according to land surveys [40]. The agricultural area is characterized as a mosaic patchwork of fields, with a dominance of dry-land fields (22% of the total area) and rice paddy fields (8%) as the remaining. Rice paddies (Oryza sativa L., cv. Odae) are cultivated at less than 500 m a.s.l. in the catchment [41]. LAI was measured by plot harvests and a leaf area meter (LI-3000A, LI-COR Inc., Lincoln, NE, USA) throughout the measurement period (from April to October 2010) at the rice paddy field.

2.1.2. Mase, Japan (MSE)

The Mase site is located in the rural area (36°3′14″N, 140°1′38″E, 15 m a.s.l.) of Tsukuba City in Central Japan. The rice paddy (Oryza sativa L.) was ca. 2 ${\mathrm{km}}^2$ and was managed as a single rice cropping field following practices common in the area [42]. LAI was measured from 2002 to 2005 with an optical area meter (AAM-7, Hayashi Denkoh, Tokyo, Japan). Further information of LAI measurements are described in Saito et al. (2005) [42].

2.1.3. Shizukuish, Japan

The Shizukuish site is located near Morioka, Iwate Prefecture in northern Japan (39°38′N, 140°57′E) where grows a large portion of Japan’s rice crop [43]. Seedlings of rice (Oryza sativa L.) were hand-transplanted into the field in late May, and the field was flooded throughout the growing season except several days of summer drainage and prior to harvest. LAI was measured through destructive samplings of biomass and further information of LAI measurement can be found in Kim et al. (2003) [43].

2.1.4. Aso, Japan

The Aso site is a commercial rice paddy in Aso Basin in Kyushu, Japan (32°56.9′N, 131°3.3′E). The rice paddy (Oryza sativa L., cv. Koshihikari) was cultivated following the common management practices. LAI was measured at 15–20 day intervals during the growing season using a leaf area meter (LI-3050A, LI-COR Inc., Lincoln, NE, USA). Further details about the site and measurement are described in Maruyama et al. (2007, 2008, and 2010) [44,45,46].

2.1.5. El Saler-Sueca, Spain (ESES2)

The El Saler Sueca site is located in the protected wetland area of La Albufera Natural Park in the Valencia region of Spain (39°16′32″N, 0°18′55″E, 10 m a.s.l.). El Saler, which is in a sub-arid Mediterranean climate, experiences hot summers with almost no rain and cold winters with substantial rainfall. The rice paddy was ca. 15 ${\mathrm{km}}^2$, and had been managed in same way for 200 years [47,48]. LAI was measured by plant samples from plot harvests. Further details about the site, agricultural management, and measurements are described in [47,48].

2.2. Vegetation Index from Remote Sensing

NDVI and EVI from 2002–2013 were obtained from daily gridded L3G (level-3) composite data at 250 m resolution, which is embedded in the MODIS Terra surface reflectance products (MOD13Q1) from the Warehouse Inventory Search Tool [49]. MOD13Q1 provides 16-day composite vegetation condition derived from blue (459–479 nm), red (620–670 nm), and near-infra-red (NIR, 841–876 nm) surface reflectance in the sinusoidal projection. NDVI and EVI are calculated as follows:

$$\mathrm{NDVI}=\left(({\rho }_{NIR}-{\rho }_{red})\right)/\left(({\rho }_{NIR}+{\rho }_{red})\right)$$

(1)

$$\mathrm{EVI}=G({\rho }_{NIR}-{\rho }_{red})/({\rho }_{NIR}+C1/rh{o}_{red}-C2{\rho }_{blue}+L)$$

(2)

where $G=2.5,C1=6,C2=7.5$, and $L=1$ (Huete et al. [18]). Raw values of VI include noise components due to clouds, water, snow, shadow, bidirectional effects, high solar or scan angles and transmission errors. Since about 80% of the rice was grown in the field that was flooded for long time period [50], background effect on VI was minimal in this study. The TIMESAT program, developed by Jonsson and Eklundh [51], was applied to smooth VI and minimize the influence of the associated errors. VI was smoothed by the adaptive Savitzky-Golay filtering method in the TIMESAT program to produce final VI (Figure 1). The adaptive Savitzky-Golay filtering is able to follow complex fluctuations that occur with rapid increase and decreases in the reflectance data [51,52] and uses local least-squares polynomial approximation. This character of the filtering fit best for agricultural crops changing their phenological and physiological conditions over short time intervals.

2.3. Leaf Area Index Estimates

Estimates of LAI are tightly correlated with spectral reflectance, which changes with different growth stages [24,25,53,54]. For instance, regulation of protein and pigment synthesis is most likely altered as the growth of rice grains compete for plant resources [55,56], leading to changes in the reflectance of rice. To account for these differences in the reflectance along the development stages${\mathrm{LAI}}_{\max }$, we divided the entire growing season into two growth phases: the vegetative phase with increasing LAI (the beginning to the maximum leaf expansion; indicated as Before in figures) and the reproductive phase with decreasing LAI (the maximum leaf expansion to the harvest; indicated as After in figures).

In this study, the relationship between LAI and NDVI was derived from the general regression method, i.e., linear and exponential regression method and the consistent development curve method. The regression method was applied to the data from the entire growing season and the two growth phases. According to the results from the two regression approaches, the relationship extracted from an exponential regression illustrated higher correlation than from a linear regression. Due to better explicability of time-dependent change in LAI of the exponential model than of the linear model (e.g., [3,24,53,57]), we only presented the results from the exponential model.

The consistent development curve method is based on the general additive models (GAMs; [58]) and explains the variation of LAI in relation to ${\mathrm{EVI}}_{\max }$ and ${\mathrm{EVI}}_{\max }$ respectively. The observed LAI was normalized as a ratio of the observed LAI at each site to the average observed maximum LAI (${\mathrm{LAI}}_{\max }$, 5.3 from all sites). Each site had one value of ${\mathrm{VI}}_{\max }$ per a growing season except MSE in 2003 with the double peaks (DOY 183 and 215) in VIs. We choose the value of VI on DOY 215 as ${\mathrm{VI}}_{\max }$, which occurred in vicinity of the observed ${\mathrm{LAI}}_{\max }$. At each site, the number of departed days from the day of each ${\mathrm{VI}}_{\max }$ was counted at each point of the normalized LAI. Then, the GAMs were applied to this relationship in order to draw a consistent development curve of rice. Negative values indicate the observation occurred before each ${\mathrm{VI}}_{\max }$, whereas positive values indicate the observation occurred after each ${\mathrm{VI}}_{\max }$.

GAMs are non-parametric models that examine the data before a particular response function (e.g., linear and quadratic functions) is selected, and finds the response function through a smoothing process. They has been used to understand non-linear ecological responses to a wide range of environmental variables [58,59,60]. GAMs can be more general and flexible to allow a wide range of response curves to the non-linear relationships and provide an alternative approach for limitation of linear functions that fails to detect complex dependency between the predictor and response variables. When X is an independent variable, a linear function of GLMs (${\displaystyle {\sum }_{j=1}^p}{\beta }_j{X}_j$ in Equation (3)) is replaced to a sum of smooth function (${\displaystyle {\sum }_{j=1}^p}{f}_j\left(X_j\right)$ in Equation (4); [59,60,61]).

$$g\left(\mu \right)={\beta }_0+{\displaystyle \sum _{j=1}^p}{\beta }_j{X}_j$$

(3)

$$g\left(\mu \right)=f_0+{\displaystyle \sum _{j=1}^p}{f}_j\left(X_j\right)$$

(4)

where g$\left(\mu \right)$ is the link function, ${\beta }_0$ and ${\beta }_j$ are the intercept and slope of a linear regression, and $f_0$ is the intercept of a non-parametric function. After aligning the day of each ${\mathrm{VI}}_{\max }$ from the consistent development curve with the day of each ${\mathrm{VI}}_{\max }$ at specific sites, the normalized LAI was converted into normal LAI by multiplying daily values of the curve by the average observed ${\mathrm{LAI}}_{\max }$.

3. Results

3.1. Exponential Model for Estimating LAI

LAI was positively correlated with NDVI and EVI when the regressions were applied to the entire growing season (Figure 2). The regression analysis showed good correlation between LAI and NDVI (0.63 < ${\mathrm{R}}^2$ < 0.83 and 0.82 < RMSE < 1.56) and between LAI and EVI (0.62 < ${\mathrm{R}}^2$ < 0.72 and 1.06 < RMSE < 4.41;

Table 2. Summary of regression analyses between LAI and NDVI and between LAI and EVI. The analysis was conducted using the data from different growth phases: “Entire” indicates the data from the entire growing season. “Before” indicates the data collected before the maximum VI, whereas “After” indicates the data collected after the maximum VI. ${\mathrm{R}}^2$ is the adjusted coefficient of determination, and RMSE is the root mean square error. HK = Haean (S. Korea), MSE = Mase (Japan), and ESES2 = El Saler-Sueca (Spain).

Site	Year	Growth Phase	NDVI					EVI
			Regression Equation	${\mathbf{R}}^2$	RMSE	p		Regression Equation	${\mathbf{R}}^2$	RMSE	p
HK	2010	Entire	y = ${\mathrm{e}}^{(16.53x-12.24)}$	0.83	0.82	0.003		y = ${\mathrm{e}}^{(18.28x-9.93)}$	0.72	1.06	0.01
MSE	2002–2005	Entire	y = ${\mathrm{e}}^{(9.54x-6.41)}$	0.63	1.56	<0.001		y = ${\mathrm{e}}^{(8.01x-3.49)}$	0.62	4.41	<0.005
ESES2	2007–2008	Entire	y = ${\mathrm{e}}^{(6.39x-3.87)}$	0.67	1.27	<0.001		y = ${\mathrm{e}}^{(6.72x-2.7)}$	0.72	1.98	<0.005
HK	2010	Before	y = ${\mathrm{e}}^{(16.52x-12.51)}$	0.90	1.21	0.009		y = ${\mathrm{e}}^{(17.81x-10.31)}$	0.96	0.32	0.01
		After	-	-	-	-		y = ${\mathrm{e}}^{(10.72x-4.93)}$	0.97	0.3	0.08
MSE	2002–2005	Before	y = ${\mathrm{e}}^{(11.28x-8.04)}$	0.82	0.94	<0.005		y = ${\mathrm{e}}^{(9.73x-4.37)}$	0.72	6.63	<0.005
		After	y = ${\mathrm{e}}^{(1.98x-0.2)}$	0.42	0.78	0.007		y = ${\mathrm{e}}^{(0.74x+0.83)}$	0.01	1.13	0.4
ESES2	2007–2008	Before	y = ${\mathrm{e}}^{(7.30x-4.82)}$	0.84	0.84	<0.005		y = ${\mathrm{e}}^{(7.68x-3.42)}$	0.88	1.22	<0.005
		After	y = ${\mathrm{e}}^{(1.81x+0.08)}$	0.93	0.24	<0.005		y = ${\mathrm{e}}^{(1.70x+0.5)}$	0.73	0.53	0.009
Asia	2002–2010	Before	Y = ${\mathrm{e}}^{(11.57x-8.38)}$	0.80	0.97	<0.005		y = ${\mathrm{e}}^{(9.13x-4.4)}$	0.61	3.43	<0.005
		After	y = ${\mathrm{e}}^{(2.12x-0.38)}$	0.27	0.93	0.02		y = ${\mathrm{e}}^{(1.23x+0.51)}$	0.05	1.3	0.2
Asia & Europe	2002–2010	Before	y = ${\mathrm{e}}^{(10.11x-7.21)}$	0.76	1.03	<0.005		y = ${\mathrm{e}}^{(8.81x-4.18)}$	0.66	2.88	<0.005
		After	y = ${\mathrm{e}}^{(2.11x-0.31)}$	0.36	0.85	0.001		y = ${\mathrm{e}}^{(1.41x+0.49)}$	0.14	1.2	0.05

). At MSE and ESES2, where multiple years of data were included, greater scatter was found in the relationship with both cases than when individual years were analyzed (cf. HK for 2010). The greater scatter at MSE and ESES2 may result from the seasonality of phenological development of rice varying from year to year.

The regressions were applied to the data from the two phases separately. The increasing phase was associated with the increase in LAI (from 0 to 6), NDVI (from 0.3 to 0.9), and EVI (from 0.1 to 0.8), whereas the decreasing phase was associated with the decrease in LAI (from 6 to 2), NDVI (from 0.9 to 0.6), and EVI (from 0.8 to 0.2). The regressions at each phase were quite different from each other at all sites. The relationship exhibited higher ${\mathrm{R}}^2$ (0.84–0.96) in the increasing phase (Before in Figure 3 and Table 2) than the decreasing phase (${\mathrm{R}}^2$ = −0.01–0.97; After in Figure 3 and Table 2) for NDVI and EVI. An increasing rate of LAI along the change in NDVI and EVI was greater in the increasing phase than in the decreasing phase. Due to a large scatter in the data, the relationship between LAI and EVI was weaker than that between LAI and NDVI particularly at MSE. Compared to the correlations between LAI and VIs from the entire growing season, they were mostly higher from the increasing and decreasing phases, demonstrating the importance of appropriate consideration of phenological development at different phases in estimating seasonal LAI.

In order to establish relationships between LAI and VIs that can be generally applicable and “universal” rather than site-specific, we pooled the data for Asian sites (HK and MSE) and all sites (HK, MSE, and ESES2), respectively. During the summer monsoon in the Asian region, the remotely-sensed vegetation indices were influenced by frequent cloud cover and high water vapor density in the air. For this reason, Asian rice paddy sites were treated separately. Using the data obtained over the entire growing season, the LAI-VI relationship for Asian sites and all sites was examined in Figure 4. The LAI-NDVI relationship was similar with the LAI-EVI relationship (${\mathrm{R}}^2$ of 0.60) for the Asian sites and all sites. As for the two phases, the LAI-VI relationship showed much better agreement in the increasing phase (Before, 0.61 < ${\mathrm{R}}^2$ < 0.80 and 0.97 < RMSE < 3.43) than in the decreasing phase (After, 0.05 < ${\mathrm{R}}^2$ < 0.36 and 0.85 < RMSE < 1.3; Figure 5 and Table 2). As indicated in the analysis for individual sites (Figure 3), separation of the seasonal phases provided a better explanation of LAI changes with VIs in the pooled data. The correlation between LAI and NDVI was better than that between LAI and EVI for the two growth phases.

Figure 6 presents seasonal dynamics in the LAI estimated by exponential regressions produced from the entire growing season (Figure 4) and the two phases (Figure 5) for the pooled sites (Asia and Europe). The estimated LAI from both cases ranged from 0 to 8.6 with NDVI and from 0 to 24 with EVI, showing either overestimation or underestimation at most sites. LAI estimated from both NDVI and EVI was less satisfactory. Separating the season into two phases led to an abrupt change in seasonal LAI, where the transition from Before to After occurred in the mid-season (near Day of Year (DOY) 200), resulting in no better estimation than the estimation obtained with the entire growing season. LAI estimates from both cases showed similar variation and magnitude with the measure LAI in the increasing phase, but most of the discrepancy in the estimates occurred in the decreasing phase. The timing of the peak in LAI was acceptable in comparison to the peak of the measured LAI (within ±3 to 6 days) except MSE in 2003 showing double peaks (DOY 183 and 215) associated with the summer monsoon. However, the largest difference between the estimated and the measured LAI and between the two cases occurred during the peak.

3.2. Estimation of LAI from the Consistent Development Curve Method

A consistent development curve for rice across the study sites was presented in Figure 7. The scaled LAI slowly increased for about 80 days before the ${\mathrm{VI}}_{\max }$ (Day 0). The maximum in the scaled LAI occurred between Day 0 and 10 days after Day 0. A relatively slow decrease occurred after this peak. The estimated LAI from the consistent development curve method ranged from 0 to 4.7 (Figure 8). Magnitudes and seasonal patterns in LAI from NDVI and EVI were well reproduced across all sites except a few points of noticeable deviations from the measured LAI (e.g., HK in 2010). Compared to the results from the exponential regressions of the entire growing season and the two phases (Figure 6), the consistent development curve method had better estimates of LAI in the peak and the decreasing phase at all sites (

Table 3. Statistics of the correlation between the measured and the estimated LAI. Estimates of LAI were calculated from consistent development curve method (Figure 7) and exponential regression method (Figure 4). ${\mathrm{R}}^2$ is the adjusted coefficient of determination, RMSE is the root mean square error, and CV is coefficient of variation of RMSE. HK = Haean (S. Korea), MSE = Mase (Japan), and ESES2 = El Saler-Sueca (Spain).

VI	Site	Year	Consistent Development Curve					Exponential Regression
			${\mathbf{R}}^2$	RMSE	CV (%)	p		${\mathbf{R}}^2$	RMSE	CV (%)	p
NDVI	HK	2010	0.78	0.89	59.08	<0.005		0.75	0.93	57.12	<0.005
	MSE	2002	0.90	0.38	23.18	<0.005		0.65	1.19	59.25	<0.005
		2003	0.90	0.33	22.41	<0.005		0.24	0.84	56.44	<0.005
		2004	0.92	0.21	11.49	<0.005		0.75	1.30	56.37	<0.005
		2005	0.95	0.77	41.92	<0.005		0.89	0.37	28.00	<0.005
		2006	0.92	0.57	25.11	<0.005
	ESES2	2007	0.92	0.77	27.46	<0.005		0.68	0.65	54.22	<0.005
		2008	0.93	0.72	23.77	<0.005		0.60	0.37	57.85	<0.005
	Shizukuishi	2000	0.98	0.54	19.42	<0.005
	Aso	2003	0.98	0.77	41.79	<0.005
EVI	HK	2010	0.73	1.19	73.38	0.008		0.56	1.11	69.45	0.03
	MSE	2002	0.93	0.39	21.00	<0.005		0.86	1.57	69.53	<0.005
		2003	0.96	0.20	14.48	<0.005		0.74	1.54	57.30	<0.005
		2004	0.91	0.23	11.90	<0.005		0.85	1.22	43.22	<0.005
		2005	0.88	0.64	42.93	<0.005		0.68	5.76	342.96	<0.005
		2006	0.86	0.69	31.70	<0.005
	ESES2	2007	0.89	0.82	29.25	<0.005		0.64	1.63	90.44	0.01
		2008	0.88	0.82	27.43	<0.005		0.69	1.04	40.86	<0.005
	Shizukuishi	2000	0.97	0.57	20.41	0.01
	Aso	2003	0.98	0.72	39.73	<0.005

). Although the consistent development curve method had a limited scale due to the use of the average observed ${\mathrm{LAI}}_{\max }$ (5.3) being set as ${\mathrm{LAI}}_{\max }$, it improved applicability of the seasonal course of LAI much better than the exponential model did.

In order to further validate the consistent development curve method, we applied the method to an independent dataset (MSE, 2006, Shizukuishi, 2000, and Aso, 2003 in Japan; [43,45]), which were not used to draw the consistent development curve in Figure 7. The measured LAI and VIs were collected and processed with the same methods described in the method section. The estimated LAI from VIs showed a good comparison with the measured LAI at both sites (${\mathrm{R}}^2$ = 0.92 for NDVI and ${\mathrm{R}}^2$ = 0.86 for EVI at MSE, ${\mathrm{R}}^2$ = 0.98 for NDVI and ${\mathrm{R}}^2$ = 0.97 for EVI at Shizukuishi, and ${\mathrm{R}}^2$ = 0.98 for NDVI and ${\mathrm{R}}^2$ = 0.98 for EVI at Aso; Figure 9), capturing well the increasing and decreasing phases of and the magnitudes of LAI. The estimated LAI using EVI showed a slight shift of the peak compared to that of the measured LAI at MSE in 2006.

4. Discussion

4.1. Estimation of LAI using Different Growth Phases

We hypothesized that the phenological and physiological features of rice are different at the beginning of the growing period compared to those at the senescence period. When the exponential regression method was applied to estimate LAI, LAI estimates from the two growth phases generally fitted better to the LAI measurements than those from the entire growing season did. LAI estimates from the latter were frequently underestimated during the decreasing phase, confirming the necessity of considering different developments of seasonal phenology. Large deviations of LAI estimates from LAI measurements were observed in the middle of the growing season (e.g., MSE, 2004 and ESES2, 2007 for NDVI and EVI, Figure 6) when saturation of VIs occurred at high LAI. Due to hysteresis effects centering on ${\mathrm{LAI}}_{\max }$ in the relationship between LAI and VIs, the abrupt change occurred in transition between the two phases in the middle of the growing season (Figure 6). No better seasonal estimates in LAI from the two phases was produced than those from the entire growing season. This hysteresis effect may result from changes in leaf pigment and chlorophyll content at different vegetative growth stages [62,63,64]. Therefore, attention must be paid to seasonal estimates of LAI using the relationship between LAI and VIs particularly in the second half of the growing season when VIs were saturated.

In relation to the LAI estimates of rice, neither NDVI nor EVI was preferable with the general regression approach. Over- and under-estimation of LAI were frequently observed from the LAI-NDVI and the LAI-EVI relationships. When NDVI was greater than 0.85 at specific sites and years (e.g., MSE, 2005 and ESES2, 2007–2008) in the middle of the growing season, LAI tended to be overestimated. When EVI was greater than 0.8 at MSE (2005), an extreme overestimation of LAI was observed. On the contrary, LAI was underestimated when EVI was low (<0.6) at MSE (2002). A change rate of LAI along the change in NDVI and EVI was greater in the exponential regression as VIs were saturated. Estimation of LAI became very sensitive to a small change in VIs, resulting in the difference ranging from −3 to 20 between the estimated and the measured ${\mathrm{LAI}}_{\max }$. These results led to a caution in using the exponential regression approach due to saturation effect on the relationship between LAI and VIs.

Double peak in the estimated LAI at MSE (2003) reflects the limitation of TIMSAT in recovering the seasonal variation of VIs (Figure 1 and Figure 6). VIs were more scattered in mid-summer due to the influence of summer monsoon than other sites. This scatter affected subsequent analyses to define a smoothed curve and estimate daily VIs. VIs, which were smoothed by TIMESAT, could not capture the specific signal of LAI [65]. We found that statistical methods to obtain seasonal VIs and LAI did not necessarily provide correct data for the rice.

4.2. Performance of the Consistent Development Curve Method

We attempted to identify a better approach based on consistencies in the developmental processes of the rice crops, and provided a good reference case in that the seasonal change in the LAI of rice paddies can be determined with both NDVI and EVI using the consistent development curve method. This method showed significant improvement in estimating the seasonal course for LAI, resulting in a better estimate in LAI than an exponential model did (Figure 8). Reliance on consistent biological regulation of phenological development through the detection of the maximum in LAI and VIs allows the identification of seasonal dynamics in the LAI of the rice paddies. An exponential tendency in the relationship between LAI and VIs is extracted mostly from single or more localized sites [23,24,27]. In cases where study locations are geographically spread across different climate regions and where year-to-year climate variation occurs, such simple relationships may not be sufficient to estimate LAI universally. This study, however, has proved that the consistent development curve method can identify seasonal course of LAI spreading different geographical locations and climates.

Less satisfactory agreement between the measured ${\mathrm{LAI}}_{\max }$ and the estimated ${\mathrm{LAI}}_{\max }$ was generated from a failure in detecting a local ${\mathrm{LAI}}_{\max }$. Differences between a local ${\mathrm{LAI}}_{\max }$ at a specific site and the average observed ${\mathrm{LAI}}_{\max }$ (5.3) was carried into the conversion of the scaled LAI from curve to normal LAI by multiplying the daily values of the curve by the average observed ${\mathrm{LAI}}_{\max }$. For instance, a local ${\mathrm{LAI}}_{\max }$ (6.0) at MSE in 2006 was higher than the average observed ${\mathrm{LAI}}_{\max }$, resulting in lower estimate of ${\mathrm{LAI}}_{\max }$. On the contrary, that the local ${\mathrm{LAI}}_{\max }$ (3.5) at Aso in 2003 was lower, leading to an overestimation of ${\mathrm{VI}}_{\max }$. The consistent development curve method is strong in detecting the seasonal phenological development of the relationship between LAI and VIs. However, the need still remains to improve accuracy in estimating the maximum value.

NDVI and EVI appropriately estimated seasonal pattern and magnitude of LAI (Figure 8 and Figure 9, Table 3). The LAI-NDVI relationship was slightly better in finding the timing of ${\mathrm{LAI}}_{\max }$ than the LAI-EVI relationship did. For instance, the timing of ${\mathrm{LAI}}_{\max }$ from NDVI matched well with that of the measured ${\mathrm{LAI}}_{\max }$, but ${\mathrm{LAI}}_{\max }$ from EVI occurred in several days later (e.g., 12 days in 2005 and 7 days in 2006) due to the delayed peak of EVI (Figure 1). Because the number of departed days from the peak day of EVI was accounted in the estimated LAI, this results in a slight shift of the entire seasonal pattern of the estimated LAI toward a later season. We speculate that EVI was not strongly convinced to capture the maturation of rice in these years, hence relatively weaker explanation was observed in estimating LAI with EVI than with NDVI.

5. Conclusions

Improving seasonal estimation in LAI can provide valuable information to aid accurate estimating of biophysical processes by process-based models and interpretation of vegetation dynamics [28]. Estimation of LAI in crops is still a challenge due to crops’ distinctive phenological changes before and after ${\mathrm{VI}}_{\max }$ [66]. Varying management practices at individual sites hinder the generally applicable and the universal relationship between LAI and VIs to be adopted to specific species. Consistent development curves with VIs in this study show improvement in detecting phenological development and estimating the seasonal LAI of rice. There is still a shortcoming in the use of VIs to estimate LAI due to the limitation of the smoothing process (e.g., TIMESAT) in recovering the seasonal variation of VIs.

Despite a few limitations, we utilized the VIs as a detector of the seasonal phenological development. It also illustrates its applicability across different climatic conditions and management cycles. We expect that the consistent development curve method with VIs can improve the seasonal estimation of LAI for the various crop types across space. With the consistent development curve approach, which performs more accurately than traditional regression approaches, estimates of critical parameters such as crop growth and productivity in the process-based models can be refined and applied to other crops to determine their phenological stages. Further testing will demonstrate whether this “biophysically-based” approach, as in the case of rice, can effectively estimate the critical LAI values needed in regional process-based models.