Confounding effects of snow cover on remotely sensed vegetation indices of evergreen and deciduous trees: An experimental study

Located at northern latitudes and subject to large seasonal temperature fluctuations, boreal forests are sensitive to the changing climate, with evidence for both increasing and decreasing productivity, depending upon conditions. Optical remote sensing of vegetation indices based on spectral reflectance offers a means of monitoring vegetation photosynthetic activity and provides a powerful tool for observing how boreal forests respond to changing environmental conditions. Reflectance‐based remotely sensed optical signals at northern latitude or high‐altitude regions are readily confounded by snow coverage, hampering applications of satellite‐based vegetation indices in tracking vegetation productivity at large scales. Unraveling the effects of snow can be challenging from satellite data, particularly when validation data are lacking. In this study, we established an experimental system in Alberta, Canada including six boreal tree species, both evergreen and deciduous, to evaluate the confounding effects of snow on three vegetation indices: the normalized difference vegetation index (NDVI), the photochemical reflectance index (PRI), and the chlorophyll/carotenoid index (CCI), all used in tracking vegetation productivity for boreal forests. Our results revealed substantial impacts of snow on canopy reflectance and vegetation indices, expressed as increased albedo, decreased NDVI values and increased PRI and CCI values. These effects varied among species and functional groups (evergreen and deciduous) and different vegetation indices were affected differently, indicating contradictory, confounding effects of snow on these indices. In addition to snow effects, we evaluated the contribution of deciduous trees to vegetation indices in mixed stands of evergreen and deciduous species, which contribute to the observed relationship between greenness‐based indices and ecosystem productivity of many evergreen‐dominated forests that contain a deciduous component. Our results demonstrate confounding and interacting effects of snow and vegetation type on vegetation indices and illustrate the importance of explicitly considering snow effects in any global‐scale photosynthesis monitoring efforts using remotely sensed vegetation indices.


| INTRODUC TI ON
The boreal zone covers about 1.89 billion ha in the northern hemisphere and is defined as the high northern latitude broad vegetation zone that is covered principally with forests, wetlands, lakes, and rivers (Brandt et al., 2013).Boreal forest function and health, including photosynthesis and respiration, are closely linked to temperature and precipitation patterns (Bonan, 2008).Altered growing season temperature and moisture availability with changing climate will lead to further effects on plant productivity (Bunn et al., 2007).Warmer temperatures are likely to cause an earlier spring onset and a delayed fall decline in photosynthesis.Thus, forests may yield greater carbon uptake and productivity as a result of a longer period with green foliage and an extended period of favorable conditions (Estiarte & Peñuelas, 2015).While boreal species can show positive photosynthetic and growth responses to warming (Way et al., 2015;Way & Oren, 2010), warmer temperatures also lead to higher respiration (Dunn et al., 2007) and influence the frequency and severity of disturbances such as drought, fire, or insect infestation, which directly impact vegetation succession and reduce ecosystem productivity (Dial et al., 2022;Moreau et al., 2020;Walker & Johnstone, 2014;White et al., 2017).Changing climate also alters snow cover (IPCC, 2007), potentially affecting the ability of remote sensing to assess changing vegetation productivity (Myers-Smith et al., 2020).
Together, these effects lead to considerable uncertainty in the prediction of ecosystem productivity in a changing world.
Remote sensing provides a means of monitoring plant physiological and phenological processes at multiple spatial and temporal scales.Vegetation indices derived from optical remote sensing have been widely used to monitor photosynthetic phenology of different terrestrial ecosystems (DeFries & Townshend, 1994;Running et al., 2004;Schimel et al., 2015;Schimel & Schneider, 2019).For example, utilizing reflectance in the red and near-infrared bands, the normalized difference vegetation index (NDVI; Tucker, 1979) provides an estimate of vegetation greenness.NDVI can be used as an indicator for seasonal phenology of green biomass, leaf display, and photosynthesis for annual and deciduous species, but the NDVIphotosynthesis relationships are less clear for evergreen species (Gamon et al., 1995).Studies have used NDVI or similar indices as direct productivity indicators (Myneni et al., 1997) or as an input parameter to drive productivity estimation models (Goetz et al., 2005;Running et al., 2004).Such greenness indices, including NDVI, the enhanced vegetation index (Huete et al., 2002) and the newly developed near-infrared reflectance of vegetation (NIRv; Badgley et al., 2017Badgley et al., , 2019)), largely fail to detect changes in photosynthetic activity in evergreen species that undergo seasonal downregulation and upregulation without apparent changes in canopy structure (Gamon et al., 1995(Gamon et al., , 2015(Gamon et al., , 2016;;Springer et al., 2017;Wang et al., 2023;Zeng et al., 2022).This limitation has led to the development of additional indices that are more directly related to the actual photosynthetic regulatory processes invisible to NDVI and similar greenness indices.
The photochemical reflectance index (PRI) is capable of detecting these subtle changes in foliage relating to the regulation of photosynthetic activity (Gamon et al., 1992(Gamon et al., , 1997)).On shorter (diurnal) time scales, PRI changes are driven by changes in xanthophyll cycle pigments, including violaxanthin, antheraxanthin, and zeaxanthin.
Xanthophyll cycle activity leads to dissipation of excess light energy as heat, and the quenching of fluorescence via non-photochemical quenching, thereby protecting the photosynthetic apparatus (Demmig-Adams & Adams, 1996;Jahns & Holzwarth, 2012;Niyogi et al., 1997).For evergreens, zeaxanthin and other carotenoid levels remain high during cold period, maximizing dissipation of light energy to protect leaves during winter (Adams III et al., 2002;Öquist & Huner, 2003).In this case, PRI variation over long time periods is influenced by the seasonal change in leaf pigment pool sizes (Gamon & Berry, 2012;Garbulsky et al., 2011;Hmimina et al., 2014;Wong & Gamon, 2015a, 2015b), which are closely tied to seasonal gross primary productivity (GPP) patterns (Wong et al., 2020).Consequently, PRI can indicate both short-term (e.g., diurnal) downregulation due to xanthophyll cycle pigment activity, and longer-term (e.g., seasonal) adjustments in photosynthetic activity due to adjustments in photosynthetic and photoprotective pigment pools.The sensitivity to subtle changes in pigments makes PRI a powerful tool to monitor photosynthetic phenology of evergreen species (Wong et al., 2019(Wong et al., , 2020;;Wong & Gamon, 2015b).However, the application of PRI in tracking boreal forest phenology at global scales is currently limited by the lack of suitable bands in most current satellite sensors.
The chlorophyll/carotenoid index (CCI), which utilizes MODIS bands 1 and 11, has been proposed as an alternative indicator of photosynthetic activity, particularly for evergreen species (Gamon et al., 2016).Like the PRI, CCI is sensitive to changes in pigment pool sizes and can accurately track seasonally changing chlorophyll/carotenoid levels and photosynthetic activity at both the leaf and stand level over seasonal courses (Gamon et al., 2016;Pierrat, Magney, et al., 2022;Wong et al., 2019Wong et al., , 2020Wong et al., , 2022;;Yang et al., 2020Yang et al., , 2022)).
Consequently, CCI can be a useful indicator of altered photosynthetic activity tied to seasonally changing chlorophyll and carotenoid any global-scale photosynthesis monitoring efforts using remotely sensed vegetation indices.
Together, the combination of vegetation indices, such as NDVI, PRI, and CCI, can provide complementary information about monitoring ecosystem productivity phenology of evergreen and deciduous forests (Gamon, 2015).For deciduous species, NDVI is closely related to the fraction of absorbed photosynthetically active radiation and allows monitoring of changes in canopy greenness, such as leaf development in the spring and senescence and abscission in the fall (Balzarolo et al., 2016).In contrast, evergreen species maintains their needles through the year and show little variation in canopy greenness (Hmimina et al., 2013;Wu et al., 2014).In this case, the carotenoid sensitive indices (PRI and CCI) can be directly linked to evergreen productivity and used to estimate photosynthetic light use efficiency (Springer et al., 2017;Wong et al., 2020Wong et al., , 2022)).The integration of vegetation indices enables us to detect contrasting, complementary physiological, and structural controls of productivity for deciduous and evergreen species.
Despite the widespread use of satellite remote sensing techniques for studies of changing global productivity, few of these studies explicitly consider snow effects on estimates of changing productivity.Changing patterns of winter snow cover in the northern latitude and high elevation regions can influence ecosystem productivity estimation using satellite optical measurements due to the high visible reflectance of snow (Gamon et al., 2013;Jin et al., 2017;Jin & Eklundh, 2014;Springer et al., 2017), which alters different spectral regions to different degrees (Negi et al., 2009;Saha et al., 2019;Singh et al., 2010).These effects can readily alter the values of reflectance indices and confound interpretation of changing function in northern latitude ecosystems (Myers-Smith et al., 2020), yet this issue is rarely clearly addressed in most satellite-based studies that use reflectance-based vegetation indices.Similarly, the interaction of vegetation type with snow cover is not entirely clear, although considerable evidence exists to suggest that snow cover effects vary between different vegetation types (Bokhorst et al., 2016;Marsh et al., 2010).Given the complexities of varying vegetation type and snow cover, it is possible that varying snow cover can explain at least part of the widely reported failure for models to reliably depict the seasonal course of GPP (e.g., Keenan et al., 2012;Richardson et al., 2012;Rogers et al., 2022;Schaefer et al., 2012), particularly during transitional spring conditions.
The normalized difference snow index (NDSI) derived from the MODIS data utilizes the green and short-wave infrared bands and can be used to estimate snow coverage at continental and global scales (Riggs et al., 2016).However, the application of NDSI is limited by the malfunction of Aqua SWIR detectors (Gladkova et al., 2012) and its sensitivity to low solar irradiance caused by large solar zenith angle or at complex landscapes (Lv & Pomeroy, 2019;Wang et al., 2018).Despite the efforts in developing vegetation indices that may be less sensitive to snow cover (Camps-Valls et al., 2021;Jin et al., 2017), there are few explicit studies of snow effects on reflectance indices.These effects of snow cover have been primarily considered for some more common indices (e.g., NDVI; Gamon et al., 2013;Myers-Smith et al., 2020), but less studied for others (e.g., PRI and CCI), even though there is emerging evidence for such effects (Pierrat, Magney, et al., 2022).
Several recent studies have reported clear confounding effects of snow on vegetation indices from both tower-based measurements (Pierrat et al., 2021;Pierrat, Magney, et al., 2022) and satellite observations (Cheng et al., 2022;Gamon et al., 2013;Myers-Smith et al., 2020;Wang et al., 2023).The snow effects on vegetation indices vary among vegetation types (Pierrat, Magney, et al., 2022) and among sites representing different biomes (Wang et al., 2023) due to varying species composition, canopy structure, and environmental conditions.However, the overall influence of snow on remote observations of vegetation photosynthesis has not always been clearly addressed and a systematic experimental investigation of how snow influences different vegetation indices across vegetation types has been lacking.
In this study, we designed an experiment using leaf-level gas exchange and reflectance, canopy reflectance, and digital images of potted trees at the University of Alberta, Canada to understand how snow cover influences canopy reflectance and vegetation indices of several evergreen and deciduous boreal forest tree species.We also tested the complementarity hypothesis (Gamon, 2015;Gamon et al., 2016)

| Experimental setup
We used six tree species that are commonly found in North American boreal forests, including three evergreen species (Picea mariana, Picea glauca, and Pinus banksiana) and three deciduous species (Populus tremuloides, Populus balsamifera, and Larix laricina).
We grew six monocultures and three mixed stands-M-1 (L.laricina, P. mariana), M-2 (P.tremuloides, P. banksiana), and M-3 (P.balsamifera, P. glauca)-on the rooftop of the Biological Sciences Building at the University of Alberta, Canada (Latitude: 53.529°, Longitude: −113.526°; Figure 1).Trees were potted in 2.83 L deep pots in spring 2015 and repotted into 6.23 L pots in April 2016 (during a period of leaf flush) for adequate moisture and nutrient availability and to avoid potential root restriction.During growing seasons, trees were watered daily and fertilized periodically to avoid water and nutrient deficits (Springer et al., 2017).An automated weather station was set up on the same rooftop where the trees were planted.Air temperature (S-THB-M002; Onset) and photosynthetic photon flux density (S-LIA-M003; Onset) were collected every minute and aggregated to 15-min averages.

| Gas exchange
Photosynthetic rate, expressed as net CO 2 assimilation, was measured using a portable gas exchange system (LI-6400; LI-COR).Leaves, or bundles of leaves in the case of conifers, were placed inside a 6 cm 2 leaf gas exchange chamber (6400-02B; LI-COR).The chamber monitored CO 2 assimilation rates under 1500 μmol photons m −2 s −1 to determine light-saturated photosynthetic rate (μmol CO 2 m −2 s −1 ).Light saturation was confirmed through light-response curves that were performed on each species in July 2016 (Springer, 2018).The reference CO 2 was set to 400 μmol mol −1 to match atmospheric concentrations.The chamber air flow was set to 400 μmol s −1 , with temperature and humidity set to match ambient conditions.Following a 1-3 min acclimation period after the leaf clip was set on a plant, five consecutive measurements were taken from a single branch or leaf on each plant from a total of five plants of each species.The photosynthetic rate of each species at the plot level was determined by averaging the measurements (N = 25) from all individuals of a given species at a single sampling interval.Leaf area for broadleaf species (P.tremuloides and P. balsamifera) was determined by the 6 cm 2 leaf chamber area; for needle-leaved species, leaf area was determined from the size (length and width measured by a caliper) and number of needles present in the gas chamber during sampling.
For each species, mature leaf tissues were sampled whenever possible; new tissues were sampled for P. mariana and P. glauca during the spring of 2016 when the elongation of new branches did not allow for sampling of mature sun-lit tissues.For evergreen species, gas exchange data were collected from December 2015 to February 2017 approximately every 2 weeks.For deciduous species, gas exchange data were sampled between April 2016 (noticeable leaf-out) and October 2016 (full leaf senescence) at the same frequency.Gas exchange measurements were only collected for plants in the six monoculture plots but not for the same species in the mixed stands due to the time required for multiple gas exchange measurements.

| Leaf reflectance
Leaf reflectance was taken using a single-channel spectrometer (Unispec; PP System) coupled with a needle leaf clip (UNI501; PP System) with an internal halogen light source.The needle leaf clip allowed a narrow field of view (0.6-mm-diameter), which enabled sampling of small, narrow needle leaves.We randomly sampled five leaves for each plant and five plants for each species, for a total of 25 samples per species.For each plant, leaf measurements were preceded with a dark and a white reference scan (Spectralon; Labsphere).Leaf reflectance was calculated by dividing each leaf measurement by a white reference scan after subtracting a dark spectrum from each measurement.The spectrometer has a nominal spectral range from 350 to 1100 nm with 2-3 nm sampling intervals and 10 nm full width at half maximum (FWHM).A linear interpolation was used to estimate reflectance at 1-nm intervals.Leaf reflectance was measured on the same dates when gas exchange data were collected.Leaf reflectance was only collected for plants in the six monoculture plots but not plants in the mixed stands.This leaf reflectance dataset (Wang et al., 2016a) is available at the EcoSIS Spectral Library (ecosis.org).

| Canopy reflectance
Canopy reflectance was taken using a dual-channel spectrometer (Unispec DC; PP System).A downward looking fiber (Field of View approximately 20°) was used to measure the reflected radiance from the target and an upward looking fiber that connected to a cosine header (UNI435; PP Systems) was used to measure the incoming irradiance.The detector measured irradiance and radiance from 350 to 1100 nm with 10 nm spectral resolution (FWHM) and 2-3 nm sampling intervals.Canopy reflectance were collected for all the plots, including monocultures and mixed stands.Five measurements were taken for each of the nine treatment plots at different locations (center, NW, NE, SW, and SE; Springer, 2018) ca. 1 m above the top of the canopy and averaged as the reflectance of the plot.The reflectance data were cross-calibrated using measurements of a white reference panel (Spectralon; LabSphere) prior to sampling of each plot.This dual-channel design allowed us to compensate for any slight changes in sky conditions, during the reflectance sampling period (Gamon et al., 2006).Canopy reflectance was measured between December 2015 and February 2017 on an approximately 2-week basis and was taken under sunny days and within 1 h of solar noon on each sampling day.We typically measured the canopy reflectance within 3 days of the leaf-level measurements.This canopy reflectance dataset (Wang et al., 2016b) is available at the EcoSIS Spectral Library (ecosis.org).
We calculated three vegetation indices that are commonly used to monitor photosynthetic phenology of boreal forests, including NDVI, PRI, and CCI at both leaf and canopy levels using directly measured leaf and canopy reflectance, respectively, as:

| Quantifying snow effects on seasonal changes in VIs
We quantified the snow effects on vegetation indices by calculating the seasonal change in vegetation indices with and without snow present (Figure 2).During snow-covered periods, temperatures were typically below or near 0°C and there was no measurable photosynthesis (Springer et al., 2017;Wong & Gamon, 2015a).Thus, we used two snow-free dates in the winter (January 29, 2016 and November 7, 2016) as "zero" snow cover points that provided "base line" values for vegetation indices without snow contamination or detectable photosynthetic activity.We then calculated the seasonal deltas (Δ seasonal ) of

| Quantifying snow effects on VIs by calculating snow percentage cover using digital images
To independently quantify the snow percentage cover, we took photographs of each monoculture plot from a height of approximately 1.5 m above the canopy on the same dates that we sampled the canopy reflectance using a digital camera (Nikon Coolpix 520; Nikon Canada Inc).We trained a binary support vector machine (SVM) classifier using the fitcsvm function in MATLAB (2022) to identify snow pixels from each image, because supervised classification can yield more accurate snow cover estimations than methods based on applying thresholds on RGB bands (Fedorov et al., 2016).To train the binary SVM model, we selected four images that had different levels of snow coverage on both deciduous and evergreen species and visually identified snow and non-snow cover pixels.We cropped each image by using 30% of the shorter edge (height or width) at each side from the center of the image to only use the center portion of each image to estimate the snow percent cover to avoid distortions at the edges (Liu & Pattey, 2010).We tested the SVM classification accuracy using 10-fold cross-validation.
The overall accuracy of snow identification using the binary SVM classification was 0.8.Most of the misclassified pixels occurred to pixels in shadows.We then applied the binary SVM model to the cropped images to calculate the snow percent cover of each image.

| Linking gas exchange to VIs at leaf and canopy levels
We used the leaf-level gas exchange data to establish the "biological baseline" for plants' seasonal photosynthetic rate variation.We then

| Snow effects on canopy reflectance and vegetation indices
Clear seasonal changes occurred in spectral reflectance due to both biological effects and snow (Figure 3).For evergreen species, the pattern of winter canopy reflectance measured during the snow free period exhibited a clear difference from the summer reflectance, shown as lower reflectance values in the green wavelengths and higher reflectance values in the red wavelengths (Figure 3).
Evergreen reflectance changes were also visible in the near-infrared (NIR), with higher NIR reflectance in the winter, particularly for P.
banksiana (Figure 3b).As expected, deciduous species showed more striking seasonal contrasts in reflectance, with higher visible reflectance in the red and blue regions, and lower NIR reflectance in winter (Figure 3d-f).For Populus spp., winter reflectance consisted of energy scattered from branches and soil background.L. laricina, a deciduous conifer, had some dry yellow needles remained on the branches over the winter period, leading to higher NIR reflectance than Populus spp., but also higher reflectance in the visible (VIS) wavelengths (particularly in the red region) than the evergreen species (Figure 3).Snow strongly increased canopy reflectance across all spectral bands, and this increase varied with wavelength, decreasing the difference between VIS and NIR bands (Figure 3).
Along with these reflectance changes, canopy-level vegetation indices of different species varied across season and were clearly influenced by snow (Figure 4).Snow shifted the vegetation indices from measurements taken during winter dormancy, with sharp decreases in NDVI values and abrupt increases in PRI and CCI values during periods of snow (Figure 4).For evergreens, snow was the predominant influence on NDVI (Figure 4a).In deciduous species, snow led to a "two-step" shape in the winter-spring NDVI time series for the deciduous species, clearly showing the two causes (snow melt and green up) of changes in the spring NDVI values (Figure 4b).
Compared to deciduous species, limited seasonal variations in NDVI occurred in evergreen species potentially due to biological effects (changing pigmentation or leaf area index).The slight drop in canopylevel vegetation indices in April 2016 was caused by repotting trees during a period of leaf flush (Figure 4).
Using "delta values" (defined in Figure 2), we explored the relative contribution of snow and biological effects on the reflectance indices (Figure 5).For all three vegetation indices, wintertime changes due to the variable presence of snow ("snow" delta values, N = 8) were much greater than normal winter variability ("winter" delta values, N = 3 but in different months) without snow (i.e., due to biological effects alone, which were minimal in winter).For evergreen species, snow caused a variation in NDVI that was four to six times larger than the typical seasonal variation ("seasonal" delta values) calculated when snow-affected period was excluded (Figure 5).
In deciduous and mixed stands, the effect of snow on NDVI was smaller than the seasonal variation due to biological effects ("seasonal" delta values), but still yielded changes that were 40%-60% as large as the biological effects.For mixed stands, the snow effect on NDVI (expressed as "snow" delta values) was roughly equivalent to that of the biological effects ("seasonal" delta values; Figure 5).
The variation in PRI caused by snow was at a level that was comparable to the seasonal PRI variation for evergreens and mixed stands but was slightly smaller than biological PRI variation ("seasonal" deltas) for deciduous stands (Figure 5).In the case of CCI, snow effects were approximately half those of seasonal variation due to biological effects ("seasonal" deltas), indicating that while this index was least affected by snow, the effect was still substantial.
Only minor changes were noticed in the vegetation indices during the winter snow-free period ("winter" deltas, Figure 5), indicating limited biological changes due to physiological or canopy structural changes during the winter.

| Correcting for snow effects on canopy level vegetation indices
The snow percent cover data obtained from the digital images during for deciduous species than evergreen species, but the slopes were not significantly different between the two vegetation types (p = .991).

| Linking gas exchange to VIs at leaf and canopy levels
For evergreen species, photosynthetic activation began in March, with photosynthetic rates increasing to summer maxima by mid-June (Figure 7a).Following the growing season maxima, photosynthetic rates gradually decreased, reaching near zero again by November.
During the spring transition, a large drop in photosynthetic rate was observed for the spruce species (P.mariana and P. glauca) that coincided with the sampling of newly emerged branches following budburst; this drop was not seen in P. banksiana as only mature needles of this species were sampled through the transition.
For deciduous species, rapid photosynthetic changes occurred in both spring and fall, with winter periods lacking foliage for sampling (Figure 7b).Spring transition of deciduous species started in April.L. laricina developed needles at the beginning of April and the Populus spp.unfolded their leaves within a relatively short period of time (between April 8 and April 18).In deciduous species, we missed sampling photosynthesis during the early spring transition following bud burst and early leaf development due to the relatively coarse sampling intervals and the difficulty of accurate measurements on small, emerging leaves (Figure 7b,d,f,h).Fall senescence of broadleaf species started at the end of August, indicated as a decrease in photosynthetic rate and leaf yellowing due to loss of chlorophyll (Figure 7b).At the end of September, few green deciduous leaves remained, which can also be confirmed by near-zero photosynthetic rates at this time (Figure 7b).Broadleaf species exhibited higher peak  For evergreen species, stronger relationships emerged between photosynthetic rate and PRI and CCI than between photosynthetic rate and NDVI (Figure 8; Table 1).The low photosynthetic rate obtained for the new needles of the two Picea species in late April to early May affected the overall relationship between photosynthetic rate and vegetation indices, leading to outliers (Figure 8).Removing measurements of the Picea species for four outliers on April 20 and May 4 enhanced the relationship between photosynthetic rate and PRI and CCI (Figure 8; Table 1).
Relationship between photosynthetic rate and NDVI changed slightly (Figure 8; Table 1).
For deciduous species, and especially for the two Populus species, limited variation in vegetation indices occurred during the summer growing season, but all vegetation indices yielded significant correlations with photosynthetic rate (Figure 7).For deciduous species, NDVI exhibited the closest relationship (largest R 2 ) with photosynthetic rate among the three vegetation indices tested (Table 1).This was the reverse of the pattern in the evergreens, where NDVI had markedly less variation with changing photosynthetic rates than PRI and CCI (Figure 8), illustrating the complementary nature of these indices in assessing photosynthetic activity for different vegetation types (Gamon, 2015;Gamon et al., 2016;Garbulsky et al., 2011).photosynthetic rate and PRI (expressed as higher R 2 ), while weakening the relationship between photosynthetic rate and NDVI, both of which represented artifacts of snow rather than actual biological effects.
Snow cover had relatively little effect on the CCI-photosynthesis relationship (Figure 9; Table 1).We note that because most gas exchange measurements were limited to the growing season (Figure 7), a period with minimal snow cover, the effects of snow on these relationships were minimized in this figure .Because no leaf-level measurements were collected for deciduous species during the snow-affected period (before bud burst and early leaf development), we evaluated the snow correction on canopy-level vegetation indices only for evergreen species in the following analysis.
Similar to the relationships between photosynthetic rate and vegetation indices at leaf level (Figure 8), the low photosynthetic rate of the two spruce species in late April and early May (Figure 7) also weakened the relationships between photosynthetic rate and vegetation indices at canopy scale (Figure 9).enhanced the relationship between photosynthetic rate and canopy PRI and CCI (Figure 9; Table 1).
To evaluate snow effects on scaling, we further compared the leaf-level vegetation indices to canopy-level vegetation indices directly (Figure 10).In all cases, snow correction improved the relationship between leaf-level and canopy-level vegetation indices (increased R 2 , Figure 10).The agreement between leaf and canopy measurements (and indicator of "scalability") was strongest for PRI and CCI and weakest for NDVI (Figure 10).Leaf-level NDVI exhibited a much larger range than canopy-level NDVI, causing the NDVI plot to deviate markedly from the 1:1 line (Figure 10).

| Snow effects on vegetation indices
The magnitude and patterns of changes in reflectance spectra in winter are largely influenced by vegetation canopy structure that impacts the retention and detection of snow in canopies (Pomeroy et al., 2002), and this varied with species and functional type.Among evergreens, different taxa (e.g., Picea spp. vs. Pinus spp.) vary in their ability to retain Relationships between photosynthetic rate and leaf-level vegetation indices for evergreen (a, c, e) and deciduous (b, d, f) species.Regressions (solid lines) were applied to combined evergreen and deciduous species, respectively.The outliers (Picea spp. on April 20 and May 4) are labeled as open symbols (panels a, c, e).Removing these outliers enhanced the relationship (dash line) between photosynthetic rate and PRI and CCI, as indicated in parentheses on the figures.Regression results are summarized as Table 1.4).In natural landscape, deciduous branches may collapse and be buried beneath the snow cover, particularly for small trees or shrubs, influencing the overall optical signal of deciduous forests (Marsh et al., 2010;Ménard et al., 2014;Ray & Bret-Harte, 2022).Additionally, shrub branches can extend above the snowpack and alter albedo, surface temperature, and soil thermal regimes, which further accelerate snow melt (Myers-Smith & Hik, 2013;Pomeroy et al., 2006).These varying snow effects among evergreen and deciduous species can lead to variations in optical properties at broader scales (e.g., Cheng et al., 2022;Pierrat et al., 2021;Pierrat, Magney, et al., 2022), which would vary by species and vegetation type, in part due to varying snow retention.tude of changes that coincide with bud burst and leaf flush for deciduous species and can far exceed that due to biological effects (due to canopy development or altered pigmentation) in evergreens (Figure 11).
Together, these effects of snow on NDVI can lead to substantial errors in phenology estimation.Being inversely impacted by snow (relative to NDVI), PRI and CCI increases due to snow in wintertime to levels that are comparable to growing season maxima would suggest a sudden "turning on" of photosynthesis in the middle of winter for evergreens, and even deciduous species with no foliage.These results reveal clear and distinct artifacts of snow melt on vegetation indices apart from any direct effect of vegetation structure (e.g., leaf area index) or physiology (e.g., pigmentation and associated photosynthetic activity).

| Snow effects on estimation of vegetation photosynthetic activity and productivity
Remote sensing data can be used to quantify carbon fluxes through direct empirical relationships between vegetation indices and GPP, light use efficiency concepts, terrestrial biosphere models,  1.
or machine learning methods (reviewed in Ryu et al., 2019;Xiao et al., 2019).However, all reflectance-based vegetation indices commonly used for monitoring terrestrial ecosystem photosynthesis are sensitive to snow cover, especially for the high latitude and altitude areas, which endure long periods of winter snow (Eklundh et al., 2011).Compared to snow-affected data, NDVI exhibited limited seasonal variation for evergreen species, as expected based on previous studies quantifying NDVI phenology for evergreens in the absence of snow (e.g., Gamon et al., 1995).The small seasonal variation in NDVI for the evergreen species (Figure 5) is undoubtedly due to NDVI's sensitivity to green canopy structure (which has little seasonal change in evergreens).However, in northern latitude evergreen forests, greenness-based vegetation indices such as NDVI may exhibit strong relationships with GPP derived from eddy covariance measurements.Our results suggest that the strong relationship between GPP and greenness-based vegetation indices appears in part to be a "false signal" due to the confounding effects of snow and the coincidence of snow melt and snow fall with periods of photosynthetic activation and deactivation (Wang et al., 2023).The large errors in reflectance-based estimates of GPP during seasonal transitions (Keenan et al., 2012;Richardson et al., 2012;Schaefer et al., 2012), are likely partly attributable to such confounding snow effects.
The apparent seasonal response of greenness-based indices in evergreen forests can also be due to the significant contribution of   rooftop experiment, indicating that snow is likely to have even greater confounding effects in satellite observations, as has recently been observed in a parallel study of MODIS satellite data (Wang et al., 2023).

| Caveats and future work
In addition to the possible differences between our experimental conditions and natural landscapes, several other issues are worth noting when considering the application of snow corrections.Despite the linear relationships between snow cover and vegetation indices reported here (Figure 8), additional attention would likely be required when exploring these precise snow cover effects in natural settings.
Since snow cover estimates from imagery will likely vary with view angle, angular effects should also be considered in studies designed to assess snow impacts on remotely sensed signals.For example, images from the PhenoCam Network, are typically taken from oblique angles, which yields a different snow cover estimate than nadir views due to enhanced visibility of vegetation at oblique angles (Brown et al., 2016;Toomey et al., 2015).Snow estimation using images taken at oblique angles typically requires a transformation to orthophotos.The accuracy of snow cover estimation is subject to errors in the transformation and illumination conditions (Fedorov et al., 2016;Hinkler et al., 2002).Combining imagery collected by drone (Belmonte et al., 2021) or CubeSat (Cannistra et al., 2021;John et al., 2022) with machine learning methods may offer new approaches to deriving high spatial resolution snow coverage products.The accuracy of these methods is sensitive to forest canopy structure and landscape topography (Belmonte et al., 2021;Cannistra et al., 2021).In any case, the advent of global network of digital cameras and CubeSat observations hold promising potential for snow estimation, which could further advance the forest productivity estimation using remote sensing by eliminating snow's confounding effects.

| CON CLUS IONS
Vegetation indices are widely used in monitoring terrestrial plant photosynthesis and GPP across spatial and temporal scales.We For accurate GPP estimation, these large snow cover effects on vegetation indices and implications for relationships between vegetation indices and productivity should be addressed in any large-scale satellite studies of GPP using these vegetation indices.
In particular, the large and contrasting effects of snow on different vegetation indices and for different vegetation types suggest that many models of GPP driven from these indices could have considerable errors, particularly during the critical period of spring snow melt and photosynthetic activation, as has previously been noted (e.g., Richardson et al., 2012;Rogers et al., 2022).Failure to remove the snow effects could lead to incorrect estimations of GPP using either empirical methods between GPP and vegetation indices (Wang et al., 2023) or using machine learning methods (Cheng et al., 2022;Pierrat, Bortnik, et al., 2022;Pierrat, Magney, et al., 2022).Rather than simple linear effects, the influence of , which argues that different vegetation indices (e.g., NDVI and PRI) having contrasting sensitivity to canopy structure and physiology.Thus, different vegetation indices provide complementary information about plant photosynthesis.In particular, we hypothesize that (1) snow has substantial effects on vegetation indices and confounds the estimations of photosynthetic properties based on optical properties (further illustrated in Figure 2 below) and (2) effects of snow on vegetation indices vary among vegetation indices and between deciduous and evergreen species.Because the combined effects of snow and vegetation type are difficult to assess from satellite measurements at large spatial scales, an initial clarification of these effects in an experimental setting is a necessary first step toward improved global assessment of primary productivity based on vegetation indices.
vegetation indices as difference between the maximum and minimum extremes for the snow-free data collected in 2016, which captured the seasonal biological effect on vegetation indices apart from any snow influence.The snow deltas (Δ snow ), captured the direct effects of snow on vegetation indices, were calculated between the maximum and minimum extremes for all the data collected in winter (November 2015-March 2016).We also calculated winter deltas (Δ winter ) between the maximum and minimum extremes for all the snow-free data collected in winter (November 2015-March 2016), which captured the winter variation in vegetation indices apart from any snow effects and was typically close to zero.(1) NDVI = 800 − 630 800 + 630 .Conceptual figure indicating how snow influences the vegetation indices for evergreen (a, c, e) and deciduous (b, d, f) tree canopies.The blue line denotes the hypothetical seasonal vegetation index trend without snow, while the red line denotes the snowaffected vegetation index in the winter period.The calculations of seasonal (Δ seasonal ), winter (Δ winter ), and snow (Δ snow ) deltas are illustrated in panel (b), which overlays hypothetical data points on the hypothetical trends.
related vegetation indices at both leaf level and canopy level to the gas exchange data to investigate how different vegetation indices provide complementary information about plant photosynthetic activities of deciduous and evergreen species at different scales.The comparison between gas exchange data and vegetation indices at leaf and canopy levels also enabled us to evaluate how snow influenced the canopylevel vegetation indices for both evergreen and deciduous species.We compared leaf-level vegetation indices to canopy-level vegetation indices to test how the different vegetation indices can be scaled from leaf to canopy scale.To further examine the snow effects on phenology of the vegetation indices, we used asymmetric Gaussian functions(Jönsson & Eklundh, 2004) to fit the snow-corrected data to estimate the theoretical seasonal snow-free vegetation index variations.
winter period also allowed us to quantify how snow-affected vegetation indices differently among species and functional types apart from any seasonal effects due to changing photosynthetic activity (Figure6).Two snow free dates in the winter (January 29, 2016 and November 07, 2016) were used as "zero" snow cover point that provided baseline values for vegetation indices without snow contamination.Overall, stronger relationship (larger R 2 ) between vegetation indices and snow cover was found in deciduous than evergreen, while sharper changes (steeper slopes, p < .001) in NDVI and PRI with increasing snow percent cover were found in evergreen species than deciduous species (Figure6a,b).The different slopes for PRI and NDVI indicated strong interactive effects of snow and vegetation type on these indices.For CCI, a larger intercept (p < .001)occurred photosynthetic rates in the summer growing season and faster decreases in photosynthetic rate during the fall transition than L. laricina, a deciduous conifer.The ability to use vegetation indices to track the seasonal photosynthetic rate varied among species and functional types (Figure 7).For evergreen species, clear seasonal patterns matching the trajectory of photosynthetic rates were exhibited by leaflevel PRI and CCI (Figure 7e,g).A slight seasonal variation of NDVI occurred in P. banksiana (pine) but not in the two Picea (spruce) species.During the wintertime, limited changes occurred in all three vegetation indices of the evergreens (Figure 7c,e,g).For deciduous species, clear seasonal patterns were found for all three vegetation indices tested, expressed as high values in the summer growing season and abrupt drops during the fall senescence (Figure 7d,f,h).L. laricina, a deciduous conifer, exhibited a more gradual change in vegetation indices during senescence than the broadleaf deciduous species.
Removing measurements of two spruce species on April 20 and May 4 (during leaf flush) F I G U R E 5 Delta vegetation index values, including NDVI (a), PRI (b), and CCI (c), for evergreen (left), deciduous (center), and mixed stands (right) on different temporal scales.Seasonal deltas (green) were taken from the maximum and minimum extremes from 2016, and represent the changes in index values due to biological effects during the growing season.Winter deltas (blue) were taken from the maximum and minimum extremes from November 2015 to March 2016 and represent the change in index values due to biological effects in the winter.Seasonal and winter deltas excluded data collected with snow present on canopies.Snow deltas (yellow) represent the changes in index values due to snow alone, and were taken from the maximum and minimum extremes from November 2015 to March 2016, and compare index values before and after snowfall.Calculation of these delta values are further explained in Figure 2.F I G U R E 6 Vegetation index (NDVI: a, PRI: b, and CCI: c) variation with snow percent coverage by canopy type (deciduous, evergreen, or mixed, represented as different colors).Snow percent coverage was estimated using digital images taken on the same days when canopy reflectance was sampled.The mixed plot data (blue points) are illustrated but not included in regressions.p-value < .001for all the regressions.To isolate snow effects from biological effects, only winter data were used in this figure, assuming little change in physiology or structure during winter.
Figure11).For deciduous species, snow correction removed the initial increase exhibited in spring NDVI values caused by snow melt but not the increase due to biological effects, because snow melt happened prior to

TA B L E 1
photosynthetic activity for Picea spp. to a greater degree.On the other hand, the absence of winter foliage for deciduous species leads to the combination of snow-covered soil with bare stems and branches as influences on reflectance and vegetation indices (Figures1 and 4).In natu- Due largely to their different formulation using dissimilar spectral bands, contrasting vegetation indices were differentially affected by snow cover, indicating contradictory and confounding effects of snow on these indices (Figure 4; Table 2).These contradictory effects include exaggerated seasonal change (in the case of NDVI) and reduced seasonal change (in the case of PRI and CCI) and indicate artifacts in the sense they can easily lead to false conclusions regarding vegetation greenness, pigmentation, and photosynthetic activity.In accordance with previous studies (Eklundh et al., 2011; Gamon et al., 2013; Springer et al., 2017) the fastest NDVI increases happened during the snow melt period in spring (Figure 4) when photosynthetic activity remained low.Additionally, without accounting for snow cover, NDVI can indicate a sharp decline in apparent "greenness" for deciduous species well after fall senescence and might incorrectly suggest a rapid loss or sudden emergence of green foliage in the winter for evergreens.Changes in NDVI due to differential snow cover in the early transition from winter to spring rival the magni- deciduous plants (including forest understory species) to surface reflectance.It has been reported that the forest understory, often including a large deciduous component, can account for about 61% total GPP in an Alaska black spruce forest(Ikawa et al., 2015).In such cases, the seasonal change in deciduous components can provide a considerable contribution to the overall remotely sensed optical signals (as illustrated by the mixed stands in this study), leading to strong relationships between NDVI-based indices and GPP in the boreal forests, despite snow artifacts.However, the confounding effects of snow on vegetation indices, which are clearly different for evergreen and deciduous species, are likely to complicate assessment of seasonal change from remote sensing in such mixed stands.
This experimental study allowed us to quantify snow effects on vegetation indices of deciduous and evergreen species and to separate snow artifacts from actual biological effects.While this study demonstrates the principles of snow effects on canopy reflectance and vegetation indices, the details would vary between this study and natural landscapes due in part to difference in species composition, canopy structure, and landscape effects (e.g., variation in topography and forest density).Our photosynthesis measurements were limited to periods of the growing season F I G U R E 1 0 Snow correction enhanced the relationships between leaf-level and canopy-level vegetation indices (NDVI: a & b, PRI: c & d, CCI: e & f) for evergreen species.Canopy data with snow visible in the scene are labeled using stars.Solid line denotes the regression line and dash line indicates the 1:1 reference line.Regressions were fitted combining data from all three evergreen species.when leaves were available, which tended to minimize the confounding effects of snow on GPP estimation due to the lack of sampling in the winter and early spring for deciduous trees.In this experiment, the tree canopies tended to have less snow cover than comparable stands in natural landscapes, in part due to the sun-exposed and slightly warmer rooftop conditions, reducing the snow influence relative to natural landscapes.The relatively tight spacing of plants in synthetic plots likely influenced the retention of snow and likely altered the confounding effect of snow relative to that of natural forests.Larger vegetation has different canopy structure compared to younger trees, which could further influence the retention of snow and alter the exact snow cover impacts in natural landscapes.For these reasons, vegetation indices in this small-scale experiment were likely less sensitive to snow cover than would be apparent in most natural landscapes indicating that this canopy-level experiment probably presented a conservative estimate of the confounding effects of snow.Natural landscapes tend to have more persistent snow cover than observed in this F I G U R E 11 Effect of snow and snow removal on phenology of remotely sensed vegetation indices.Examples of one evergreen (Picea mariana, a, d, g), one deciduous (Populus tremuloides, b, e, h), and one mixed stand (Populus balsamifera & Picea glauca, c, f, i) are shown.Open points at the beginning and end season indicate snow-affected measurements, identified from digital images.Asterisks indicate snowcorrected values, and solid points indicate snow-free values.Instead of using data from calendar years, we expressed data as "days from December 15, 2015" to include more snow-affected period, due to a relatively short snow season in 2016.TA B L E 2 Snow effects on vegetation index values, season changes, and photosynthetic rate-vegetation index relationships.
found that vegetation indices from reflectance were all affected by snow cover change, but in different directions and to different degrees.These effects of snow cover on vegetation indices also varied with functional types (evergreen and deciduous) and species (e.g., Picea sp. vs. Pinus sp.).Snow cover decreased NDVI-based vegetation indices and snow melt in the spring, exaggerating the correlation between NDVI-based vegetation indices and GPP, while the opposite effects were seen for CCI and PRI.These results indicate that large errors can result in GPP estimation for large regions of the world exposed to seasonally varying snow cover, including the boreal and high-altitude regions.
snow on vegetation indices are complicated by species, vegetation type, and view angle effects, making simple and accurate snow corrections challenging.These findings of snow and vegetation type effects on vegetation indices have large implications for GPP estimation in boreal forests, where snow cover, species composition, and photosynthetic activity are all changing.While we demonstrate that empirical corrections using snow coverage estimation based on digital imagery can reduce snow artifacts to better reveal the underlying biological effects, the confounding snow effects on reflectance-based vegetation indices and relationships between GPP and these vegetation indices deserve further study if we are to develop and apply reliable methods of accounting for snow cover.This is a particular need for satellite studies of boreal forests, considering the critical role the boreal forests play in regulating the global carbon cycle and its sensitivity to global climate change.Proper correction for snow effects could lead to improved accuracy of GPP assessments for large regions of the world prone to changing snow cover.
Note: N denotes the sample size.p value < .001for all regressions.At canopy level, regressions were only calculated for evergreen species to illustrate snow effects (before and after snow correction), because no photosynthesis measurements were collected for deciduous species during the snow-affected period.R 2 with Picea spp.outliers on April 20 and May 4 removed are shown in parentheses.