Inﬂuence of Varying Solar Zenith Angles on Land Surface Phenology Derived from Vegetation Indices: A Case Study in the Harvard Forest

: Vegetation indices are widely used to derive land surface phenology (LSP). However, due to inconsistent illumination geometries, reﬂectance varies with solar zenith angles (SZA), which in turn affects the vegetation indices, and thus the derived LSP. To examine the SZA effect on LSP, the MODIS bidirectional reﬂectance distribution function (BRDF) product and a BRDF model were employed to derive LSPs under several constant SZAs (i.e., 0 ◦ , 15 ◦ , 30 ◦ , 45 ◦ , and 60 ◦ ) in the Harvard Forest, Massachusetts, USA. The LSPs derived under varying SZAs from the MODIS nadir BRDF-adjusted reﬂectance (NBAR) and MODIS vegetation index products were used as baselines. The results show that with increasing SZA, NDVI increases but EVI decreases. The magnitude of SZA-induced NDVI/EVI changes suggests that EVI is more sensitive to varying SZAs than NDVI. NDVI and EVI are comparable in deriving the start of season (SOS), but EVI is more accurate when deriving the end of season (EOS). Speciﬁcally, NDVI/EVI-derived SOSs are relatively close to those derived from ground measurements, with an absolute mean difference of 8.01 days for NDVI-derived SOSs and 9.07 days for EVI-derived SOSs over ten years. However, a considerable lag exists for EOSs derived from vegetation indices, especially from the NDVI time series, with an absolute mean difference of 14.67 days relative to that derived from ground measurements. The SOSs derived from NDVI time series are generally earlier, while those from EVI time series are delayed. In contrast, the EOSs derived from NDVI time series are delayed; those derived from the simulated EVI time series under a ﬁxed illumination geometry are also delayed, but those derived from the products with varying illumination geometries (i.e., MODIS NBAR product and MODIS vegetation index product) are advanced. LSPs derived from varying illumination geometries could lead to a difference spanning from a few days to a month in this case study, which highlights the importance of normalizing the illumination geometry when deriving LSP from NDVI/EVI time series. effect on the retrieved LSP. BRDF parameters provided by the MCD43A1 (BRDF) product were used as inputs in the RTLSR BRDF model to simulate reﬂectance at a constant SZA conﬁguration. The MCD43A4 (NBAR) and MOD13A1 (VI) products were used as baselines. Our results illustrate that both NDVI and EVI time series are sensitive to varying SZAs but in an opposite manner: NDVI increases with increasing SZAs, while EVI decreases. We also compared the retrieved LSP to ground measurements, ﬁnding limited strength of remote sensing data in retrieving EOS compared to retrieving SOS. Although NDVI and EVI have equivalent accuracy in deriving SOS, EVI demonstrates advantages in predicting EOS. The discrepancies between these two vegetation indices in terms of deriving LSP are also reﬂected in whether the retrieved phenology is advanced or delayed. Speciﬁcally, SOSs derived from NDVI time series are advanced, yet from EVI time series are delayed. However, EOSs derived from NDVI time series are delayed, and from simulated EVI time series under consistent illumination geometry are also delayed, but EOSs derived from baseline products (i.e., MCD43A4 (NBAR) and MOD13A1 (VI)) are generally advanced. Additionally, the LSP derived from varying illumination geometries could lead to a difference spanning from a few days to a month in this case study, which has important implications for normalizing the illumination geometries in deriving LSP from remote sensing data. End-users of multi-angle remote sensing products, therefore, are suggested to take the SZA effect into consideration.


Introduction
The plant growth cycle, including the germination, growth, flowering, fruit-bearing, aging, and death stages, is accompanied by changes in plant physiology and appearance [1]. The physical, chemical, and biological properties of plants exhibit seasonal variations that can be described using the spectral characteristics of vegetation [2][3][4]. These changes are considered as a plant's phenology and can be effectively monitored with remote sensing technologies [5]. Vegetation phenology can be used to detect the footprint of climate change, guide farming activities, and alert the appropriate parties of the potential outbreak of allergic diseases [6][7][8].
Phenological information is commonly derived from spectral vegetation indices (VIs), which combine signals of multiple spectral bands. Recent studies also focus on the angular vegetation indices that have been successfully used to retrieve various vegetation structural parameters [9][10][11]. However, the angular vegetation indices characterize three-dimensional vegetation structures [12], which have limited help in deriving land surface phenology (LSP) because phenological phases are closely related to variations in leaf pigments, and these variations are efficiently indicated by the spectral vegetation indices. Among the spectral vegetation indices, the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI) are broadly used in regional or global phenology studies [13][14][15][16][17]. By analyzing NDVI and EVI time series, the patterns of vegetation growth, their seasonal characteristics, and their interannual variations can be derived.
To accurately capture phenological information from vegetation indices, previous works focus on improving phenology detection methods [5,[18][19][20][21], evaluating the performance of different vegetation indices [22][23][24], investigating the spatial and temporal scaling effect on LSP metrics [25][26][27][28], and exploring the multilevel validation approach [29][30][31]. However, less attention has been paid to the influence of varying illumination geometries, which might make the retrieved LSP metrics deviated from their true dates. The inconsistency of illumination geometry results from several factors, especially when long-timespan and multiangle data are accessible at present [32]. First, the solar zenith angle (SZA) at local solar noon (LSN) is not constant throughout the year but exhibits seasonal amplitudes [33]. Second, SZA not only changes with seasons but also changes greatly throughout the day. The problem of inconsistent SZAs is exacerbated by the improved temporal resolution of remote sensing images from geostationary or equatorial-orbit satellites that capture several images per day [34,35]. Third, the assumption of a constant illumination geometry is jeopardized by the effect of satellite orbit drift [36]. For example, the drifts associated with the Landsat and EO-1 satellites have modified the time of image acquisition and, therefore, the acquisition of SZA [37,38]. Fourth, high-relief topography and changes in terrain slopes and aspects intensify inconsistent solar illumination, thereby affecting the retrieved LSP [39][40][41]. Fifth, multi-sensor data integration raises issues associated with incorporating the potential effects of inconsistent SZAs in the time series because different products provide inconsistent angle configurations [42][43][44][45][46].
Ignoring changes in the SZA affects the accuracy of vegetation indices, which inevitably leads to an inaccurate assessment of the growth conditions [47]. Bhandari et al. [47] verified that varying illumination geometries had impacts on spectral vegetation indices and further on phenology and disturbance detection, pointing out that NDVI time series was more effective in capturing such changes. Ji and Brown [48] testified how orbital-driftinduced SZA change impacted LSP metrics, noting that these two factors were significantly correlated in 4.1-20.4% of the contiguous United States area. Vegetation indices also show different sensitivities to varying SZA. Breunig et al. [49] studied how SZA impacted the spectral vegetation indices individually and synergistically with leaf area index, revealing the different dependence of NDVI and EVI on sun-sensor geometry.
Despite the progress made so far, fewer studies have directly explored to what extent varying SZAs impacted the derivation of LSP metrics, especially at high latitudes where SZA at LSN shows prominent variations during a year. However, the influence of the changing illumination geometries should be taken into consideration before retrieving Remote Sens. 2021, 13, 4126 3 of 22 phenological information from vegetation indices time series, especially when end-users tend to ignore the uncertainties coming from the inconsistent illumination geometries that are inherently incorporated in remotely sensed observations. The specific objectives of this study are to investigate the impact of changing illumination geometries on NDVI and EVI profiles as well as the cascading effect on retrieved LSP metrics, and to assess the associated uncertainties. The semi-empirical RossThick-LiSparse Reciprocal (RTLSR) bidirectional reflectance distribution function (BRDF) model was adopted to simulate surface reflectance under a set of constant SZAs (i.e., 0 • , 15 • , 30 • , 45 • , and 60 • ), with the MCD43A1 (BRDF) product providing kernel parameters as inputs. Then LSP metrics under the five constant SZAs were retrieved based on the simulated surface reflectance. LSP metrics retrieved from the MOD13A1 (VI) and the MCD43A4 (NBAR) products were adopted as baselines because those two products are broadly used in phenology studies yet without correcting angle effects for the former, and with view angle corrected to the nadir but solar angle varies at local solar noon for the latter.

Study Area
In this study, we would like to explore whether a large variation of SZA throughout the year impacts the vegetation indices and exerts further influences on derived LSP, so we need a research area located at a high latitude where SZA at LSN shows prominent variations during a year. The average annual precipitation is 1100 mm, and precipitation remains at approximately the same level throughout the year. We used an image chip (3 by 3 MODIS 500 m pixel window) as the study site to match the area of ground observations in the Harvard Forest. We chose the image chip centered at 42.535 • N and 72.185 • W to guarantee that the ground measurements (i.e., Data Archive HF-003) are located within the image chip and thus can be used to validate the LSP metrics derived from remote sensing data. In addition, choosing image chips instead of pixels can minimize the pixel-to-pixel difference, and therefore has been adopted in several angle effect studies [50,51].

Ground Measurements
Three types of spring phenology are recorded in the Harvard Forest [52][53][54]: the percentage of buds on the trees that have broken open (BBRK), the percentage of leaves on the tree that is at least 75% of their final size (L75), and the percentage of leaves on the tree that is greater than or equal to 95% of their final size (L95). Two types of autumn phenology data are also provided: one is the proportion of leaves remaining on the tree that have changed color (LCOLOR), another is the proportion of defoliation (LFALL). In this study, BBRK30 (days when 30% percent of buds on the tree have broken open) and LFALL50 (50% leaves have fallen from a given tree) are used to verify the accuracy of NDVI/EVI-derived LSPs.
The Harvard Forest site is dominated by deciduous species, including red maple (Acer Rubrum, accounts for 22% basal area) and red oak (Quercus rubra, accounts for 52% basal area) [55], and the averaged phenology dates of those two species are widely used in previous studies as ground phenology references [56][57][58][59][60]. Good relative consistency exists among species and among individuals within species in the Harvard Forest, so we consider the mean phenology dates of five red maple trees and four red oak trees qualified to validate the LSP derived from remote sensing data.  [61,62] and its quality control dataset MCD43A2 [63] were used as inputs to simulate reflectance under a set of constant SZAs. Although MODIS instruments are onboard the Terra and Aqua satellites (with Terra overpasses at 10:30 A.M. and 10:30 P.M. local time and Aqua overpasses at 1:30 A.M. and 1:30 P.M. local time), BRDF parameters from the MCD43A1 product allow simulating surface reflectance under any angle configurations based on a BRDF model, which will be introduced in Section 3.1. We also used the MCD43A4 (NBAR) product (16-day composition with an 8-day interval) as a baseline because it has been broadly used as source data in phenology studies [5,33,64,65]. Although its view zenith angle is normalized to nadir, its solar zenith angle is adjusted to the local solar noon of the day, suggesting that the impact of varying view zenith angles is eliminated yet the impact of the seasonal variation in the solar zenith angle still exists. We included the MOD13A1 (VI) product as well because it has been broadly used in phenology studies [24,66] but the varying angle configurations have not been corrected. The MOD13A1 (VI) product uses the constraint view angle-maximum value composite (CV-MVC) method as the primary compositing algorithm [67,68], which selects the observation with a smaller view angle between the two highest vegetation index values to minimize the impact of large view zenith angle, but the view zenith angles are still inconsistent over time. Besides, due to the lack of restrictions on solar zenith angles, the solar direction of the chosen observation can be either from the zenith or deviated. Therefore, the performance of the MOD13A1 (VI) product in deriving LSP metrics should be examined. To minimize the impacts of clouds and snow on the reflectance in the Harvard Forest, we only kept the pixels with high-quality full inversion and snow-free BRDF parameters retrieved. Full inversion means that sufficient high-quality observations are available when performing matrix inversion to compute the kernel parameters in the MODIS BRDF/Albedo product (i.e., MCD43A1) [69][70][71]. Although the MCD43 V006 product is retrieved on a daily basis using 16-day observations while the MCD43 V005 product is retrieved every eight days, we chose Version 005 instead of Version 006 because we considered eight days as a reasonable time interval for generating the 10-year NDVI/EVI time series that is used to extract LSP metrics. In addition, the angle effect is consistent from Version 005 to Version 006 because the latest product does not improve upon the elimination of the influence of inconsistent illumination geometries, indicating that our study with Version 005 is also enlightening for studies using Version 006.

Methods
In this study, to compare LSPs derived from simulations with constant SZAs and from baseline products with inconsistent SZAs, four major processes are required, including NDVI/EVI time series simulation based on a kernel-driven model, NDVI/EVI time series reconstruction, land surface phenology derivation, and accuracy assessment ( Figure 1). series reconstruction, land surface phenology derivation, and accuracy assessment (Figure 1).

NDVI/EVI Time Series Simulation Based on the Kernel-Driven Model
We used the semiempirical RossThick-LiSparseReciprocal (RTLSR) BRDF model [12,61,72,73] to simulate reflectance under fixed illumination geometries. The RTLSR BRDF model simplifies the directional distribution of land surface reflectance as a linear combination of three components, i.e., the geometric optical scattering, volumetric optical scattering, and isotropic scattering [74][75][76][77], denoted by where ( , , , ) is the reflectance in waveband under a designated viewing and illumination geometry. as anisotropic parameters. and are the geometric optical scattering kernel and the volumetric optical scattering kernel, respectively [74,78,79].
, , and are the view zenith angle, solar zenith angle, and relative azimuth angle, respectively. Given that the MCD43A1 (BRDF) product provides the three model weighting parameters ( ), ( ), and ( ), RTLSR can be used to standardize the illumination geometry and generate comparable phenology metrics. The performance of the RTLSR BRDF model has been fully evaluated in previous studies [12,72,73], and it has been adopted as the standard model in operational MODIS BRDF/Albedo/NBAR products [61,69,70,80].  Figure 1. The process chain of deriving LSP from simulations with constant SZAs and from baseline products with inconsistent SZAs. The green rectangles from top to the bottom represent NDVI/EVI time series simulation based on the kernel-driven model, NDVI/EVI time series reconstruction, land surface phenology derivation, and accuracy assessment, respectively. Blue solid rectangles are used products or derived results, while blue dashed rectangles are used models or integrated algorithms.

NDVI/EVI Time Series Simulation Based on the Kernel-Driven Model
We used the semiempirical RossThick-LiSparseReciprocal (RTLSR) BRDF model [12,61,72,73] to simulate reflectance under fixed illumination geometries. The RTLSR BRDF model simplifies the directional distribution of land surface reflectance as a linear combination of three components, i.e., the geometric optical scattering, volumetric optical scattering, and isotropic scattering [74][75][76][77], denoted by where ρ(θ v , θ s , ϕ, λ) is the reflectance in waveband λ under a designated viewing and illumination geometry. f iso (λ), f geo (λ), and f vol (λ) are spectrally dependent model parameters, with f iso (λ) denoting an isotropic parameter, and f geo (λ), and f vol (λ) serving as anisotropic parameters. K geo and K vol are the geometric optical scattering kernel and the volumetric optical scattering kernel, respectively [74,78,79]. θ v , θ s , and ϕ are the view zenith angle, solar zenith angle, and relative azimuth angle, respectively. Given that the MCD43A1 (BRDF) product provides the three model weighting parameters f iso (λ), f geo (λ), and f vol (λ), RTLSR can be used to standardize the illumination geometry and generate comparable phenology metrics. The performance of the RTLSR BRDF model has been fully evaluated in previous studies [12,72,73], and it has been adopted as the standard model in operational MODIS BRDF/Albedo/NBAR products [61,69,70,80]. By verifying the longitude and latitude of the research site, we read the BRDF model parameters from the MCD43A1 (BRDF) product from 2001 to 2010. We also used the MCD43A2 (QA) product to remove the data contaminated by clouds, cloud shadows, and snow. The eligible BRDF model parameters were then put into the RTLSR BRDF model to simulate the surface reflectance with both the view zenith angle and the relative azimuth angle set to 0 • , but the SZA was set to 0 • , 15 • , 30 • , 45 • , or 60 • to test the solar angle effect. Although the real SZA at LSN in the Harvard Forest varies from 19 • to 65 • , variations in the illumination geometry come not only from the seasonal amplitude of the SZA at LSN [33] but also from the improved temporal resolution [34,35], satellite orbit drift [37,38], different data source integration [42][43][44], and terrain illumination [39][40][41], etc. Because of the abovementioned factors, the coupled impact could lead to a larger SZA range. Therefore, this study focused on SZAs spanning from 0 • to 60 • , with an interval of 15 • . We start this range from 0 • because we also want to analyze the impact of the hot spot, where the directional reflectance at nadir view is near to the illumination direction and therefore results in a spike in the reflectance profile [69,81]. Besides, we consider that the measured reflectance under SZAs larger than 60 • is less accurate. We set a constant SZA throughout the year to eliminate the confounding effect of varying SZAs on deriving LSPs. To demonstrate that the impact of varying SZAs on the derived LSPs was not a special case in a certain year, we used a ten-year time series as a repeated annual experiment to strengthen the conclusion.
We selected NDVI and EVI in our study to derive LSP metrics because these two vegetation indices have been shown to be efficient indicators of vegetation growth conditions and have been largely used in phenology studies [13][14][15][16]. Based on Equations (2) and (3), the ten-year NDVI/EVI time series are obtained.
R NIR , R RED and R BLUE are the reflectance in the near-infrared band, red band, and blue band, respectively.

Reconstructing NDVI/EVI Time Series
LSP metrics were derived from continuous NDVI/EVI time series. However, the removal of unqualified values resulting from clouds, cloud shadows, and snow led to incomplete NDVI/EVI time series [82]. Therefore, we reconstructed the NDVI/EVI profiles in TIMESAT [83]. TIMESAT embeds three methods in its time-series-smoothing module, including adaptive Savitzky-Golay filter, asymmetric Gaussians, and double logistic functions [83], from which we chose both the SG filter and the asymmetric Gaussians. For NDVI/EVI time series with sufficient signals after removing unqualified values, the SG filter has the highest fidelity and is able to retain more details in the NDVI/EVI time series [84]; for NDVI/EVI time series with more than 3 consecutive months without signals, we used the asymmetric Gaussian due to its better performance in capturing the overall seasonal structure and in demonstrating less sensitivity to incomplete time-series data [85,86].

Deriving Land Surface Phenology
The reconstructed NDVI/EVI time series were then used to derive LSP metrics, including the start of season (SOS) and the end of season (EOS) [87], with SOS and EOS referring to the point at which the curve intersected a proportion of the seasonal amplitude measured from the left minimum and from the right minimum, respectively [51,[88][89][90][91][92][93].
TIMESAT has been recognized as an efficient and robust software in phenology studies [94]. With its graphical user interface (GUI), TIMESAT is user-friendly by setting key parameters, e.g., the size of the moving window and the number of iterations [57].
To determine the SOS and EOS in the Harvard Forest from remote sensing data, we set the threshold at 0.5 to find the point at which the curve intersected 50% of the seasonal amplitude. For the SOS, the exposure of the understory reduces the NDVI at the beginning of each year, which extends the amplitude of the NDVI/EVI time series from the left side, Remote Sens. 2021, 13, 4126 7 of 22 so 50% of the amplitude closely matches the ground validations determined by BBRK30, as we introduced in Section 2.2.1. For the EOS, 50% of the amplitude from the right base level corresponds to the ground phenology metric of LFALL50.

Accuracy Assessment
The accuracy of the LSP derived from sensor data is measured by RMSE (root mean square error) and Spearman's rank correlation coefficient. RMSE is used to measure the difference between the VI-derived LSP and the ground measurements, as shown in Equation (4).
where (DOY) sensor and (DOY) ground denote the LSP metrics (SOS or EOS) derived from sensor data and ground measurements, respectively. N indicates the total year numbers in the time series, which is 10 in this study.
In contrast to RMSE, which measures how close the VI-derived LSPs are to ground measurements, Spearman's rank correlation coefficient describes whether the trend of the 10-year time series in the VI-derived LSP aligns well with the trend in the ground-measured LSP. A large positive rank correlation coefficient indicates that the trend of the VI-derived LSP is consistent with the trend of the ground-measured LSP. Spearman's rank correlation coefficient is calculated using Equation (5).
Here, ρ denotes the Spearman's rank correlation coefficient; N is 10 in this study. d i is computed by sorting the VI-derived LSPs and ground-measured LSPs, calculating the difference between those two sorted arrays; the difference between corresponding elements in the two arrays is d i . NDVI/EVI time series from the MCD43A4 (NBAR) product (cyan line) is close to the time series simulated at 15° SZA (orange line); when the SZA is approximately 60°, the cyan line is close to the time series simulated at 60° SZA (green line) (Figure 2). This pattern not only demonstrates that the MODIS BRDF product has a high degree of accuracy in simulating the surface reflectance under a certain angle configuration but also shows the impact of varying SZAs on NDVI/EVI time series.   Figure 3 suggests that the NDVI-derived SOS tends to underestimate the SOS compared to ground measurements [95,96]. In Table 2, when SZA is 0 • and 60 • , the RMSE is high (11.4 and 10.8, respectively), which indicates discrepancies between NDVI-derived SOSs and ground-measured SOSs. The rank correlation between NDVI-derived SOSs and ground-measured SOSs is only 0.285 when SZA equals 60 • , suggesting a large difference between the ground-measured SOS and the NDVI-derived SOS under 60 • SZA in terms of the trend of SOS across the 10-year period. When the SZA equals 30 • , the retrieved SOS has the best performance, exhibiting a lower RMSE (8.1) and a higher correlation (0.539). This could be explained by the fact that the actual SZA at the date of SOS in each year is close to 30 • . Specifically, when averaging the corresponding SZA of each SOS, the 10-year averaged actual SZA is 28.23 • , which is close to the 30 • SZA we set in this angle configuration.

SOS Derived from NDVI Time Series
has the best performance, exhibiting a lower RMSE (8.1) and a higher correlation This could be explained by the fact that the actual SZA at the date of SOS in each close to 30°. Specifically, when averaging the corresponding SZA of each SOS, the 1 averaged actual SZA is 28.23°, which is close to the 30° SZA we set in this angle co ration.   Figure 4 reveals that EVI-derived SOSs are delayed relative to the ground-me SOSs, which is inverse to the pattern observed from NDVI-derived SOSs. EVI-d SOSs show fewer discrepancies across different SZA configurations, with a ma RMSE of 11.1 days, and a minimum RMSE of 10.0 days ( Table 3). The Spearman   Figure 4 reveals that EVI-derived SOSs are delayed relative to the ground-measured SOSs, which is inverse to the pattern observed from NDVI-derived SOSs. EVI-derived SOSs show fewer discrepancies across different SZA configurations, with a maximum RMSE of 11.1 days, and a minimum RMSE of 10.0 days ( Table 3). The Spearman's rank correlation coefficient (ρ) also demonstrates a good rank correlation between EVI-derived SOSs and ground-measured SOSs, excluding SOSs retrieved under 60 • SZA and retrieved directly from the MOD13A1 (VI) product.
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 23 correlation coefficient (ρ) also demonstrates a good rank correlation between EVI-derived SOSs and ground-measured SOSs, excluding SOSs retrieved under 60° SZA and retrieved directly from the MOD13A1 (VI) product.

Comparison of SOS Derived from NDVI and EVI Time Series
Comparisons of the SOSs derived from the NDVI and EVI time series from 2001 to 2010 are shown in Figure 5. The x-axis represents the NDVI-derived SOS, and the y-axis indicates the EVI-derived SOS. The SOSs derived from MCD43A1 (BRDF) product are distributed along the diagonal and exhibit a higher correlation (0.57 < R 2 < 0.88), indicating that the NDVI-derived SOSs and EVI-derived SOSs under a constant SZA are relatively close; the SOSs derived from the MCD43A4 (NBAR) and MOD13A1 (VI) products show a relatively low correlation, with an R 2 of 0.49 for the former and only 0.09 for the latter. In addition, the distribution of points occurs mainly above the diagonal line, suggesting that EVI-derived SOSs generally occur later than NDVI-derived SOSs.

Comparison of SOS Derived from NDVI and EVI Time Series
Comparisons of the SOSs derived from the NDVI and EVI time series from 2001 to 2010 are shown in Figure 5. The x-axis represents the NDVI-derived SOS, and the y-axis indicates the EVI-derived SOS. The SOSs derived from MCD43A1 (BRDF) product are distributed along the diagonal and exhibit a higher correlation (0.57 < R 2 < 0.88), indicating that the NDVI-derived SOSs and EVI-derived SOSs under a constant SZA are relatively close; the SOSs derived from the MCD43A4 (NBAR) and MOD13A1 (VI) products show a relatively low correlation, with an R 2 of 0.49 for the former and only 0.09 for the latter. In addition, the distribution of points occurs mainly above the diagonal line, suggesting that EVI-derived SOSs generally occur later than NDVI-derived SOSs.

EOS Derived from NDVI Time Series
In contrast with NDVI-derived SOS, NDVI-derived EOS generally overestimates the EOS relative to that of the ground measurements, excluding 2006. From Table 4, the NDVI-derived EOS shows a large RMSE when the SZA equals 60°, indicating a poor performance of the simulated NDVI time series in deriving EOS under a large SZA. Compared to the rank correlation between NDVI-derived SOSs and ground-measured SOSs (Table 2), the rank correlation between NDVI-derived EOSs and ground-measured EOSs is low (Table 4), aligning well with the pattern shown in Figure 6, where the trend of ground-measured EOS over ten years shows less variation, yet the NDVI-derived EOS demonstrates drastic variation throughout the evaluated years.

EOS Derived from NDVI Time Series
In contrast with NDVI-derived SOS, NDVI-derived EOS generally overestimates the EOS relative to that of the ground measurements, excluding 2006. From Table 4, the NDVIderived EOS shows a large RMSE when the SZA equals 60 • , indicating a poor performance of the simulated NDVI time series in deriving EOS under a large SZA. Compared to the rank correlation between NDVI-derived SOSs and ground-measured SOSs (Table 2), the rank correlation between NDVI-derived EOSs and ground-measured EOSs is low (Table 4), aligning well with the pattern shown in Figure 6, where the trend of ground-measured EOS over ten years shows less variation, yet the NDVI-derived EOS demonstrates drastic variation throughout the evaluated years.

EOS Derived from EVI Time Series
The EOS derived from the MCD43A4 (NBAR) and MOD13A1 (VI) products are generally underestimated in contrast to the EOS derived from the simulated EVI time series under constant illumination geometry (Figure 7). Similar to NDVI-derived EOSs, EVI-derived EOSs are substantially advanced in 2006 relative to the ground-measured EOS, with an average underestimation of 13.4 days. From Tables 4 and 5, the RMSE and Spearman's rank correlation coefficients are generally improved in comparison with NDVI-derived EOSs, suggesting the better performance of the EVI time series to derive EOSs.

EOS Derived from EVI Time Series
The EOS derived from the MCD43A4 (NBAR) and MOD13A1 (VI) products are generally underestimated in contrast to the EOS derived from the simulated EVI time series under constant illumination geometry (Figure 7). Similar to NDVI-derived EOSs, EVI-derived EOSs are substantially advanced in 2006 relative to the ground-measured EOS, with an average underestimation of 13.4 days. From Tables 4 and 5, the RMSE and Spearman's rank correlation coefficients are generally improved in comparison with NDVI-derived EOSs, suggesting the better performance of the EVI time series to derive EOSs.

EOS Derived from EVI Time Series
The EOS derived from the MCD43A4 (NBAR) and MOD13A1 (VI) products are generally underestimated in contrast to the EOS derived from the simulated EVI time series under constant illumination geometry (Figure 7). Similar to NDVI-derived EOSs, EVI-derived EOSs are substantially advanced in 2006 relative to the ground-measured EOS, with an average underestimation of 13.4 days. From Tables 4 and 5, the RMSE and Spearman's rank correlation coefficients are generally improved in comparison with NDVI-derived EOSs, suggesting the better performance of the EVI time series to derive EOSs.

Comparison of EOS Derived from NDVI and EVI Time Series
A comparison between NDVI-derived EOSs and EVI-derived EOSs is shown in Figure 8, with a poor correlation compared to VI-derived SOSs. When the SZA is within the range from 15 • to 45 • , 0.45 < R 2 < 0.58, and when the SZA is 0 • and 60 • , R 2 is 0.07 and 0.21, respectively. The R 2 between retrieved EOSs from the MOD13 (VI) NDVI and the EVI time series is only 0.0004, indicating a compromised consistency. Additionally, the distribution of points is located below the diagonal, which indicates that EOSs derived from EVI are generally earlier than those derived from NDVI.

NDVI and EVI Time Series
The NDVI and EVI time series depict similar seasonal variations in the deciduous Harvard Forest, but discrepancies still exist. First, EVI profiles are lower than NDVI profiles. Huete et al. [67] indicated that the reflectance of the red and blue bands in a phenological period were relatively stable, but the near-infrared band exhibited clear seasonal characteristics: when the growing season reached its peak, the near-infrared curve also reached its highest value, so NDVI would be larger than EVI under the same vegetation condition and geometric structure. Second, the EVI profile shows less fluctuation at the mature stage relative to the NDVI profile due to its improvement in reducing atmospheric impacts and the influence of soil background. Third, in each summer season, except for the NDVI/EVI profiles derived from the MOD13A1 (VI) product, other NDVI/EVI profiles show clear drops, especially lines at 0° because of the impact of the hot spot. This phenomenon appears every year, presumably due to the concentration of precipitation because the presence of clouds in the summer leads to a sharp decrease in the NDVI.

Comparison of VI-Derived Land Surface Phenology Metrics
We find a limited strength of remote sensing data in retrieving EOS compared to SOS, which is consistent with previous studies [22,24,95,97,98], irrespective of sensors or plant functional types [96]. Testa et al. [24] compared the NDVI and EVI-derived LSP metrics of French deciduous forests, pointing out that SOSs are more accurate (±10 days) than EOSs (±23 days). For the SOS, the results from both the NDVI and EVI time series are relatively close to the ground measurements, with absolute mean differences of 8.01 days for the former (Figure 3) and 9.07 days for the latter (Figure 4). The RMSE of the difference between VI-derived SOSs and ground-measured SOSs is less than 12 (Tables 2 and 3). However, a considerable lag exists when deriving the EOS from NDVI time series, with 14.67 days as the absolute mean difference from the ground measurements across ten years ( Figure 6). Table 4 also shows that the NDVI-derived EOS drastically overestimates the EOS compared to the ground-measured EOS, with a maximum RMSE of 33.7. Multiple reasons contribute to the uncertainty in retrieving EOS from remote sensing data, including longer and slower changes in canopy greenness in fall [97,98], the temporal resolution

NDVI and EVI Time Series
The NDVI and EVI time series depict similar seasonal variations in the deciduous Harvard Forest, but discrepancies still exist. First, EVI profiles are lower than NDVI profiles. Huete et al. [67] indicated that the reflectance of the red and blue bands in a phenological period were relatively stable, but the near-infrared band exhibited clear seasonal characteristics: when the growing season reached its peak, the near-infrared curve also reached its highest value, so NDVI would be larger than EVI under the same vegetation condition and geometric structure. Second, the EVI profile shows less fluctuation at the mature stage relative to the NDVI profile due to its improvement in reducing atmospheric impacts and the influence of soil background. Third, in each summer season, except for the NDVI/EVI profiles derived from the MOD13A1 (VI) product, other NDVI/EVI profiles show clear drops, especially lines at 0 • because of the impact of the hot spot. This phenomenon appears every year, presumably due to the concentration of precipitation because the presence of clouds in the summer leads to a sharp decrease in the NDVI.

Comparison of VI-Derived Land Surface Phenology Metrics
We find a limited strength of remote sensing data in retrieving EOS compared to SOS, which is consistent with previous studies [22,24,95,97,98], irrespective of sensors or plant functional types [96]. Testa et al. [24] compared the NDVI and EVI-derived LSP metrics of French deciduous forests, pointing out that SOSs are more accurate (±10 days) than EOSs (±23 days). For the SOS, the results from both the NDVI and EVI time series are relatively close to the ground measurements, with absolute mean differences of 8.01 days for the former (Figure 3) and 9.07 days for the latter (Figure 4). The RMSE of the difference between VI-derived SOSs and ground-measured SOSs is less than 12 (Tables 2 and 3). However, a considerable lag exists when deriving the EOS from NDVI time series, with 14.67 days as the absolute mean difference from the ground measurements across ten years ( Figure 6). Table 4 also shows that the NDVI-derived EOS drastically overestimates the EOS compared to the ground-measured EOS, with a maximum RMSE of 33.7. Multiple reasons contribute to the uncertainty in retrieving EOS from remote sensing data, including longer and slower changes in canopy greenness in fall [97,98], the temporal resolution of sensor data [96], species and topography [99], and different sensitivities that different vegetation indices show to vegetate growth [100].
The LSP metrics derived from NDVI and EVI time series also show discrepancies. Specifically, NDVI-derived SOSs advance the green-up dates compared to the ground-measured SOSs (Figure 3), but delayed SOSs are retrieved from the EVI time series (Figure 4). This pattern is further verified in Figure 5, with the distribution of points essentially located above the diagonal, illustrating that the SOSs retrieved from EVI time series are later than those retrieved from NDVI time series [101,102]. For the VI-derived EOSs, the results derived from the EVI time series are generally earlier than those from the NDVI time series (Figure 8). Specifically, the EOSs derived from the NDVI time series overestimate the ground-measured EOSs (Figure 6), the EOSs derived from the simulated EVI time series under constant illumination geometries also overestimate, but the EOSs derived from EVI time series generated from baseline products (i.e., MCD43A4 (NBAR) and MOD13A1 (VI)) generally underestimate the ground-measured EOSs. Compared to the NDVI-derived EOS, the EVI-derived EOS has a better performance with a smaller RMSE and a higher Spearman's rank correlation coefficient (Tables 4 and 5), which echoes the previous finding that the LSP derived from different vegetation indices varies, suggesting different sensitivities of vegetation indices to vegetation growth [24,100].

Impacts from Illumination Geometry on Deriving Land Surface Phenology Metrics
The NDVI and EVI time series show an opposite pattern with increasing SZA: the larger the SZA, the higher the NDVI profiles, but the lower the EVI profiles. This result aligns well with previous studies [34,38]. The increase of SZA results in decreased reflectance in both red and near-infrared bands, but NDVI increases at large SZA because of a greater decrease in the red than in the near-infrared [103][104][105], while EVI decreases at large SZA because of greater NIR dependence [106][107][108]. In addition, the sensitivities of NDVI and EVI to varying SZAs are different. Specifically, EVI is more sensitive to varying SZAs than NDVI, especially in growing seasons, where the NDVI profiles are less scattered than the EVI profiles. The stronger sensitivity of EVI to varying SZAs echoes previous research conducted by Galvão et al. (2011). They compared 10 narrowband vegetation indices in terms of their sensitivity to view-illumination effects and found EVI was more anisotropic than NDVI. They also found shade fraction was negatively related to EVI. In this study, both view zenith angle and relative azimuth angle are set to 0, so shade fraction increases with the increase of SZA. Therefore, EVI decreases with the increase of SZA (Figure 2b).
The impact of varying SZAs on the NDVI/EVI time series leads to discrepancies in the retrieved LSP. For example, when using the NDVI time series to extract the SOS in the Harvard Forest in 2010, the results from the NDVI time series under a 0 • SZA and from the MCD43A4 (NBAR) product are DOY 86 and DOY 116, respectively, exhibiting a difference of 30 days (Figure 3). In 2007, the EOS derived from the MOD13A1 (NDVI) product is DOY324, whereas the mean EOS of other products that take the influence of the SZA into account is DOY284, which is closer to the ground measurement of DOY 278 ( Figure 6). Although the above evidence demonstrates that the SZA effect should be corrected uniformly when deriving LSP from NDVI/EVI time series, the different responses of the retrieved LSP to varying SZAs suggest the difficulty of applying a standard SZA correction method to all VI-derived LSPs. For example, the NDVI-derived SOS under a 30 • SZA has the best performance of all the VI-derived SOSs, with an RMSE of 8.1 and a Spearman's rank correlation coefficient of 0.539 (Table 2); however, the NDVI-derived EOS under a 30 • SZA shows a large discrepancy relative to the corresponding ground-measured EOS, with an RMSE of 21.3 and a Spearman's rank correlation coefficient of 0.321 (Table 4). This highlights the need of exploring how site-specific characteristics impact the correction of the SZA effect in future studies.
It is worth noting that we use MODIS products in this study because the MODIS BRDF/Albedo product provides a possibility to normalize the angle configurations based on the RTLSR BRDF model. However, for any dataset that has inconsistent angle configurations, the impact of inconsistent SZAs on the accuracy of retrieved LSPs should be considered.

Uncertainties
The VI-derived LSP under a set of constant SZAs is evaluated using ground observations, which raises the issue of scaling mismatch. This scaling issue has existed since the very beginning of remote sensing techniques, and satellite-based phenology studies were not immune to the spatial mismatch with the field data. However, the defects cannot obscure the virtues, satellite remote sensing products provide long-term, continuous, local-to global-scale, retrospective observations, which expands the horizon of traditional phenology [89,109,110]. Ground observations have also been widely used to validate satellite-derived LSP metrics [26,110,111]. In addition, most scaling issues come from the heterogeneity of the Earth's surface [26,27,89]; however, a good relative consistency exists among species and among individuals within a species in this study area. Therefore, we consider that the data gathered in the Harvard Forest could be used to validate our results, especially as this dataset is mature and qualified [53,54].
A discrepancy also exists between ground LSP measured from the overstory and from the understory, separating the over-and understory would enable more accurate analysis of LSP [112,113]. But this study focuses on the discrepancy among LSP derived under different SZAs, and this discrepancy exists due to inherent biases in angle configurations, even if without comparing to ground measurements. Compared to ground validation data helps to quantify how LSP derived from different remote sensing products or derived under different angle configurations deviate from each other. Although we focus more on differences in the LSP metrics resulting from SZA effects instead of optimizing the method to derive the LSP metrics that are closest to the ground measurements in this study, the uncertainties due to potentially mismatched satellite-derived LSP and ground-measured LSP should be mentioned.

Conclusions
In this study, we demonstrat that changing SZAs not only influence the NDVI/EVI time series but also exerts a cascading effect on the retrieved LSP. BRDF parameters provided by the MCD43A1 (BRDF) product were used as inputs in the RTLSR BRDF model to simulate reflectance at a constant SZA configuration. The MCD43A4 (NBAR) and MOD13A1 (VI) products were used as baselines. Our results illustrate that both NDVI and EVI time series are sensitive to varying SZAs but in an opposite manner: NDVI increases with increasing SZAs, while EVI decreases. We also compared the retrieved LSP to ground measurements, finding limited strength of remote sensing data in retrieving EOS compared to retrieving SOS. Although NDVI and EVI have equivalent accuracy in deriving SOS, EVI demonstrates advantages in predicting EOS. The discrepancies between these two vegetation indices in terms of deriving LSP are also reflected in whether the retrieved phenology is advanced or delayed. Specifically, SOSs derived from NDVI time series are advanced, yet from EVI time series are delayed. However, EOSs derived from NDVI time series are delayed, and from simulated EVI time series under consistent illumination geometry are also delayed, but EOSs derived from baseline products (i.e., MCD43A4 (NBAR) and MOD13A1 (VI)) are generally advanced. Additionally, the LSP derived from varying illumination geometries could lead to a difference spanning from a few days to a month in this case study, which has important implications for normalizing the illumination geometries in deriving LSP from remote sensing data. End-users of multi-angle remote sensing products, therefore, are suggested to take the SZA effect into consideration.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The remote sensing data that support the findings of this study are openly available on Google Earth Engine platform. The ground measurements are from the Harvard Forest Data Archive: https://harvardforest1.fas.harvard.edu/exist/apps/datasets/showData.html? id=HF003 (accessed on 1 October 2021). The dataset that supports the findings of this study is available from the corresponding author upon request.

Conflicts of Interest:
The authors declare no conflict of interest.