Decomposing reflectance spectra to track gross primary production in a subalpine evergreen forest

Photosynthesis by terrestrial plants represents the majority of CO2 uptake on Earth, yet it is difficult to measure directly from space. Estimation of gross primary production (GPP) from remote sensing indices represents a primary source of uncertainty, in particular for observing seasonal variations in evergreen forests. Recent vegetation remote sensing techniques have highlighted spectral regions sensitive to dynamic changes in leaf/needle carotenoid composition, showing promise for tracking seasonal changes in photosynthesis of evergreen forests. However, these have mostly been investigated with intermittent field campaigns or with narrow-band spectrometers in these ecosystems. To investigate this potential, we continuously measured vegetation reflectance (400–900 nm) using a canopy spectrometer system, PhotoSpec, mounted on top of an eddy-covariance flux tower in a subalpine evergreen forest at Niwot Ridge, Colorado, USA. We analyzed driving spectral components in the measured canopy reflectance using both statistical and processbased approaches. The decomposed spectral components covaried with carotenoid content and GPP, supporting the interpretation of the photochemical reflectance index (PRI) and the chlorophyll/carotenoid index (CCI). Although the entire 400–900 nm range showed additional spectral changes near the red edge, it did not provide significant improvements in GPP predictions. We found little seasonal variation in both normalized difference vegetation index (NDVI) and the nearinfrared vegetation index (NIRv) in this ecosystem. In addition, we quantitatively determined needle-scale chlorophyllto-carotenoid ratios as well as anthocyanin contents using full-spectrum inversions, both of which were tightly correlated with seasonal GPP changes. Reconstructing GPP from vegetation reflectance using partial least-squares regression (PLSR) explained approximately 87 % of the variability in observed GPP. Our results linked the seasonal variation in reflectance to the pool size of photoprotective pigments, highPublished by Copernicus Publications on behalf of the European Geosciences Union. 4524 R. Cheng et al.: Decomposing reflectance spectra to track gross primary production lighting all spectral locations within 400–900 nm associated with GPP seasonality in evergreen forests.

Abstract. Photosynthesis by terrestrial plants represents the majority of CO 2 uptake on Earth, yet it is difficult to measure directly from space. Estimation of gross primary production (GPP) from remote sensing indices represents a primary source of uncertainty, in particular for observing seasonal variations in evergreen forests. Recent vegetation remote sensing techniques have highlighted spectral regions sensitive to dynamic changes in leaf/needle carotenoid composition, showing promise for tracking seasonal changes in photosynthesis of evergreen forests. However, these have mostly been investigated with intermittent field campaigns or with narrow-band spectrometers in these ecosystems. To investigate this potential, we continuously measured vegetation reflectance (400-900 nm) using a canopy spectrometer system, PhotoSpec, mounted on top of an eddy-covariance flux tower in a subalpine evergreen forest at Niwot Ridge, Colorado, USA. We analyzed driving spectral components in the measured canopy reflectance using both statistical and process-based approaches. The decomposed spectral components covaried with carotenoid content and GPP, supporting the interpretation of the photochemical reflectance index (PRI) and the chlorophyll/carotenoid index (CCI). Although the entire 400-900 nm range showed additional spectral changes near the red edge, it did not provide significant improvements in GPP predictions. We found little seasonal variation in both normalized difference vegetation index (NDVI) and the nearinfrared vegetation index (NIRv) in this ecosystem. In addition, we quantitatively determined needle-scale chlorophyllto-carotenoid ratios as well as anthocyanin contents using full-spectrum inversions, both of which were tightly correlated with seasonal GPP changes. Reconstructing GPP from vegetation reflectance using partial least-squares regression (PLSR) explained approximately 87 % of the variability in observed GPP. Our results linked the seasonal variation in reflectance to the pool size of photoprotective pigments, high-

Introduction
Terrestrial gross primary production (GPP), the gross CO 2 uptake through photosynthesis, is the largest uptake of atmospheric CO 2 (Ciais et al., 2013), yet the uncertainties are large, hampering our ability to monitor and predict the response of the terrestrial biosphere to climate change (Ahlström et al., 2012). Hence, accurately mapping GPP globally is critical. In contrast to unevenly distributed ground-level measurements such as Fluxnet (Baldocchi et al., 2001), satellites can infer GPP globally and uniformly. Remote sensing techniques are based on the optical response of vegetation to incoming sunlight, which can track photosynthesis via the absorption features of photosynthetic and photoprotective pigments (Rouse et al., 1974;Liu and Huete, 1995;Gamon et al., 1992Gamon et al., , 2016. Progress is particularly important for evergreen forests, which can have large seasonal dynamics in photosynthesis but low variability in canopy structure and color. However, these promising techniques still lack a comprehensive evaluation/validation using both continuous in situ measurements and process-based simulations. GPP can be expressed as a function of photosynthetically active radiation (PAR), the fraction of PAR absorbed by the canopy (fPAR), and light-use efficiency (LUE): with LUE representing the efficiency of plants to fix carbon using absorbed light (Monteith, 1972;Monteith and Moss, 1977). The accuracy of remote-sensing-derived GPP is limited by the estimation of LUE, which is more dynamic and difficult to measure remotely than PAR and fPAR, particularly in evergreen ecosystems. There have been many studies inferring the light absorbed by canopies (i.e., fPAR) from vegetation indices (VIs) that estimate the "greenness" of canopies (Running et al., 2004;Zhao et al., 2005;Robinson et al., 2018;Glenn et al., 2008), such as the normalized difference vegetation index (NDVI; Rouse et al., 1974;Tucker, 1979), the enhanced vegetation index (EVI; Liu and Huete, 1995;Huete et al., 1997), and the near-infrared vegetation index (NIRv; Badgley et al., 2017). Current GPP data products derived from Eq. (1) rely on the modulation of abiotic conditions to estimate LUE (Xiao et al., 2004). LUE is derived empirically by defining a general timing of dormancy for all evergreen forests with the same plant functional type (e.g., Krinner et al., 2005) or the same meteorological thresholds (e.g., Running et al., 2004). However, within the same climate region or plant functional type, forests are not identical -leading to uncertainties in estimated LUE (Stylinski et al., 2002;Gamon et al., 2016;Zuromski et al., 2018), which propagate to the estimation of GPP. Because evergreen trees retain most of their needles and chlorophyll throughout the entire year , LUE in evergreens is regulated by needle biochemistry. As LUE falls with the onset of winter due to unfavorable environmental conditions and seasonal downregulation of photosynthetic capacity, evergreen needles quench excess absorbed light via thermal energy dissipation that involves the xanthophyll cycle and other pigments (Adams and Demmig-Adams, 1994;Demmig-Adams and Adams, 1996;Verhoeven et al., 1996;Zarter et al., 2006). Thermal energy dissipation is a primary de-excitation pathway measured by pulse-amplitude fluorescence as non-photochemical quenching (NPQ; Schreiber et al., 1986). At the same time, a small amount of radiation, solar-induced fluorescence (SIF), via the de-excitation of absorbed photons is emitted by photosystem II (Genty et al., 1989;Krause and Weis, 1991).
Some vegetation indices are sensitive to photoprotective pigments (e.g., carotenoids) and can characterize the seasonality of evergreen LUE with some success. For instance, the photochemical reflectance index (PRI; Gamon et al., 1992Gamon et al., , 1997 and chlorophyll/carotenoid index (CCI; Gamon et al., 2016) both use wavelength regions that represent carotenoid absorption features around 531 nm at the leaf level (Wong et al., 2019;Wong and Gamon, 2015a, b) and show great promise for estimating photosynthetic seasonality (Hall et al., 2008;Hilker et al., 2011a). Due to the relatively invariant canopy structure in evergreen forests, CCI and PRI have been applied at the canopy level as well (Gamon et al., 2016;Garbulsky et al., 2011;Middleton et al., 2016). In addition, the green chromatic coordinate (GCC; Richardson et al., 2009Richardson et al., , 2018Sonnentag et al., 2012), an index derived from the brightness levels of RGB canopy images, is also capable of tracking the seasonality of evergreen GPP . However, the full potential of spectrally resolved reflectance measurements to explore the photosynthetic phenology of evergreens has not been comprehensively explored at the canopy scale. The evaluation of pigment-driven spectral changes in evergreen forests over the course of a season is necessary to determine where, when, and why certain wavelength regions could advance our mechanistic understanding of canopy photosynthetic and photoprotective pigments. However, this has not been done with both empirical and process-based methods using continuously measured canopy hyperspectral reflectance and in situ pigment samples.
In addition to seasonal changes in pigment concentrations, canopy SIF was found to correlate significantly with the seasonality of photoprotective pigment content in a subalpine coniferous forest (Magney et al., 2019a). Steady-state SIF is regulated by NPQ and photochemistry (Porcar-Castell et al., 2014), and it provides complementary information on canopy GPP. Yang and van der Tol (2018) justified that the relative SIF, SIF normalized by the reflected near-infrared radiation, is more representative of the physiological variations in SIF as it is comparable to a SIF yield (Guanter et al., 2014;Genty et al., 1989). Our continuous optical measurements make it possible to differentiate mechanisms undergoing seasonal changes by comparing the decomposed reflectance spectrum against relative far-red SIF. Additionally, using relative SIF can effectively correct for incoming irradiance and account for the sunlit and shade fractions within the observation field of view (FOV) of PhotoSpec (Magney et al., 2019a).
In the present study, we analyzed continuous canopy reflectance data from PhotoSpec at a subalpine evergreen forest at the Niwot Ridge AmeriFlux site (US-NR1) in Colorado, US, and sought to understand the mechanisms controlling the seasonality of photosynthesis using continuous hyperspectral remote sensing. We first explored empirical techniques to study all seasonal variations in reflectance spectra, identified specific spectral regions that best explained the seasonal changes in GPP, and then linked these spectral features to pigment absorption features that impacted both biochemical and biophysical traits. We also used full-spectral inversions using a canopy RTM to infer quantitative estimates of leaf pigment pool sizes. Finally, we compared the spring onset of photosynthesis captured by different methods, VIs, and relative SIF to determine the underlying mechanisms that contributed to photosynthetic phenology.

Study site
The high-altitude (3050 m above sea level) subalpine evergreen forest near Niwot Ridge, Colorado, US, is an active AmeriFlux site (US-NR1, lat: 40.0329 • N, long: 105.5464 • W; tower height: 26 m; Monson et al., 2002;Burns et al., 2015Burns et al., , 2016Blanken et al., 2019). Three species dominate: subalpine fir (Abies lasiocarpa var. bifolia), Engelmann spruce (Picea engelmannii), and lodgepole pine (Pinus contorta) with an average height of 11.5 m, a leaf area index of 4.2 (Burns et al., 2016), and minimal understory. The annual mean precipitation and air temperature are 800 mm and 1.5 • C, respectively (Monson et al., 2002). The high elevation creates an environment with cold winters (with snow present more than half the year), while the relatively low latitude (40 • N) allows for year-round high solar irradiation (Monson et al., 2002). Thus, trees have to dissipate a considerable amount of excess sunlight during winter dormancy, which makes this forest an ideal site for studying seasonal variation in NPQ including the sustained component of it during dormancy Magney et al., 2019a).

Continuous tower-based measurements of canopy reflectance
PhotoSpec (Grossmann et al., 2018) is a 2D scanning telescope spectrometer unit originally designed to measure SIF. It also features a broadband Flame-S spectrometer (Ocean Optics, Inc., Florida, USA), used to measure reflectance from 400 to 900 nm at a moderate (full width at half maximum = 1.2 nm) spectral resolution with a FOV of 0.7 • (more details in Grossmann et al., 2018;Magney et al., 2019a). In the summer of 2017, we installed a PhotoSpec system on the top of the US-NR1 eddy-covariance tower, from where we can scan the canopy by changing both viewing azimuth angle and zenith angles. On every other summer day and every winter day, PhotoSpec scans the canopy by changing the view zenith angle with small increments at fixed view azimuth angles, i.e., elevation scans. Only one azimuth position is kept after 18 October 2017 to protect the mechanism from potentially damaging winter conditions at the site. Spectrally resolved reflectance was calculated using direct solar irradiance measurements via a cosine diffuser mounted in the upward nadir direction (Grossmann et al., 2018) as well as reflected radiance from the canopy. The reflectance data used in this study are from 16 June 2017 to 15 June 2018.
Here, we integrated all elevation scans to daily-averaged reflectance (every other day before 18 October 2017) by using all scanning viewing directions with vegetation in the field of view over the course of a day, filtering for both lowlight conditions and thick clouds by requiring PAR to be both at least 100 µmol m −2 s −1 and 60 % of theoretical clear-sky PAR. A detailed description of data processing can be found in Appendix B. To further test whether bidirectional reflectance effects impacted our daily averages, we compared the NDVI and NIRv at various canopy positions given a range of solar zenith and azimuth angles (Figs. A1-A3). Neither of the daily averaged VIs was substantially impacted by the solar geometry supporting the robustness of daily averaged canopy reflectance. An additional analysis (Fig. A4) has also shown the variation in phase angle at a daily time step is not a critical factor for the change in reflectance.
About 49 winter days exhibited significantly higher reflectances, attributable to snow within the field of view, which we corroborated with canopy RGB imagery from the tower. After removing data strongly affected by snow and excluding the days of instrument outages, 211 valid sample days remained, among which 96 valid sample days were between DOY 100 and 300. The daily-averaged reflectance was computed as the median reflectance from all selected scans for a single day, which was then smoothed by a 10-point (3.7 nm) box-car filter over the spectral dimension (400-900 nm) to remove the noise in the spectra. Figure 1a shows the seasonally averaged and spectrally resolved canopy reflectances measured by PhotoSpec.
To further emphasize the change in reflectance as a result of changes in pigment contents, we transformed the reflectance (shown as R λ ) using the negative logarithm (Eq. 1), as light intensity diminishes exponentially with pigment contents (Horler et al., 1983).
Here σ is the absorption cross section of pigments. Therefore, the log-transformed reflectance (Fig. 1b) should correlate more linearly with pigment contents (shown as C). We also considered a variety of typical VIs using the reflectance data from PhotoSpec. (Richardson et al., 2009). (3e) In order to calculate GCC, we convolved the reflectance using the instrumental spectral response function (  . For comparison, we normalized the reflectance by the value at 800 nm on each day. Here, we referred to 13 November-18 April as dormancy, and 2 June-21 August as the main growing season. The seasonal averaged canopy reflectance is composed of 39 daily-average reflectance in the growing season and 113 daily-averaged reflectance in the dormancy. In addition to the reflectance measurements, we also included relative SIF, far-red SIF normalized by the reflected near-infrared radiance at 755 nm. The far-red SIF (745-758 nm, Grossmann et al., 2018) was measured simultaneously with reflectance with a QEPro spectrometer (Ocean Optics, Inc., Florida, USA). The daily relative SIF was processed in the same fashion as the reflectance.

Eddy covariance measurements and LUE
Observations of net ecosystem exchange (net flux of CO 2 , NEE), PAR, and meteorological variables made at the US-NR1 tower are part of the official AmeriFlux Network data (Burns et al., 2016). GPP was estimated in half-hourly intervals (Reichstein et al., 2005) using the REddyProc package (Wutzler et al., 2018), allowing us to compute LUE (Goulden et al., 1996;Gamon et al., 2016) at half-hourly intervals.
According to the light response curves, GPP is a nonlinear function of PAR ( Fig. 2; Harbinson, 2012). Magney et al. (2019a) showed that fPAR does not significantly vary with seasons. We started to observe a photosynthetic saturation between 500 and 1000 µmol m −2 s −1 of PAR (Fig. 2), when the carboxylation rate, driven by maximum carboxylation rate (V cmax ), became the limiting factor (Farquhar et al., 1980). Thus, we defined the light-saturated GPP (GPP max ), as the mean half-hourly GPP at PAR levels between 1000 and 1500 µmol m −2 s −1 , a range which was covered throughout the year (Fig. 2), even in winter. Therefore, GPP max was less susceptible to short-term changes in PAR. Yet, due to the lower light intensity during storms, GPP max was not always available. As suggested by the low PAR value at which light saturation happened, plants remained in a light-saturated condition for most of the daytime. A higher GPP max indicates a greater V cmax and maximum electron transport rate (J max ) when the variation in GPP max is independent of stomatal conductance and intercellular CO 2 concentration (Leuning, 1995). Therefore, GPP max was closely correlated with daily LUE driven by physiology (see Sect. S2.4 in the Supplement).
We refrained from normalizing GPP max by absorbed photosynthetically active radiation (APAR) due to some of the APAR measurements (see Sect. S2.1 in the Supplement) not available in the beginning of growing season. GPP max was significantly linearly correlated with normalized GPP max by APAR (Fig. S2c).
We also included air temperature (T air ) and vapor pressure deficit (VPD) provided from the AmeriFlux network data. Daytime daily mean T air and VPD were computed from averaging the half-hourly T air and VPD when PAR was greater than 100 µmol m −2 s −1 .

Pigment measurements
To link canopy reflectance with variations in pigment contents, we used pigment data Bowling and Logan, 2019;Magney et al., 2019a) at monthly intervals over the course of the sampling period. Here, we focused on the xanthophyll cycle pool size (violaxanthin + antheraxanthin + zeaxanthin, V + A + Z), total carotenoid content (car), and total chlorophyll content (chl) measured on Pinus contorta and Picea engelmannii needles with units of moles per unit fresh mass. Car includes V + A + Z, lutein, neoxanthin, and beta-carotene. We also computed the ratio of chlorophyll to carotenoid contents (chl : car), because CCI derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) can track chl : car (Gamon et al., 2016). Overall, we can match 10 individual leaf-level sampling days for both pine and spruce samples with reflectance measured within ±2 d. Among these 10 valid sample days, 6 sample days are between DOY 100 and 300.

Data-driven spectral decomposition
We assumed that the spectrally resolved reflectance is a result of mixed absorption processes by different pigments. This allowed us to apply an independent component analysis (ICA; Hyvärinen and Oja, 2000) to decompose the logtransformed reflectance matrix (day of the year in rows and spectral dimension in columns) into its independent components. An advantage of the ICA is that it can separate a multivariate signal into additive subcomponents that are maximally independent, without the condition of orthogonality (Comon, 1994). We extracted three independent components, which explained more than 99.99 % of the variance, using the ICA algorithm (FastICA, Python package scikitlearn v0.21.0; Sect. S4 in the Supplement), such as where i is the ith component in spectral space. The decomposed spectral components revealed characteristic features that explain most of the variance in the reflectance matrix, which dictated the time-independent spectral shapes of pigment absorption features based on Eq. (1). The corresponding temporal loadings showed temporal variations in these spectral features, i.e., the variations in pigment contents. We will introduce the method of extracting pigment absorption features in a quantitative model-driven approach in Sect. 2.6.
In addition to analyzing the transformed reflectance alone, we empirically correlated the reflectance with GPP max using partial least-squares regression (PLSR, Python package scikit-learn v0.21.0). PLSR is a predictive regression model which solves for a coefficient that can maximally explain the linear covariance of the predictor with multiple variables (Wold et al., 1984;Geladi and Kowalski, 1986). PLSR has been used to successfully predict photosynthetic properties using reflectance matrices in previous studies from the leaf to canopy scales (e.g., Serbin et al., 2012Serbin et al., , 2015Barnes et al., 2017;Silva-Perez et al., 2018;Woodgate et al., 2019). Applying the PLSR to the hyperspectral canopy reflectance and GPP max resulted in a time-independent coefficient that emphasizes the key wavelength regions which contribute to the covariation of reflectance and GPP max , such as We implemented another set of PLSR analyses on the reflectance with individual pigment measurement as the target variable, such as the mean values of V + A + Z, car, and chl : car, such as We did not include chl as one of the target variables in this PLSR analysis since Bowling et al. (2018) and Magney et al. (2019a) have already shown chl did not vary seasonally in our study site. Fitting the minimal variance in chl will lead to overfitting the PLSR model. Comparing the PLSR coefficient of pigment measurements at the leaf level with the PLSR coefficient of GPP max connected the changes in GPP max to the pool size of photoprotective pigments, because the reflectance is regulated by the absorption of pigments.
We used PROSAIL (with PROSPECT-D; Féret et al., 2017) to compute the derivative of the daily-averaged negative logarithm transformed reflectance with respect to individual pigment contents, namely chlorophyll content (chlorophyll Jacobian, ∂−log(R) ∂C chl ) and carotenoid content (carotenoid Jacobian, ∂−log(R) ∂C car ) (Dutta et al., 2019). This helped explain the decomposed spectral components from the empirical analysis.
We also used PROSAIL to infer pigment contents (i.e., C chl , C car , C ant ) by optimizing the agreement between PROSAIL-modeled reflectance and measured canopy daily-mean reflectance from PhotoSpec. We fixed canopy structural parameters (e.g., the leaf area index (LAI) to 4.2, as reported in Burns et al., 2015) and fitted leaf pigment compositions as well as a low-order polynomial for soil reflectance (Appendix C), similar to Vilfan et al. (2018) and Féret et al. (2017). The cost function J in Eq. (7) represents a leastsquares approach, whereR is the modeled reflectance.
We used the spectral range between 450 and 800 nm, which encompasses most pigment absorption features.
3 Results and discussion

Seasonal cycle of GPP max and environmental conditions
As can be seen in Fig. 3, the subalpine evergreen forest at Niwot Ridge exhibits strong seasonal variation in GPP, T air , VPD, GPP max , and PAR. GPP and GPP max dropped to zero while sufficient PAR, required for photosynthesis, was still available in the dormancy, which suggests that the abiotic environmental factors impact photosynthesis seasonality nonlinearly and jointly. Abiotic factors played a strong role in regulating GPP max in this subalpine evergreen forest over the course of the season. For instance, there was a strong dependence of GPP max with T air . However, photosynthesis completely shut down during dormancy, even when the T air exceeded 5 • C (Fig. 3). During the onset and cessation periods of photosynthesis, GPP max rapidly increased with temperature ( Fig. S3a  left panel), potentially because needle temperature co-varied with T air , and needle temperature controls the activity of photosynthetic enzymes which affect V cmax . Spring warming approaches the optimal temperature for photosynthetic enzymes, leading to activation of photosynthesis, while cooling in the early winter inhibits these enzymes (Rook, 1969). Warming in spring melted frozen boles and made them available for water uptake , and thus caused the recovery of GPP max (Monson et al., 2005). Once the temperature was around the optimum (in the growing season), T air was no longer the determining factor for photosynthesis. Higher VPD caused by rising T air can stress the plants such that stomata closed, intercellular CO 2 reduced, and photosynthesis decreased (Fig. S3a right panel). When intercellular CO 2 concentration was not a limiting factor, GPP max was more representative of V cmax and did not vary T significantly.

Seasonal cycle of reflectance
In Fig. 4, the Jacobians show the maximum sensitivity of the reflectance spectral shape to carotenoid content at 524 nm, and near 566 and 700 nm for chlorophyll. The first peak of the chlorophyll Jacobian covers a wide spectral range in the visible range, while the second peak around the red edge is narrower.
It can be seen that the first spectral ICA component has a similar shape as the chlorophyll Jacobian. The corresponding temporal loading has a range between −0.2 and 0.2 without any obvious seasonal variation, consistent with a negligible seasonal cycle in chlorophyll content as shown in the pigment analysis. However, there is a gradual increase before DOY 50 in the first temporal loading, which appears to be anti-correlated with the temporal loading of the second ICA structure.
Two major features in the second spectral component can be observed. One is a negative peak centered around 530 nm, which aligns with the carotenoid Jacobian. At the negative logarithm scale, the negative values resulting from the negative ICA spectral peak multiplied by the positive ICA temporal loadings (growing season in Fig. 4 middle plots) indicate there were fewer carotenoids during the growing season ( Eqs. 1 and 4). Conversely, positive values resulting from a negative spectral peak multiplied by the negative temporal loadings (dormancy in Fig. 4 middle plots) indicate there were more carotenoids during dormancy (i.e., sustained photoprotection via the xanthophyll pigments; Bowling et al., 2018). Another feature is the valley-trough shape, which is co-located with the chlorophyll Jacobian center at the longer wavelength in the red-edge region. The center of this feature occurs at the shorter-wavelength edge of the chlorophyll Jacobian but does not easily explain changes in total chlorophyll content, which should show equal changes around 600 nm. The corresponding temporal loading apparently varied seasonally with GPP max .
The second temporal loading transitioned more gradually from dormancy to the peak growing season than GPP max . Unfortunately, we were missing data to evaluate the relative timing of GPP max cessation.
The third spectral component is similar to the mean shape of reflectance spectra. Its temporal loading remained around zero throughout the year.
Overall, the second ICA spectral component is more representative of the seasonal variation in the magnitude of total canopy reflectance than the other spectral components. The spectral changes around the red edge in the second component are interesting and might be related to structural needle changes in chlorophyll-a and chlorophyll-b contributions (de Tomás Marín et al., 2016;Rautiainen et al., 2018), which are not separated in PROSPECT.
CCI and PRI (Fig. 5a-b) followed the seasonal cycle of GPP max closely. CCI and PRI use reflectance near the center of the 530 nm valley feature (Eqs. 3c-3d), the spectral range that is most sensitive to the change of carotenoid content, so that they matched changes in GPP max very well. PRI was the smoothest throughout the year, without any significant fluctuations within the growing season, as opposite to what was observed in GPP max , which co-varied with T air and VPD ( Fig. S3a and b). This performance is intriguing given that PRI was originally developed to track short-term variations in LUE (Gamon et al., 1992), such as day-to-day and subseasonal scales. GCC (Fig. 5c) also correlated well with GPP max , but less than CCI and PRI. As can be seen in Fig. S1, the peak of the green channel used for GCC is close to the carotenoid Jacobian peak, while the red channel feature covers a part of the chlorophyll Jacobian feature. This explained the sensitivity of the GCC to changes in both carotenoid content and chlorophyll. The bands used in GCC are broader than the ones used by PRI and CCI; however it still captured these variations and can be computed using RGB imagery. Gentine and Alemohammad (2018) found that the green band helps to reconstruct variations in SIF using reflectances from MODIS. While they speculated that most variations in SIF are related to variations in PAR · fPAR (Gentine and Alemohammad, 2018), we suggest here that the green band indeed captures variations in LUE as well. NDVI (Fig. 5e) and NIRv (Fig. 5f) did not show an obvious seasonal variability.
Similar to the ICA components, all VIs were quite noisy during dormancy, especially prior to DOY 50. This noise may be due to snow because we only removed the reflectance when the canopy was snow covered. Scattered photons possibly still reached the telescope when there was snow on the ground, which is true for our study site as snowpack exists in winter .

PLSR coefficients of reflectance with GPP max and pigment measurements
The spectral shape of the PLSR coefficient with GPP max highlighted a peak (centering at 532 nm) near that of the carotenoid Jacobian with the same valley-trough feature observed near the second peak of the chlorophyll Jacobian (Fig. 6a). The reconstructed GPP max captured the onset and cessation of growth, while the day-to-day noise in reflectance during dormancy propagated to the reconstructed GPP max (−2 to 5 µmol m −2 s −1 ). During the growing season, the day-today variations in GPP max were not captured by any of the methods using pigment absorption features (Figs. 5a-c and 6b), which indicates those variations were not related to pigment content, but rather changes in environmental conditions that lead to day-to-day changes in photosynthesis (Fig. S3a). Overall, the observed GPP max was significantly correlated with the PLSR reconstruction (Pearson r 2 = 0.87), but very similar compared to CCI and PRI.
A similar PLSR model of reflectance but with pigment measurements (Fig. 7) showed a direct link between pigment contents and reflectance. It can be seen that the PLSR coefficients of reflectance are very similar, irrespective of the target variable. They feature a valley near the peak of the carotenoid Jacobian and a valley-trough feature near the peak at the longer wavelength of the chlorophyll Jacobian. This spectral shape is also very similar to the second ICA spectral component and PLSR coefficients of GPP max . V + A + Z, chl : car, and car were all nicely reconstructed by using the PLSR coefficients and reflectance (Fig. 7b). The reconstructed V + A + Z, car, and chl:car are correlated with the measured ones with Pearson r 2 values of 0.84, 0.71, and 0.93, respectively.
The second ICA component and PLSR empirically showed the seasonality of reflectance using two different empirical frameworks. ICA only used the reflectance, while the PLSR model accounts for variations in both reflectance and GPP max or pigment content. Yet, both ICA and PLSR agreed on similar spectral features that co-varied seasonally with GPP max . This indicates that the resulting spectral features were primarily responsible for representing this seasonal cycle. The overlap of these features with the chlorophyll/carotenoid absorption features showed that the seasonality of GPP max was related to variation in pigment content at the canopy scale, which was directly validated with a similar PLSR coefficient of reflectance and pigment contents. These results are consistent with leaf-level measurements of a higher ratio of chlorophyll to carotenoid content during the growing season in this forest (Fig. 7).
The highlighted spectral feature around 530 nm from ICA and PLSR closely overlaps with one of the bands used in CCI, PRI, and GCC (Eqs. 3a-3e), which provides a justification that these VIs can remarkably capture the LUE seasonality. The comparable Pearson r 2 values of PLSR, CCI, and PRI with GPP max suggest the pigment-driven seasonal cycle of GPP max is sufficiently represented by CCI and PRI. The spectral feature around the red edge does not make PLSR significantly more correlated with GPP max than CCI or PRI, which implies the feature is not driven by total chlorophyll or carotenoid contents.

Process-based estimation of pigment content
PROSAIL inversion results further supported the link between canopy reflectance, pigment contents, and GPP max . Figure 8 shows a continuous time series of C chl , C car , anthocyanin content (C ant ), and C chl C car derived from the PROSAIL canopy RTM inversion model. Examples of simulated and measured reflectance spectra shown are in Fig. C1. Anthocyanins are another type of photoprotective pigment (Pietrini et al., 2002;Lee and Gould, 2002;Gould, 2004) that protects the plants from high light intensity (Hughes, 2011). The pigment inversions closely matched the seasonality of GPP max .
C chl C car showed the greatest sensitivity in capturing the seasonal cycle, with the strongest correlation to leaf level measurements (Fig. 8c). The inverted C chl had the weakest empirical relationship with the measured one ( Fig. 8a right panel). Apparently, some of the inversion errors of individual C car and C chl contents canceled out in the ratio, making the ratio more stable. C ant performed similarly to C car , since they both are photoprotective, and the anthocyanins absorb at 550 nm (Sims and Gamon, 2002), which is close to the center of carotenoid absorption feature. Even though we lacked field measurements of anthocyanins to validate anthocyanin retrievals, the inversions showed that more than just carotenoid content can be obtained from full-spectral inversions.
Strictly speaking, the complex canopy structure of evergreens makes the application of 1D canopy RTMs such as PROSAIL difficult Zarco-Tejada et al., 2019). Yet, Moorthy et al. (2008), Ali et al. (2016), andZarco-Tejada et al. (2019) reasonably discussed the pigment retrieval in conifer forests with careful applications. In our study, the reflectance was collected from needles with a very small FOV, and our study site has a very stable canopy structure throughout the year (Burns et al., 2016). Thus, the inversion results are meaningful for discussing the seasonality of pigment contents. In the future, radiative transfer models that properly describe conifer forests, such as LIBERTY (Dawson et al., 1998), could be used.

Comparison across methods
Although decomposing the hyperspectral canopy reflectance and using relative SIF (Fig. 5d) both successfully tracked the seasonal cycle of evergreen LUE, they underlie different de-excitation processes. During the growing season, environmental conditions primarily drove the day-to-day variations in GPP max . Relative SIF responded to such environmental stresses  so that it appeared to track sub-seasonal variations better than reflectance, particularly during the growing season (Fig. S5f). Yet, reflectance decompositions and VIs were less sensitive to such day-to-day variations (Figs. 6, S3b).
There was also some variability between reflectance-based methods and relative SIF during the transition periods between the growing season and dormancy. We focused on the growing season onset since the reflectance measurements were not available during the cessation period. The onset (DOY 60 to 166) described by all the methods mentioned above as well as the relative SIF are compared in Fig. 9, using a sigmoid fit to available data (Fig. D1). The observed GPP max had the most rapid yet latest growing onset. The methods and VIs derived from or related to the pigment contents increased earlier than GPP max -such as the ICA component, PLSR coefficient, PROSAIL C chl C car , and CCI. However, they built up slowly to reach the maximum, which sug-gests that reduction of the carotenoid content is a slower process than the recovery of LUE. Reflectance-based VIs (Fig. 5) and decomposing methods (Figs. 4 and 8b, c) had a slower growing season onset than GPP max , as found in Bowling et al. (2018) as well. On the other hand, relative SIF started the onset at almost the same time as the GPP max , and it quickly reached the maximum. Therefore, using both SIF and reflectance to constrain the LUE prediction  can further improve the prediction accuracy.

Conclusion and future work
In this study, we analyzed seasonal co-variation in GPP and the spectrally resolved visible and near-infrared reflectance signal, as well as several commonly used VIs. The main spectral feature centered around 530 nm is most important for inferring the seasonal cycle of reflectance (400-900 nm) 4534 R. Cheng et al.: Decomposing reflectance spectra to track gross primary production Figure 8. The left panels are the estimations of (a) C chl , (b) C car , C ant , and (c) C chl C car from the PROSAIL overlaid with the GPP max . We normalized two metrics because they report the pigment contents in different units. The vertical dashed line divides the observations from DOY for the years 2017 and 2018. The plots on the right compare the pigment contents from leaf-level measurements and using PROSAIL: (a) chl vs. C chl , (b) car vs. C car , and (c) chl : car vs. C chl C car . The correlations are statistically significant except C chl . Figure 9. Temporal evolution of the growing season onset using sigmoid fits (scaled) of PLSR, ICA, CCI, chlorophyll-to-carotenoid ratio, and relative SIF. and LUE, which corresponds to changes in carotenoid content. This explains why CCI, PRI, and GCC track GPP seasonality so well, as most variations are driven by carotenoid pool changes. Our analysis included RTM simulation and in situ pigment measurements throughout the season, confirming the link between reflectance/VIs and pigment con-tents. The comparison of reflectance/VIs and relative SIF reveals differences in the timing of the growing season onset, pigment changes, and SIF, indicating the potential of using both reflectance and SIF to track the seasonality of photosynthesis. However, the close correspondence between both SIF and reflectance suggests that hyperspectral reflectance alone provides mechanistic evidence for a robust approach to track photosynthetic phenology of evergreen systems. Because seasonal variation in pigment concentration plays a strong role in regulating the seasonality of photosynthesis in evergreen systems, our work will help to inform future studies using hyperspectral reflectance to achieve accurate monitoring of these ecosystems. While indices like PRI and CCI are performing sufficiently as our methods which use the full-spectrum analysis at the canopy scale, the application of the full spectrum might be more robust for space-based measurements. In addition, we found seasonal changes of canopy reflectance near the red-edge region, which could be related to leaf structural changes or changes in chlorophyll a and b. Our PLSR coefficients are good references for customizing VIs to infer the photosynthetic seasonality in evergreen forest when there are restrictions to use the specific bands from currently existing VIs (such as PRI and CCI). While our current study is limited to a subalpine evergreen forest and canopyscale measurements, applications to other regions, vegetation types, and observational platforms will be a focus for future research.
R. Cheng et al.: Decomposing reflectance spectra to track gross primary production 4535 Appendix A: Bidirectional reflectance effect

A1 NDVI and NIRv
The impact of geometry and small FOV is relatively negligible. First, our method only used the scans when FOV is on the needles by setting a NDVI threshold. Second, we plotted the NDVI and NIRv against the solar geometry at each individual tree target throughout a year. NDVI and NIRv are quite homogeneous regardless of various solar geometries as shown in the following figures.

4536
R. Cheng et al.: Decomposing reflectance spectra to track gross primary production A2 PLSR on phase angle and reflectance We did a PLSR analysis on individual measurements of phase angle and reflectance for 3 summer days (1 to 3 July 2017). The results are the same from other sample days. Indeed, the reflectance has different sensitivities to the phase angle. However, the poor correlation of PLSR reconstructed phase angle and the measured one suggests the variations in phase angle should not be the critical factor for the change in reflectance. In our study, we primarily removed the bidirectional impact by averaging all the individual reflectance that was measured at different solar geometry and viewing geometry. Appendix B: Detailed processes on integrating daily-averaged canopy reflectance First, we chose scans targeting vegetation only by requiring an NDVI greater than 0.6. Second, it is important to ensure that the solar irradiation did not change between the acquisition of the solar irradiance and the reflected radiance measurement. To achieve this, we matched the timestamps of a PAR sensor (LI-COR LI-190SA, LI-COR Environmental, Lincoln, Nebraska, US) to the timestamps of PhotoSpec, and we compared the PAR value from the PAR sensor during the PhotoSpec irradiance acquisition with PAR during the actual target scan of the reflected radiance from vegetation. We only used the scans when the ratio of the two was 1.0 ± 0.1, ensuring stable PAR conditions. Third, in order to avoid unstable PAR because of clouds (Dye, 2004), we also removed cloudy scenes by requiring PAR to be at least 60 % of a theoretical maximum driven by solar geometry (Fig. B1). Further, only data when PAR was greater than 100 µmol m −2 s −1 were considered to eliminate the impact of low solar angles on reflectance data. The VIs shown in Fig. 5 were extracted in the same fashion as above. Figure B1. The distribution of the ratio of the measured PAR to the PAR at theoretical maximum from all individual scans.

4538
R. Cheng et al.: Decomposing reflectance spectra to track gross primary production Appendix C: PROSAIL fits We used the following range constraints for variables included in the state vector of the PROSAIL inversion.  . Individual sigmoid fits of the onset of growth from different methods and more VIs. The fitted curve has been expressed as the derivation as above. The Pearson r 2 and p values listed in each subplot were calculated from the correlation of observed and fitted variables. The residual was calculated as the average L2 norm of the difference between observed (y) and fitted variables (ŷ) normalized by the observation, i.e., 1 n i ( y−ŷ y ) 2 . The fittings are overall good. Because the ICA loading lacks a clear sigmoid shape, ICA has a larger residual.
The first derivative of y is At the half maximum point (x = x half ), y = a+b 2 . Therefore, we need to solve Hence, x half = d.