Draft Remote sensing evaluation of High Arctic wetland depletion following permafrost disturbance by thermo-erosion gullying processes

Complete List of Authors: Perreault, Naïm; Universite du Quebec a Trois-Rivieres, Sciences de l'Environnement; Centres d'etudes nordiques, Universite Laval Levesque, Esther; Université du Quebec à Trois-Rivières, Sciences de l'Environnement; Centre d'etudes nordiques, Universite Laval Fortier, Daniel; Universite de Montreal, Geographie; Centres d'etudes nordiques, Universite Laval Gratton, Denis; Universite du Quebec a Trois-Rivieres, Sciences de l'Environnement Lamarque, Laurent; Universite du Quebec a Trois-Rivieres, Sciences de l'Environnement; Centres d'etudes nordiques, Universite Laval


Introduction
Wetlands, accounting for up to 80% of the total surface area in some Arctic regions (Pienitz et al. 2008;Rautio et al. 2011), represent prominent ecosystems for tundra biodiversity and livelihood of local communities.They offer critical staging and nesting homes to many birds (Hughes et al. 1994;Lepage et al. 1998;CAFF 2013), and provide many herbivores refuges from predators and preferred grazing grounds with abundant graminoid species (Henry 1998;Kristensen et al. 2011;Doiron et al. 2014).Northern wetlands also play a significant role in global gas exchanges by storing large quantities of carbon due to slow organic matter decomposition rates (Tarnocai et al. 2009;Hugelius et al. 2014;Bouchard et al. 2015) and by producing 5 to 20% of total natural methane emissions to the atmosphere (Olefeldt et al. 2013;McEwing et al. 2015).
Northern wetlands are highly sensitive to current climate warming.Their sustainability relies on the integrity of perennially frozen ground that provides permanent or seasonal water supply and prevents vertical percolation (Woo 2012;Natali et al. 2015).However, observations over the last decades indicate that permafrost thawing and active layer deepening, closely related to the ratio of ground ice content and ground thermal conditions (Burn and Kokelj 2009;Kanevskiy et al. 2013;Rudy et al. 2017), have intensified with increasing temperatures (Romanovsky et al. 2010;Rowland et al. 2010;Liljedahl et al. 2016;Vincent et al. in review).The resulting changes in hydrological conditions can have profound consequences for biochemical cycles and ecological communities (see Wrona et al. 2016 for a review).The rapid development of new drainage channels may for instance reduce the areal extent of wetlands globally (Klein et al. 2005;Fortier et al. 2007;Avis et al. 2011;Godin et al. 2012a).
Continuous permafrost zones with widespread ice-wedge networks are particularly prone to instability given the large amount of ice found near the top of permafrost (Jorgenson et Kanevskiy et al. 2013;Liljedahl et al. 2016), and wetlands in these regions are highly susceptible to thermo-erosion gullying (Kokelj and Jorgenson 2013;Kokelj et al. 2014).This permafrost geosystem feature results from thermal and mechanical ice-wedge degradation subsequent to ice melt through convective heat flow between running water and ice/frozen ground (Fortier et al. 2007;Morgenstern et al. 2013;Godin 2016).Thermo-erosion gullying has been reported from many sites throughout the circumpolar North and the Antarctic (Levy et al. 2008;Grosse et al. 2011;Günther et al. 2013;Liljedahl et al. 2016).On Bylot Island, Nunavut, thermo-erosion processes led to substantial permafrost tunneling, gullying and subsequent mass wasting with an eroded area estimated at nearly 158 000 m 2 within 10 to 15 years after gully inception (Fortier et al. 2007;Godin and Fortier 2012a;Godin et al. 2014;Veillette et al. 2015).In consequence, adjacent wetlands were severely affected by rapid drainage.The decrease in soil moisture and thaw-front depth in disturbed low-centered polygons culminated in a relatively rapid transition from wet to mesic plant communities (Godin et al. 2016;Perreault et al. 2016).
Satellite data have been used increasingly in the Arctic for various geomorphological purposes such as mapping polygonal tundra landscape components (Skurikhin et al. 2013), detecting ice-wedge degradation (Liljedahl et al. 2016), assessing permafrost coast erosion (Obu et al. 2016), understanding effects of flooding following storm (Lantz et al. 2015), and monitoring permafrost thaw dynamics (Ulrich et al. 2014;Beck et al. 2015;Günther et al. 2015;Liu et al. 2015, Saruulzaya et al. 2016).Techniques such as InSAR (Interferometric synthetic aperture radar) as well as airbone and terrestrial laser scanning have been used to document detailed topographical changes (Chen et al. 2013;Hubbard et al. 2013;Liu et al. 2014;Wolfe et al. 2014), thermokarst (Barnhart and 2013), active layer thickness (Gangodagamage et al. 2014;Schaefer et al. 2015;Widhalm et al. 2017) and freeze-thaw cycles (Daout et al. 2017).While remote sensing studies have also characterized the D r a f t distribution of thermo-erosion landforms (Morgenstern 2012;Belshe et al. 2013;Godin et al. 2014), they have yet to examine the impacts of these processes on surrounding vegetation.
Using the satellite-derived normalized difference vegetation index (NDVI) coupled with field surveys of plant community composition, our objective was to test the feasibility of remotelysensed data for quantitatively assessing the extent of landscape-scale permafrost disturbance on the depletion of adjacent wetlands.

Study area
This work was conducted in the Qarlikturvik valley, on Bylot Island, Nunavut, Canada (73°09'N, 79°57'W; Fig. 1).Located on the southern plain of the island, which is part of the Sirmilik National Park, the valley has a surface area of 75 km 2 and is ~ 17 km long and 4-5 km wide.It is bounded to the north and south by plateaus (~500 m a.s.l.) and to the east by the C-79 and C-93 glaciers (Inland Waters Branch 1969).A proglacial river runs in a glaciofluvial outwash plain through braided channels which prograde into the Navy Board Inlet to the west forming a wide delta.Streams and rills coming from sub-perpendicular valleys and from alluvial fans can be temporarily connected to the river.The plain is bordered on both sides by low gradient terraces which have developed as the accumulation of ~ 4-5 m of ice-rich peat mixed with aeolian and alluvial silts and sands following glacier retreat during the Holocene (Fortier and Allard 2004;Fortier et al. 2006).Networks of syngenetic ice wedges formed during the Late Holocene and resulted in the development of a polygonalpatterned ground which features high-and low-centered polygons, polygon ponds and thaw lakes (Hughes et al. 1994;Fortier and Allard 2004;Bouchard et al. 2015).In this area, mean air temperature averaged -14.7 °C at 24 m a.s.l. for the 1995-2015period (CEN 2016)).Spring and fall temperatures have risen by 2.8 and 4.3 °C, respectively, over the past 35 years D r a f t 6 (Gauthier et al. 2011), while the annual cumulative number of thawing degree-days has increased from 381 in 1989 to 521 in 2011 (Gauthier et al. 2013).
Field surveys were specifically conducted around three gullies that were selected among the 36 reported by Godin and Fortier (2012a) in the low-centered polygon landscape (Fig. 1b and 2).These gully networks, closely monitored since 1999, originated from the collapse of underground tunnels which followed the infiltration of snowmelt water into cavities formed in the frozen active layer and permafrost (Fortier et al. 2007;Godin and Fortier 2010).The gullies R08p and R06, respectively 805 and 717 m long, were characterized by ongoing thermal-erosion while the gully RN08, 180 m long, has been very weakly active in recent years (Godin andFortier 2012a, Veillette et al. 2015).Vegetation surveys were carried out during the 2009 and 2010 plant growing seasons at 197 sampling sites (n = 95, 81 and 21 around gully R08p, R06 and RN08, respectively) which were selected using a stratified D r a f t random sampling design to encompass intact wet and mesic sites as well as disturbed polygons (n = 62, 48 and 87 for wet polygons, mesic environments and disturbed polygons, respectively; see Perreault (2012) and Perreault et al. (2016) for a detailed description of the methods and habitats).The position of each sampling site was recorded using a 5-m accuracy GPS receiver (Garmin eTrex Venture HC).Plant species richness was highest in disturbed polygons where species of wet polygons and mesic environments were present (total species richness = 36, 47 and 54 in wet polygons, mesic environments and disturbed polygons, respectively).Species abundance in disturbed polygons represented a middle-range state between wet and mesic conditions, and were thus covered by both hydrophilic species such as Carex aquatilis and Dupontia fisheri and mesic species such as Salix arctica and Arctagrostis latifolia (Table 1).The moss carpet was mainly made of living Drepanocladus spp. in wet polygons, living Aulocomnium spp. in mesic environments and dried (i.e., dead) species in disturbed polygons (Table 1).

Remote sensing data analysis
We used a high-resolution satellite image of the Qarlikturvik valley acquired by the GeoEye sensor (GeoEye Inc., Herndon, VA, USA) on September 2, 2010 at 17:40 GMT.This sensor covered an area of 14 km long x 19 km wide with a 0.5 m spatial resolution panchromatic band and four multispectral bands between 450 and 920 nm of 2 m spatial resolution (GeoEye Elevating Insight, GeoEye Inc.).The sun elevation and azimuth angles at the time of image acquisition were 24.37 and 185.86 degrees, respectively, whilst the sensor look angle and elevation were 81.40 and 220.18 degrees.Cloud cover represented less than 5% of the scene.
The image was delivered geometrically corrected in UTM coordinates (Zone 17 N in the world geodetic system, WGS 84).

Image analysis
Image processing and spatial data analyses were conducted using the ENVI image analysis (Exelis Geospatial Systems, Rochester, NY, USA) and the ArcGIS version 9.2 (Esri, Redlands, CA, USA) software packages.We first pre-processed the satellite image by applying a geographical correction based on the differential global positioning system (DGPS) coordinates of the studied gullies (Godin andFortier 2010, 2012a,b) and by converting the raw digital numbers of red and near-infrared bands to Top Of Atmosphere (TOA) reflectance values (Goward et al. 2003).Following cloud masking, a subset scene of the valley encompassing areas of interest was obtained by manually digitizing boundaries.
Spectral bands were merged with the panchromatic image to create 0.5 m standard pansharpened multispectral images, i.e., true and false color composite images.This step ensured the accurate positioning of sampling sites, as gully and sampling site GPS positions were projected onto the GeoEye satellite image.Analysis areas -one per sampling site -were identified following the thorough description of plant communities and environmental conditions from Perreault et al. (2016).Each analysis area, which excluded polygon rims, incorporated as many pixels as possible (2 m x 2 m per pixel) in order to have a good representation of its corresponding sampling site.A minimum of 6 pixels was integrated per analysis area, for a total of 5,406 selected pixels (1,832 pixels for wet polygons, 931 pixels for mesic environments and 2,643 pixels for disturbed polygons).Analysis areas were finally digitized as Regions Of Interest (ROI) in ENVI and used to create NDVI values for each habitat type (i.e., for wet polygons, mesic environments and disturbed polygons; Fig. 3).

NDVI calculation
The NDVI was determined as follows: NDVI = (P nir -P red ) / (P nir + P red ), where P nir and P red respectively represent the spectral reflectance in the near-infrared (from the plant canopy) and D r a f t 9 in the red (where chlorophyll absorbance is maximal; Raynolds et al. 2008).This global vegetation index, which measures the vegetation relative greenness, yields values between -1 and 1, with those close to -1 representing areas devoid of vegetation and those close to 1 referring to dense vegetation areas (Rouse et al. 1973;Tucker and Sellers 1986;Pettorelli 2013).A proportionate stratified random sampling method (by vegetation class) was used to select 80% of the analyzed pixels.NDVI statistics (means, medians, quartiles and percentiles) were calculated for each type of habitat based on this selection.

Image classification and mapping of affected areas
A simple threshold-based classification technique was performed on the NDVI image to determine class thresholds, using the 80% of the analysis area pixels which were selected beforehand.The residual 20% were used to assess the accuracy of the classification.A standard confusion matrix was used to assess the accuracy between class thresholds and habitat types.Before the majority filter was applied, we calculated the overall accuracy (i.e., the proportion of correctly classified pixels) for the classified image and the Kappa coefficient, which indicates the percentage of correctly allocated pixels (Thomas et al. 2003;Congalton & Green 2009).Because data following pixel-based classification show so-called salt-pepper noise (Rongqun and Daolin 2011), a 3 x 3 pixel majority filter was applied to the classified image to improve type unification.Low NDVI areas located around the gullies were vectorized considering eight adjacent pixels for the connectivity and a minimum of 20 pixels per region.Vectorized images were exported to ArcGIS and these areas surrounding gullies were retained for the determination of disturbed areas.For this study, the contours of the gullies were mapped using the gully cartography from Godin andFortier (2010, 2012a).

Habitat NDVI
The NDVI analysis emphasized a marked difference between undisturbed sites -from either wet or mesic environments -and disturbed areas, with greater NDVI values for wet polygons and mesic environments compared to disturbed polygons (Fig. 4).The NDVI value of 0.27, which corresponded to both the lower quartile of mesic environments and the upper quartile of the disturbed polygons, separated the greater vegetation cover of wet and mesic sites from the lower vegetation cover of disturbed areas.There was a strong overlap in NDVI values between wet and mesic sites, and the NDVI value of 0.33, corresponding to the median of wet polygons and the upper quartile of mesic environments, was used to differentiate these two baseline habitats.

Image classification
The NDVI values ranged from -0.73 to 0.74 for the processed scene (Fig. 5).The final satellite-derived NDVI image was classified into five classes with the last three related to the habitat types of interest: (i) class 1 (-0.73 to -0.1) represented water bodies; (ii) class 2 (-0.1 to 0.15) represented bare grounds such as proglacial braided river bars and aeolian sand covers; (iii) class 3 (0.15 to 0.27) represented areas poorly covered by vegetation surrounding alluvial fans and drainage systems; (iv) classes 4 (0.27 to 0.33) and 5 (0.33 to 0.74) were associated with mesic and wet environments, respectively (Table 2).The classification was supported by an overall accuracy of 61.7% and a Kappa coefficient of 0.41.When wet polygons and mesic environments were pooled together, the overall accuracy increased to 79.5% with a Kappa coefficient of 0.59, providing another evidence of the similarity in NDVI trends between these two types of habitats.

Mapping of affected areas
In the three studied cases of thermo-erosion gullying and subsequent permafrost mass wasting and wetland drainage, disturbed areas were directly adjacent to new drainage systems formed by the gullies and distributed all along their length on both sides of gully channels (Fig. 6).
Nearby sites had higher NDVI values, highlighting the dominance of wet polygons and mesic environments (Fig. 5).The overall wetland area affected by the three gullies was approximately 95,430 m 2 , corresponding to 56 m 2 of disturbed area per linear metre of gullying.Considering that wetlands cover 34% of the Qarlikturvik valley total area (19 km 2 out of 56 km 2 ; Masse 1998), thermo-erosion gullying along these three gullies only has drained 0.5% of the total wetland area within a decade.The extent of disturbed areas due to drainage strongly differed among gullies with the greatest disturbances found around the active R08p and R06 gullies (Table 3, Fig. 2).

Discussion
Having previously reported that the initiation of thermo-erosion gullying has led to the relatively rapid replacement of hydrophilic plant species by those indigenous to mesic environments (Perreault et al. 2016), we sought here to assess and quantify the extent of disturbed wetlands by permafrost disturbance using fine-scale remotely-sensed data.Climatedriven changes to surface moisture and wind regimes at the century timescale have been documented for the study site over the Late Holocene (Fortier et al. 2006).However, drastic and rapid changes to surface wetness that are linked to the development of thermo-erosion gullies are unprecedented (Fortier et al. 2007;Godin et al. 2014) While NDVI values for wet and mesic sites overlapped significantly due to similar total live plant cover present in these habitats at the time the image was taken, the analysis we carried out from GeoEye data allowed for a good discrimination between disturbed and undisturbed wetland and mesic vegetation types.The areas affected by thermo-erosion processes were restricted to the gully edges, illustrating the strong control of ice-wedge networks in shaping polygon landscape hydrology (Liljedahl et al. 2012;Vonk et al. 2013;Liljedahl et al. 2016).Water drainage subsequent to thermo-erosion gullying has induced a significant decrease in soil moisture of disturbed polygons where wetland vegetation has not been able to persist and is being replaced by mesic plant species (Godin et al. 2016;Perreault et al. 2016).Considering that the stable state of mesic vegetation has not yet been reached, disturbed polygons currently encompass both wet and mesic plant species, with a total vascular plant cover similar to that in undisturbed sites (approximately 25%).The low NDVI characterizing these zones is therefore due to the higher proportion of dried mosses, as increases in NDVI values are tightly linked to the green moss layer thickness and moss water content (Douma et al. 2007;Engstrom et al. 2013).
The three studied gullies have already drained within 10 years ca.95,000 m 2 of wetland areas, distributed along 1,702 m of gullies.Such disturbance, which takes into consideration both the wetland drainage per se and the gully areas, that also represent a net loss of habitat since no stable vegetation has established to date, amounts to a loss of 0.5% of the total wetland area of the valley.With an average of 56 m 2 of affected area per linear meter of gullying, thermo-erosion-induced drainage of wetlands is worrisome for the near-future sustainability of these systems, especially since (i) R08p and R06 gullies have shown rapid thermo-erosion progression (i.e., 23 m y -1 ) (Godin and Fortier 2012a); and (ii) a total of 36 D r a f t 13 gullies have been characterized in the Qarlikturvik valley since 1999 for an overall eroded area already estimated at 158,000 m 2 (Godin andFortier 2010, 2012a;Godin et al. 2014), corresponding to a net habitat loss of 0.8%.Thermo-erosion processes thus represent a major threat to wetland carrying capacity, which will be even more pronounced when accounting for the wetland areas depleted by drainage around all the 36 gullies in the valley.Indeed, these wetlands support graminoids such as Dupontia fisheri and Eriophorum scheuchzeri that are preferred food resources of greater snow geese (Cadieux et al. 2008;Doiron et al. 2014).
Between 1990 and 2012, snow geese, which can number up to 650 breeding pairs in the valley, consumed on average 30% of the above-ground biomass that was annually produced in the wetlands (Gauthier et al. 2013).Their nesting success may also be impacted by wetland depletion, as a decrease in water availability implies increased distances to water sources and therefore reduced ability to defend nests against predators such as the Arctic fox (Lecomte et al. 2008).The rapid decrease of wetland habitats could also affect other tundra species (Berteaux et al. 2016); for example, the community structure of Arctic arthropods is strongly dependent, even at small spatial scales, on vegetation composition changes (Hansen et al. 2016).Although wetlands are common throughout the valley, it is worth noting that their distribution is heterogeneous, hence, the consequences of their depletion due to permafrost disturbance may be exacerbated in place where they are scarcer.Finally, thermo-erosion gullying is likely to alter greenhouse gas emissions at large scales.Methane emissions which are strongly regulated by wetland habitats (Olefeldt et al. 2013;Tveit et al. 2015), may decrease as a result of adjacent wetland drainage.On the contrary, carbon release may increase given wetlands generally act as carbon sinks (Bouchard et al. 2015;Treat et al. 2015), whilst the export of sediments via the gully channels would provide greater fluxes of particulate, dissolved carbon and nutrients (Godin et al. 2014).

D r a f t 14
The impact of thermo-erosion processes strongly varied between gullies.The greatest disturbances found along the R08p and R06 gullies can be explained by the substantial lengths of the main channels and the presence of secondary channels (Godin and Fortier 2012b).The R08p gully has been expanding simultaneously on six distinct secondary channels while the R06 gully main channel progression has been among the fastest thermoerosion development rates in the valley (Godin and Fortier 2012a).Comparatively, the RN08 gully, which has not been active recently, is significantly shorter and does not have secondary axes (Godin and Fortier 2012a;Godin 2016).Given that these three gullies have all developed on peaty sandy to silty deposits, differences in erosion activity features among them may not be related to the sedimentary environments, although Veillette et al. (2015) found that erosion patterns of gullies are indirectly related to the type of soils and sizes of ice wedges where they are formed.
Satellite remote sensing data have played a prominent role in the observation and interpretation of spatial and temporal land cover changesover large spatial extents in the Arctic (Stow et al. 2004;Muster et al. 2013;Jorgenson and Grosse 2016;Park et al. 2016).As a result, NDVI analyses have been applied increasingly in studies of Arctic ecology for various purposes such as reporting climate-induced plant community composition changes (Rudy et al. 2013;Beck et al. 2015;Pattison et al. 2015), quantifying permafrost degradation (Fraser et al. 2011;Belshe et al. 2013), determining deglaciation effects on tundra vegetation distribution (Raynolds and Walker 2009) or predicting peaks of above-ground plant biomass (Doiron et al. 2013).In this study, undisturbed baseline wet and mesic communities were adequately separated from disturbed areas; hence we suggest that finer-spatial resolution analyses of NDVI trends can be effective for monitoring the drying or recovery of wetlands following inception of thermo-erosion gullying.This approach is likely to yield higher accuracy provided the satellite data are acquired earlier in the growing season, e.g., when wet D r a f t 15 polygons are more water-saturated, to better differentiate these two baseline vegetation types.
Differences between vegetation types are indeed greater when NDVI is calculated around the peak growing season when there is maximum green leaf area (Riedel et al. 2005;Bratsch et al. 2016), however, in our case satellite data were acquired at the end of summer when wet polygons and mesic environments show similar values of total above-ground biomass (50.5 g.m -2 ± 2.8 in wetlands and 44.2 g.m -2 ± 6.8 in mesic tundra; Legagneux et al. 2012).
Additional strategies to increase differentiation between non-disturbed wet and mesic vegetation communities may include considering plant structure and variations in soil moisture (Laidler and Treitz 2003;Laidler et al. 2008;Atkinson and Treitz 2012) as well as using other high-resolution remote sensing platforms.In this regard, Unmanned Aerial Vehicles, which are relatively inexpensive and have successfully been tested for monitoring temporal dynamics of landslide and tundra vegetation (Turner et al. 2015;Fraser et al. 2016), as well as air borne laser scanning and Interferometric Synthetic Aperture Radar used to monitor topographic changes and active layer thickness (Gangodagamage et al. 2014;Liu et al. 2014;Schaefer et al. 2015;Windham et al. 2017), represent complementary techniques full of potential.

Conclusion
We found that the transformation of northern polygon landscapes and shifts in plant communities induced by thermo-erosion gullying processes can be accurately estimated by remotely-sensed data.We presented a novel quantification of wetland disturbances that are induced by thermo-erosion gullying processes using fine-scale NDVI analyses.These results open up new opportunities to examine the causes of spatial heterogeneity in tundra response to warming.They will allow for the monitoring of the rapid changes associated with wetland disruption in the High Arctic as well as for the evaluation of the contributions of permafrost D r a f t 16 degradation to global greenhouse gas emissions.This approach may also be useful for studying the long-term effects of slower hydrological regime changes on the transition between wet and mesic environments, as well as for quantifying the feedbacks of vegetation succession on the stabilization of disturbed permafrost terrains.1. Mean cover (%) of the most commun taxa present in the three types of habitat (wet polygons, mesic environments and disturbed polygons) studied along thermo-erosion gullies in the Qarlikturvik valley, Bylot Island, Nunavut (modified from Perreault et al. 2016).<<: cover < 0.01%; <: cover < 0.1%.Total covers (%) of both common and rare taxa are also indicated.
Table 2. Confusion matrix showing results of the classification accuracy.Proportion (%) of residual pixels per NDVI class for each of the three habitat types (wet polygons, mesic environments and disturbed polygons) encountered around thermo-erosion gullies in the Qarlikturvik valley, Bylot Island, Nunavut.NDVI classes 1 and 2 for which the analysis was not suitable are not shown.
Table 3. Extent of affected areas around the three studied gullies of the Qarlikturvik valley, Bylot Island, Nunavut.Length and area of gullies R08p and R06 were retrieved from Godin and Fortier (2012) while length of gully RN08 was measured using the GeoEye-1 panchromatic band and the measurement tool in ENVI.The asterisk refers to the ratio of affected areas and gully lengths.Classes 1 and 2 represent water bodies and bare ground, while classes 3 to 5 are related to the three studied habitat types, i.e., disturbed polygons (in yellow), mesic environments (in light green) and wet polygons (in dark green).Disturbed areas were estimated via NDVI values derived from pseudo-colour images.Hatched areas represent the gully topographic contours extracted from Godin andFortier (2010, 2012a).Red lines outlining the yellow colour of class 3 indicate the extent of the thermo-erosion-induced disturbed areas using both gully cartography from Godin and Fortier and the NDVI analyses.The black arrows indicate the general water flow direction.

D r a f t
Table 1.Mean cover (%) of the most commun taxa present in the three types of habitat (wet polygons, mesic environments and disturbed polygons) studied along thermo-erosion gullies in the Qarlikturvik valley, Bylot Island, Nunavut (modified from Perreault et al. 2016).<<: cover < 0.01%; <: cover < 0.1%.Total covers (%) of both common and rare taxa are also indicated.

Habitats
, and may imply a fundamental change in the hydrological cycle in response to changing snowmelt water run-off high-resolution NDVI may therefore help to understand the magnitude of change and to follow vegetation changes in time as the climate warms.

Figure 1 .
Figure 1.Location of the study area.(a) The Qarlikturvik valley (black circle), Bylot Island (black rectangle), Nunavut, Canada.(b) False color GeoEye-1 image of the valley (taken on September 2, 2010 at 17:40 GMT) and the three studied thermo-erosion gullies, located on syngenetic ice-wedge polygon terraces bordering a glacio-fluvial outwash plain, terminating in a delta prograding into Navy Board Inlet.IWPT refers to the ice-wedge polygon terraces.

Figure 2 .
Figure 2. Oblique aerial views of the three thermo-erosion gullies studied in the Qarlikturvik valley, Bylot Island, Nunavut (photos taken on July 23, 2013 by D. Fortier).(a) Gully R08p (805 m long), (b) gully R06 (717 m long) and (c) gully RN08 (180 m long).Each ice wedge polygon located along the gullies is ca.15-20 m wide.The gullies follow closely the network of ice wedges in the ground, which explains the angularity of the gully channels.The widths of the gully cross-sections vary along their length but generally range from 2 to 20 m.The depths of the gully cross-sections are a function of the base level of the gully channel at the gully outlet.The depths increase upstream according to the topographic gradient of the terrace to reach over 4 m at the gully head.The black arrows indicate the general water flow direction.

Figure 3 .
Figure 3. Location of the 197 analysis areas around (a) gully R08p, (b) gully R06 and (c) gully RN08 (wet polygons (n = 62); mesic environments (n = 48); disturbed polygons (n = 87)).The types of habitat were characterized by geomorphological and vegetation field surveys, while the size of each site was estimated via image analysis.

Figure 4 .
Figure 4.For each type of habitat (WP = wet polygons, ME = mesic environments, DP = disturbed polygons), (a) NDVI values calculated based on 80% of the analysis area pixels (n = 1,465 pixels from 62 sites for wet polygons; n = 745 pixels from 48 sites for mesic environments; n = 2,127 pixels from 87 sites for disturbed polygons); the 10 th percentile,

Figure 5 .
Figure5.GeoEye-1 NDVI image for the Qarlikturvik valley, Bylot Island, Nunavut.Classes 1 and 2 represent water bodies and bare ground, respectively, while classes 3 to 5 are related to the three studied habitat types, i.e., disturbed polygons (in yellow), mesic environments (in light green) and wet polygons (in dark green).IWPT refers to the ice-wedge polygon terraces.

Figure 6 .
Figure6.Extent of thermo-erosion impact along (a) gully R08p, (b) gully R06 and (c) gully RN08.Classes 1 and 2 represent water bodies and bare ground, while classes 3 to 5 are related to the three studied habitat types, i.e., disturbed polygons (in yellow), mesic environments (in light green) and wet polygons (in dark green).Disturbed areas were estimated via NDVI values derived from pseudo-colour images.Hatched areas represent the gully topographic contours extracted fromGodin and Fortier (2010, 2012a).Red lines outlining the yellow colour of class 3 indicate the extent of the thermo-erosion-induced disturbed areas using both gully cartography from Godin and Fortier and the NDVI analyses.The black arrows indicate the general water flow direction.

Table 2 .
Confusion matrix showing results of the classification accuracy.Proportion (%) of residual pixels per NDVI class for each of the three habitat types (wet polygons, mesic environments and disturbed polygons) encountered around thermo-erosion gullies in the Qarlikturvik valley, Bylot Island, Nunavut.NDVI classes 1 and 2 for which the analysis was not suitable are not shown.

Table 3 .
Extent of affected areas around the three studied gullies of the Qarlikturvik valley, Bylot Island, Nunavut.Length and area of gullies R08p and R06 were retrieved fromGodin and Fortier (2012)while length of gully RN08 was measured using the GeoEye-1 panchromatic band and the measurement tool in ENVI.The asterisk refers to the ratio of affected areas and gully lengths.