## 1 Introduction  
Phenology is a key and sensitive indicator of a terrestrial ecosystem's response to climate change through its contributions to global carbon, water, and energy cycles. Phenological estimates and predictions are thus crucial for understanding the interactions between climate and the biosphere [1; 2; 3]. Phenology is often measured at different scales, from species-scale to the ecosystem-scale [4]. Among the different scales, Land surface phenology (LSP), which was defined as the study of the timing of recurring seasonal pattern of variation in vegetated land surfaces observed from synoptic sensors [5], is measurable by remote sensing methods [6; 7; 8]. However, the scale disparity between the spatial extent of the organisms and the spatial resolution of the sensor results in an ambiguous mixture of target and background or signal and noise [9; 10]. LSP metrics are often compared with direct observations of physical changes in the canopy structure. Eddy covariance (EC) measurements of carbon exchange offer a spatially and temporally perspective for extracting key phenological dates by assessing seasonal cycles of gross primary productivity (GPP) [8; 11]. The long-term durability of CO\({}_{2}\) flux measurements in EC tower sites worldwide allows comparisons with remote sensing-based estimates at spatial footprints similar to coarse and medium resolution satellite pixels [4; 12; 13]. Productive efforts have been made to improve the estimates of phenological transitions based on remote sensing observations [12; 14; 15; 16], using new methods such as digital cameras [8; 17], developing improved vegetation indices, and applying remotely sensed solar-induced fluorescence (SIF) from the Global Ozone Monitoring Experiment-2 (GOME-2) and the Orbiting Carbon Observatory-2 (OCO-2) [18; 19; 20]. Optical remote sensing observations are widely applied to estimate key phenological dates, often using vegetation indices (VIs) based on a combination of reflectances in different spectral bands. The most commonly applied VI is the normalized difference vegetation index (NDVI), which involves red and near-infrared reflectances. The NDVI is sufficiently stable to permit meaningful comparisons of seasonal dynamics and inter-annual variations of vegetation, and its ratio concept reduces many sources of noise. However, it also has significant drawbacks, including saturation over high-biomass areas [21]. The enhanced vegetation index (EVI) was developed to optimize the vegetation signal with improved sensitivity in high-biomass regions and to further reduce atmospheric influences. EVI is more sensitive than NDVI to structural canopy variations. However, uncertainties still exist over the use of traditional VIs to extract phenology metrics [8; 14; 22; 23; 24]. Unlike traditional VIs, the phenology index (PI) was developed specifically to monitor vegetation phenology. It combines the merits of NDVI and the normalized difference infrared index (NDII), and aims to remove the confounding effects of brightness and wetness from greenness signals [4; 5]. Furthermore, a physical-based plant phenology index (PPI) was proposed to improve the satellite-based observations of vegetation phenology in boreal regions. Compared with traditional VIs, PPI was reported to be less influenced by snow and soil brightness, thus making it better suited for phenology retrieval in boreal forests. However, for the pan-Arctic regions, the use of VIs to estimate land surface phenology were rarely validated [25; 26]. It is therefore necessary to compare the performances of different VIs for extracting phenological metrics across different landscapes. Phenological metrics, i.e., the start of season (SOS) and end of season (EOS), can be extracted from seasonal dynamics using fitting functions and models [14; 27]. Although a simple logistic function has been used extensively in phenological studies, it has limitations for representing changes in greenness from extreme events [28].