Sun-Angle Effects on Remote-Sensing Phenology Observed and Modelled Using Himawari-8

Satellite remote sensing of vegetation at regional to global scales is undertaken at considerable variations in solar zenith angle (SZA) across space and time, yet the extent to which these SZA variations matter for the retrieval of phenology remains largely unknown. Here we examined the effect of seasonal and spatial variations in SZA on retrieving vegetation phenology from time series of the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) across a study area in southeastern Australia encompassing forest, woodland, and grassland sites. The vegetation indices (VI) data span two years and are from the Advanced Himawari Imager (AHI), which is onboard the Japanese Himawari-8 geostationary satellite. The semi-empirical RossThick-LiSparse-Reciprocal (RTLSR) bidirectional reflectance distribution function (BRDF) model was inverted for each spectral band on a daily basis using 10-minute reflectances acquired by H-8 AHI at different sun-view geometries for each site. The inverted RTLSR model was then used to forward calculate surface reflectance at three constant SZAs (20◦, 40◦, 60◦) and one seasonally varying SZA (local solar noon), all normalised to nadir view. Time series of NDVI and EVI adjusted to different SZAs at nadir view were then computed, from which phenological metrics such as start and end of growing season were retrieved. Results showed that NDVI sensitivity to SZA was on average nearly five times greater than EVI sensitivity. VI sensitivity to SZA also varied among sites (biome types) and phenological stages, with NDVI sensitivity being higher during the minimum greenness period than during the peak greenness period. Seasonal SZA variations altered the temporal profiles of both NDVI and EVI, with more pronounced differences in magnitude among NDVI time series normalised to different SZAs. When using VI time series that allowed SZA to vary at local solar noon, the uncertainties in estimating start, peak, end, and length of growing season introduced by local solar noon varying SZA VI time series, were 7.5, 3.7, 6.5, and 11.3 days for NDVI, and 10.4, 11.9, 6.5, and 8.4 days for EVI respectively, compared to VI time series normalised to a constant SZA. Furthermore, the stronger SZA dependency of NDVI compared with EVI, resulted in up to two times higher uncertainty in estimating annual integrated VI, a commonly used remote-sensing proxy for vegetation productivity. Since commonly used satellite products are not generally normalised to a constant sun-angle across space and time, future studies to assess the sun-angle effects on satellite applications in agriculture, ecology, environment, and carbon science are urgently needed. Measurements taken by new-generation geostationary (GEO) satellites offer an important opportunity to refine this assessment at finer temporal scales. In addition, studies are needed to evaluate the suitability of different BRDF models for normalising sun-angle across a broad spectrum of vegetation Remote Sens. 2020, 12, 1339; doi:10.3390/rs12081339 www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 1339 2 of 23 structure, phenological stages and geographic locations. Only through continuous investigations on how sun-angle variations affect spatiotemporal vegetation dynamics and what is the best strategy to deal with it, can we achieve a more quantitative remote sensing of true signals of vegetation change across the entire globe and through time.

Within the SE Australia study area, three well-characterised local sites were selected, representing different biome types and rainfall regimes: Tumbarumba (Eucalyptus forest, mean annual precipitation, MAP = 1924 mm), Cumberland Plains (Eucalyptus woodland, MAP = 806 mm), and Yanco (pasture, MAP = 472 mm) ( Figure 1 and Table 1). The spatial extent of each site was 3 km × 3 km corresponding to the 3 × 3 pixels window to extract H-8 AHI time series for each site. These sites were used to investigate the sun-angle influence on vegetation indices and phenology at the site level.

Himawari-8 Advanced Himawari Imager Data
Himawari-8 is a Japanese geostationary satellite launched on 7 October 2014 and is positioned above 140.7°E and 0.02°S [47]. The AHI on board the H-8 satellite scans the Asia-Pacific region every 10 minutes at spatial resolution 500 m for the Red band (channel 3, 640 nm) and 1000 m for the Blue (channel 1, 470 nm) and near infrared (NIR, channel 4, 860 nm) bands, respectively. The data from Red band were averaged over 2 × 2 pixels to match the resolution of the other bands. Nearly two years of 10-minutes H-8 AHI surface reflectance data were processed from March 2016 to December 2017 by the Australian Bureau of Meteorology (BoM) using the multi-angle implementation of Within the SE Australia study area, three well-characterised local sites were selected, representing different biome types and rainfall regimes: Tumbarumba (Eucalyptus forest, mean annual precipitation, MAP = 1924 mm), Cumberland Plains (Eucalyptus woodland, MAP = 806 mm), and Yanco (pasture, MAP = 472 mm) ( Figure 1 and Table 1). The spatial extent of each site was 3 km × 3 km corresponding to the 3 × 3 pixels window to extract H-8 AHI time series for each site. These sites were used to investigate the sun-angle influence on vegetation indices and phenology at the site level.

Himawari-8 Advanced Himawari Imager Data
Himawari-8 is a Japanese geostationary satellite launched on 7 October 2014 and is positioned above 140.7 • E and 0.02 • S [47]. The AHI on board the H-8 satellite scans the Asia-Pacific region every 10 minutes at spatial resolution 500 m for the Red band (channel 3, 640 nm) and 1000 m for the Blue (channel 1, 470 nm) and near infrared (NIR, channel 4, 860 nm) bands, respectively. The data from Red band were averaged over 2 × 2 pixels to match the resolution of the other bands. Nearly two years of Remote Sens. 2020, 12, 1339 5 of 23 10-minutes H-8 AHI surface reflectance data were processed from March 2016 to December 2017 by the Australian Bureau of Meteorology (BoM) using the multi-angle implementation of atmospheric correction (MAIAC) algorithm [48]. The MAIAC was originally developed for MODIS to retrieve aerosol and perform atmospheric correction over the land surface [49] and has shown to be suitable for the processing of H-8 AHI data as well [42]. It should be noted that the MAIAC code we used in this study for the H-8 AHI data is still under development in order to optimise its performance, although we do not expect that future code adjustment should have significant impacts on the results of this study. For the selected local sites, 10-minute measurements within a 3×3 pixels window (9 km 2 ) of H-8 AHI data were extracted across the two years of study period. As the three study sites were chosen with the consideration of homogeneity in vegetation cover types, the 3 km × 3 km window should sufficiently account for the geolocation error of H-8 AHI, which was reported to be 150 m (north-south) × 500 m (east-west) [50]. Figure 2 shows a workflow diagram indicating major steps involved in H-8 data preprocessing, BRDF modelling (Sec. 2.3), and phenology retrieval (Sec. 2.5).
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 23 atmospheric correction (MAIAC) algorithm [48]. The MAIAC was originally developed for MODIS to retrieve aerosol and perform atmospheric correction over the land surface [49] and has shown to be suitable for the processing of H-8 AHI data as well [42]. It should be noted that the MAIAC code we used in this study for the H-8 AHI data is still under development in order to optimise its performance, although we do not expect that future code adjustment should have significant impacts on the results of this study. For the selected local sites, 10-minute measurements within a 3×3 pixels window (9 km 2 ) of H-8 AHI data were extracted across the two years of study period. As the three study sites were chosen with the consideration of homogeneity in vegetation cover types, the 3 km × 3 km window should sufficiently account for the geolocation error of H-8 AHI, which was reported to be 150 m (north-south) × 500 m (east-west) [50]. Figure 2 shows a workflow diagram indicating major steps involved in H-8 data preprocessing, BRDF modelling (Sec. 2.3), and phenology retrieval (Sec. 2.5).

Figure 2.
Workflow diagram of Japanese Himawari-8 (H-8) Advanced Himawari Imager (AHI) data preprocessing, bidirectional reflectance distribution function (BRDF) modelling, and phenological metrics retrieval. Figure 3 shows an example H-8 AHI sun-view geometry at one of the three local sites. H-8 AHI observes the same target always from a fixed angle, resulting in fixed view zenith angle and view azimuth angles for each pixel. For the three sites we selected, the view zenith angles (VZAs) are 42.28° for Tumbarumba, 40.72° for Cumberland Plains, and 41.13° for Yanco. In contrast to the fixed view angles, the solar zenith angle (SZA) and solar azimuth angle (SAA) of each individual H-8 observation can vary diurnally from sunrise to sunset ( Figure 3B and Figure 3C).

Figure 2.
Workflow diagram of Japanese Himawari-8 (H-8) Advanced Himawari Imager (AHI) data preprocessing, bidirectional reflectance distribution function (BRDF) modelling, and phenological metrics retrieval. Figure 3 shows an example H-8 AHI sun-view geometry at one of the three local sites. H-8 AHI observes the same target always from a fixed angle, resulting in fixed view zenith angle and view azimuth angles for each pixel. For the three sites we selected, the view zenith angles (VZAs) are 42.28 • for Tumbarumba, 40.72 • for Cumberland Plains, and 41.13 • for Yanco. In contrast to the fixed view angles, the solar zenith angle (SZA) and solar azimuth angle (SAA) of each individual H-8 observation can vary diurnally from sunrise to sunset ( Figure 3B,C).

Bidirectional Reflectance Distribution Function (BRDF) Modelling
The surface reflectance can be described using a bidirectional reflectance distribution function (BRDF), which is a function of solar zenith angle, view zenith angle, and both solar azimuth angle and view azimuth angle, with respect to a reference direction [51]. BRDF can be used to standardise surface reflectance taken from varying sun-view geometries to a common geometry to facilitate quantitative comparisons [52]. BRDF is determined by land surface structure and optical properties such as shadow-casting, soil condition, mutual view shadowing, and the spatial distribution of vegetation elements [53].

Bidirectional Reflectance Distribution Function (BRDF) Modelling
The surface reflectance can be described using a bidirectional reflectance distribution function (BRDF), which is a function of solar zenith angle, view zenith angle, and both solar azimuth angle and view azimuth angle, with respect to a reference direction [51]. BRDF can be used to standardise surface reflectance taken from varying sun-view geometries to a common geometry to facilitate quantitative comparisons [52]. BRDF is determined by land surface structure and optical properties such as shadow-casting, soil condition, mutual view shadowing, and the spatial distribution of vegetation elements [53].
The inset of Figure 2 shows the major steps taken for running the RTLSR model first in "inversion-mode" and subsequently in "forward mode". The inversion mode takes the 10-min H-8 surface reflectance and observational Sun-view geometry for any given pixel to invert the three model parameters, fiso, fvol, and fgeo, on a daily basis using multiple linear regression. The BRDF of Equation 1 was assumed to be the same as the measured reflectance from the H-8 AHI sensor. We imposed a lower limit for the model parameters so that all inversions that resulted in negative model parameters (fiso, fvol, and fgeo) were excluded from the analysis. We did not impose an upper limit on the three parameters but, instead, assessed the goodness of fit for model inversion to ensure that the observations can be predicted using the fitted model parameters within a predefined error threshold. In addition, the sufficient amount of input observations from H-8 AHI to fit the model allowed a high degree of freedom during the multiple linear regression and the three free parameters, therefore, can be inverted successfully on a daily basis. Before running the inversion, quality flags of H-8 AHI data were used to preserve only good quality observations as the model input. Observations taken when In this study, the semi-empirical RossThick-LiSparse Reciprocal (RTLSR) BRDF model [53,54] was used to normalise H-8 AHI surface reflectance to nadir view and different solar zenith angles. The RTLSR model has been used as the standard model for the MCD43 BRDF/NBAR/Albedo product [55]. The utility of the RTLSR model for normalising surface reflectance anisotropy of H-8 geostationary data has also been tested [42,56]. The RTLSR model takes the form given by [57], where R(θ s , θ v , φ s , φ s , Λ) is the BRDF in waveband Λ; θ s is solar zenith angle (SZA); θ v is VZA; φ s is solar azimuth angle (SAA); φ v is view azimuth angle (VAA); Λ is waveband of width ∆λ; f iso (Λ) is isotropic parameter of BRDF at waveband Λ; K vol is the RossThick volume scattering kernel; K geo is the Li-Sparse-Reciprocal geometric scattering kernel. Full mathematical expressions of K vol and K geo can be found in [53,57,58].
The inset of Figure 2 shows the major steps taken for running the RTLSR model first in "inversion-mode" and subsequently in "forward mode". The inversion mode takes the 10-min H-8 surface reflectance and observational Sun-view geometry for any given pixel to invert the three model parameters, f iso , f vol , and f geo , on a daily basis using multiple linear regression. The BRDF of Equation (1) was assumed to be the same as the measured reflectance from the H-8 AHI sensor. We imposed a lower limit for the model parameters so that all inversions that resulted in negative model parameters (f iso , f vol , and f geo ) were excluded from the analysis. We did not impose an upper limit on the three parameters but, instead, assessed the goodness of fit for model inversion to ensure that the observations can be predicted using the fitted model parameters within a predefined error threshold. In addition, the sufficient amount of input observations from H-8 AHI to fit the model allowed a high degree of freedom during the multiple linear regression and the three free parameters, therefore, can be inverted successfully on a daily basis. Before running the inversion, quality flags of H-8 AHI data were used to preserve only good quality observations as the model input. Observations taken when SZA is greater than 70 • (very low sun), which usually occurred in early morning or late afternoon, were also abandoned as larger SZA means a longer path of sunlight travelling through the atmosphere and causes the atmospheric correction model to be less effective.
The goodness of fit for model inversion was evaluated by the coefficient of determination adjusted for the degree of freedom (adjusted R 2 ), while the accuracy of the model inversion was assessed by the root-mean-square error (RMSE) associated with inversion. Low R 2 and high RMSE signal low model performance and hence low reliability in using the model parameters to adjust reflectance values to the targeting sun-view geometry. Only inversions that resulted in RMSEs lower than a pre-defined threshold value were used later in the "forward-mode". For the three local sites, RMSE for the Red band ranges from 0 to 0.11 (mean = 0.006 ± 0.007), from 0 to 0.038 (mean = 0.004 ± 0.003) for the Blue band, and from 0 to 0.09 (mean = 0.017 ± 0.015) for NIR band. Given the fact that different spectral bands have much different ranges of reflectance values and RMSEs, we used a relative threshold for each band based on normalised RMSE (nRMSE), computed as RMSE normalised by mean reflectance of each band, with nRMSE lower than 20% considered as acceptable quality inversion. Preliminary analysis showed that after excluding the inversions with nRMSE greater than 20%, the mean nRMSEs across three local sites and across two years were 6.7 ± 5.0%, 8.1 ± 5.4%, and 4.9 ± 4.5%, for Red, Blue, and NIR bands respectively.
Once the three RTLSR model parameters were reliably inverted for each day, the model was then operated in the "forward mode" to calculate daily reflectance at the desired Sun-view geometry ( Figure 3). Since the focus of this study is to assess the effect of sun-angle variations on vegetation indices and phenology retrievals, VZA was, therefore, normalised to zero (i.e., nadir view) and SZA was normalised to different scenarios. Four SZA scenarios were used, including three constant SZAs: 20 • , 40 • , 60 • , which means that SZA is constant across time for any given site, as well as constant across sites; one seasonally varying SZA: local solar noon (LSN), which means that SZA is set as the value as acquired at local solar noon of each day for each site, and hence it varies across time for any given site and is also different across latitude for any given time ( Figure 3C). The LSN scenario is equivalent to the MODIS MCD43A4 NBAR-Nadir BRDF adjusted reflectance product [59], which also fixes VZA to nadir but allows SZA to vary across time at local solar noon. The RTLSR model inversion and forward calculation were carried out for Red, Blue, and NIR reflectances, from which NDVI and EVI at different SZA configurations can be computed:

Vegetation Indices
Phenology retrieval from satellite remote sensing commonly uses time series of vegetation indices (VIs) [60][61][62][63][64]. VIs are remote-sensing proxies of canopy "greenness", integrating canopy properties such as canopy structure, green leaf area, and canopy leaf chlorophyll content [65,66]. The two most widely used VIs for phenology retrievals are the NDVI and EVI. The EVI was proposed as an optimised version of NDVI that effectively reduces influences from varying soil background reflectance and atmospheric conditions [13]. The equations defining NDVI [67] and EVI are: where ρ nir , ρ red and ρ blue are reflectances of the near infrared, red, and blue bands, corresponding to channel 4 (860 nm), channel 3 (640 nm), and channel 1 (470 nm) of the H-8 AHI sensor, respectively [46].

Phenology Metrics Retrieval Method
The SSA-Pheno (singular spectrum analysis for phenology) algorithm, graphically depicted by Figure 4, was used to retrieve the phenological metrics from VI time series [68,69]. Six phenological metrics, including the start of growing season (SGS), peak of growing season (PGS), end of growing season (EGS), length of growing season (LGS), seasonal maximum VI (VI max ) and annual integrated VI (IntVI), were retrieved from the time series of H-8 AHI NDVI and EVI time series. Phenological metrics retrieved from time series of VIs normalised to different SZAs will be indicated using the subscripts, e.g., SGS S-40 and IntNDVI S-20 means the SGS and IntNDVI retrieved from time series of NDVI S-40 .
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 23 metrics, including the start of growing season (SGS), peak of growing season (PGS), end of growing season (EGS), length of growing season (LGS), seasonal maximum VI (VImax) and annual integrated VI (IntVI), were retrieved from the time series of H-8 AHI NDVI and EVI time series. Phenological metrics retrieved from time series of VIs normalised to different SZAs will be indicated using the subscripts, e.g., SGSS-40 and IntNDVIS-20 means the SGS and IntNDVI retrieved from time series of NDVIS-40. The SSA-Pheno algorithm has been tested over Australia across a wide-range of vegetation types and has been demonstrated to be a robust and reliable method for extracting phenological metrics from noisy VI time series [68,69]. The SSA (Singular Spectrum Analysis) is a data-adaptive method that has been found to be well-suited to the analysis of nonlinear dynamics in geophysical datasets [70][71][72]. SSA-Pheno detects SGS when VI reached the value equal to the pre-season minimum VI value LG S The SSA-Pheno algorithm has been tested over Australia across a wide-range of vegetation types and has been demonstrated to be a robust and reliable method for extracting phenological metrics from noisy VI time series [68,69]. The SSA (Singular Spectrum Analysis) is a data-adaptive method that has been found to be well-suited to the analysis of nonlinear dynamics in geophysical datasets [70][71][72]. SSA-Pheno detects SGS when VI reached the value equal to the pre-season minimum VI value (or baseline) plus 10% of the seasonal VI amplitude (seasonal maximum minus pre-season minimum). Similarly, EGS is detected when VI reaches the values equal to the minimum value after the growing season plus 10% of seasonal amplitude during the browndown phase ( Figure 4). The difference between EGS and SGS is the LGS and the timing when VI reaches its seasonal maximum is PGS. VImax Remote Sens. 2020, 12, 1339 9 of 23 is the seasonal maximum VI and IntVI is the annual integrated VI (with soil background VI subtracted), both have been used as a remote sensing proxy for vegetation productivity [69,73,74]. Soil background VI was empirically assigned as 0.15 for NDVI and 0.08 for EVI, respectively.

Land-Cover Map
We used the Dynamic Land Cover Dataset (DLCD) from Geoscience Australia and Bureau of Agricultural and Resources Economics and Sciences (https://ecat.ga.gov.au/geonetwork/srv/eng/catalog. search#/metadata/83868) [75]. The raw DLCD data is in 250 m spatial resolution, which was resampled to 1km to analyze with the comparative H-8 AHI data.

Statistics
The RMSE was calculated to assess the sun-angle variation caused differences in phenological metrics retrieved from time series of VIs normalised to different SZA scenarios, where x 0 is the phenological metrics for any given site from time series of VI S-LSN , and x is the phenological metrics for any given site from time series of VI S-40 , and the n is the number of samples.
To assess the uncertainty relative to sample mean, a normalized RMSE (nRMSE, %), defined as the ratio between RMSE and the sample mean, was also computed. The open source R programming language (version 3.6.2) was used for data processing, statistical analysis, and data visualisation [76], with R packages contributed by the user community (http: //cran.r-project.org).

Seasonal Profiles of Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI)
Normalised to Different Solar Zenith Angles (SZAs) Figure 5 shows the time series of NDVI and EVI normalised to different SZAs. Smaller SZA resulted in lower NDVI, with the magnitude of NDVI S-LSN was close to that of NDVI S-20 due to the selection of the lowest SZA within each day to generate the NDVI S-LSN ( Figure 5). In comparison, sun-angle had a much smaller effect on the profiles of EVI ( Figure 5). Greater SZA resulted in higher EVI over the minimum greenness period (the period when EVI declines from highest to the lowest within a growing season), but this effect was reversed in the peak greenness period (the period when EVI increases form the lowest to the highest within a growing season) ( Figure 5). Figure 6 presents the cross-site comparison of the variations of NDVI and EVI in relation to SZAs over two phenological stages, the peak greenness period and the minimum greenness period. Sensitivity of VI to change in SZA (i.e., change in VI per degree SZA change) was also computed. Since NDVI and EVI can have different sensitivities merely because they have different dynamic range (NDVI tends to have higher value than EVI), and different sites and phenological stages can also have different VI sensitivity to SZA simply because of the difference in the amount of green foliage. Therefore, to facilitate the comparisons of VI sensitivity to SZA variations between VIs, among sites and between phenological stages, normalised VIs were computed by dividing all VIs by their corresponding VI S-20 at any given site/date. NDVI increased with the increase of SZA, with the lowest NDVI found when the SZA is the smallest ( Figure 6A). By contrast, EVI did not show a single response to SZA variations between phenological stages, with negative response to the increase in SZA during the peak greenness period and positive response during the minimum greenness period ( Figure 6B). respectively. Data from this period will be used to assess the sensitivities of VIs to SZA variations in relation to phenological stages. Figure 6 presents the cross-site comparison of the variations of NDVI and EVI in relation to SZAs over two phenological stages, the peak greenness period and the minimum greenness period. Sensitivity of VI to change in SZA (i.e., change in VI per degree SZA change) was also computed. Since NDVI and EVI can have different sensitivities merely because they have different dynamic range (NDVI tends to have higher value than EVI), and different sites and phenological stages can also have different VI sensitivity to SZA simply because of the difference in the amount of green foliage. Therefore, to facilitate the comparisons of VI sensitivity to SZA variations between VIs, among sites and between phenological stages, normalised VIs were computed by dividing all VIs by their corresponding VIS-20 at any given site/date. NDVI increased with the increase of SZA, with the lowest NDVI found when the SZA is the smallest ( Figure 6A). By contrast, EVI did not show a single response to SZA variations between phenological stages, with negative response to the increase in SZA during the peak greenness period and positive response during the minimum greenness period ( Figure 6B). During both phenological stages and across all sites (biome types), NDVI exhibited much stronger dependency to SZA variations than EVI ( Figure 6C,D). On average, NDVI sensitivity to SZA was 4.6 (peak greenness period) and 5.1(minimum greenness period) times greater than EVI sensitivity ( Figure 6C,D; Table 2). For NDVI, there was a 75% increase in SZA sensitivity from the peak greenness to the minimum greenness period, with the highest (130%) increase observed over the Yanco grassland site ( Figure 6C,D). There was also an increasing trend in the difference between NDVI and EVI sensitivities to SZA from forest (NDVI sensitivity is 2.2 times of EVI sensitivity), woodland (3.6 times), to pasture (5.4 times) sites ( Figure 6C,D; Table 2).  Figure 7 presents the SSA-reconstructed VI time series normalised to different SZAs. The phenological transitional dates retrieved from each time series were also labelled on the plot. Figure 7 clearly shows that the sun-angle effect on individual observations was propagated across time due to the interaction of seasonal sun-angle variations and vegetation structural change, leading to alternations in both the shape and magnitude of VI temporal profiles (Figure 7). The effect not only caused difference between VI S-40 or VI S-60 (seasonally constant SZA cases) and VI S-LSN (seasonally varying SZA cases), but also between the two constant SZA cases VI S-40 and VI S-60 (Figure 7). Therefore, VIs with SZA normalised to local solar noon, which is the common input to many global and regional phenology products, do not necessarily align with phenology profiles generated from VIs that were normalised to constant SZA (Figure 7). Error bars on the plot indicate 95% confidence interval of the mean and the statistical significance of the slope is indicated as follows: "***" (p < 0.0001), "*" (p < 0.01), "ns" (non-significant, p > 0.05).\.

Sun-Angle Effect on Vegetation Phenology at Site Level
During both phenological stages and across all sites (biome types), NDVI exhibited much stronger dependency to SZA variations than EVI ( Figure 6C and Figure 6D). On average, NDVI sensitivity to SZA was 4.6 (peak greenness period) and 5.1(minimum greenness period) times greater than EVI sensitivity ( Figure 6C, 6D; Table 2). For NDVI, there was a 75% increase in SZA sensitivity from the peak greenness to the minimum greenness period, with the highest (130%) increase observed over the Yanco grassland site ( Figure 6C and Figure 6D). There was also an increasing trend in the difference between NDVI and EVI sensitivities to SZA from forest (NDVI sensitivity is 2.2 times of EVI sensitivity), woodland (3.6 times), to pasture (5.4 times) sites ( Figure 6C, 6D; Table 2). 60 • (by every 5 • ) during the peak greenness period and the minimum greenness period respectively, at three local sites; (C,D) Comparison of the sensitivity of NDVI and EVI to changes in SZA during different phenological stages. For comparison among sites, VI values have been normalised by VI SZA20 for each site. Sensitivity is defined as change in VI per degree change in SZA. Error bars on the plot indicate 95% confidence interval of the mean and the statistical significance of the slope is indicated as follows: "***" (p < 0.0001), "*" (p < 0.01), "ns" (non-significant, p > 0.05).\. Figure 8 presents the cross-site relationships between phenological transitional dates retrieved from VI S-40 time series (x-axis, used as a reference) and those retrieved from VI S-60 or VI S-LSN time series (y-axis). Overall, seasonal sun-angle variations introduced non-negligible uncertainties on phenological transitional dates retrieved from NDVI and EVI. For NDVI, the uncertainties (RMSE) of using time series of NDVI S-LSN for retrieving SGS, PGS, EGS, and LGS, as compared to the use of NDVI S-40 , were 7.5, 3.7, 6.5, and 11.3 days respectively ( Figure 8A-D). The equivalent uncertainties for Remote Sens. 2020, 12, 1339 12 of 23 EVI were 10.4, 11.9, 6.5, and 8.4 days ( Figure 8E-H). Difference in phenology retrievals between two constant SZA cases can be also observed ( Figure 8). Figure 7 presents the SSA-reconstructed VI time series normalised to different SZAs. The phenological transitional dates retrieved from each time series were also labelled on the plot. Figure  7 clearly shows that the sun-angle effect on individual observations was propagated across time due to the interaction of seasonal sun-angle variations and vegetation structural change, leading to alternations in both the shape and magnitude of VI temporal profiles (Figure 7). The effect not only caused difference between VIS-40 or VIS-60 (seasonally constant SZA cases) and VIS-LSN (seasonally varying SZA cases), but also between the two constant SZA cases VIS-40 and VIS-60 ( Figure 7). Therefore, VIs with SZA normalised to local solar noon, which is the common input to many global and regional phenology products, do not necessarily align with phenology profiles generated from VIs that were normalised to constant SZA (Figure 7).  for EVI were 10.4, 11.9, 6.5, and 8.4 days ( Figure 8E-H). Difference in phenology retrievals between two constant SZA cases can be also observed ( Figure 8).  Figure 9 shows the cross-site comparisons between VImax (seasonal maximum VI) and IntVI (annual integrated VI), retrieved from time series of VIS-40 (x-axis, used as a reference) and those retrieved from time series of VIS-60 or VIS-LSN (y-axis). As expected from a positive NDVI response to the increase in SZA, NDVImaxS-LSN and IntNDVIS-LSN were on average 7% and 12% lower than NDVImaxS- 40 and IntNDVIS-40, respectively ( Figure 9A). These uncertainties were smaller for EVI, with EVImaxS-LSN and IntEVIS-LSN were on average 4% and 2% higher, than EVImaxS-LSN and IntEVIS-LSN respectively  Figure 9 shows the cross-site comparisons between VI max (seasonal maximum VI) and IntVI (annual integrated VI), retrieved from time series of VI S-40 (x-axis, used as a reference) and those retrieved from time series of VI S-60 or VI S-LSN (y-axis). As expected from a positive NDVI response to the increase in SZA, NDVI maxS-LSN and IntNDVI S-LSN were on average 7% and 12% lower than NDVI maxS-40 and IntNDVI S-40 , respectively ( Figure 9A). These uncertainties were smaller for EVI, with EVI maxS-LSN and IntEVI S-LSN were on average 4% and 2% higher, than EVI maxS-LSN and IntEVI S-LSN respectively ( Figure 9C). As a result of the differential sensitivity of VIs to SZA among sites (biome types), the betweensites relative differences in VImax or IntVI were altered by the sun-angle effect. For instance, IntNDVIS- 40 for Tumbarumba forest site was on average 151% higher than that of the Yanco pasture site. However, such differences increased to 174% if using IntNDVIS-LSN (Table 3). This change could be ignored for EVI, with IntEVIS-40 and IntEVIS-LSN for Tumbarumba forest site were on average 77% and 77% higher than those of Yanco pasture site respectively. Taking the Tumbarumba and Yanco site pair as an example, the estimated distortions in proportion of vegetation productivity between any two sites (or any two pixels), as caused by sun-angle effect alone, was 15% and 0% for IntNDVI and IntEVI, respectively (Table 3).  As a result of the differential sensitivity of VIs to SZA among sites (biome types), the between-sites relative differences in VI max or IntVI were altered by the sun-angle effect. For instance, IntNDVI S-40 for Tumbarumba forest site was on average 151% higher than that of the Yanco pasture site. However, such differences increased to 174% if using IntNDVI S-LSN (Table 3). This change could be ignored for EVI, with IntEVI S-40 and IntEVI S-LSN for Tumbarumba forest site were on average 77% and 77% higher than those of Yanco pasture site respectively. Taking the Tumbarumba and Yanco site pair as an example, the estimated distortions in proportion of vegetation productivity between any two sites (or any two pixels), as caused by sun-angle effect alone, was 15% and 0% for IntNDVI and IntEVI, respectively (Table 3).

Discussion
In this study we investigated the influence of seasonal sun-angle variations on the temporal VI profiles and phenology retrievals by coupling measurements taken by H-8 geostationary satellite and BRDF modelling. Our results showed that NDVI was more sensitive to SZA than EVI, and such sensitivity varies between biome types and phenological stages. The sun-angle effect that propagated into the VI time series not only introduced large uncertainties in retrieving phenological transitional dates, but also led to errors in estimating vegetation productivity. These results call for an urgent need to take into account sun-angle effects on monitoring vegetation dynamics and phenology using satellite observations. Our results also demonstrated that the sensors onboard the new generation GEO satellites, with refined radiometric and spatial resolutions and augmented by a much denser temporal (angular) sampling than LEO satellites, offer important opportunities to enhance global vegetation monitoring applications.

Sun-Angle Dependency of Vegetation Indices (VIs)
Differential sensitivities to sun-angle variations between VIs, over phenological stages, and across biome types were observed, which highlighted the complexity of the sun-angle effect in the remote sensing of vegetation. Our findings of a strong sun-angle NDVI dependency was largely consistent to what was reported in previous studies [16,28,77]. The large SZA sensitivity of NDVI can be attributed to the fact that NDVI is functionally more related to red band reflectance. Due to its high canopy absorption, the red reflectance has a low canopy transmittance, resulting in a more prominent shadow effect, rendering the red reflectance being more sensitive to the changes in Sun (illumination) angle than NIR reflectance which has a high canopy transmittance.
Overall, NDVI was found to be several times more sensitive to sun-angle variations than EVI. This result is consistent with previous findings in Australian tropical forests [25] and savannas [28], but different from what was reported over Amazon tropical forests [26]. Galvão et al. found that NDVI was not as responsive to sun-angle variations as EVI over Amazon forests [26]. These apparently conflicting conclusions might be reconciled with the following reasoning. First, it is well-known that NDVI tends to respond to changes in vegetation biomass in a nonlinear manner and saturate over high biomass areas [13]. Previous studies have found that there was a decreasing trend in NDVI sensitivity to sun-angle with the increase in leaf area index [78] or site-average greenness [28], implies that the muted response of NDVI to sun-angle variations as was observed by [26] was likely due to the fact that NDVI was simply saturated over these high biomass areas. Though more careful investigations are needed to confirm our speculation. Second, the behaviour of NDVI with sun-angle variations was found to be dependent on the reflectance properties of underlying soil [17], and such soil-induced effect would be minimal over dense canopies, additionally explained the reduced sensitivity of NDVI to sun-angle variations over dense Amazonian tropical forests. Nonetheless, these studies highlighted the need to assess the sun-angle effect on VIs more systematically across a broad spectrum of vegetation structural classes and geographic regions to draw a general conclusion and more importantly, to reveal the many underlying factors that dictate the sun-angle dependency of VIs and their relative importance at different occasions.

Sun-Angle Effect on Retrievals of Vegetation Phenology and Productivity
It is well-known that sun-angle variations can alter reflectances and VIs, however, it is much less well known the extent to which such effects would be propagated to the retrievals of vegetation phenology parameters, primary productivity, among other biophysical and biochemical properties of vegetation. This is likely because the commonly used LEO satellites can only take 1-2 measurements within a day (or even less frequent if spatial resolution is higher), limiting their ability to sample a wide range of sun-angles. Common practice for studying or correcting sun-angle effect on VIs from LEO involved the use of BRDF models and measurements taken over 8-to 16-day periods, assuming that surface anisotropy would remain unchanged within that interval of time [11,28,[78][79][80]. An apparent shortcoming of this approach is that an intrinsic uncertainty is imposed by the model assumptions, especially for time-sensitive metrics such as phenological transitional dates. Such shortcomings can be circumvented by the use of measurements from GEO satellites, from which a BRDF model can be inverted on a daily basis and the assumption of a constant surface anisotropy within several days is no longer needed. This enabled the assessment of sun-angle effect on phenology retrieval at finer temporal scales, a major benefit of using GEO.
The sun-angle effect alone caused more than one week of uncertainty in retrieving most phenological metrics. Our results using H-8 are largely consistent to a recent study using MODIS over tropical savannas [28] and add more biome types. Furthermore, we showed that the sun-angle effect can introduce 7% and 10% uncertainty in estimating NDVI max and IntNDVI, and a reduced uncertainty for EVI. This result has important implications as both the VI max and IntVI have been used extensively for estimating vegetation primary productivity (or crop yield and rangeland forage production) [81][82][83][84][85], for land carbon uptake modelling [86,87], and for food security assessment for famine early-warning systems [88,89]. Our findings therefore stress the need to consider proper corrections of the sun-angle effect to achieve more reliable use of vegetation indices in a variety of applications.

Limitations and Future Perspectives
The results from this study confirmed that GEO satellites with their advanced temporal and angular sampling offer unprecedented opportunities to better quantify the sun-angle effect on the remote sensing of vegetation dynamics. Nonetheless, several limitations of the existing approach can still be identified. These limitations call for innovations in satellite data collection strategies and sun-angle correction algorithms from the remote-sensing community on one hand, and on the other hand call for user community attention to more carefully examine the relevance of the sun-angle effect on their specific applications.
The first limitation of using GEO satellites is the under-sampling of view angles. By contrast with LEO satellites such as EOS-Terra/Aqua which are limited in both view/sun-angle sampling at any given day, GEO satellites can scan the entire disk more frequently within a day and hence offer much improved sub-daily sampling of sun-angle. However, this sampling comes with the tradeoff of having only a single view angle for any given location. A few potential technical pathways to overcome this limitation can be considered. First, with the assumption of BRDF reciprocity, the RTLSR model can be inverted by interchanging Sun and view geometry [33], allowing varying Sun angles from GEO to be used to fill view angle gaps. Another potential pathway, although it has yet to be tested, is the combined use of measurements taken by two or more GEO platforms located at different longitudes which gives at least two view angles for any given location, or even combining the measurements from GEO and LEO platforms. Of course, the second pathway demands many efforts in inter-calibrating measurements taken by multiple sensors [90], which is currently a rapidly developing field [91][92][93].
The second limitation of the existing approach for normalising sun-angle, not limited to this study, is the reliance on a BRDF model and hence the results obtained have an intrinsic dependency on the model's assumptions. Indeed, it is possible to select measurements taken by GEO to form a time series of reflectance or VI of constant SZA. However, uncertainty would still be introduced into the time series due to seasonally varying solar azimuth angle, in addition to the spatial variations in view-geometry. Only for regions that are located directly underneath the GEO platform (e.g., for H-8, the Papua New Guinea), where time series of VIs with a constant SZA at nadir view angle can be empirically formed without the reliance on a model. Theoretically speaking, unless there are sufficient near-nadir measurements for any given location/date, possibly from multiple narrow FOV (field of view) sensors onboard satellites such as Landsat, from which both solar zenith and azimuth angles can be fixed across time, otherwise it would be very difficult, if not impossible, to empirically correct the sun-angle effect. Even with this kind of data, higher latitude would still present a great challenge due to the very low maximum daily sun-angle in wintertime, forcing the selection of observations made at large SZA deemed to be less reliable.
Two technical directions can be taken to address the second limitation. First, from an observation point of view, studies can take advantage of increasingly available GEO data to assess the relative sensitivity of reflectances and VIs to SZA and RAA variations in relation to biome types and phenological stages. Conclusions drawn from these studies can be used to make decisions whether constant SZA or constant RAA is to be targeted (as the two cannot be fixed in the same time). A conceptually similar approach was proposed recently to match GEO (H-8 AHI) and LEO (Terra MODIS) based on the criteria of equal SZA or equal RAA [93]. Second, from a modelling point of view, the inversion of more complicated 3-D radiative-transfer models may become feasible with the combined use of high-resolution spaceborne light-detection and ranging (LiDAR, e.g., those from NASA's Global Ecosystem Dynamics Investigation, or GEDI mission), multi-spectral / multi-angle measurements, from which a more reliable surface reflectance correction for BRDF effect can be achieved.
The third limitation is the direct ground validation of the BRDF model simulated reflectance. Indeed, it is a great challenge to obtain ground measurements at different sun-view geometries at the field size comparable to the size of satellite pixels (e.g., 1km for H-8 AHI). Besides, normalising SZA to a constant value, e.g., 45 • , implies that for mid-/high-latitude there would essentially be no available real observations at such a SZA during the winter time (as the sun never reaches those SZAs in the winter for mid-/high-latitude). Ground-truth measurements would be relatively easier to make by comparing to high-resolution satellite measurements, but the dilemma here is that a high-resolution satellite means narrower FOV and hence even less frequent revisits and measurements for BRDF model inversion.
Given these dilemmas, alternative and indirect validation strategies may be taken to get around the validation challenge. The first strategy is to take multi-angle measurements using unmanned aerial vehicles (UAVs) over the relatively homogeneous area to obtain vegetation BRDF [94]. These measurements can be used to study the sun-view geometry effect on VIs and their temporal dynamics at small scale directly, but can also be used to compare to VIs normalised to defined sun-view geometry based on satellite measurements. The second strategy is to use the measurements from rapid developing global eddy-covariance flux tower networks, from which canopy structure (e.g., photosynthetic capacity) and function (e.g., Gross Primary Productivity, GPP) can be derived at the footprint size comparable to the size of moderate spaceborne sensors (e.g., 100 m-1000 m) [95,96]. Satellite and flux tower measurements have been frequently used for monitoring ecosystem dynamics [96][97][98][99][100][101][102][103][104]. The attempts of verifying the temporal vegetation dynamics as observed from satellites have also been made [32,67,100,105]. Additionally, measurements taken by more cost-effective time-lapse multispectral phenology-cameras (or phenocam), that now are widely installed on flux towers and other ecological monitoring sites can also serve as ground verification data [106][107][108][109][110].
Indeed, all the aforementioned challenges in normalising sun-angle variations with satellite observations originated by the inability of satellites to maneuver to a variety of positions relative to the sun. Much efforts are urgently needed to refine our understanding of the interactions between sun-angle variations, vegetation phenology, canopy structure, among many other external factors. Aside from the methodological development and the increasingly available large amount of data from space, it is perhaps equally important to reflect back on an early question asked by Middleton in 1992: "which solar zenith angle was best for acquiring surface reflectance measurements in order to estimate canopy variables?" [18].

Conclusions
Here we evaluated the sun-angle effect on vegetation indices and phenology retrievals from time series of VIs using H-8 AHI geostationary satellite measurements. Our results revealed that spatiotemporal sun-angle variations can pose many challenges to resolving vegetation dynamics by introducing considerable uncertainties in phenological metrics retrievals. Future studies are urgently needed to assess the relevance of the sun-angle effect on retrievals of biophysical and ecosystem function variables from remote-sensing measurements. Only through continuous investigations on how sun-angle variations affect spatiotemporal vegetation dynamics and what is the best strategy to deal with it, can we achieve a more quantitative remote sensing of true signals of vegetation change across the entire globe and through time.