The environmental drivers of annual variation in forest greenness are variable in the northern Intermountain West, USA

Peak snowpack in the western USA has decreased 21% since the mid-20th century, and these trends are forecasted to continue over the next century. In water-limited systems, forest productivity during the growing season is assumed to be linked with snowpack during the previous winter, but this linkage has proven not to be universal in studies across varying ecosystems. We compared peak Normalized Difference Vegetation Index (NDVI, a proxy for vegetation productivity) at relatively fine scales (30-m pixel) to point measurements of snow water equivalent (SWE) at 169 locations scattered among the northern Intermountain West for the 25-yr period between 1993 and 2017. We further compared NDVI to additional environmental variables including topographical variables (elevation, aspect, and slope), and climatic variables (Palmer Drought Severity Index, snowmelt date, and snowmelt length). We hypothesized that peak SWE would be a strong predictor of peak NDVI, but the strength of the SWE–NDVI relationships would vary spatially in our topographically complex study area. We observed weak relationships between snowpack and vegetation productivity, with only 10% of the locations displaying statistically significant correlations between annual peak SWE and peak NDVI. Furthermore, we determined that the response of NDVI to snowpack and environmental variables differed among subregions of our study area. For example, we observed strong SWE–NDVI associations in the continental-most portion of our study area and weak correlations in other subregions of our study area. Our results highlight that in regions characterized by complex topography, (1) environmental drivers of NDVI are highly variable and (2) fine-resolution studies more accurately capture the spatial heterogeneity of NDVI compared to courseresolution studies. However, one consequence of fine-resolution studies may be that microclimate effects at each site swamp the signal of predictor-ecological response relationships, whereas these relationships may be more discernable in larger scale remote sensing studies.


INTRODUCTION
Since the mid-20th century, peak snowpack (often quantified as snow water equivalent on April 1) in the western USA has decreased 21%, totaling a 36-km 3 loss in snowpack (Mote et al. 2018). These declines have been linked with regional climate change (Gutzler and Robbins 2011, IPCC 2014, Mote et al. 2018, including increases in wintertime temperatures and decreases in the fraction of winter precipitation as snowfall (Stewart et al. 2004, Knowles et al. 2006, Mote 2006. These trends are forecasted to continue over the next century (Cox et al. 2000, Nogués-Bravo et al. 2007, IPCC 2014. In water-limited systems, forest productivity during the growing season has been linked with wintertime snowpack (Grippa et al. 2005, Peng et al. 2010, Trujillo et al. 2012, Trujillo and Molotch 2014, Knowles et al. 2018, Wang et al. 2018, because the amount and duration of available water during the growing season is related to winter snowpack. However, increasing winter and springtime temperatures may also affect snowmelt dynamics by shifting snowmelt seasonality, reducing high snowmelt rates, and thus altering soil moisture availability (Barnett et al. 2005, Barnhart et al. 2016, Musselman et al. 2017. Greater winter snowpack can increase vegetation productivity during the summer due to greater amount and/or longer duration of available soil moisture during the growth season (Trujillo et al. 2012, Winchell et al. 2016, Knowles et al. 2018. For example, greater net ecosystem exchange (carbon-dioxide uptake, g CO 2 Ám −2 Ás −1 ) occurred in years with greater peak snowpacks with later melt dates in a subalpine forest (Winchell et al. 2016). However, a large and/or late-melting snowpack also has the potential to reduce growing season length, which can reduce overall summer productivity and produce negative relationships between productivity and snowpack. Day et al. (1989) reported that photosynthesis was suppressed in several subalpine conifer species by snowpack that persisted through the early summer, due to cold soil temperature limitations on root growth and nutrient absorption. It is even possible that there will be little to no statistically significant relationships between summertime vegetation productivity and wintertime snowpack exist (Peng et al. 2010, Kariyeva and van Leeuwen 2011, Wang et al. 2018. This seems to be especially so in regions with complex terrain, where topographical variables such as elevation and slope strongly modify the influences of moisture availability, radiation, and temperature on NDVI at small spatial scales, swamping the influence of wintertime snowpack on summertime productivity (Kariyeva and van Leeuwen 2011). Clearly, more understanding of how vegetation responds to changes in snowpack and other environmental factors is needed in order to predict how ecosystems may react to changing snow conditions in future climate scenarios. The wide variability in the strength and directionality (positive/negative) of the relationship between snowpack and plant productivity in the existing body of observational literature (based on satellite-derived observations) could be due to several reasons. One is that these studies occurred in different geographic regions, with varying climates, vegetation composition, and plant canopy characteristics (summarized in Table 1). Another explanation for the variability in snowpack-productivity relationships could be due to differences in methodology. Most of the existing studies have used relatively coarse resolution imagery (e.g., 1 km, 8 km, or 1°pixels) to investigate snowpack-vegetation relationships over macroscales (10-1000 km 2 ) and mesoscales (100 m-10 km 2 ). At local scales (nominally, <100 m 2 ), snowpack and other environmental parameters may vary due to micro-scale variation in elevation, aspect, topography, vegetation structure, and prevailing wind patterns (Pomeroy et al. 2002). For example, canopy structure can have micro-scale effects on wind speeds and solar irradiance, affecting snowfall interception, wind scour, and sublimation (Stegman 1997, Pomeroy et al. 2002, Winkler et al. 2005, Woods et al. 2006, Molotch et al. 2007, Veatch et al. 2009, Tennant et al. 2017. Finer scale studies that utilize satellite imagery to quantify snowpack-NDVI relationships are lacking in the literature, and such studies may be able to improve our understanding of snowpack-productivity relationships, especially in regions with complex terrain. For example, in three regions in Central Asia of contrasting topography, NDVI was related to different sets of environmental drivers (Kariyeva and van Leeuwen 2011). Also, forest ecosystem processes can operate differently at different spatial scales. For example, photosynthetic light response (CO 2 -flux vs. sunlight intensity), and environmental regulators of this relationship, varies from leaf to crown to canopy scales (Arkebauer et al. 2009).
Herein, we present a study in which we compared vegetation productivity at relatively fine (30-m pixel) scales to point values of SWE data within study plots scattered among the northern Intermountain West (nIMW) of the USA. We made comparisons during the 25-yr period between 1993 and 2017 using data from satellite platforms for NDVI (Normalized Difference Vegetation Index), which we used as a proxy for forest productivity (Rouse et al. 1974, Townshend and Justice 1986, Grippa et al. 2005. We then compared these data to annual winter/spring peak SWE measurements (as well as other environmental variables) collected at 169 Natural Resources Conservation Service (NRCS) SNOw TELemetry (SNOTEL) stations across the region. Based on previous studies in regions near our study region (Trujillo et al. 2012, Knowles et al. 2018, we hypothesized (H1) that peak SWE would be a strong predictor of peak NDVI. Alternatively, we hypothesized (H1a) that peak NDVI would be associated with environmental parameters other than peak snowpack and that SWE-parameter relationships would vary across subregions (e.g., Kariyeva and van Leeuwen 2011) and elevation (e.g., Trujillo et al. 2012) in our topographically complex study area.

Study area
The area of study was the northern Rocky Mountain/Intermountain West region, approximately 42°N, 109°W to 49°N, 117°W (Fig. 1). This region is characterized by high variability in latitude (42°-49°N), topography and elevation (1536-3078 m), climate (MAP 20-250 cm, MAT 10°-18°C), and vegetation composition and cover (open canopy scrub forest to high-density conifer forest). Such patchiness of environmental and ecological variables results in high variability of snow deposition and redeposition (Nayak et al. 2010, Tennant et al. 2017. Because of this heterogeneity, we subdivided our study area into three terrestrial ecoregions (U.S. EPA; http://www.e pa.gov/wed/pages/ecoregions.htm) for analysis -the Northern Rockies/Columbia Mountains (hereafter referred to as Northern Rockies), the Middle Rockies, and the Idaho Batholith. We confined the Northern Rockies ecoregion to the U.S. border, because this was the extent of the snow station network used in our study. For some analyses, we subdivided our study area into 13 climate divisions (https://www.ncdc. noaa.gov/monitoring-references/maps/us-clima te-divisions.php) to further constrain heterogeneity of environmental characteristics among stations.

Snow metrics
Within the study area, there are 169 Natural Resources Conservation Service (NRCS) SNO-TEL stations, which have been recording data on snow depth and water content since the late 1970s. SNOTEL stations are located in areas thought to represent snowpack conditions specific to that elevation and aspect, and previous studies show strong correlations with surrounding snow course data (r 2 = 0.90; Dressler et al. 2006). We downloaded daily snow water equivalent (SWE) values from the SNOTEL network (https://wcc.sc.egov.usda.gov/reportGenerator/) for every station for all days, for each year SWE data were available (years 1993-2017). Some stations did not have a snow pillow installed until several years into our study. We determined peak SWE for each station for each year manually, and a lag SWE using a running average for two-, three-, and four-year SWE. It has been shown that lag responses to severe drought are not detectable beyond four years in forests in waterlimited systems (Anderegg et al. 2015). For any given year, in the instances of unreliable SWE data, we did not use those data for those stations that year. This never exceeded two stations for any given year, and never occurred for an individual station more than two times. Additionally, we estimated the snow aridity index (SAI) for each year using the sum of 1 April -31 August potential evaporation (PE) divided by annual max SWE from each station (Knowles et al. 2017(Knowles et al. , 2018. We calculated PE using the Blaney-Criddle formula (Blaney and Criddle 1962) to preserve site-specific PE, because sufficient meteorological data from each station were not available for the Penman-Monteith equation. Other variables determined from the SNOTEL data set included date of peak SWE, date of snowmelt (i.e., ablation date, defined as the last continuous five days in which daily SWE > 1 mm is observed; Wang et al. 2018), air temperature (which was further used to calculate growing degree days; GDD), and the geographical location of each station (latitude, longitude, and elevation).

Normalized Difference Vegetation Index (NDVI) data
We retrieved satellite imagery from the United States Geological Survey (USGS) Landsat 5, 7, and 8 satellites as analysis ready data (ARD; available at https://earthexplorer.usgs.gov/), and downloaded all images overlapping our study area spanning May through August, from the years 1993-2017. The ARD dataset comes atmospherically corrected and was downloaded as surface reflectance values. This dataset has a spatial resolution of 30 m for all platforms, and a return interval of 16 days. Using a 10 × 10 grid of 30-m cells (to simulate Landsat platform resolution) around each of the 169 SNOTEL stations within the study area (see Fig. 1), we extracted Normalized Difference Vegetation Index (NDVI) data in the local area around each station by performing a zonal statistics function for each of the grids using Python 3.6.6. This resulted in 29,160 v www.esajournals.org satellite images used for analysis. We excluded areas of our grid with NDVI values less than 0.2 to eliminate the influence of rock, bare soil, water, and snow (Bannari et al. 1995). Additionally, we removed water years (1 October -30 September) containing less than four Landsat images with sufficient pixel quality for the entire 10 × 10 grid for each station (which might occur due to cloud cover, e.g., Trujillo et al. 2012). We then determined peak NDVI values for all 100 grid cells (which might have occurred at slightly different dates for each grid cell) around each station and averaged them as a single peak NDVI value for that year (see example code in Appendix S1).
Due to a narrower range in the red and nearinfrared bands in Landsat 8 (relative to previous Landsat platforms), we converted surface NDVI for all Landsat 8 data (years 2013-2017) using an ordinary least squares equation using Eq. 1 (Roy et al. 2016): where ETM + is the output value equivalent to Landsat 7 Enhanced Thematic Mapper (ETM+) and OLI is the Landsat 8 Operational Land Imager (OLI) input NDVI values. Without this adjustment, we found that Landsat 8 NDVI values were noticeably greater (5-6%) than previous Landsat satellites (Appendix S1: Fig. S1). This transformation allowed for the continuity of NDVI values across multiple satellite platforms throughout our time series analysis.

Environmental and vegetative characteristics of the sites
We extracted the aspect and degree of slope for each station from aspect and slope layers created in ArcGIS Pro, respectively. We created these layers by stitching together 63 Digital Elevation Models to span our study region from the Advanced Spaceborne Thermal Emission and Reflectance Radiometer (ASTER, 30-m resolution; available at https://earthexplorer.usgs.gov/). These variables were of particular importance to us, as heterogeneity in topography can explain inconsistencies in snow depth (Tennant et al. 2017. Previous studies in the western USA have suggested that vegetation may have differing sensitivities to snowpack based on elevation (Trujillo et al. 2012). To explore this possibility in our study, we separated each of the three terrestrial ecoregions into three elevation classes, the ranges of which varied by ecoregion. This is because the elevation ranges of the mountains varied greatly among the ecoregions. To do this, we performed a variety of breaking methods based on various gradients in ecological (tree species composition), meteorological (precipitation), and statistical variables. In the end, we found that the simplest approach was a natural breaks (jenks) test in Arc-GIS, which we also determined generally agreed with perceived breaks in ecological and meteorological variables. The resulting elevation zones (low, medium, and high) of each ecoregion can be found in Table 2.
We downloaded land-cover classes from the National Gap Analysis Project (GAP) provided and distributed by the US Geological Survey (available at https://gapanalysis.usgs.gov/gapla ndcover/). Land-cover classes are often utilized to compare the productivity response of dissimilar vegetation types (Peng et al. 2010, Trujillo et al. 2012). In our case, the 169 SNOTEL stations were located within 16 vegetation cover classes. We lumped these 16 classes into five broader vegetation categories for easier analysis. These were as follows: Deciduous Forest, Grassland, Shrub/Sage, Mixed Conifer Forest, and Riparian Forest. By analyzing multiple vegetation types, we aimed to gain a broader understanding of montane ecosystem productivity, whereas other similar studies in mountain systems typically focused on conifer forested areas.
We downloaded Palmer Drought Severity Index (PDSI) data as a monthly value (May--August) for all WYs (1993WYs ( -2017 from the National Climatic Data Center (NCDC; https:// www7.ncdc.noaa.gov/CDO/CDODivisionalSelec Note: The ranges varied among each terrestrial ecoregion and were determined using natural breaks methods, which generally corresponded to ecological transitions. v www.esajournals.org t.jsp#). PDSI is a useful drought index because it accounts for water balance components such as evapotranspiration and runoff, rather than just precipitation inputs alone (Trenberth et al. 2014). It does this by using a assessing a combination of precipitation, air temperature, and soil moisture measurements (Dai 2013, Dai et al. 2018.

Statistical analyses
We first tested peak NDVI and peak SWE for spatial autocorrelation using Moran's I test in ArcGIS Pro. Because we found significant autocorrelation, we conducted some analyses using the three EPA ecoregions (to maximize generalizability when possible) and other analyses (independently from the first set of analyses) using the smaller-sized 13 NCDC climate divisions (to increase precision when needed). We also tested for serial autocorrelation within each grouping (ecoregion or climate division) and for each station individually using the Box-Jenkins method using JMP Pro v 14 Statistical Software (SAS Institute, Cary, North Carolina, USA). No serial autocorrelation was detected at any scale.
To evaluate snowpack-NDVI relationships, we performed ordinary linear regression (OLR) analysis between annual peak NDVI against annual peak SWE, lag SWE, and SAI using both the annual 25-yr means for all 25 yr at each of the 169 SNOTEL sites and the 25-yr means for all combined sites. The relationships between NDVI, SWE, lag SWE, and SAI were then stratified by elevation zones (low, mid, and high) within each of the three ecoregions because a previous study showed that NDVI was more responsive to changes in SWE in middle elevation zones (Trujillo et al. 2012). Additionally, we investigated the effect of snowpack on forest productivity using OLR analysis stratified among different vegetation classes, slope, and aspect.
To determine what environmental factors other than peak SWE were associated with NDVI on an annual basis, we implemented an exploratory approach to investigate different geographic regions, ecoregions, climate divisions, elevations, latitudinal gradients, and community compositions. To increase replication and empirical rigor, we performed multiple regression analysis within each of the 13 climate divisions of our study region using JMP. To determine which of the 14 original predictor variables were best suited for analysis, we first performed a data reduction approach by stepwise removal of all possible multi-correlated predictors by calculating the variance inflation factors (VIF) using R CRAN package (v 3.4.2, R Core Team 2017), and excluding variables exhibiting VIF values >5 (Fox and Monette 1992). We did not include geographic variables (e.g., latitude), as these were accounted for by our use of the climate zones. After stepwise removal of multi-correlated parameters, we were left with five predictor variables: SWE, snowmelt date (SMD), degree of slope, elevation, and PDSI. Within each of the 13 climate divisions, we then used these five predictor variables to determine (1) the significance (based on individual P-values) of each of the individual variables to NDVI, and (2) the best model for predicting NDVI, which was chosen based on the lowest root mean square error (RMSE; this model may not have included all five variables).

Temporal and spatial patterns of NDVI and SWE
Averaged over the 25-yr study period (WY 1993(WY -2017, peak NDVI ranged from 0.36 to 0.86 and peak SWE ranged from 13 to 150 cm among the sites in the area of study (Fig. 2). There was an obvious geospatial trend in NDVI and SWE, with clustering of similar values throughout the study area (Fig. 2). This spatial clustering was statistically significant for both peak NDVI (z = 6.13, P < 0.0001) and peak SWE (z = 4.34, P < 0.0001; Appendix S1: Fig. S2). Across the 25yr study period, annual peak NDVI averaged across all 169 stations ranged from 0.59 to 0.65 and annual peak SWE ranged from 38 to 84 cm (Fig. 3). There was no obvious temporal trend in either (Fig. 3) and no statistical relationship between the two variables on an annual basis (r 2 = 0.04; P = 0.34; Fig. 4). While greater snow accumulation occurred with elevated forest productivity in some years (e.g., 1999 and 2017; Fig. 3), this was not the case for all years with high snowpack (e.g., 1997). Similarly, years with below average snow accretion occurred with below average NDVI in some years (e.g., 1994 and 2007; Fig. 3), but this was not a consistent trend across the dataset (e.g., 2005 and 2015; Fig. 3). There was high variability in SWE and v www.esajournals.org 6 August 2020 v Volume 11(8) v Article e03212 NDVI on an annual basis, and no temporal autocorrelation was detected at any spatial scale ( Fig. 4; the only significant model was AR (0) at all spatial scales).

Correlations of NDVI with snowpack metrics, elevation, land cover, and aspect
Averaged across all stations over the entire study period, mean annual NDVI was not statistically related to mean annual SWE (Fig. 4a). Comparing annual means for each station individually over the entire study period, neither peak SWE nor snowmelt date was found to be strong predictors of NDVI (Fig. 4b, c). Furthermore, we found no significant relationships between NDVI and SAI across individual stations (r 2 = 0.009, P = 0.20), averaged by year (r 2 = 0.001, P = 0.68), or stratified by ecoregion (Northern Rockies r 2 = 0.03, P = 0.47; Idaho Batholith r 2 = 0.001, P = 0.48; Middle Rockies r 2 = 0.002, P = 0.60; data not shown). NDVI was significantly related to snow ablation period length, though the strength of this association was low (r 2 = 0.04; Fig. 4d).
NDVI was negatively correlated with elevation, and this relationship was monotonic when comparing across all sites, or when stratified by ecoregion (Fig. 5). When we further subdivided our analysis by ecoregion and elevation zone, we observed statistically significant correlations in only two of the nine possible comparisons (only in high-elevation zones in the Northern Rockies and Middle Rockies ecoregions; Appendix S1: Fig. S3). In addition, we found very few significant relationships between lag SWE and SAI with NDVI; of all 48 combinations of lag SWE and SAI with NDVI among the various elevation classes and ecoregions, on two (Idaho Batholith 2-yr SWE at mid-elevation and Middle Rockies 2-yr SWE at low elevations) were significant (data not shown). Like the significant spatial autocorrelation with NDVI and SWE, there was Fig. 2. Peak NDVI (a) and peak SWE (b) for each SNOTEL station within the study area averaged across the 25-year study period. Note, not all stations contain 25 yr of data, as some years with unreliable NDVI and SWE data were removed. The timing (date) of peak NDVI and peak SWE was manually determined for each individual station for each year. Fig. 3. Annual variation in mean peak NDVI (blue) and peak SWE (box plots) averaged across all 169 stations for each year for the period 1993-2017. Outliers indicate the 5th and 95th percentiles; vertical lines indicate the 10th and 90th percentiles of peak SWE; boxes indicate the 25th and 75th percentiles; horizontal black and red lines designate the median and mean peak SWE, respectively. obvious spatial clustering in the strengths of the NDVI-SWE correlations (Fig. 6).
Of the land-cover types present in our study area, significant NDVI-SWE relationships were found in only the mixed conifer land-cover class, though the strength of this association was very low (r 2 = 0.06; Appendix S1: Table S1). However, only shrub ecosystems were found to be correlated with lag SWE (3-year, r 2 = 0.17, P = 0.02) while none of the land-cover types were significantly correlated with SAI (data not shown). Comparing NDVI-SWE relationships among aspect, NDVI was significantly related to SWE only on south-and west-facing slopes, though with low-to-moderate correlation coefficients (Appendix S1: Fig. S4).

Modeling environmental driving factors of NDVI
Using five non-auto-correlated environmental variables, we evaluated the strength of the significance of each of these to NDVI. This varied among climate divisions (Fig. 7). NDVI was most related (lowest P-value of the five predictor variables) to SWE (four divisions), slope (four divisions), and elevation (five divisions). Our longitudinal model building resulted in statistically significant models for each of the 13 climate divisions, with low-to-high predictability (r 2 values ranging from 0.13 to 0.93; P-values < 0.001; Table 3). Evaluating the sensitivity of NDVI to the variables used in each model (based on the highest coefficient for all model parameters) again revealed differences among climate divisions. NDVI was most sensitive to SWE (three divisions), PDSI (four divisions), slope (four divisions), and elevation (one division).

DISCUSSION
There is lack of consensus with respect to how summer forest productivity is related to wintertime snowpack, especially in water-limited systems. Some studies have found positive relationships between NDVI and SWE (Trujillo et al. 2012, Knowles et al. 2018, while others have found little to no relationship between NDVI and wintertime snow metrics (Peng et al. 2010, Kariyeva and van Leeuwen 2011, Wang et al. 2018. For example, a previous investigation throughout China documented that NDVI was related to changes in snowpack only in arid ecosystems (Peng et al. 2010). The body of literature exploring this relationship ranges from ecosystem-scale investigations using both satellite imagery and eddy covariance data (Trujillo et al. 2012, Knowles et al. 2018, 2019 to individual tree-scale investigations using tree ring and/ or isotopes (Hu et al. 2010, Martin et al. 2017. However, all previous studies using satellite imagery use large-scale (<1 km 2 ) grain size to address this relationship, which may not be appropriate for regions characterized by high heterogeneity of topography, forest community composition, and climate (e.g., Jiang et al. 2006, Kariyeva andvan Leeuwen 2011). To address current knowledge gaps, we quantified the variation in NDVI, SWE, and NDVI-SWE relationships, across 25 yr for 169 small-scale sites in the nIMW. By using smaller grain size in our study, we were able to increase replication, empirical rigor, and the testability of our hypotheses.
We hypothesized that peak NDVI would be correlated with peak SWE on an annual basis because a high proportion of mountain precipitation in the nIMW arrives as snowfall. Our results generally did not support our hypothesis. Within   6. The spatial patterns of the linear regression correlation coefficients (r 2 ) between peak annual NDVI and peak annual SWE between 1993 and 2017 for each of the 169 SNOTEL stations within the study area. Larger circles indicate a greater relationship between NDVI and snowpack. Significant relationships between NDVI and SWE are indicated with pink (P < 0.10) and yellow (P < 0.05) pins, where appropriate (relationships were not significant at all stations). Fig. 7. Comparison of the relationships between various environmental predictor variables and NDVI across 13 climate divisions. LogWorths over 2 indicate a P-value less than 0.01 for any parameter. The effect of snow water equivalent (SWE, blue bars), snowmelt date (SMD, yellow bars), degree of slope (green bars), elevation (orange bars), and Palmer Drought Severity Index (PDSI, red bars) on NDVI are shown. An * indicates statistical significance at the P = 0.01 level for any parameter. v www.esajournals.org 10 August 2020 v Volume 11(8) v Article e03212 our study area, significant correlations between winter snow accumulation and satellite-based vegetation greenness were found at only 10% (16/169) of our sites (Fig. 7). There was, however, a conspicuous clustering of sites with significant NDVI-SWE relationships in the high-elevation mountains in the eastern portions of our study area, where climate is more arid and continental. This is consistent with Winkler et al. (2005) who showed that only in arid environments (deserts and grasslands) was NDVI related to changes in snowpack. In more mesic environments, environmental conditions during the growing season might regulate NDVI more strongly than winter snow (Winkler et al. 2005, Hu et al. 2010, Knowles et al. 2019). This concept was further supported by Berkelhammer et al. (2017), who showed that mountain environments consistently have a twofold stronger response of GPP (as measured by solar-induced florescence) to summer rain relative to winter snowpack, in more mesic environments. We found that peak SWE and NDVI, and other metrics of snowpack, varied geographically, as did the strength of univariate relationships between them (Figs. 2, 4, and 7). In general, peak NDVI and SWE increased with latitude, with highest values in the northern mountains of our study area. This is consistent with these northern mountains residing in a wetter climate typical of the Pacific Northwest, including a denser and greener forest. Besides peak SWE, we found that NDVI was not strongly associated with snowmelt date, ablation period length (r 2 < 0.05 in both cases; Fig. 4), lag SWE, SAI (not shown), or land-cover class (Appendix S1: Table S1), although other studies have reported significant relationships between these variables (Trujillo et al. 2012, Winchell et al. 2016. For example, Winchell et al. (2016) found that earlier snowmelt resulted in decreased vegetation productivity during early spring in subalpine forests at Niwot Ridge, Colorado, USA. The weak univariate associations found in our study may be due to the spatial scales used in our study. Widely varying microclimates and topographies at the 30-m scale may have hampered our ability to detect broad, generalizable relationships (discussed more below).
Investigating elevation relationships, NDVI, but not SWE (data not shown), was strongly and negatively related to elevation, and this was true  within all ecoregions (Fig. 5). This may reflect the more northerly location of our study area relative to previous studies in the southern Rocky Mountains and Sierra Nevada Mountains, as higher latitudes tend to have greater energy constraints on productivity. Previous studies have shown that NDVI is especially sensitive to variation in SWE in mid-elevation zones in the Sierra Nevada Mountains (Trujillo et al. 2012), but our analysis did not detect this (Appendix S1: Fig. S3). This lack of relationship may be due to differences in geographical location. Berkelhammer et al. (2017) showed that summer rains in the Rocky Mountains can greatly influence forest productivity, while the Sierra Nevada Mountains experience much less summer rain. Stratified within ecoregion, we found that NDVI and SWE were significantly related in only two of the nine elevation combinations evaluated (Appendix S1: Fig. S3).
Our alternative hypothesis was that NDVI would be best predicted by environmental variables other than SWE, or a combination of variables which included SWE. We further predicted that environmental parameters, and their relationship with NDVI, would differ across our study area, which has highly variable elevation, topography, and climate. This alternative hypothesis was supported, as the predictors of NDVI varied among the 13 climate divisions. For example, NDVI in some of the highest elevation, most interior (most continental) climate divisions (Central Montana 6, Snake River Drainage 10, and Big Horn River Drainage 11) was most significantly controlled by SWE relative to the other environmental predictors (Fig. 7). In nearly all of southern-most climate divisions (Eastern Highlands 3, Northeastern Valleys 4, Snake River Drainage 10, and the Wind River Drainage 13), NDVI was strongly controlled by slope. It is important to note that in climate divisions adjacent to each other, predictors of NDVI were rarely the same, because topography, elevation, and/or vegetation cover varied so much between climate divisions. Kariyeva and van Leeuwen (2011) also found that phenological variables including NDVI were controlled by distinctly different suites of predictor variables in highly contrasting regional landscapes, which they attributed to the complex heterogeneity in topography among their subregions. While regional climate phenomena such as the Pacific Decadal Oscillation and El Nino Southern Oscillation largely determine vegetation community composition and NDVI, topographical variables such as elevation and slope also strongly modify regional climate patterns via the influences of moisture availability, radiation, and temperature, which can differ greatly across short distances in complex terrain.
There are several considerations that constrained our interpretations. The first is spatial scaling. Physiological processes such as photosynthesis have been shown to have varying responses to environmental drivers across differing spatial scales (e.g., Arkebauer et al. 2009), and similar results have been observed for NDVI (Jiang et al. 2006). Most previous satellite investigations quantifying NDVI-SWE relationships used 1-km resolution. Our study was conducted at a much finer resolution and thus most likely better captured the heterogeneity of the landscapes, and NDVI, in our study area. However, this large heterogeneity caused high variance of predictor and response variables and may have hampered our ability to find strong relationships (Weins 1989, Scholtz et al. 2018) and detect generalizable patterns sometimes found in studies that used lower resolution imagery (e.g., Grippa et al. 2005, White et al. 2005, Trujillo et al. 2012, Table 1). Topography has been demonstrated to greatly influence hydrometeorological processes in topographically complex areas such as western Montana , underscoring that even meter-scale differences in topographic position can impact ecosystem productivity. This highlights a key message of our study, which is that environmental drivers of NDVI can vary widely in regions characterized by very patchy environmental and ecological characteristics, at meso-to micro-scales. Furthermore, this underscores a trade-off that future remote sensing studies should consider: Remote sensing studies of complex regions use lower resolution imagery may over-generalize ecological patterns occurring in heterogeneous areas, while approaches using higher-resolution imagery may encounter issues with statistical power and being able to generalize patterns. The second consideration was that previous studies have shown that environmental drivers such as soil moisture availability, soil organic carbon, and seasonally averaged precipitation and temperature are important predictors of NDVI (e.g., Kariyeva and van Leeuwen 2011). We did not specifically include these variables in our analysis, although some were indirectly covered through synthetic variables that we did use such as SAI, PDSI, and GDD. Instead, we constrained our predictor variables of interest to compare to those used in other studies which focused on wintertime precipitation and snowpack. This is because precipitation during plantdormant seasons is thought to impact productivity during the growing season through soil moisture recharge and/or soil biogeochemical processes, especially in water-limited systems (LaMalfa and Ryle 2008). Therefore, our hypotheses were mostly related to snowpack metrics since hydrological regimes in the mountain forests of the nIMW are snow dominated. Lastly, we encourage further research into the controlling effect of snowpack on vegetation productivity using different remote sensing instruments, and/or studies comparing responses of NDVI to environmental variables at multiple spatial scales. While research has been conducted globally (Table 1), we have yet to see repeated satellite studies in the same location to compare methodologies and the effect of image resolution on vegetation indices.
Our study adds another system to the body of literature quantifying the impact of snowpack on summertime productivity. We showed that environmental controls on NDVI to be widely varying and hard to predict from one subregion to the next. This highlights that remotely sensed studies of this nature may not be appropriate in all regions of the world and/or that studies performed at different resolutions/scales can affect results and interpretations.