Using MODIS LAI Data to Monitor Spatio-Temporal Changes of Winter Wheat Phenology in Response to Climate Warming

: Understanding spatio-temporal changes in winter wheat ( Triticum aestivum L) phenology and its response to temperature will be vital for adapting to climate change in the coming years. For this purpose, the heading date (HD), maturity date (MD), and length of the reproductive growth period (LRGP) were detected from the remotely sensed leaf area index (LAI) data by a threshold-based method during the harvest year 2003 to 2018 across the North China Plain. The results show that there was high spatial heterogeneity of winter wheat phenology in pixel scale across the whole area, which could not be detected in previous site-based studies. The results also verified that climate warming could explain part of the change in the HD. However, for the LRGP, the potential impact of non-climate effects should be further investigated. This study presents the spatiotemporal changes both in winter wheat phenology and corresponding mean temperature and then analyzes their relationships in pixel scale. Additionally, this study further discusses the potential impact of non-climate effects on the LRGP.


Introduction
Over the past several decades, the annual mean of daily mean temperatures in the North China Plain (NCP) has experienced an increasing trend of approximately 0.2 ℃ per decade [1,2]. Climate change has influenced the terrestrial ecosystems [1][2][3][4][5]. In the future, a warmer growth environment could further affect the vegetation and the whole ecosystem. Phenological shifts can be an important indicator for us to better understand some of the changes to ecosystems [6][7][8][9]. For native and perennial plant species, increasing spring temperatures could generally boost vegetation growth, lengthen the growing season, and enhance the net ecosystem productivity [10][11][12]. However, for the winter wheat in the NCP, a warmer climate could accelerate senescence, shorten the growing period, and therefore lead to yield loss [13][14][15][16][17].
Phenological shifts can respond to climate change effectively [18,19]. For instance, due to global warming, the advancement of green-up, heading, and maturity dates in winter wheat (Triticum

Study Area
The NCP is the largest agricultural region in China (approximately 400,000 km 2 ) and is the main wheat-producing region in East Asia. Our study area is a wheat-maize rotation cropping system with relatively adequate irrigation and fertilizer ( Figure 1) [46]. In this area, winter wheat is sowed in October and harvested in June of the following year. Although most studies have concluded that climate warming has reduced global wheat production, the actual yields in the NCP have in fact increased over the past 50 decades [47,48]. The adaptation options and techniques could offset the negative impact of climate warming [31,49]. Due to the interaction between climate warming and non-climate effects, the phenological change in this region will be more uncertain in the future.

Phenological Metrics from Remotely Sensed LAI Data
The LAI product (MCD15A3H) was a four-day composite dataset with 500 m pixel size, including the best pixel available from all the acquisitions of the Moderate Resolution Imaging Spectroradiometer from within the four-day period (https://lpdaac.usgs.gov/products/mcd15a3hv006/). The LAI data during winter wheat growing seasons from 2003 to 2018 were obtained from the website of the Earth Observing System Data and Information System (EOSDIS) (https://search.earthdata.nasa.gov/).
In this study, the time series of remotely sensed LAI data were used to estimate winter wheat phenological metrics. Considering the characteristics of winter wheat canopy structure, LAI time series could be consistent with field phenological observations [50]. Due to cloud and haze effects, a filtering algorithm was applied to remove noise [51]. In this study, the Harmonic Analysis of Time Series (HANTS) filtering algorithm was used to fit the curve on a pixel-by-pixel basis using the ENVI image analysis software (version 5.3) and Interactive Data Language (IDL, version 8.5). Moreover, it is certain that the HANTS filtering algorithm can be used to detect vegetation phenology [52,53].
This study used the threshold-based method to estimate phenological metrics. The thresholdbased method has been proven to be an acceptable method for estimating phenological information [24,37,54]. The phenological node was defined as the fitted curve reaching a fixed threshold. In this study, the heading date (HD) was defined as the corresponding date of the maximum LAI value for each year, and the maturity date (MD) of winter wheat was defined as the date when the fitted curve of the LAI time series reached 20% of its maximum amplitude for each year [24,37]. The length of the reproductive growth period was calculated as the number of days between the HD and MD for each year.

Winter Wheat Planting Area
Most of the previous studies analyzed the change in wheat phenology with the assumption that the wheat planting area was unchanged during the research period [24,36]. However, this assumption could lead to a high uncertainty of estimated wheat phenology in each pixel. In this study, the winter wheat planting area for each year was extracted based on a simple decision tree method. Decision tree has been widely used in land cover classification and can improve our results [55,56]. A binary decision tree was built to distinguish the pure pixels of winter wheat from other pixels (the non-winter-wheat and mixed pixels). The estimated heading and maturity dates (HD and MD) have been specified as two phenological nodes ( Figure 2). Based on observational ground phenology data and previous studies, the ranges of the estimated HD and MD were set between day of year (DOY) 93 and 145, and DOY 125 and 181, respectively [23,26,46]. The pixel value should be within both of these ranges, otherwise it would be set to nodata.

Climate and Observational Ground Phenology Data
Monthly gridded air temperature (T) was obtained from the Global Land Data Assimilation System (GLDAS), with 0.25 arc degrees spatial resolution (https://ldas.gsfc.nasa.gov/) [57]. The heading and maturity dates during the harvest years 2003 to 2018 were obtained from agrometeorological observation stations (http://data.cma.cn/en/). The performance of satellite-derived phenological metrics was evaluated by using the mean absolute error (MAE): where i y and ˆi y are the observed and estimated values, respectively, and m is the number of observations.

Analysis
A Theil-Sen slope estimator was used to calculate the trends of phenology and temperature for each pixel, respectively [58,59]. Next, the nonparametric Mann-Kendall (M-K) test was used to quantify the statistical significance of the results in pixel scale [60,61]. According to the Z-values from the M-K test, if |Z-value| > 1.65, 1.96, and 2.32, the trend is at the 99%, 95%, and 90% significance level (p < 0.01, 0.05, and 0.1), respectively. Additionally, the Pearson correlation coefficient was calculated to determine the correlation between phenological metrics (the HD and LRGP) and the corresponding mean temperature in pixel scale.
It should be noted that this study used the monthly mean temperature instead of the daily mean temperature in the analysis. Although the daily mean temperature would be more accurate to obtain the real temperature of the corresponding growing season, we had to consider how to match with the satellite-based phenology in each pixel. Here, this study assumed that the mean temperature in March-April and May-June have primary effects on the HD and LRGP, respectively [24].

Spatio-temporal Changes in Satellite-based Winter Wheat Phenology
The estimated phenological metrics were validated and compared with ground-observed records ( Figure 3). The satellite-based heading date (HD, R = 0.364, p < 0.01) and maturity date (MD, R = 0.851, p < 0.01) both correlated significantly with the ground-observed dates. The mean absolute errors (MAE) of the HD and MD were 8.2 days and 4.5 days, respectively. Overall, threshold-based phenological metrics could be used to produce a robust spatio-temporal distribution of winter wheat phenology.    (Figure 5), respectively. In pixel scale, the advancement of the HD and MD and the shortening of the LRGP occurred significantly in 4.7%, 1.9%, and 4.3% of the pixels, respectively, which were mostly distributed in the north of the NCP (p < 0.1). On the contrary, the delay trend of the HD and MD and the prolongation of the LRGP were significantly in 8.1%, 14.4%, and 16.4% of the pixels, respectively, and mostly distributed in the south of the NCP (p < 0.1). The spatio-temporal changes in winter wheat phenology showed high spatial heterogeneity. Our findings suggest that the phenological shifts could be responding to some effect factors.

Spatio-temporal Changes in Corresponding Mean Temperature
The results showed the spatial distributions of the two-month mean temperature in March-April (two months before the HD) and May-June (two months before the MD), averaged from the harvest year 2003-2018 ( Figure 6). Considering the impact of temperature on winter wheat phenology, we further calculated the change in temperature that would further affect the HD and LRGP (Figure 7). Our findings show that the two-month mean temperature in March-April was increased by 0.1 ℃ per decade during 2003-2018, which could be significant in the north and east part of the NCP. Meanwhile, the two-month mean temperature in May-June increased by 0.6 ℃ per decade. Overall, the results indicated that the warming trend in spring would be more obvious than in summer across the NCP.

The Responses of Winter Wheat Phenology to Temperature
Pearson's correlation coefficient was calculated to further demonstrate whether the spatial patterns of HD/LRGP were correlated with temperature. The results show that the HD was negatively correlated with the March and April mean temperature in most parts of this study area, and 15.7% of the pixels were significant (p < 0.05) (Figure 8a,c); but no significant correlations were found between the LRGP and the May and June mean temperature (Figure 8b,d). Overall, the advanced trend of the HD could be attributed to climate warming in March and April. However, there was no significant correlation between the LRGP and mean temperature, which suggests that some non-climate effects could have an impact on the change in the LRGP.

Improved Monitoring on Wheat Phenology Using MODIS LAI Data
Using the satellite-based phenological metrics and gridded temperature data, this study assessed the spatio-temporal changes in winter wheat phenology and climate warming across the NCP during the harvest years 2003-2018. The results show better temporal and spatial resolution than previous studies [24,36,46,62]. Moreover, our study could detect more details of spatial heterogeneity in wheat phenology than some site-based studies [14,15,29]. Based on the remotely sensed LAI product, we can better conduct the monitoring of phenology. In our study, we derived the spatio-temporal changes in the HD, the MD, and the LRGP by using the MODIS LAI product. Our findings suggest that LAI data could be used to detect phenological metrics with high and stable spatio-temporal resolution [44,45].
In addition, unlike previous studies that rarely consider land use types and land use changes, this study extracted the winter wheat planting area each year by using the decision tree method, which can reduce the uncertainty of inter-annual changes in the crop planting area [24,36]. For the spatio-temporal changes in temperature, our study was consistent with several previous studies at the site-and region-scales [2,26,33]. Therefore, compared with previous studies, this study presented a viable approach for mapping the spatio-temporal changes in winter wheat phenology and detecting more detailed spatial heterogeneity in crop phenology at a regional scale.
Previous site-based studies showed that there are not consistent trends of winter wheat phenology across the NCP. Tao et al. [29] reported that the HD advanced about 1.28 days per year at 40 stations in China. However, our results showed that the mean HD delayed 0.4 days per year across the NCP. Similarly, the delay trend of HD was observed in North America [24]. Different scales of studies would lead to these diverse results. Additionally, the delayed trend of HD was mainly distributed in the south of the NCP, where late-flowering winter wheat cultivars were used to avoid late frost. Moreover, the difference among pixels increased with the diverse non-climate effects. For the MD, Xiao et al. [21] found that the MD was advanced with an average of 0.17 days per year at most stations of the NCP, which is lower than the 0.37 days per year in our study. For the LRGP, Tao et al. [29] reported that the LRGP of winter wheat was prolonged at 58 stations but shortened at 40 stations in China. Our study found that the LRGP was shortened with a mean changing rate of −0.32 days per year, which was consistent with Ren et al. (−0.36 days per year) [24], but contrary with Xiao et al. (0.13 days per year) [21]. The inconsistent studies suggest that different data would lead to different results in crop phenological trends. Further research is required in order to contrast the obtained results across different data.

The Potential Impact of Non-climate Effects on the LRGP
Using site-based observational records, previous studies have found that non-climate effects impact winter wheat phenology [26,[30][31][32][33]. To explain the non-significant correlation between winter wheat phenology and climate warming in this study, we further discuss the potential impact of nonclimate effects on the LRGP. We used an indirect method to assess the potential impact of non-climate effects on the LRGP according to the following equation: v obs C C T     (2) where v C (days/year) represents the change in the LRGP only under the impact of non-climate effects, obs C (days/year) represents the satellite-observed phenological trend of the LRGP,  (days/℃) represents the mean temperature sensitivity of the LRGP determined by crop models in previous studies, which could be stable and representative across the NCP [26,33], and T  (℃/year) represents the change trends of corresponding mean temperature. The result showed that the impact of climate warming on the LRGP had changed by -0.1 ± 0.27 days/year during 2003-2018 with large spatial differences across the NCP (Figure 9a). Furthermore, the impact of non-climate effects on the LRGP had a mean changing rate of -0.2 ± 1.46 days/year (Figure 9b). Although the results ( Figure 9) are not convincing enough to prove the impact of non-climate effects on the LRGP, they improve our understanding of the non-significant correlations between climate warming and winter wheat phenology (Figure 8b,d). The impact of non-climate effects could lead to more spatial differences, especially between parts of the north and the south. The results suggest that spatial heterogeneity could be induced by non-climate effects, which may be unable to be validated at present. Contrasting with the spatial interpolation used in site-based studies, this indirect method can be used to map the spatio-temporal changes of non-climate effects in continuous space [29,30]. We expect that our findings can be further evaluated and extended with different crops and in different regions.

Uncertainties and Limitations
Although the NCP is a typical irrigated agriculture region, with relatively enough water and fertilizer, no additional data can be used to validate this in pixel scale [49,63,64]. The complicated interactions between climate warming and non-climate effects were also not considered in this study. In addition, the rough monthly temperature data we used in this study, similarly to other studies, could lead to some uncertain results [17,24,37]. As a remote sensing study, another limitation that we have to acknowledge is that the results can monitor and map spatio-temporal changes in winter wheat phenology but may not give a full explanation for mechanisms. Although this study analyzed the relationships between winter wheat phenology and climate warming and discussed the potential impact of non-climate effects, our results cannot separate these impacts without using more parallel field experiments. Moreover, studies on the assimilation of remote sensing data into crop growth models could also improve our understanding of how crops respond to climate warming [65]. Consequently, the impact of non-climate effects, such as crop management and cultivar improvement, should be further analyzed in future research. In addition, more studies will be needed to reduce the uncertainties induced by different study scales and data.

Conclusions
Based on the four-day composite LAI product, the spatio-temporal changes in satellite-based winter wheat phenology were derived with more detail. The results show that there were spatial differences and heterogeneity across the NCP in the HD, MD, and LRGP, with a mean changing rate of 0.4 days per year, −0.37 days per year, and −0.32 days per year (p < 0.1), respectively. Meanwhile, the spatio-temporal changes in temperature showed that the warming trend could be obvious during the winter wheat growth period. The relationships between temperature and the HD showed that 15.7% of pixels in the NCP were significantly negative (p < 0.05). However, there was no significant correlation between the LRGP and temperature. The non-climate effects could weaken the impact of climate warming and lead to high spatial heterogeneity. More data (e.g., estimated pixel values regarding cultivar and agronomic measures) should be obtained to assess and verify the non-climate effects on crop phenology. Moreover, the impact of interaction between climate warming and nonclimate effects should be considered in future studies.