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

Photosynthesis by terrestrial plants represents the majority of CO₂ 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 co-varied 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 near-infrared vegetation index (NIRv) in this ecosystem. In addition, we quantitatively determined needle-scale chlorophyll-to-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, highlighting all spectral locations within 400–900 nm associated with GPP seasonality in evergreen forests.

1 Introduction

Terrestrial gross primary production (GPP), the gross CO₂ uptake through photosynthesis, is the largest uptake of atmospheric CO₂ (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., 1992, 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):

$$\begin{array}{}\text{(1)}& \mathrm{GPP}=\mathrm{PAR}\cdot \mathrm{fPAR}\cdot \mathrm{LUE},\end{array}$$

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 (Bowling et al., 2018), 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., 1992, 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., 2009, 2018; Sonnentag et al., 2012), an index derived from the brightness levels of RGB canopy images, is also capable of tracking the seasonality of evergreen GPP (Bowling et al., 2018). 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.

Here, we used continuous measurements in both spectral space (full spectrum between 400 and 900 nm) and time (daily over an entire year) to evaluate the potential of hyperspectral canopy reflectance for better understanding the sensitivity of VIs to pigment changes that regulate GPP in evergreen forests. Continuous measurements of spectrally resolved reflectance at the canopy scale have so far been sparse at evergreen forest sites (Gamon et al., 2006; Hilker et al., 2011b; Porcar-Castell et al., 2015; Rautiainen et al., 2018; Wong et al., 2020). There are only a few empirical studies on hyperspectral canopy reflectance in evergreen forests (Smith et al., 2002; Singh et al., 2015). Yet, empirically decomposed canopy spectral reflectance has been used as a predictor of maximum photosynthetic capacity (Serbin et al., 2012; Barnes et al., 2017; Dechant et al., 2017; Silva-Perez et al., 2018; Meacham-Hensold et al., 2019), GPP (Matthes et al., 2015; Huemmrich et al., 2017; DuBois et al., 2018; Huemmrich et al., 2019; Dechant et al., 2019), and other physiological properties (Ustin et al., 2004, 2009; Asner et al., 2011; Serbin et al., 2014).

In contrast to empirical methods, process-based approaches, such as canopy radiative transfer models (RTMs) can help to quantitatively link canopy photosynthesis with leaf-level contents of photosynthetic/photoprotective pigments (Feret et al., 2008; Jacquemoud et al., 2009). With RTMs, we can use spectrally resolved reflectance to directly derive leaf pigment contents (Féret et al., 2017; Jacquemoud et al., 1995) and plant traits (Féret et al., 2019).

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.

2 Material and methods

2.1 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., 2015, 2016; Blanken 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 (Bowling et al., 2018; Magney et al., 2019a).

2.2 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 low-light conditions and thick clouds by requiring PAR to be both at least 100 µmol m⁻² s⁻¹ 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.

Figure 1 (a) Seasonally averaged canopy reflectance in winter dormancy (red) and the growing season (black) from PhotoSpec. (b) Seasonally averaged negative logarithm transformation of reflectance (400–900 nm). 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.

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. 2), as light intensity diminishes exponentially with pigment contents (Horler et al., 1983).

$$\begin{array}{}\text{(2a)}& {\displaystyle }& {\displaystyle }{R}_{\mathit{\lambda }}\propto \exp (-C\cdot \mathit{\sigma }(\mathit{\lambda }\left)\right)\text{(2b)}& {\displaystyle }& {\displaystyle }C\cdot \mathit{\sigma }\left(\mathit{\lambda }\right)\propto -\log \left(R_{\mathit{\lambda }}\right)\end{array}$$

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.

$$\begin{array}{}\text{(3a)}& {\displaystyle }\mathrm{NDVI}={\displaystyle \frac{R_{\mathrm{800}}-R_{\mathrm{670}}}{R_{\mathrm{800}}+R_{\mathrm{670}}}}\text{(Rouse et al., 1974)}\text{(3b)}& {\displaystyle }\mathrm{NIRv}=\mathrm{NDVI}\cdot R_{\mathrm{800}}\text{(Badgley et al., 2017)}\text{(3c)}& {\displaystyle }\mathrm{PRI}={\displaystyle \frac{R_{\mathrm{531}}-R_{\mathrm{570}}}{R_{\mathrm{531}}+R_{\mathrm{570}}}}\text{(Gamon et al., 1992)}\text{(3d)}& {\displaystyle }\mathrm{CCI}={\displaystyle \frac{R_{\mathrm{526}\text{–}\mathrm{536}}-R_{\mathrm{620}\text{–}\mathrm{670}}}{R_{\mathrm{526}\text{–}\mathrm{536}}+R_{\mathrm{620}\text{–}\mathrm{670}}}}\text{(Gamon et al., 2016)}\text{(3e)}& {\displaystyle }\begin{array}{rl}\mathrm{GCC}& ={\displaystyle \frac{R_{\mathrm{Green}}}{R_{\mathrm{Red}}+R_{\mathrm{Green}}+R_{\mathrm{Blue}}}}\\ & \text{(Richardson et al., 2009)}.\end{array}\end{array}$$

In order to calculate GCC, we convolved the reflectance using the instrumental spectral response function (Fig. S1 in the Supplement; Wingate et al., 2015) of the StarDot NetCam SC 5 MP IR (StarDot Technologies, Buena Park, CA, USA), which is the standard camera model used by the PhenoCam Network protocol (Sonnentag et al., 2012).

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.

2.3 Eddy covariance measurements and LUE

Observations of net ecosystem exchange (net flux of CO₂, 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⁻² s⁻¹ 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⁻² s⁻¹, 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₂ 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⁻² s⁻¹.

Figure 2 Half-hourly GPP as a function of PAR during the measurement period. Points were colored by month. Bold points were the median GPP when PAR was binned every 100 µmol m⁻² s⁻¹ approximately. The solid lines represent the canopy light response curve.

2.4 Pigment measurements

To link canopy reflectance with variations in pigment contents, we used pigment data (Bowling et al., 2018; 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.

2.5 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 log-transformed 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 scikit-learn v0.21.0; Sect. S4 in the Supplement), such as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}-\log \left(R_{\mathit{\lambda },\mathrm{DOY}}\right)=\sum _{i=\mathrm{1},\mathrm{2},\mathrm{3}}(& \mathrm{spectral}{\mathrm{component}}_{\mathit{\lambda }}^i\\ & \cdot \mathrm{temporal}{\mathrm{loading}}_{\mathrm{DOY}}^i),\end{array}\end{array}$$

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. (2). 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., 2012, 2015; Barnes 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

$$\begin{array}{}\text{(5)}& {\mathrm{GPP}}_{\max ,\mathrm{DOY}}=-\log \left(R_{\mathit{\lambda },\mathrm{DOY}}\right)\times \mathrm{PLSR}{\mathrm{coefficient}}_{\mathit{\lambda }}^{{\mathrm{GPP}}_{max}}.\end{array}$$

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

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathrm{pigment}& \mathrm{measurement}=-\log \left(R_{\mathit{\lambda },\mathrm{DOY}}\right)\\ & \times \mathrm{PLSR}{\mathrm{coefficient}}_{\mathit{\lambda }}^{\mathrm{pigment}\mathrm{measurement}}.\end{array}\end{array}$$

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.

2.6 Process-based methods

PROSPECT + SAIL (PROSAIL; Jacquemoud et al., 2009) is a process-based 1-D canopy RTM that models canopy reflectance, given canopy structure information (SAIL) as well as leaf pigment contents (PROSPECT) (Jacquemoud and Baret, 1990; Vilfan et al., 2018).

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, $\frac{\partial -\log \left(R\right)}{\partial C_{\mathrm{chl}}}$) and carotenoid content (carotenoid Jacobian, $\frac{\partial -\log \left(R\right)}{\partial C_{\mathrm{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 least-squares approach, where $\widehat{R}$ is the modeled reflectance.

$$\begin{array}{}\text{(7)}& J=\sum _{\mathit{\lambda }=\mathrm{450}\mathrm{nm}}^{\mathrm{800}\mathrm{nm}}(R_{\mathit{\lambda }}-{\widehat{R}}_{\mathit{\lambda }}{)}^{\mathrm{2}}.\end{array}$$

We used the spectral range between 450 and 800 nm, which encompasses most pigment absorption features.

3 Results and discussion

3.1 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 (Bowling et al., 2018), 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₂ reduced, and photosynthesis decreased (Fig. S3a right panel). When intercellular CO₂ concentration was not a limiting factor, GPP_max was more representative of V_cmax and did not vary T significantly.

Figure 3 Daily-averaged time series of air temperature (T_air), vapor pressure deficit (VPD), photosynthetically active radiation (PAR), gross primary production (GPP) from half-hourly data when PAR was greater than 100 µmol m⁻² s⁻¹, and time series of GPP_max. DOY 166 (2017) was the first day of observation. The vertical dashed line divides the observations from day of year (DOY) for the years 2017 and 2018.

3.2 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.

Figure 4 A set of three spectral components (top, colored) and corresponding temporal loadings (bottom, colored) from ICA decomposition. The first spectral component is overlaid with the chlorophyll Jacobian ($\frac{\partial -\log \left(R\right)}{\partial C_{\mathrm{chl}}}$, dashed–dotted), and the second spectral component is overlaid with the carotenoid Jacobian ($\frac{\partial -\log \left(R\right)}{\partial C_{\mathrm{car}}}$, dotted). The third spectral component is overlaid on the annual mean shape of transformed reflectance spectra. Temporal loadings are overlaid with GPP_max (grey line). The axis of Jacobians is not shown because its magnitude is arbitrary here. The vertical dashed line divides the observations from DOY for the years 2017 and 2018.

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. 2 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 sub-seasonal 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 (Bowling et al., 2018).

Figure 5 Magenta points are time series of VIs: (a) CCI, (b) PRI, (c) GCC, (d) relative SIF, (e) NDVI, (f) NIRv. The grey points in the background show GPP_max. The Pearson r² values of regressing VIs and GPP_max are noted in each plot. The p values of all correlations in this figure are less than 0.005. The vertical dashed line divides the observations from DOY for the years 2017 and 2018.

3.3 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⁻² s⁻¹). During the growing season, the day-to-day 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²=0.87), but very similar compared to CCI and PRI.

Figure 6 (a) The PLSR coefficient of reflectance with GPP_max is the blue line. The overlaid dashed–dotted and dotted lines are chlorophyll and carotenoid Jacobians, respectively. The overlaid orange solid line is the second ICA spectral component, which was scaled to fit to the plot without a y axis. (b) The reconstructed GPP_max (blue) by PLSR is overlaid with the observed GPP_max (red). The vertical dashed line divides the observations from DOY for the years 2017 and 2018.

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² values of 0.84, 0.71, and 0.93, respectively.

Figure 7 (a) PLSR coefficients of reflectance and three pigment measurements. The overlaid dashed–dotted and dotted lines are chlorophyll and carotenoid Jacobians, respectively. The overlaid solid grey line is the second ICA spectral component, which is scaled to fit to the plot without a y axis. (b) The reconstructed pigment measurements (blue) by PLSR are overlaid with the measured mean pigment measurements (red). The error bar is 1 standard deviation of the measurements. The vertical dashed line divides the observations from DOY for the years 2017 and 2018.

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² 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.

3.4 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 $\frac{C_{\mathrm{chl}}}{C_{\mathrm{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. $\frac{C_{\mathrm{chl}}}{C_{\mathrm{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 (Jacquemoud et al., 2009; Zarco-Tejada et al., 2019). Yet, Moorthy et al. (2008), Ali et al. (2016), and Zarco-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.

Figure 8 The left panels are the estimations of (a) C_chl, (b) C_car, C_ant, and (c) $\frac{C_{\mathrm{chl}}}{C_{\mathrm{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. $\frac{C_{\mathrm{chl}}}{C_{\mathrm{car}}}$. The correlations are statistically significant except C_chl.

3.5 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 (van der Tol et al., 2014) 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 $\frac{C_{\mathrm{chl}}}{C_{\mathrm{car}}}$, and CCI. However, they built up slowly to reach the maximum, which suggests 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 (van der Tol et al., 2014) can further improve the prediction accuracy.

Figure 9 Temporal evolution of the growing season onset using sigmoid fits (scaled) of PLSR, ICA, CCI, chlorophyll-to-carotenoid ratio, and relative SIF.

4 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) 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 contents. 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 canopy-scale measurements, applications to other regions, vegetation types, and observational platforms will be a focus for future research.

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.

Figure A1 NDVI and NIRv of all scans targeting a pine at different solar azimuth angles and solar zenith angles throughout a year.

Figure A2 NDVI and NIRv of all scans targeting a fir at different solar azimuth angles and solar zenith angles throughout a year.

Figure A3 NDVI and NIRv of all scans targeting a spruce at different solar azimuth angles and solar zenith angles throughout a year.

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.

Figure A4 PLSR analysis on phase angle and reflectance.

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⁻² s⁻¹ 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.

Appendix C: PROSAIL fits

We used the following range constraints for variables included in the state vector of the PROSAIL inversion.

• Leaf mesophyll structure (N): 0.9–1.1

• Chlorophyll content (C_chl): 0–120 µmol cm⁻²

• Carotenoid content (C_car): 0–70 µmol cm⁻²

• Anthocyanin content (C_ant): 0–10 µmol cm⁻²

• Brown pigments (C_brown): 0–0.6

• Water content (C_w): 0–0.2 cm

• Dry matter content (C_m): 0–0.2 g cm⁻²

• Xanthophyll cycle status (C_x) 0–1

• Leaf area index (LAI): fixed to 4.2

Figure C1 The observed and fitted reflectance spectra at low (a) and high (b) $\frac{C_{\mathrm{chl}}}{C_{\mathrm{car}}}$.

Appendix D: Sigmoid fit

The sigmoid equation is

$$y=b+{\displaystyle \frac{a-b}{\mathrm{1}+\exp \left(\frac{d-x}{c}\right)}}.$$

In this form, a and b represent the maximum and minimum values of the sigmoid fit. And d is the half maximum of the fit. We obtained the optimal values of these parameters.

Proof:

If $x\to +\mathrm{\infty }$, $\exp \left(\frac{d-x}{c}\right)\to \mathrm{0}$. So,

$$\underset{x\to +\mathrm{\infty }}{lim}y=a.$$

If $x\to -\mathrm{\infty }$, $\exp \left(\frac{d-x}{c}\right)\to +\mathrm{\infty }$. So,

$$\underset{x\to -\mathrm{\infty }}{lim}y=b.$$

Figure D1 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² 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 ($\widehat{y}$) normalized by the observation, i.e., $\frac{\mathrm{1}}{n}{\sum }_i(\frac{y-\widehat{y}}{y}{)}^{\mathrm{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

$${\displaystyle \frac{\mathrm{d}y}{\mathrm{d}x}}={\displaystyle \frac{a-b}{(\mathrm{1}+\exp \left(\frac{d-x}{c}\right){)}^{\mathrm{2}}}}\exp \left({\displaystyle \frac{d-x}{c}}\right){\displaystyle \frac{\mathrm{1}}{c}}.$$

At the half maximum point (x=x_half), $y=\frac{a+b}{\mathrm{2}}$. Therefore, we need to solve

$${\displaystyle \frac{a+b}{\mathrm{2}}}=b+{\displaystyle \frac{a-b}{\mathrm{1}+\exp \left(\frac{d-x_{\mathrm{half}}}{c}\right)}}.$$

Hence, x_half=d.