Quantifying snow controls on vegetation greenness

Snow is a key driver for biotic processes in Arctic ecosystems. Yet, quantifying relationships between snow metrics and biological components is challenging due to lack of temporally and spatially distributed observations at ecologically relevant scales and resolutions. In this study, we quantified relationships between snow, air temperature, and vegetation greenness (using annual maximum normalized difference vegetation index [MaxNDVI] and its timing [MaxNDVI_DOY]) from ground-based and remotesensing observations, in combination with physically based models, across a heterogeneous landscape in a high-Arctic, northeast Greenland region. Across the 98-km distance from the Greenland Ice Sheet (GrIS) to the coast, we quantified significant inland–coast gradients of air temperature, winter precipitation (using pre-melt snow-water-equivalent [SWE]), and snowmelt timing (using snow-free day of year [SnowFree_ DOY]). Near the coast, the mean annual air temperature was 4.5°C lower, the mean SWE was 0.3 m greater, and the mean SnowFree_DOY was 37 d later, than near the GrIS. The regional continentality gradient was eight times stronger than the south-to-north air–temperature gradient along the Greenland east coast. Across this strong gradient, the mean vegetation greening-up period (SnowFree_DOY-MaxNDVI_ DOY) varied spatially by 24–57 d. We quantified significant non-linear relationships between the vegetation characteristics of MaxNDVI and MaxNDVI_DOY, and SWE, SnowFree_DOY, and growing degree-days-sums during greening-up (Greening_GDD) across the 16-yr study period (2000–2015). These demonstrated that the snow metrics, both SWE and SnowFree_DOY, were more important drivers of MaxNDVI and MaxNDVI_DOY than Greening_GDD within this seasonally snow-covered region. The methodologies that provided temporally and spatially distributed snow, air temperature, and vegetation greenness data are applicable to any snowand vegetation-covered area on Earth.


INTRODUCTION
Seasonal snow is a key driver of climate change effects, and it controls a range of Arctic ecosystem processes (Jones 1999, Post et al. 2009, Brooks et al. 2011, Bokhorst et al. 2016. Snow affects vegetation growth both directly and indirectly (Cooper 2014). For instance, the presence and absence of snow cover influence the surface energy balance (Marks andDozier 1992, Stiegler et al. 2016), which in turn affects the below-ground surface thermal regime. The snow cover acts as an efficient insulator during winter with its low thermal conductivity (Goodrich 1982, Sturm et al. 1997) and high thermal resistance . This insulating effect keeps soil thermal conditions relatively stable during snow-covered periods (Taras et al. 2002, Zhang 2005 and protects the vegetation from frost damages (Bokhorst et al. 2011). In turn, soil thermal conditions regulate decomposition and sequestration rates in the soil (Schimel et al. 2004, Johansson et al. 2013 as well as microbial activity, respiration rate, and amount of soil organic carbon produced during winter (Elberling 2007). Hence, the snow-depth evolution through autumn and winter governs the amount and timing of plantavailable nutrients at the end of winter and the following spring in tundra ecosystems (Schimel et al. 2004, Sturm et al. 2005, Buckeridge and Grogan 2008. Snow cover limits the direct light availability to the below-snow vegetation, and the timing of spring snowmelt initiates the potential vegetation-growing season in the Arctic. Furthermore, precipitation accumulated in the snowpack during winter is released as meltwater in spring (Jones 1999) and provides moisture for vegetation growth at the start and during the growing season , Ellebjerg et al. 2008.
To quantify these multiple, snow-related, direct and indirect effects on vegetation, temporally and spatially distributed observations are needed; this is particularly true at ecologically appropriate resolutions and scales in the Arctic and other remote areas, where ground observations are often sparsely distributed in both space and time. The Arctic remoteness and the polar night, in particular, make winter and spring fieldwork a challenge. Year-round, ground-based investigations and monitoring of interseasonal relationships between snow and vegetation are sporadic in these areas. Furthermore, the spatial extent and resolution of available snow and vegetation observations range from established monitoring sites representing few landscape types, to global remote-sensing observations of 10s of km grid resolutions. Hence, in ecological impact studies, a mismatch in resolution and extent often arise due to missing snow observations at ecologically relevant scales (Bokhorst et al. 2016, Pedersen 2017. Hence, quantifying relationships between snow metrics and vegetation growth is challenging, but highly needed.
To increase the relevance and impact of a given Arctic data set (or any data set), observations can be collected across environmental gradients, which may constitute changes in specific variables, for example, temperature or precipitation, during a defined period of time in the past or future (Henry and Molau 1997). The gradients may be latitudinal (e.g., Kottek et al. 2006, Legagneux et al. 2014 or along an elevation increase (K€ orner 2007). For example, since winter precipitation is a key driver of Arctic vegetation phenology (e.g., Cooper 2014), the investigation of vegetation phenology and biomass responses can be conducted at a plot-scale and along experimental winter-precipitation gradients using snow fences (e.g., Leffler andWelker 2013, Semenchuk et al. 2013). Additionally, responses of vegetation growth to temperature gradients have been observed across multiple Arctic sites using standardized plots, where the different latitude and climate zones of each site act as surrogates for a gradient in temperature and snow conditions (Bienau et al. 2014, Myers-Smith et al. 2015. Such plot-based studies are often used in space-for-time substitution approaches, where vegetation differences along environmental gradients are assumed to represent conditions in the past or future . Other driving factors of Arctic vegetation phenology include snowmelt timing and winter air temperature (Pudas et al. 2008, Hollesen et al. 2015. A delay in timing of snowmelt can even counteract the expected phenological responses to increased temperatures (Bjorkman et al. 2015). Such uncertainties encourage studies along naturally varying environmental gradients using spatially continuous data such as remote-sensing imagery, spanning multiple latitude degrees that represent temperature, irradiation, and moisture gradients (Zeng et al. 2011, Zeng andJia 2013). Seasonal peak vegetation greenness, as inferred from normalized difference vegetation index (NDVI) data (Riedel et al. 2005, Raynolds et al. 2012, shows a positive relationship with growing-season temperature across the Arctic. However, the strength of this relationship has declined between 1982 and 2011 (Piao et al. 2014). Spatially distributed vegetation phenology and timing of peak vegetation greenness show a significant positive correlation with timing of snowmelt in parts of the Arctic (Zeng and Jia 2013), while the direction of this relationship varies with latitude (Grippa et al. 2005).
To address these scale-and data-related challenges, we used a combination of ground-based, remote-sensing, and modeling data, to quantify snow-vegetation relationships across local-(100s of m), landscape-(1000s of m), and regional-scale (10s of km) environmental distributions and gradients in high-Arctic Greenland. First, we calculated and mapped the inland-coast distributions and gradients of air temperature, pre-melt snowwater-equivalent (SWE), and snow-free day of year (DOY), within an ice-sheet-free region in northeast Greenland. These key environmental drivers of Arctic vegetation growth and phenology were modeled using the MicroMet and SnowModel (Liston and Elder 2006a, b) temporally and spatially distributed modeling tools. In situ observations of meteorological and snow variables were used to validate and verify the modeled variables. Second, we investigated the importance of snow and temperature as drivers for the spatial distribution and temporal evolution of maximum vegetation greenness in the region. For this, we quantified the relationships among vegetation greenness, snow, and greening-up temperature sums. For the vegetation greenness variables, we used the annual maximum normalized difference vegetation index (MaxNDVI) and its timing (MaxNDVI_DOY); for snow, we used the pre-melt SWE and Snow-Free_DOY; and for the greening-up temperature sums (Greening_GDD), we calculated the number of growing degree-days (GDD) between the SnowFree_DOY and MaxNDVI_DOY.

Study area
The northeast Greenland study region is centered at 74°27 0 N, 20°34 0 W and includes the Zackenberg valley area. It extends from the Greenland Ice Sheet (GrIS) margin to the Greenland east coast near the Greenland Sea. The fjord system within the study area is covered with sea ice up to eight months of the year. The study region covers an area 98 km west-east by 75 km south-north (3°15 0 longitude by 0°40 0 latitude), and the topographic relief varies between sea level and 1570 m a.s.l. (Fig. 1). This high-Arctic location (Kottek et al. 2006) has a tundra climate (Bliss and Matveyeva 1992) defined by bioclimate subzone C (Walker et al. 2005). The mean annual air temperature equaled À9.0°C in the center area near Zackenberg during the period 1996-2015. The Zackenberg valley consists of a spatial mosaic of vegetation types ranging from highly productive and moist fen to low-productive tundra scrubs (Dryas and Cassiope heath; Elberling et al. 2008

MicroMet and SnowModel
In addition to air temperature, we focused our study on two ecologically relevant snow metrics: the pre-melt SWE (i.e., the maximum amount of water available in the snowpack prior to spring snowmelt in meters depth) and the SnowFree_ DOY (i.e., the DOY when SWE equaled zero during spring melt). In previous studies, both snow metrics have proven to be important controls on Greenland vegetation growth and phenology (Grøndahl et al. 2006, Tamstorf et al. 2007, Lund et al. 2012, Hollesen et al. 2015. To quantify both the regional-scale distributions and landscape-scale heterogeneity of daily mean air temperature, pre-melt SWE, and SnowFree_DOY across the region, we used the MicroMet high-resolution meteorological model (Liston and Elder 2006b) coupled with the SnowModel snow-evolution modeling tool (Liston and Elder 2006a). Nine automated weather stations (AWS) distributed within the region ( Fig. 1) provided model inputs of air temperature, relative humidity, wind speed, and wind direction. Also, gridded reanalysis meteorological inputs including precipitation rate were provided as model input by the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Re-Analysis (ERA-Interim; Dee et al. 2011). We used the 12 ERA-Interim grid cells nearest the simulation domain (Fig. 1). The reanalysis data provided meteorological information (particularly precipitation) during the winter months when no other data were available.
MicroMet and SnowModel were run on a 300m resolution grid for the period 1 August 1999 through 31 July 2015. Since synoptic weather systems producing precipitation often occur on sub-daily time scales across Greenland (Schuenemann et al. 2009), we used a 3-hourly model time increment to resolve the diurnal cycle. Furthermore, we analyzed the model output data at a daily temporal resolution to resolve the timing of variables and events that may occur from one day to the next. To make the reanalysis precipitation values consistent with SWE observations, SnowModel was run while coupled to Snow-Assim (Liston and Hiemstra 2008). This terrestrial Greenland SnowModel configuration was described in detail by Pedersen et al. (2015).

Snow observations
The SnowAssim SWE assimilation used observed pre-melt SWE observations collected in snowpits during the years 2004-2015 in the Zackenberg valley by GeoBasis Zackenberg (Pedersen et al. 2016). In addition, observations of snowpack properties, including SWE, were collected across the domain at 55 sites during an April 2014 snow field campaign (Appendix S1: Table S1). The SWE from five randomly selected sites (Appendix S1: Table S1, data type "a" for ) of the study area in northeast Greenland including regional place names (black letters) and glaciers and ice caps (white shading). Gray points mark the center of 12 European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Re-Analysis (ERA-Interim) grid cells; each 0.75 degrees longitude by 0.75 degrees latitude. Red points and numbers mark the locations of nine automated weather stations (AWS; location names in red). The AWS meteorological data included in the analysis spanned these periods: 1: 2009-2014, 2: 2008-2014, 3: 2008-2014, 4: 2011-2014, 5: 2003-2010, 6: 2003-2015, 7: 1999-2015, 8: 2003-2014, and 9: 1999-2015. assimilation) were included in the SWE assimilation. The remaining 50 SWE observations (Appendix S1: Table S1, data type "v" for validation) were used to validate the spatially distributed SWE across the region for the winter 2013/2014. These spatially distributed SWE observations included local-scale hilltops and depressions that were not resolved by the 300-m SnowModel simulation. Therefore, this localscale variability was removed from the validation data points using the same interpolation methods used in SnowAssim. Snow-depth observations from a snow-depth sensor installed in the Zackenberg valley (Stn. 7, Fig. 1) were used to estimate the SnowFree_DOY, and these data were used to validate the SnowModel simulated SnowFree_DOY.

Precipitation adjustment
A gradient in pre-melt SWE across the region was identified in the snow observations made during April 2014 (Appendix S1: Table S1). This inland-coast gradient was quantified by calculating the slope of the linear regression between the SWE observations and the west-east distance across the region. The observed significant SWE gradient was +0.0059 m (5.9 mm)/km distance from inland to the coast. In contrast, an inlandcoast gradient in modeled SWE, when forced with the original ERA-Interim precipitation data, was +0.0011 m (1.1 mm)/km (Fig. 2). To correct this, we randomly selected five SWE observations across the region and included them in the SWE assimilation to generate a smoothed distribution of precipitation-correction factors across the domain using SnowAssim. The correction factors varied between 0.8 and 1.8 from inland to the coast, respectively.
The original ERA-Interim precipitation values were adjusted using this west-east correctionfactor prior to being included in the final Snow-Model simulation. Initially, this correction was applied to a SnowModel simulation for the year 2013/2014. The resulting modeled, significant SWE gradient, with corrected precipitation rates, equaled +0.0052 m (5.2 mm)/km, which matched the observed SWE gradient. Linear regression revealed that neither slope nor intercept of the observed and modeled SWE (with adjustment) data set was significantly different (P = 0.519; Fig. 2). Based on this, the precipitation-correction factors were applied to all other years in the SnowModel simulation. This approach was also supported by the fact that the original ERA-Interim precipitation data showed an inlandcoast gradient in precipitation in all years from 2000 through 2015; it was just weaker than our field observations showed. The ERA-Interim precipitation gradient may be weaker because it was estimated from the six grid cells that cover the west-east dimension of the study area; this resolution is likely too coarse to resolve the precipitation effects of the coastal mountains within the study area.

Quantifying inland-coast gradients
One goal of this study was to characterize the regional snow-related growing conditions and quantify the inland-coast gradients in air temperature, pre-melt SWE, and SnowFree_DOY for the low-altitudinal and vegetated areas across the region. Therefore, to minimize the orographic effect on the inland-coast gradient estimates, while also including vegetated valley areas, we extracted each modeled variable (air temperature, pre-melt SWE, and SnowFree_DOY) in all land-covered grid cells at elevations below 300 m a.s.l. and used those in our analyses. This comprised 66% of all vegetated grid cells within the region. Using data from these grid cells, we computed an annual average of each variable for every 300-m longitudinal increment along the west-east dimension of the domain. These 16 longitudinal averages were averaged to yield the overall gradients across the 98-km west-east distance for the variables of interest. To evaluate the relative strength of our regional gradients in a Greenland-wide context, we compared our regional inland-coast temperature gradient with a large-scale gradient along the east coast of Greenland. Using nine Danish Meteorological Institute (DMI) coastal weather stations located in east Greenland between 60°N and 80°N (Cappelen 2014(Cappelen , 2016.

NDVI
To assess the impact of the quantified environmental gradients on vegetation growth patterns across the region, we used a data set of satellitederived NDVI, which is an integrated measure of vegetation growth, greenness, and aboveground phytomass (Tucker 1979, Myneni et al. 1997. NDVI was calculated by combining Moderate Resolution Imaging Spectroradiometer (MODIS) Daily Surface Reflectance MOD09GQ (Collection 6) at 250-m resolution (reflectance band 1, 620-670 nm, and band 2, 841-876 nm) and the State Quality Assessment (QA) flag from MOD09GA (Collection 6) at 1-km resolution (Vermote et al. 2016). MOD09GA is considered the primary and recommended overall quality-assessment data source for MODIS surface reflectance products (Vermote et al. 2016); it was used to qualityscreen the 250-m reflectance data. Both data sets were regridded to the 300-m SnowModel grid using nearest-neighbor resampling, and reflectance bands 1 and 2 were filtered using the QA flags to identify cloud-free data for the analysis. NDVI was determined from the difference between the infrared radiance (band 2) and red radiance (band 1) divided by their sum (Tucker 1979). The quality-checked, daily NDVI data set spanned 1 January 2000 through 31 December 2015 and was used to estimate MaxNDVI and MaxNDVI_DOY.
To ensure that only vegetation-covered areas of the region were included in the analysis, we selected grid cells from the original NDVI data set that had a minimum of 15 NDVI values above 0.2 NDVI per year to be included in the estimation of MaxNDVI (Raynolds et al. 2006, Ding et al. 2016). This minimum of 15 valid data points per growing season was based on counts of valid grid cells in known vegetated areas, for example, the Zackenberg valley, and was chosen to ensure that these ecologically important areas were included in the analyses. The original NDVI data set included several gaps, from either MODIS data files missing from the original NASA archive or from data removed as part of the quality checks. These gaps were filled using temporal linear interpolation between the available valid original NDVI values, under the assumption that the NDVI values from one day to the next were temporally auto-correlated. MaxNDVI was estimated from these interpolated NDVI values, where MaxNDVI was defined to be the average of the 10% highest NDVI values per growing season occurring between 1 April and 31 October.
Daily, ground-based NDVI observations were available from three stations in the Zackenberg valley, each with a sensor located above a different vegetation type: Dryas heath, Cassiope heath, and fen. The MODIS MaxNDVI was validated against MaxNDVI, which was derived from these ground-based NDVI observations using the same method, that is, the average of the 10% highest NDVI values per growing season.
In these daily, ground-based NDVI data, MaxNDVI_DOY for each vegetation type was clearly visible (not shown) and the DOY was identified manually and defined to be the observed MaxNDVI_DOY. In order to develop a method to distinguish the MaxNDVI_DOY automatically, we used these ground-based NDVI observations. We cumulated the NDVI values above 0.2 NDVI during the period 1 April through 31 October and identified the DOY, where different percentage thresholds (of the total accumulated NDVI for the season) were reached (30-60%, in 10% increments). The 40%threshold resulted in the lowest mean difference between the observed MaxNDVI_DOY and estimated MaxNDVI_DOY. Therefore, this method of 40%-threshold was applied to the daily MODIS NDVI data. The annual MaxNDVI_DOY from MODIS NDVI was extracted for the grid cells covering the three sensor locations in years with available ground-based NDVI observations. These were compared against the ground-based MaxNDVI_DOY to assess the difference between the ground-based and MODIS MaxNDVI_DOY. The duration of the greening-up period, defined to be the number of days between SnowFree_ DOY and MaxNDVI_DOY, was estimated for all 16 yr for all grid cells in the study domain. The daily temporal resolution limited the temporal uncertainty of the estimated greening-up period to 1-2 d.
As a measure for the air temperature during the first part of the growing season, from Snow-Free_DOY to MaxNDVI_DOY, we summed the GDD, that is, daily mean air temperature above 0.0°C, throughout the annual greening-up per grid cell in the study area. The annual GDD sums for the greening-up (Greening_GDD) were used to investigate the importance and control of temperature on MaxNDVI and MaxNDVI_DOY. Greening_GDDs were significantly correlated with annual mean air temperature (r = 0.76, slope = 0.013, R 2 = 0.94, P < 0.001). The regional mean annual air temperature was used to evaluate the inland-coast gradient against airtemperature gradients observed along the entire east coast of Greenland (i.e., spanning 23°l atitude).

Statistical analyses
Linear regression was used to compare groundbased and MODIS MaxNDVI and MaxNDVI_ DOY. To determine the spatial association between the snow and vegetation variables, we examined the relationships between spatially distributed pre-melt SWE and SnowFree_DOY, and MaxNDVI and MaxNDVI_DOY, using Pearson's product-moment correlation coefficient, r, based on annual averages for the period 2000-2015. Furthermore, we quantified the relationships between snow metrics, NDVI metrics, and Green-ing_GDD using the longitudinal averages for each variable spanning the years 2000-2015 using generalized linear models (GLMs) with successive removal of non-significant (P > 0.05) predictor variables. Predictor variables were included as both linear and quadratic. The statistical analysis was conducted in R (www.R-project.org.) and using FORTRAN programming for computational efficiency.

Temporal linkages between snow and NDVI data sets
We found clear temporal linkages between the two spatially distributed and temporally evolving, but completely independent, data sets we developed and analyzed: SnowModel outputs and MODIS NDVI (Fig. 3). The timing of daily SnowModel SWE evolution and snow disappearance corresponded well with the greening-up evolution defined using the daily NDVI data. For instance, during the spring months, NDVI values started increasing (vegetation greening-up) when the snowpack had melted away (after the Snow-Free_DOY). Furthermore, in the autumn, when the first snow storms had occurred and a persistent snow cover was established, the NDVI values were again below 0.1. Hence, the annual cycle of snow and NDVI dynamics temporally matched despite that the analyzed data originated from different sources, procedures, and tools. This suggests the validity of our methodology for describing the abiotic and biotic dynamics of this system, and the relationships among them. An example of the spatial and daily temporal evolution of the SnowModel snow depth and the MODIS NDVI values across the region in year 2003 is provided in Video S1.

Model validation
We found a linear correlation between modeled and observed SWE in April 2014 (P < 0.001, Fig. 4a, b). The modeled SWE captured the observed spatial distribution of SWE within the heterogeneous landscape of the region, where SWE increased with distance from the GrIS margin (Fig. 4a). Furthermore, the SnowModelestimated SnowFree_DOY and the SnowFree_ DOY, derived from snow-depth sensor measurements at Zackenberg valley (Stn. 7, Fig. 1), were strongly correlated (R 2 = 0.72, P < 0.001, Fig. 4c).
❖ www.esajournals.org The estimated timing of 40% seasonally cumulated NDVI (Estimated MaxNDVI_DOY) and the observed timing of maximum NDVI per season for each NDVI sensors (Observed MaxNDVI_ DOY) in the Zackenberg valley were significantly correlated (P < 0.001, Fig. 5a). MaxNDVI levels derived from the MODIS data were not significantly correlated with MaxNDVI derived from the three Zackenberg NDVI sensors (r = 0.43, P = 0.069, Fig. 5b). MaxNDVI from the three sensor locations was significantly different (regression lines not shown, P < 0.001), but only the MaxNDVI observed at Dryas heath was significantly correlated with the MODIS MaxNDVI (r = 0.65, P < 0.05), whereas the MaxNDVI at Cassiope heath and fen was not (P > 0.05, Fig. 5b). This lack of direct correspondence of MaxNDVI in all vegetation types may originate from the higher variability in NDVI values within the MODIS 300-m (90,000 m 2 area) grid cell (i.e., subgrid variability in vegetation types) than within the NDVI station-sensor footprint of a few square meters. Yet the MODIS MaxNDVI values, ranging 0.43-0.68, were within the range in ground-based MaxNDVI of 0.39-0.76 (Fig. 5b).
The ground-based MaxNDVI_DOY and MODISderived MaxNDVI_DOY for all station locations were significantly correlated (r = 0.83, P < 0.001, Fig. 5c). The validation for both snow and NDVI variables allowed us to derive spatial representations of the heterogeneous patterns of pre-melt SWE, SnowFree_DOY, MaxNDVI, and MaxND-VI_DOY. These data were then available to quantify the inland-coast gradients of these variables.

Regional snow and NDVI distributions and gradients
The spatial distribution of pre-melt SWE showed a significant gradient across the region, with less SWE available in the inland area near Payers Land and western Clavering Island compared to the more snow-rich Wollaston Foreland area near the east coast (Fig. 6a, see place names in Fig. 1). The SnowFree_DOY followed the distribution of pre-melt SWE; areas with less SWE had earlier SnowFree_DOYs, and vice versa (Fig. 6b). Comparing winters with contrasting SWE conditions, for example, the snow-poor winter 2008/2009 and the snow-rich winter 2011/2012 revealed that these gradients in pre-melt SWE and SnowFree_DOY existed across the region during all simulation years, regardless of precipitation quantity (Appendix S1: Fig. S1a, b).
Our analysis also revealed that, superimposed on this overall inland-coast gradient, was a high degree of heterogeneity due to local-and landscape-scale topography-related patterns. For example, the pre-melt SWE increased with elevation and the SnowFree_DOY was delayed with increasing elevation, due to the simulated precipitation increase with elevation, for example, in west Clavering Island (Fig. 6a, b). Similarly, the landscape-scale distribution of MaxNDVI showed a decrease in vegetation greenness with increasing elevation (Fig. 6c), whereas the spatial distribution of MaxNDVI_DOY showed both a delay with increasing elevation and across the region from west to east (Fig. 6d).
From the spatial distributions between the GrIS margin in the west, to the coast in the east, we found horizontal gradients in air temperature of À0.05°C/km, in pre-melt SWE of +0.003 m/km, in SnowFree_DOY of +0.4 d/km, and in the MaxNDVI_DOY of +0.1 d/km (Fig. 7). These significant gradients equaled a total difference across the 98-km-wide domain of À4.5°C in average air temperature and 0.3 m average pre-melt SWE. The average SnowFree_DOY and MaxND-VI_DOY varied up to 37 and 7 d, respectively, across the region. However, no significant gradient of MaxNDVI was detected between the GrIS and the coast (P > 0.05; Figs. 6c, 7e); although the annual average MaxNDVI was pixel-wise significantly positively correlated with annual average pre-melt SWE and SnowFree_DOY, especially within the abundantly vegetated valley areas of the region (P < 0.05, r > 0.5; Fig. 6c). In the Zackenberg valley, the level of MaxNDVI was rather stable across years compared to the pronounced interannually varying pre-melt SWE (Fig. 3). Furthermore, the persistent pattern of MaxNDVI across the region was also identified in years with distinctly different SWE conditions (Appendix S1: Fig. S1a, d). In contrast, the MaxNDVI_DOY co-varied with levels of premelt SWE and the associated dynamics of Snow-Free_DOY (Appendix S1: Fig. S1a-c). In order to relate these regional gradients to the temperature difference in Greenland, we calculated a significant mean annual air-temperature gradient of À0.0058°C/km from south to north along the east Greenland coast and averaged over the period 2000-2015 (R 2 = 0.96, P < 0.001). Hence, the À0.05°C/km regional gradient, from inland to the coast in our study area, was approximately eight times stronger than this large-scale east Greenland gradient.

Non-linear relationships between NDVI, snow, and greening-up temperatures
Generally, MaxNDVI_DOY was earlier in years with early SnowFree_DOY (e.g., 2008/ 2009). Similarly, MaxNDVI_DOY was later in years with later SnowFree_DOYs, which were associated with above-average pre-melt SWE (e.g., 2011/2012; Appendix S1: Fig. S1a, b, d). Over the 16-yr NDVI period, we found an average greening-up period ranging between 25 d and 57 d in the vegetated areas across the region (Fig. 8a). Across years, the greening-up period was consistent and varied only 1 d to 5 d (i.e., standard error on the mean) for each pixel (Fig. 8b). The longest greening-up period was found in the inland area 0-45 km from the GrIS margin and averaged 54 d AE 2 d. The annual greening-up period decreased in duration with distance from GrIS margin. In some areas of the domain, the greening-up period differed from the average inland-coast gradient; at higher elevations, the greening-up period was shorter due to either late SnowFree_DOY or early MaxND-VI_DOY, and the greening-up period was longer primarily in the inland areas because of the relatively early SnowFree_DOY there (Fig. 8a).
However, these areas were generally located outside the areas where we found significant positive correlations between SnowFree_DOY and MaxNDVI_DOY (P < 0.05).

Non-linear snow-vegetation greenness relationships
Previous Arctic studies have identified snow conditions and the summer temperature regime as significant drivers of vegetation greenness (Walker et al. 1999, Dye and Tucker 2003, Bhatt et al. 2013. The majority of these studies have quantified a positive linear correlation between vegetation greenness (MaxNDVI) and summer temperature (e.g., using a monthly summer warmth index; e.g., Jia et al. 2006, Raynolds et al. 2008. Other studies find no or even a weakening of the relation between MaxNDVI and summer temperatures (Blok et al. 2011, Piao et al. 2014. In contrast, we quantified significant nonlinear relationships between MaxNDVI and MaxNDVI_DOY, and the snow conditions of the preceding winter (both pre-melt SWE and Snow-Free_DOY) and the temperature sums during the greening-up period. Previous studies that include moisture sources as a driver of vegetation growth also found non-linearity between the water resources (i.e., soil moisture, which is related to pre-melt SWE amounts) and vegetation greenness at high latitude (Jeong et al. 2013). In addition, Papagiannopoulou et al. (2017) found non-linear relations between climate variables, for example, air temperature and In the present study, the consistent effects of snow and air temperature on MaxNDVI and MaxNDVI_DOY were observed across a strong gradient in pre-melt SWE and air temperature. In the analyses, we included a wide range of snow conditions; pre-melt SWE ranging from 0.05 m to 1.06 m, and SnowFree_DOY ranging from 134 (14 May) to 212 (31 July). Similarly, the Green-ing_GDDs ranged from 3.8°C to 815.0°C within the study region. Furthermore, the vegetation greenness response within our study area included MaxNDVI levels ranging from 0.23 to 0.84; these values likely represent a broad range of vegetation types adapted to high-Arctic growth conditions. The MaxNDVI range within our study region is wider than that found locally within the vegetated Zackenberg valley (Buus-Hinkler et al. 2006). The span in SnowFree_DOY of 78 d (~2.5 months) within the study region exceeds the variation in SnowFree_DOY (~30 d) found in most plot-scale experiments investigating vegetation responses to snow-free date (Wipf and Rixen 2010). Within the study region, the significant non-linear NDVI-snow relationships are likely due, at least in part, to the large spread in premelt SWE and SnowFree_DOY. We investigated the response of MaxNDVI and MaxNDVI_DOY, to SnowFree_DOY, over a period that varied from early snow-free date in mid-May to late snow-free date in early August. Therefore, the non-linear relationships were established across a wider range and larger gradient of SnowFree_ DOY, and thereby larger differences in air temperature and incoming solar radiation, than included in previous studies.
The strength of these non-linear relationships demonstrated that a larger part of the spatial and temporal variability in both MaxNDVI and MaxNDVI_DOY was explained by pre-melt SWE and SnowFree_DOY, than by Greening_GDD. This result supports the findings of Buus-Hinkler et al. (2006) who showed that moisture available from snowmelt and the timing of snow-free ground are more important drivers of NDVI than summer temperatures in the snow-dominated ecosystems within the Zackenberg valley. Similarly, snow accumulation and timing of snowmelt have been identified as more important drivers of biomass production (measured as NDVI) during the growing season, than summer temperatures, at both a plot scale (Wipf and Rixen 2010) and a continental scale (Grippa et al. 2005).

Relationships established across strong gradients
The snow-vegetation relationships presented herein have been established for a northeast Greenland area that spanned pronounced differences in vegetation growth conditions from inland to the coast. Within the region, the difference between the first and last SnowFree_DOY, averaged 37 d and the MaxNDVI_DOY spanned a week from inland to the coast. Furthermore, the annual mean pre-melt SWE ranged more that an order of magnitude (0.05-1.06 m) across the region. The average pre-melt SWE gradient of +0.003 m/km was less steep than the observed SWE gradient in 2014 (Fig. 2). Inland-coast gradients of this magnitude have been previously identified in Greenland data sets (e.g., Ohmura and Reeh 1991, Schuenemann et al. Fig. 9. Longitudinal averages (AE standard error on the means; gray error bars) for the study region during the period 2000-2015 of (a) MaxNDVI and the pre-melt SWE (m), and (b) MaxNDVI_DOY and SnowFree_DOY. The color scale shading is growing degree-days-sums during vegetation greening-up (Greening_GDD;°C). Lines show the quadratic models for the significant non-linear relationships. Summary statistics are shown in Table 1. 2009, Lucas-Picher et al. 2012); however, none of these studies include the orographic effects and did not map the details presented in the current study. In another coastal Greenland study, the inland-coast annual precipitation gradient for west and southwest Greenland equaled +0.0029 m/km and +0.0071 m/km, respectively (Mernild et al. 2014). These are similar to the premelt SWE gradient of +0.003 m/km found in our northeast Greenland region (Fig. 7b). Greenland air-temperature gradients were investigated by Taurisano et al. (2004), who found an inlandcoast gradient for Nuuk Fjord in west Greenland in summer (March-August) of À0.04°C/km. This differed slightly from the summer (March-August) À0.065°C/km air-temperature gradient in our study area. Quantifications of inland-coast air-temperature gradients based on spatially and temporally continuous data are scarce in Greenland. Hansen et al. (2008) found an inverse, yet weak, summer signal locally between two locations 25 km apart: Zackenberg and Daneborg (Fig. 1). The distance between these two locations represents less than a third of the 90 km inland-coast distance across our study area. Hence, the gradient may originate from localscale difference in, for example, topographic aspect and distance to coast that influence air temperature. In contrast, our gradients are resolved from a spatially continuous data set of daily air temperatures that are averaged for the period March through August during several years and across several grid cells.
Other studies have shown that within the coastal, ice-sheet-free area of all of Greenland, which spans 23 latitudinal degrees, there are pronounced south-north air-temperature gradients along the west and east coasts, which are strongest in winter and spring months (November-April; Abermann et al. 2017). We found that our regional inland-coast temperature gradient was approximately eight times stronger than this large-scale gradient along the east coast of Greenland. Furthermore, when comparing the Zackenberg area regional inland-coast gradients to other Arctic regions, we found that gradients of these magnitudes have not previously been documented.
In a similar study, Kobayashi et al. (2016) estimated the latitudinal gradients in snowmelt timing (i.e., SnowFree_DOY) across the boreal forest and tundra from south to north across Alaska (61-71°N ) using satellite imagery. They found spring SnowFree_DOY gradients of +0.018 d/km in 2011 and +0.053 d/km in 2012 from south to north. In contrast, the average SnowFree_DOY gradient across our 98-km-wide Greenland study region was approximately an order of magnitude stronger (+0.4 d/km). These gradient differences are likely related to differences in land-cover types (i.e., boreal forest/tundra vs. tundra), topographic relief, and relationships between synoptic-scale atmospheric circulation patterns and geographic features such as distance from the coast or ice sheet.
We surmise that the pronounced topographic relief in our Greenland study region may affect the distribution of the inland-coast snow precipitation gradients. Within the first 35 km from the GrIS margin, the average pre-melt SWE and SnowFree_DOY remained relatively constant (Fig. 7b, c). These average gradients were extracted from grid cells in the valley parts of the area below 300 m a.s.l. in order to limit the orographic effect in the gradient estimates. In the inner part of the study area fjord system, the valleys become narrower, the mountains are higher, and their slopes are steeper than in the lower mountains near the coast. Inland, the high mountains on Payers Land, Clavering Island, and the A. P. Olsen Ice Cap (Fig. 1) may affect the precipitation rates at lower elevations downstream of the predominant moisture advection trajectories, since they may act as barriers for the snow precipitation. Hence, although we have accounted for the effects of increased precipitation rates with elevation by including only data from lowland areas in our gradient estimates, the effect of these relatively large-scale landscape features acting as precipitation barriers may still exist in the extracted snow data. If we were to exclude the westernmost 35 km of relatively constant premelt SWE (Fig. 7b), the resulting regional gradient in both pre-melt SWE and SnowFree_DOY would result in an even stronger gradient across the remaining distance between there and the coast.

Ecologically relevant resolution
By combining high-resolution SnowModel simulations with daily MODIS NDVI spatial representations, we quantified gradients in ecologically relevant variables across complex mountainous terrain within a region of ice-sheetfree northeast Greenland. These strong regional gradients were superimposed upon a high degree of landscape-scale variability in the heterogeneous terrain, which was captured in both the SnowModel outputs and the MODIS data sets (Figs. 6-8). For instance, the mapped MaxNDVI (Fig. 6c) showed that vegetation greenness varied spatially within in the Zackenberg valley. In addition, this vegetation-greenness heterogeneity was present across all valley areas of the region, displaying landscape-scale variability. This spatially distributed variability of highly productive and less productive vegetation originated from local-and landscape-scale processes and dynamics. These are associated with precipitation and air temperature lapse rates, and topographic aspect relative to wind and solar radiation exposure, factors that affect the distribution and evolution of snow and other key environmental variables. Within a small fraction of our study area, the local-scale (10s of meters) variability in vegetation greenness is quantified using daily camera imagery (Westergaard-Nielsen et al. 2017). These results suggest that ecological processes occur at finer spatial resolution than 300 m, which likely contribute to the sub-grid variability in MaxNDVI estimates in each 300-m grid cell across our study area.
Our results emphasize that the landscape-scale gradients are complex and non-linear and may differ from the overall inland-coast regional gradients defined in Fig. 7. The 300-m spatially and temporally varying gridded data sets proved useful to study these local-and landscape-scale relationships because they adequately resolved the NDVI patterns and the variations in SWE, and the temporal evolution of both. They captured the variations between valleys, hills, and mountain slopes, and they limited the accuracy of the greening-up estimates to 1-2 d because of the daily temporal resolution of the SnowFree_ DOY and MaxNDVI_DOY data sets. This high temporal resolution is important particularly in a high-Arctic climate in order to resolve year-toyear differences in the relatively short growing season and pronounced variability in the length of the snow-free period (Pedersen et al. 2016).

Future SnowModel and MODIS NDVI applications
The temporally and spatially distributed snow and vegetation greenness data used in this project were produced at ecologically relevant spatial and temporal scales. These types of data sets can be used for scaling carbon and energy exchange from local scale in situ observations in the Zackenberg valley (Lund et al. 2012, Mastepanov et al. 2013, Stiegler et al. 2016 to the larger inlandcoast region. Furthermore, the spatially and temporally distributed data sets can be used to understand the ecology of plants and animals (Kankaanp€ a€ a et al. 2018), and are particularly suited for studying mobile species such as muskoxen whose movements and foraging are dependent on dynamic snow characteristics and vegetation greenness and phenology across the region (Forchhammer et al. 2005, Mosbacher et al. 2016, Schmidt et al. 2016. The data sets were generated using methodologies, that is, MicroMet/ SnowModel and MODIS NDVI processing, that are applicable to any seasonally snow-and vegetation-covered area on Earth. MicroMet and SnowModel have been validated across the Arctic and in other seasonally snow-covered areas around the world (Liston and Sturm 1998, Greene et al. 1999, Liston et al. 2000, Hasholt et al. 2003, Hiemstra et al. 2006, Liston and Hiemstra 2008, Randin et al. 2009, Stuefer et al. 2013, Mernild et al. 2016. Daily, 250-m resolution, MODIS NDVI data are freely available and cover the entire globe since 2000. Hence, using these data sets as input, the NDVI processing method presented herein can provide MaxNDVI and its timing at a landscape scale for virtually any study area of interest inside and outside the pan-Arctic region. This will enable continued investigations of seasonal snow's importance and controlling role on vegetation greenness in other vegetated landscapes.

CONCLUSIONS
We quantified vegetation-related environmental gradients across an area of complex mountainous terrain in an ice-sheet-free region of northeast Greenland, for the period 2000-2015.
The gradients were defined for the following variables: annual average air temperature, premelt SWE, the timing of snowmelt (SnowFree_ DOY), MaxNDVI, and the timing of MaxNDVI (MaxNDVI_DOY). NDVI data were provided using NASA daily MODIS reflectance data, and the meteorological and snow variables were simulated using MicroMet and SnowModel. Both MODIS NDVI estimates and modeled snow variables were validated against ground observations. Our findings include the following: 1. Significant inland-coast gradients in premelt SWE and SnowFree_DOY were identified in all years from 2000 to 2015. These gradients were similar during contrasting snow-rich and snow-poor years, and during average snow years. 2. While the MaxNDVI level was rather stable throughout the 16-yr study period, the MaxNDVI_DOY varied strongly from year to year. This interannual variation was mainly driven by variability in SnowFree_ DOY. 3. The greening-up period, that is, number of days between SnowFree_DOY and MaxND-VI_DOY, varied spatially by 24-57 d across the region; for any given position on the landscape, it only varied by a few days throughout the 16 yr.

Significant non-linear relationships between
MaxNDVI and MaxNDVI_DOY, and the SnowFree_DOY, pre-melt SWE, and greening-up temperature sums were quantified. For this region of strong environmental gradients, the snow metrics were more important drivers of MaxNDVI and MaxNDVI_DOY than air temperature during the greening-up period. 5. The study region's inland-coast mean annual air-temperature gradient was eight times stronger than the mean annual airtemperature gradient from south to north along the entire east coast of Greenland.

ACKNOWLEDGMENTS
We wish to thank the logistics team at Zackenberg Research Station, Aarhus University, for their support during the snow fieldwork in April 2014. Furthermore, we gratefully acknowledge the logistic support of Arctic Research Centre (ARC), Aarhus University. Support was also provided by the Canada Excellence Research Chair (CERC). We thank the institutions and people responsible for gathering the data made available through the Greenland Ecosystem Monitoring Programme. We also thank the Climate Research Section, Central Institute for Meteorology and Geodynamics (ZAMG) in Austria, for providing meteorology and snow data from Freya Glacier. We are thankful to the Danish Meteorological Institute for providing meteorological data from the stations they manage, and we are indebted to NASA's Earth Science Program for the MODIS data. This study was funded by the Danish Environmental Protection Agency, the Danish Energy Agency, and United States National Science Foundation Grant 1023562, and is a contribution to the Arctic Science Partnership (ASP) asp-net.org.