## 1 Introduction  
Heavy metal contamination in agricultural soil is one of the most severe problems affecting food security and public health at global and local levels [1]. Generally, pollution is caused against a backdrop of urban industrialization and is commonly characterized as relatively stable in spatial and temporal distribution [2; 3]. A number of agronomy studies have exhibited the variations in spectral properties, pigment content, dry matter, photosynthesis, and transpiration of plants under heavy metal stress [4; 5; 6]. Based on these symptoms, researchers have made attempts to (1) extract the characteristics of multispectral and hyperspectral signatures [7; 8; 9]; (2) retrieve physiological and biochemical parameters [10; 11]; and (3) integrate remote sensing data and crop growth models [12; 13; 14; 15] to detect heavy metal contamination in agricultural soils. Since satellite remote sensing has enabled the continuous monitoring of crop growth, it provides a new insight to indicate environmental interactions in agricultural systems through variations in crop phenology [16; 17; 18; 19]. Heavy metals affect crop development mainly by interfering with physiological activity, and it has been observed that the onset of heading is delayed when rice is exposed to cadmium [4; 6]. Basically, the change of the growth stages will certainly impact the vegetation index (VI) time series, thereby reflecting in the satellite-derived phenometrics. However, the detection of heavy metal stress-induced variations in rice phenology through the remote sensing approach is still problematic. One reason is the temporal mismatch between the satellite-derived rice phenometrics and the actual phenological events. Typically, a satellite-based photometric is the certain date of a feature point on the VI curve, while an actual phenological event refers to a certain stage in rice development. Moreover, for rice agroecosystems, the absolute dates do not mean much in identifying the stress-induced change in the growth state, because they are mainly decided by the timing of transplanting. [PERSON] et al. [20] analyzed the differences in the lengths of phenophases, which are defined as the period between two specific rice phenometrics, across a heavy metal stress gradient. However, the length of a phenophase may also be influenced by meteorological factors (e.g., air temperature and solar radiation), soil conditions, and management practices [21; 22; 23]. Therefore, more effort is needed to explore the relationships between the delaying of the actual heading stage and the variations of satellite-derived phenometrics under heavy metal stress. Optical observations of crop phenology have mostly used greenness VIs, such as the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI) [24; 25; 26; 27]. However, they may not identify the greenness-independent features, such as canopy water content, which could also be sensitive to the phenophase transition. The normalized difference water index (NDWI) has also been used to characterize the phenological cycle of vegetation--however, mostly for reducing the effect of snow in boreal natural ecosystems [28; 29; 30]. In addition to a single-purpose VI, [PERSON] et al. [31] combined the advantages of NDVI and NDWI to detect vegetation phenology with the small effect of soil and snow. Likewise, [PERSON] et al. [32] derived phenometrics in a seasonally snow-covered landscape using a phase-space diagram, which was constructed by the NDVI-NDWI time series over an entire phenological year. The observed pixels in the phase-space were segmented into snow-covered and snow-free phases by K-means clustering, which made the phenological descriptors (e.g., start of the growing season) free of the influence of the onset/offset of snow. Moreover, [PERSON] et al.