Arctic greening associated with lengthening growing seasons in Northern Alaska

Many studies have reported that the Arctic is greening; however, we lack an understanding of the detailed patterns and processes that are leading to this observed greening. The normalized difference vegetation index (NDVI) is used to quantify greening, which has had largely positive trends over the last few decades using low spatial resolution satellite imagery such as AVHRR or MODIS over the pan-Arctic region. However, substantial fine scale spatial heterogeneity in the Arctic makes this large-scale investigation hard to interpret in terms of vegetation and other environmental changes. Here we focus on one area of the northern Alaskan Arctic using high spatial resolution (4 m) multispectral satellite imagery from DigitalGlobe™ to analyze the greening trend near Utqiaġvik (formerly known as Barrow) over 14 years from 2002 to 2016. We found that tundra vegetation has been greening (τ = 0.65, p = 0.01, NDVI increase of 0.01 yr−1) despite no overall change in vegetation community composition. The greening is most closely correlated to the number of thawing degree days (R2 = 0.77, F = 21.5, p < 0.001) which increased in a similar linear trend over the 14 year study period (1.79 ± 0.50 days per year, p < 0.01, τ = −0.56). This suggests that in this Arctic ecosystem, greening is occurring due to a lengthening growing season that appears to stimulate plant productivity without any significant change in vegetation communities. We found that vegetation communities in wetter locations greened about twice as fast as those found in drier conditions supporting the hypothesis that these communities respond more strongly to warming. We suggest that in Arctic environments, vegetation productivity may continue to rise, particularly in wet areas.


Introduction
Vegetation trends are important to the carbon balance of Arctic ecosystems (Joos et al 2001, Mishra andRiley 2012). Organic carbon content in Arctic soils is about 1300 Pg (Hugelius et al 2014); roughly twice the current total atmospheric carbon content (IPCC 2013). Accelerated warming is occurring in the Arctic as a result of climate change and positive feedbacks (Chapin et al 2005, Serreze andFrancis 2006, IPCC 2013) putting this large carbon pool at risk of loss to the atmosphere (Schuur et al 2013, Schuur et al 2015. Increasing permafrost temperatures (Romanovsky et al 2017), lateral flow of organic matter (Spencer et al 2015), and disturbance (Price et al 2013) are contributing to carbon losses. Carbon exchanges vary between wet and dry vegetation communities where wet communities (often standing water dominated by wetland sedges) are typically characterized by strong uptake of carbon dioxide (CO 2 ) (Sturtevant andOechel 2013, Treat et al 2018) and strong methane (CH 4 ) emissions (Davidson et al 2016b, Treat et al 2018. Dry communities (dominated by mosses, lichens and shrubs with a sub-surface water table) typically exhibit weaker uptake or sometimes net release of CO 2 and low emissions of CH 4 (Natali et al 2015, Treat et al 2018.
Changes in vegetation productivity have been measured by monitoring the normalized difference vegetation index (NDVI) calculated from satellite imagery (Jia et al 2003, Goetz et al 2005, Bhatt et al 2010. NDVI is positively correlated to vegetation biomass across Arctic tundra biomes (Jia et al 2003, Raynolds et al 2011, Epstein et al 2012 therefore NDVI increases are defined as 'greening' and NDVI decreases are defined as 'browning'. Researchers have detected greening trends across the pan-Arctic over the past several decades (Jia et al 2003, Bhatt et al 2010, Bhatt et al 2014; however, greening trends have started to weaken (Piao et al 2014) or reverse with large browning areas, specifically in the European Arctic and Seward Peninsula (Bhatt et al 2013, Phoenix and Bjerke 2016, Bhatt et al 2017, Lara et al 2018. Greening has predominantly been attributed to large-scale climate conditions including rising temperatures linked to reduced albedo over the ocean due to sea ice decline (Bhatt et al 2010, Bhatt et al 2014, Macias-Fauria et al 2017. In lower latitude Arctic regions, greening has been associated with shrub encroachment (Tape et al 2006, Forbes et al 2010, Myers-Smith et al 2011, which also has a positive feedback with warming (Chapin et al 2005). Greening has also been correlated to shifting moisture regimes (Bhatt et al 2017, Westergaard-Nielsen et al 2017 and increases in temperature and precipitation due to movements of air masses from lower latitudes (Macias-Fauria et al 2017).
Most studies on Arctic greening use coarse satellite imagery including AVHRR and MODIS with spatial resolutions of 1-8 km and 0.25-1 km, respectively. These images provide broad coverage of the pan-Arctic, are relatively continuous temporally, and are useful for determining global trends. However, due to the coarse spatial resolution, fine details of patterns and processes are indiscernible and has limited our understanding of greening and browning (Bartsch et al 2016, Myers-Smith et al 2019, particularly in ecosystems characterized by fine spatial heterogeneity such as Arctic tundra (Webber 1978, Billings andPeterson 1980).
Recent studies have taken more detailed approaches to assess greening in the Arctic by utilizing Landsat satellites (30 m spatial resolution) (Frost et al 2014, Nitze and Grosse 2016, Raynolds and Walker 2016, Lara et al 2018. Some studies found that browning can be muddled by surface water due to the strong absorption radiance by water despite greening of vegetation (Raynolds and Walker 2016). Others have shown how the heterogeneity of vegetation cover, geomorphic landscape type, and climate regimes can impact the rate and trend of greening (Lara et al 2018). While Landsat's 30 m spatial resolution is able to cover a large spatio-temporal extent, it does not fully capture the heterogeneity within Arctic ecosystems, which have fine scale polygonal landforms often less than 30 m across (Webber 1978, Billings andPeterson 1980).
Here we present a high spatial resolution (4 m) analysis of greening for an area near Utqiaġvik (formerly Barrow), Alaska. As more high spatial resolution satellite imagery is collected, these time-series approaches will become increasingly valuable (Stow et al 2004) and offer new opportunities and challenges for studying these complex, heterogeneous ecosystems. These tools and their results will likely prove beneficial for modeling carbon dynamics over multiple scales and validating trends reported by larger scales studies. In our study, we aim to (1) show if vegetation communities have changed, (2) analyze fine-scale spatial NDVI trends, and (3) assess correlates of shifts in vegetation communities and NDVI.

Study Area
Our study was conducted on the Barrow Environmental Observatory (BEO) located near Utqiaġvik, Alaska (figure 1). The BEO has a long research history, resulting in a high frequency of satellite tasking, a relatively steady meteorological data record, and accessible field sites. Located in the Arctic Coastal Plain, the site is on continuous permafrost and consists of several drained lake basins (hereby basins), differing in time since drainage (Hinkel et al 2003), and icewedge polygon formations due to a seasonal freezethaw cycle (Webber 1978). Dominant vegetation in the region is characterized as a sedge/grass moss wetland (Walker et al 2005). Further details are available in the supplementary material is available online at stacks.iop.org/ERL/14/125018/mmedia.
2.2. Imagery data acquisition and pre-processing Imagery used in this study was procured through the Polar Geospatial Center at the University of Minnesota (table 1). Acquisition dates spanned 14 years, with the earliest Ikonos image acquired 18 July 2002 and the final WorldView-3 image acquired 24 July 2016. A digital elevation model ( figure 1(b)) was used to measure the relative elevation of features (Wilson et al 2013). The largest overlapping area in all images defined the study area extent with the northwest corner at 71°19′34.36″ N, 156°41′20.62″ W and the southeast corner at 71°16′9.02″ N, 156°31′26.44″ W. The total study area was 37.42 km 2 . Supplementary materials provide details on image pre-processing.

Meteorological data
Hourly air temperature and relative humidity data from the Barrow Atmospheric Baseline Observatory (BRW, NOAA), located north of the study area (71.3230°N, 156.6114°W), were used to understand climatic conditions during growing seasons. Mean growing season (June-August) air temperature, vapor pressure deficit (VPD), growing degree days (GDD, cumulative mean daily air temperature over 0°C), and thawing degree days (TDD, cumulative number of days with a mean air temperature over 0°C) were calculated. TDD and GDD were calculated for image acquisition dates allowing NDVI comparisons.

Field data collection
Field surveys were conducted in July 2018. Eight transects were surveyed with a total of 2971 m × 1 m plots. Transects ranged from 100 to 300 meters in length and plots were surveyed every 5 m. For each plot, plant types (grass, moss, lichen, shrub, open water, and forb), thaw depth, soil moisture, and canopy height were measured (see supplementary materials for descriptions of field techniques). Medians and standard deviations are reported as summary statistics. Coordinates of plots were recorded using a Trimble 5700 differential global positioning system (Trimble ® , USA).

Random forest vegetation model
Pixels with NDVI value below 0.1 (Gandhi et al 2015) were masked to remove snow and water, which have extremely high or low reflectance. Then we used histogram matching to color correct images with the 2016 WorldView-3 image as the base (results in supplementary material). Bright and dark objects were removed by the NDVI filter, therefore these objects Basins over the study area; colors represent the age of basins according to Hinkel et al (2003). The blue is the only young age basin in the study area (basin 1). There were eight medium age basins (purple) and eight old basins (yellow). Basin numbers are labelled according to age and then generally north to south. (c) Digital elevation model coverage through the study area. There is a gentle slope from the southern extent to the north with tributaries to the Elson Lagoon and basins as the lowest parts across the landscape. The background of panels b and c are a true-color RGB image taken 24 July 2016 by WV3.
will not skew histograms and allow for a better comparison of illumination. A tasseled cap transformation was done to produce three orthogonal polynomial combinations representing brightness, wetness, and greenness of objects (Kauth andThomas 1976, Yarbrough et al 2005).
Vegetation community classes from 'wet' and 'dry' locations were chosen based on functional relationships with carbon cycling (Sturtevant and Oechel 2013, Natali et al 2015, Davidson et al 2016b, Treat et al 2018 and being spectrally separable. An additional open water class was included after vegetation classification; corresponding to pixels removed by the NDVI filter. We utilized a random forest machine learning algorithm for the vegetation classification model because of its high accuracy when predicting vegetation classes (Chapman et al 2009, van Beijma et al 2014. Four spectral bands (blue, green, red, and nearinfrared (NIR)) and three tasseled cap variables were used to train the random forest model. Welch two sample t-tests were used to test for significant differences in spectral properties between vegetation communities (see supplementary material). The random forest model was then applied to all 15 images for predicting vegetation communities. Model accuracy was assessed using 30% of ground reference points not used in model training, against the classification of the latest 24 July 2016 WorldView-3 image.

Greening assessment
Greening of vegetation was measured using NDVI following many Arctic studies (Bhatt et al 2010, Epstein et al 2012, Bhatt et al 2013. Changes in NDVI were estimated prior to applying the NDVI threshold and histogram matching, since open water can contribute to browning (Raynolds and Walker 2016) which, we wanted to capture. Since NDVI changes most rapidly during green-up (Zhang et al 2018), differs seasonally and saturates (Suzuki et al 2001), only peak season (mid-July to early August) image acquisitions were used in the time-series analysis. This approach captured NDVI shifts linked to long term changes instead of seasonal differences due to acquisition dates. We also calculated a pixel level change to locate finely detailed space-time anomalies by regressing each pixel's NDVI over time.

Statistical analyzes
Non-metric multidimensional scaling (NMDS) was used to determine ecological separability between vegetation communities using the 'vegan' package (Oksanen et al 2018). Differences in the physical environmental parameters were determined using Welch two sample t-tests; the t-statistic (t), degrees of freedom (df ), and p-value (p) are reported for these tests. Homogeneity of variance was tested for all parameters before significance testing.
Changes in vegetation coverage, NDVI, and TDD over time were calculated according to Yue et al (2002) using the 'zyp' package (Bronaugh and Werner 2019), which is sensitive to potential autocorrelation in timeseries data. These changes were analyzed over the entire image and by separating the landscape into the following age based categories of basins (Hinkel et al 2003): young basins (0-50 years old), medium basins (50-300 years old), old basins (300-2000 years old), and the 'other' category that accounts for the remaining landscape (figure 1(a)). Before analyzing controls on NDVI, data were detrended to avoid temporal autocorrelation by taking the residuals of the least squares linear model between year and the given variable.
NDVI controls were analyzed using linear mixed effects models (LMEs) using maximum likelihood in the 'nlme' package (Pinheiro et al 2018) and the 'MuMIn' package (Barton 2019) was used to calculate the coefficient of determination of these LMEs. Two components of coefficients of determination are reported for each LME, marginal (R 2 m ) and conditional (R 2 c ), representing the coefficient of fixed effects and full model respectively. LMEs were compared to a null model with only year as a temporal random effect using analysis of variance. LMEs included single variable models with each variable as fixed variables (GDD, TDD, mean air temperature, and mean VPD) and a full multivariate model with TDD, mean air temperature, and mean VPD.
To determine pixel level greening and browning trends, the Theil-Sen slope method (Sen 1968, Theil 1992) was used for the rate of change and Yue et al (2002) time-series methods were again used to determine statistical significance. This method of trend detection performs better in remote sensing data than more basic least-squares regression (Fernandes and Leblanc 2005). NDVI of each pixel over time was regressed using the 'spatialEco' package (Evans 2019), resulting in rasters of Kendall's τ statistics (ranging from −1 to 1 with values closer to 0 indicating no significant trend) and Theil-Sen slope. All statistical analyzes were conducted using R v.3.5.2 statistical software (R Core Team 2018).

Vegetation classification
The random forest vegetation classification had 82% overall accuracy. Producer's accuracy was slightly higher for wet communities (86%) than dry communities (78%) and user's accuracy was comparable for wet (83%) and dry (82%) communities. Twenty-five bootstrapped repetitions of the training data for random forest model yielded an overall accuracy of 89%.

Vegetation composition by basin age
There was no noticeable change in vegetation community composition in the young basin with overall changes less than 1% (< 1 ha). Medium age basins by contrast showed a 16% increase in wet communities (85 ha) corresponding to similar losses in dry communities. Old basins had a 4% decrease in wet communities (26 ha) and similar increase in dry communities (24 ha). The remainder of the landscape showed little change with wet communities decreasing by 3% and dry communities increasing by the same amount. The only significant trend was open water in the 'other' landscape area (τ=0.77, p=0.002, slope=0.72 ha/year). All other trends were not significant.

Vegetation change at individual basins
Some significant trends in vegetation cover were observed within specific basins (table 2). Three medium age basins (3, 4 and 5) had significant increases in wet communities which, corresponded to losses in dry communities. Conversely, old basin number 10 gained dry community area but had no clear reciprocal change in open water or wet communities. Since there was only one young basin, these results do not differ from above finding no significant change.

Greening of the landscape
We found an overall greening trend over the entire study area (τ=0.65, p=0.01, NDVI increase of 0.01 yr −1 ) and within all basins, except for the young basin, which showed a slight increase with no significant change (p=0.077, τ=0.45, table 2). We also found that the NDVI trend with respect to basin age is spatially variable with faster rates of greening in medium basins than old or 'other'. Both wet and dry vegetation communities significantly increased in NDVI as well (τ=0.65, p=0.01 for both wet and dry communities). The rate of increase in NDVI was faster for wet communities (0.013±0.003 yr −1 ) than for dry communities (0.011±0.002 yr −1 ). Pixel level analysis showed similar results to those aggregated by basins or vegetation. Wet areas showed the largest change with rates of NDVI change around 0.05-0.08 y −1 . Dry ridges show NDVI changes closer to 0.01-0.04 y −1 (figure 4).

NDVI controls
TDD had the strongest correlation with NDVI (figure 5, F=21.5, p<0.001, R 2 m =0.56, R 2 c =0.77, NDVI increase of 0.004±0.001 d −1 ) and explained most of the variability (the full multivariate model was not significantly different from the individual model of TDD (p=0.88, likelihood ratio=0.25)). GDD also strongly correlated with NDVI, with a slightly lower correlation coefficient than TDD (F=20.7, p<0.001, R 2 m =0.52, R 2 c =0.77, NDVI increase of 0.0008±0.0002°C −1 ). As GDD was strongly collinear with TDD, they were not assessed together in any model (F=74.7, p<0.001, R 2 =0.84). Mean growing season air temperature and VPD were not correlated to NDVI (F=1.2, p=0.27 for air temperature and F=0.53, p=0.45 for VPD). The first day of the year with a mean air temperature above freezing   p-values show areas with significant changes were also largely located in wet areas including basins and water paths. Open water that was included in the analysis was the least significant and changed the least giving us confidence the observed changes were not due to sensor drift or artificial interferences. Figure 5. NDVI as a function of thawing degree days (TDD) across all images showing as the growing season length increases, NDVI increases. As both NDVI and TDD increased significantly over the course of the study, both were detrended in the LME testing the significance (p<0.001, R 2 m =0.56, R 2 c =0.77, N=15).
Year was kept as a random continuous variable with NDVI and TDD as fixed variables.
(i.e. first TDD) significantly shifted towards an earlier date at a rate of 1.07 d per year (τ=−0.47, p=0.02) but we found no evidence for changes in the last day above freezing (τ=−0.01, p=1).

Vegetation community characteristics
Vegetation community groupings based on functional differences in carbon cycling characteristics (Sturtevant and Oechel 2013, Davidson et al 2016b, Treat et al 2018, proved to be spectrally separable agreeing with prior studies (Lin et al 2012, Davidson et al 2016a. We were able to obtain ecological separation between vegetation communities as well using surveys of plant types (e.g. graminoid, moss, lichen, etc).
Soil moisture content was the main environmental factor driving the separation of the vegetation communities and dictating their geographical distribution. Our results agree with findings from past studies (Lin et al 2012) that showed vegetation communities near Utqiaġvik exist across a moisture gradient. Canopy heights differed between communities, with taller vegetation in wet communities and shorter vegetation in dry communities; canopy height has been used as a correlative metric for predicting above ground biomass in Arctic ecosystems (Berner et al 2018). This suggests that wet communities typically have higher biomass than dry communities and therefore the relationship between biomass and moisture may exist as well (see supplementary material). Thaw depth did not differ between communities but this does not necessarily imply equal carbon storage in the active layer beneath the vegetation due to usually higher concentrations of organic matter in wetter lowlands, i.e. higher biomass given equal soil volume (Ping et al 2008).

Vegetation community shifts at the landscape scale
We observed no shift in the overall cover of vegetation communities over the 14 year study. Measured shifts were around 1% and could be due to differences in water table, other seasonal differences, or classification errors. These results are in line with prior long-term studies of vegetation communities in this area showing 3% change over 60 years  when using high spatial resolution photography with satellite imagery (Lin et al 2012). However, Lin et al (2012) used four single images (1948, 1955, 1979 and 2008) to assess vegetation community change in their study which, did not account for intra-annual variability in standing water, a possible source of error. The seasonal progression in our data during 2010, where four images were acquired over two weeks (21 July-3 August), showed a seasonal drying. Changes in surface water during the season can erroneously translate into error in predicted cover of vegetation communities by the random forest model. This result highlights the importance of considering timing of acquisition of high spatial resolution imagery in analyzing interannual patterns and change over long timescales.
Vegetation changes were observed in more localized areas. The young basin vegetation communities did not change and consisted mostly of wet vegetation. Old basins, comprised largely of low center polygons, are thawing where ice wedge degradation occurs over time scales ∼60 years; the stage of degradation determines the degree to which the landscape is wet or dry (Liljedahl et al 2016). We saw no overall significant changes in vegetation communities through the timeseries but that does not mean degradation and shifts are not occurring. Polygon degradation effects could balance out if multiple stages are occurring simultaneously over the entire study area. According to our results, most changes occurred in medium and old age basins where medium age basins gained wet and lost dry communities and old age basins showed the opposite process (table 2, figure 3). Polygon degradation (Liljedahl et al 2016) in old basins could increase drainage to the lowest part of the nearby landscape, medium basins, due to interconnected troughs and could explain increases in wetness in these medium basins. This is supported by the change map showing the wetting in these medium basins and drying along the edges ( figure 3(a)). This is also supported by the fact that medium basins had the lowest elevations ( figure 1(b)). Evidence of erosion or rising water levels exists where edges of water paths have gains in wet communities and edges of larger tributaries show dry vegetation to open water transitions ( figure 3(b)). These results suggest that even though total vegetation community compositions have not changed over the entire study area, high spatial resolution imagery allows us to detect localized but relatively balanced shifts.

Greening of the landscape 4.3.1. Landscape scale greening trends
We found evidence of greening across vegetation communities and basins. Mean NDVI values of the whole image, split by vegetation type, basin age, and individual basins, displayed significant increasing trends in NDVI in all areas except the young basin. However, the young basin already had a relatively high NDVI (0.6 in 2016) and is at or above observed peaks for this ecosystem (Bhatt et al 2017, May et al 2017. Increases in NDVI can be due to changes in land cover type (Elmendorf et al 2012) or increased vegetation biomass (Hudson and Henry 2009). Due to the lack of change in total vegetation community composition, we suggest that the most likely explanation for the observed greening in this ecosystem is increased plant productivity and growth.

Vegetation community specific greening rates
Our results show vegetation communities greening at different rates, in agreement with prior studies (May et al 2017, Andresen et al 2018. This trend has been locally shown to stand across the North Slope of Alaska where wet communities green faster under warming (May et al 2017). Further, at a species level, wet communities have been observed to shift more than dry communities in this region over time, with sedge species replacing bryophytes (Villarreal et al 2012), which could aid in NDVI increases. When trends for each pixel were analyzed, spatial patterns became even more evident ( figure 4). Some wet areas including many medium age basins (basins 2, 4, 6 and 9) and lower lying water paths had some of the highest rates of greening, sometimes more than double that of dry ridges ( figure 4). The p-value map also supported that these same areas have more instances significant change. The faster green up in wet areas confirms studies that show increases in moisture, expected in this region under climate predictions (Zhang et al 2012), can enhance greening. Increases in moisture do not universally increase NDVI as shown in Raynolds and Walker (2016), where standing water depressed NDVI due to differential absorption of NIR and red irradiance by water. Lara et al (2018) further showed negative correlations between precipitation and NDVI.

NDVI controls and growing season length
Greening was most correlated to increases in TDD, agreeing with previous studies (Huemmrich et al 2010, Zeng et al 2011. This further suggests that NDVI is most influenced by the thawing date (Oberbauer et al 2013, Andresen et al 2018; however, this may not always be true. If there is a freezing day after initial thaw, an early thawing can have the opposite effect on growth (Oberbauer et al 2013), which causes browning (Phoenix and Bjerke 2016). Mean temperature and VPD did not influence NDVI in our study, possibly because temperature and VPD have been shown to correlate with max NDVI ( . SWI works for coarse-scale pan-Arctic studies that use maximum or time-integrated NDVI and have one value that represents the full growing season. However, in high spatial resolution investigations such as our study, a metric like TDD, that's measured on the acquisition date will better explain NDVI. Due to the tie between TDD and NDVI, we suggest that for high Arctic communities, if thawing continues to start earlier, then continued increases in NDVI are expected. NDVI and vegetation communities differed seasonally in 2010 when four images were acquired within two weeks. This emphasizes that temporal and spatial events may greatly affect NDVI values observed. We accounted for this by only using images acquired closer to peak growing season. While intraannual variability was observed in NDVI, the relationship with TDD, which also increased over the study period, gives us confidence in the overall greening trend. Further, these multiple images were used in time series analyzes including the variability in trend analysis. Ultimately more satellite tasking and image acquisitions would be ideal to continue monitoring this trend. More meteorological data, such as precipitation and evapotranspiration would be ideal to tease out larger scale influences on NDVI (Andresen and Lougheed 2015). Unfortunately, these datasets are not always available over longer time-series or spatial domains, especially within remote areas like the Arctic.

Conclusion
We found evidence of greening across the landscape, particularly in wet areas, and a balanced shift in vegetation communities. TDD was best correlated to the positive greening signal. There was no large shift in vegetation community assemblage, but localized changes were observed showing spatially variable wetting and drying. Our study emphasizes the increased ability of high-resolution remote sensing to analyze details in change detection analyzes in the Arctic. Specifically, we were able to observe that wet communities may respond to warming at a faster rate than drier communities.
Our study analyzes a small area in the Arctic but takes a detailed approach to understanding Arctic greening. Many studies have looked at large scale changes which are important in a global context but lack detailed explanations (Myers-Smith et al 2019). Only by understanding drivers of greening can we more accurately predict future vegetation productivity in the Arctic, which is currently of great importance given that experts cannot currently agree on the direction of change of the carbon balance (Abbott et al 2016). NA16SEC4810008) to WCO, from the European Union's Horizon 2020 research and innovation program under grant agreement No. 629727890 to WCO and DZ, and from the Natural Environment Research Council (NERC) UAMS Grant (NE/P002552/1) to WCO and DZ. We thank NOAA/ESRL for the meteorological data which is provided by NOAA/ ESRL BRW from their website https://esrl.noaa.gov/ gmd/obop/brw/. The data that support the findings of this study are openly available at DOI:10.25412/ iop.10324745.v1. Geospatial support for this work provided by the Polar Geospatial Center under NSF OPP award 1204263, and 1702797. This research was conducted on land owned by the Ukpeaġvik Inupiat Corporation (UIC) whom we thank for their support. Authors also thank the R developing team (R core team, Vienna, Austria) in creating the open source R statistical software. The data that support the findings in this study are openly available.

Data availability statement
The data that support the findings of this study are openly available at https://doi.org/10.25412/iop. 10324745.v1.