## 2 Materials and Methods
### NDVI-NDWI Phase-Space Construction  
High spatiotemporal resolution remote sensing data are essential in order to detect the slight changes in rice development caused by heavy metal stress. We used the enhanced spatial and temporal adaptive reflectance fusion model (ESTARFM) [40] to obtain a time series of the 8-day surface reflectance product (30 m) by blending the selected bands of Landsat 8 and MODIS data. ESTARFM focuses on improving data fusion for the mixed pixels and has made significant improvements in the accuracy of predictions for heterogeneous landscapes, which exist widely in our study area [41]. The ESTARFM algorithm uses two pairs of Landsat-MODIS images (i.e., dates \(t_{1}\) and \(t_{3}\)) and one MODIS image (date \(t_{2}\)) to predict a Landsat-like image at date \(t_{2}\). In this study, a total of 79 fine-resolution images (including 46 predicted images and the 33 original Landsat images) were obtained (Figure 2). We applied cloud and \begin{table}
\begin{tabular}{l l c c c c c} \hline \multirow{2}{*}{**Band Name**} & \multicolumn{3}{c}{**Landsat 8**} & \multicolumn{3}{c}{**MOD09A1**} \\ \cline{2-7}  & **Band No.** & **Wavelength (\(\mu\)m)** & **Resolution (m)** & **Band No.** & **Wavelength (\(\mu\)m)** & **Resolution (m)** \\ \hline Blue band & Band 2 & 0.450–0.510 & 30 & Band 3 & 0.459–0.479 & 500 \\ Green band & Band 3 & 0.530–0.590 & 30 & Band 4 & 0.545–0.565 & 500 \\ Red band & Band 4 & 0.640–0.670 & 30 & Band 1 & 0.620–0.670 & 500 \\ NIR band & Band 5 & 0.850–0.880 & 30 & Band 2 & 0.841–0.876 & 500 \\ SWIR band & Band 6 & 1.570–1.650 & 30 & Band 6 & 1.628–1.652 & 500 \\ \hline \end{tabular}
\end{table} 
Table 2: Brief characteristics of LANDSAT 8- and MODIS-selected bands. Figure 2: Scene acquisition dates by year for the Landsat, MODIS, and enhanced spatial and temporal adaptive reflectance fusion model (ESTARFM)-predicted images. cloud shadow masks generated from the pixel quality assurance in the Landsat product to all synthetic images and removed them from the overall analysis. NDVI, which is correlated with the absorption of photosynthetically active radiation, was calculated as follows [42]:  
\[NDVI=\frac{\rho_{NIR}-\rho_{Rad}}{\rho_{NIR}+\rho_{Red}}, \tag{2}\]  
where \(\rho_{NIR}\), \(\rho_{Rad}\) and \(\rho_{Blue}\) are the reflectances of the NIR band, red band, and blue band, respectively. NDWI, which is correlated with the canopy water content, was calculated as follows [43]:  
\[NDWI=\frac{\rho_{NIR}-\rho_{SWIR}}{\rho_{NIR}+\rho_{SWIR}}, \tag{3}\]  
where \(\rho_{SWIR}\) is the reflectance of the SWIR band. The application of satellite-derived VI time series is inevitably hindered by noise arising from clouds, ozone, dust, and other aerosols. Though the cloud-polluted pixels were masked out, radiance disturbance still existed in the VI time series. To solve this problem, the VI time series were fitted using a double-logistic curve as follows [44]:  
\[v\left(t\right)=a+b\left(\frac{1}{1+e^{ct+d}}+\frac{1}{1+e^{ct+f}}\right), \tag{4}\]  
where \(v(t)\) represents the value at time \(t\) (d), and \(a\), \(b\), \(\ldots\), \(f\) are the fitting parameters. The daily NDVI and NDWI time series that covered the rice growth season (from DOY 120 to DOY 305) were then generated by the fitting function (Equation (4)). We gave an example of phase-space established by the greenness and wetness VI, of which the greenness VI takes the NDVI as an example and the wetness VI refers to the NDWI (Figure 3). In the NDVI-NDWI phase-space, the temporal trajectory of a rice pixel throughout the whole growth period first drifts towards the top right areas and then towards the bottom left areas with the progression.