Nonlinear thermal and moisture response of ice-wedge polygons to permafrost disturbance increases heterogeneity of high Arctic wetland

Low-center polygonal terrains with gentle sloping surfaces and lowlands in the high Arctic have a potential to retain water in the lower central portion of ice-wedge polygons and are considered high-latitude wetlands. Such wetlands in the continuous permafrost regions have an important ecological role in an otherwise generally arid region. In the valley of the glacier C-79 on Bylot Island (Nunavut, Canada), thermal erosion gullies were rapidly eroding the permafrost along ice wedges affecting the integrity of the polygons by breaching and collapsing the surrounding rims. Intact polygons were characterized by a relative homogeneity in terms of topography, snow cover, maximum active layer thaw depth, ground moisture content and vegetation cover (where eroded polygons responded nonlinearly to perturbations, which resulted in differing conditions in the latter elements). The heterogeneous nature of disturbed terrains impacted active layer thickness, ground ice aggradation in the upper portion of permafrost, soil moisture, vegetation dynamics and carbon storage.


Introduction
Ice-wedge polygons as non-sorted patterned ground terrain type (Ballantyne, 2007) are widespread in the continuous permafrost zone characterizing the high Arctic (Black, 1976;Mackay and MacKay, 1974).High-latitude valleys (Dostovalov and Popov, 1963) and, more generally, arctic lowlands are prone to the formation of low-center polygon fields, which often typify a poor drainage and predominantly wet terrain (French, 2007;Zoltai and Tarnocai, 1975).Lowcenter polygons in arctic wetlands often develop ponds in their centers (Black, 1976) and are considered to have an important ecological role enabling suitable habitats for various macrospecies, either plants or animals (Gauthier et al., 1996(Gauthier et al., , 2005(Gauthier et al., , 2011(Gauthier et al., , 2013;;Jia et al., 2003;Massé et al., 2001;Myers-Smith et al., 2011;Woo and Young, 2012;Walker et al., 2004).Ponds and lakes that can act as carbon sinks or greenhouse gas sources (Tranvik et al., 2009;Bouchard et al., 2015) commonly form in such polygonal wetlands, therefore playing an important role in the global carbon cycle.
The main input of water in the polygons is snowmelt, while evapotranspiration contributes to water loss.The integrity of the polygon rims, the depth of the active layer and the lateral flow within the active layer further contribute to the dynamics of the storage (Helbig et al., 2013;Liljedahl et al., 2012).Distinct terrain units in a permafrost landscape (e.g., polygons, ponds, hummocky terrain) can develop and coexist over very short distances (a few meters between each terrain unit) with perceptible differences among physical characteristics such as active layer depth, ground temperature and hydrological conditions (Boudreau and Rouse, 1995).A survey at the polygon scale demonstrated significant differences in the water balance, active layer depth and plant distribution among single polygons located in similar lowland mires but a few kilometers away from each other (Minke et al., 2009), further stressing intra-site variability for almost identical terrain units.Another study demonstrated the vari-E.Godin et al.: Nonlinear thermal and moisture response of ice-wedge polygons ability in active layer depths in a nondisturbed and uniform 5 km 2 grid with high-and low-center polygons (Gangodagamage et al., 2014).
Generally, permafrost disturbance and degradation exert a range of impacts on the affected area, such as, but not limited to, ecosystem shifts and consequent changes in topography and mass transfer (Jorgenson and Osterkamp, 2005;Godin et al., 2014).When disturbance occurs in wet polygons, water could be transferred from centers to collapsed troughs when polygons evolve from low center to high center.The presence of dead tussock in the high-center polygon provides evidence of previous, more moist conditions (Jorgenson et al., 2006).Furthermore, subsidence in the scale of a few centimeters in a very low gradient landscape can exert a major impact on the hydrology of a watershed (Liljedahl et al., 2012).Disturbances in polygon fields, such as by thermo-erosion gullying of ice wedges, can occur very rapidly, with severe and immediate impacts on the terrain hydrology and ecological integrity (Fortier et al., 2007).On Bylot Island, one single gully eroded hundreds of polygon ridges over a period of 14 years and clear changes were observed in polygon moisture and vegetation conditions (Godin et al., 2014;Perreault, 2012;Perreault et al., 2015).Changes in cover and surface aspects are obvious near the gullied area and vary between eroded and intact polygons.Physical differences, between a closed-rim polygon (intact) and an open one (due to gullying) located only a few meters away, can induce very different plant communities in their respective center in response to changing moisture and active layer conditions (Perreault et al., 2015).Vegetation changes in wetlands have direct implications on the food web; for instance, in the Baffin area (and on Bylot Island), large avian herbivore populations rely on graminoids for support and the most adequate land unit (i.e., wetland) for this type of vegetation is restricted (Gauthier et al., 1996).Few detailed studies were conducted at the polygon scale with a specific focus on their thermal or moisture evolution through time while considering their microtopography (Helbig et al., 2013;Minke et al., 2009).
This paper discusses the implications of the polygon ridges' erosion by gullying and how this changed the microtopography, near-surface moisture conditions, active layer dynamics and vegetation distribution.This will be examined by comparing the aforementioned factors between undisturbed and eroded polygons.

Study site
The study site is located in the valley of the glacier C-79 (locally known as Qalikturvik valley) in the western part of Bylot Island in the eastern Canadian Arctic archipelago (73 • 09 N, 79 • 57 W) (Fig. 1a).The 4 km wide by 17 km long glacierized valley (Fig. 1b) is drained via a proglacial braided river connected to the Navy Board Inlet (Inland Water Branch, 1969).The pro-glacial river is bordered by a syngenetic ice-wedge polygon terrace composed of layers of loess and poorly decomposed peat (Fig. 2) (Fortier and Allard, 2004).Polygons in the valley have either a high center or low center and thousands of lakes and ponds are scattered over the terrace.Groups of high-center polygons are often surrounded by linear ponds that formed over ice wedges following collapse of polygon ridges (Bouchard et al., 2015).Several thermo-erosion gullies developed through the valley and increased the hydrological connectivity between the valley walls and the river (Godin and Fortier, 2012a;Godin et al., 2014).(CEN, 2014).BY-LOTPD is a reference ground temperature monitoring site installed in a low-center polygon (Allard et al., 2014).The gully R08 is located on the right side of the image, NE of the lake.Stations EP-A, EP-B, EP-C and IP-A are located near the gully margin.
Climate normals  were recorded at the Pond Inlet airport meteorological station at 62 m a.s.l.(Environment Canada, 2014) on Baffin Island, 85 km southeast from the study site.Daily average air temperature was −14.6 • C and normal yearly average precipitation was 189 mm (91 mm of rain).A cross-correlation analysis revealed close correspondence between observations at the Pond Inlet airport and those in the Qalikturvik valley (Fortier and Allard, 2005;Gauthier et al., 2013).Daily average air temperature for the 1994-2013 interval at the study site was −14.5 • C (CEN, 2014).
Four polygons located near a gully (labeled R08) experiencing active erosional activity (Godin et al., 2014) were utilized (Fig. 2).Discharge in the gully reached 1 m 3 s −1 due to the freshet and diminished rapidly to a total cessation during July, except when substantial precipitation occurred (Godin et al., 2014).
Three of those polygons were partially eroded during the development of the gully, and the fourth remained intact.These three polygons were initially eroded in 1999 and the early 2000s (Godin andFortier, 2010, 2012b) and have since stabilized.The polygons' geometry and instruments' specifications are detailed in Table 1.

Air and ground temperature
Polygons were drilled in 2012 using a portable permafrost core-drill system and the boreholes (BH), approximately 1 m deep, were equipped with a string of temperature sensors in polygons EP-A, EP-B, EP-C and IP-A (Fig. 3).Temperatures were recorded between summer 2012 and summer 2014 for these four polygons at depths of 5, 20, 50 and 90 cm (Table 1).A string of thermistors, installed in a borehole located in a low-center polygon 2 km away from the gully and connected to a logger (BYLOTPD, Fig. 2), recorded ground temperature between 2010 and the summer of 2013 at 10, 20, 30, 40 and 80 cm (Allard et al., 2014).BYLOTPD was drilled and initially utilized with thermistors in 2001, and had a sensor at the surface until 2007; this long-term series was useful when missing data occurred in the time series.Thaw depth was obtained for each site by interpolating the temperature in each borehole between sensors.Air temperature was recorded on site between 2010 and 2014 by two stations (BYLOSIL and BYLCAMP) from the SILA network (CEN, 2014).

Ground moisture and cover
An array of five moisture sensors (Table 1) were deployed in each studied polygon center during the summer of 2013 to monitor the near-surface moisture regime (Fig. 3a and b).Raw volumetric water content (VWC, m 3 m −3 ) was calculated using time domain reflectometry (TDR) sensors.Calibration of the moisture data recorded by the TDR probes was necessary (Czarnomski et al., 2005) due to the high organic content of the soils at the study site.The calibration was performed by saturating a utilized (TDR) soil sample of known mass and volume, while drying at air temperature and weighing frequently, thus building a relation between the weight of the water in the sample and the signal in the mois-  Godin and Fortier (2012b) (Fig. 3).Size was obtained by measuring the length and the width of each polygon by calculating the distance from the highest point of surrounding ridges, near their troughs.The Campbell Scientific temperature sensors (model 107-L) and Omega Type-T thermocouples were connected to a CR1000 datalogger.The Onset S-TMB-M006 temperature sensor was linked to a H21-002 Micro Station.Moisture data acquisition systems (Em5b) from Decagon were used with TDR sensors (EC-5).BYLOTPD is equipped with a 3 m long custom-made thermistor cable (Allard et al., 2014).Temperature sensors were installed in boreholes identified as BH in Fig. 3a and b. ture sensor.The equation, derived from an exponential curve obtained from the logged data (R 2 = 0.97, n = 9), was applied to all moisture readings.For the sake of data comparison between studies (Hinzman et al., 1991;Liljedahl et al., 2011), the moisture data will be presented as percent saturation.This relative metric was determined based on the expo-nential curve generated during calibration and applied to all volumetric water content values.

ID
Precipitation during the summer of 2013 was obtained daily using a Hellmann rain gauge, compact version (CEN, 2014).The snow-water equivalent relation was used to account for solid precipitations.
Maximum snow depth within and near the gully was obtained by time-lapse photography (daily, at noon during 1 year) between the summer of 2012 and summer of 2013.The camera used was a Reconyx model PC-800 Hyperfire fixed on a tripod, with graduated poles deployed in the field of view within and near the gully for reference.

Plant characterization
As a component of a larger study (Perreault, 2012;Perreault et al., 2015), each polygon was evaluated for vegetation cover using three randomly positioned 70 cm × 70 cm quadrats in July 2009 or 2010.Vertical photographs of the quadrats were taken at approximately 1.3 m from the ground (see detailed protocols in Chen et al., 2010).Vascular plants, bryophyte and lichen cover were evaluated as the projection on the ground with an abundance scale modified from Daubenmire (1959) on all pictures overlain with a 7 cm grid.Mean plant cover were calculated for each polygon.Details on the methods are available in Perreault (2012).

Active layer depth, degree days and n factors
Air temperature has an important influence on ground temperature.Thawing degree days (DDT) and freezing degree days (DDF) are simple indexes used to calculate heat induction and extraction in soil.A degree day is the difference between the mean daily temperature and 0 • C (Jumikis, 1977).The sum of the daily averages is then computed for a season, with the negative index being DDF (freezing season) and positive index being DDT (thawing season).
Active layer depth at the end of the thawing season for a given site can be estimated using Stefan's equation modified for permafrost conditions (Brown et al., 2000;Jumikis, 1977): where Z represents the active layer depth, E the edaphic factors (or physical properties of the ground and its cover) and DDT air the sum of the average daily air temperature above 0 • C at the site (Brown et al., 2000;Shiklomanov et al., 2010;Woo et al., 2007).DDT air is provided by the nearby SILA stations (CEN, 2014).
The ratio between the DDT soil (sum of thawing degree day at the soil surface) and DDT air provides the term n t , known as the thaw season n factor (Klene et al., 2001;Lunardini, 1978), as follows: DDT soil DDT air .
Similarly to n t , the term n f is used for freezing season n factor: where DDF soil and DDF air stand for freezing degree day at the ground surface and freezing degree day of air, respectively.DDF soil and DDT soil were provided by near-surface sensors buried at 0.05 m in polygons EP-B and IP-A (BH location in Figs. 2 and 3).Near surface thermal gradients were calculated for sites BYLOTPD (10-20 cm), EP-B and IP-A (5-20 cm) during August 2012 for summer and January 2013 for winter.Thermal gradient can be computed as: where i is the thermal gradient ( • C m −1 ), or the temperature change between two points in a medium (T s and T b , which are the temperature at the surface and bottom) and x the depth (thickness) of the layer considered in the measurements (Jumikis, 1977).

Statistics and landscape modeling
Quantitative analyses were performed with R (R Core Team, 2014).The graphics were prepared with the ggplot2 module.
Site microtopography was obtained using a Trimble VX station in survey mode (prism), where 92-210 spatially referenced points were recorded for each polygon and their proximal area (polygon center, ridge, troughs, nearby retrogressive thaw slumps or gully network branches).Surveyed points had an accuracy of ±1cm.In polygon centers, biological soil crust and mosses were either relatively thin (at most 5 cm thick) or sparse enough to deploy the survey pole on the ground.Ridges and gully slopes were sparsely vegetated and, consequently, the survey pole was in direct contact with the ground.Surveyed spatial data were loaded and processed in ESRI's ArcGIS v10.An interpolation (Gaussian process regression) was applied to the points to create a 3-D surface.A GeoEYE (1 pixel = 0.5 m) satellite image was used as background with an infra-red, red and green false color rendering in Fig. 2 (2 September 2010).

Winter dynamics
Near the gully, a continuous snow cover overlaid the polygons (IP-A, EP-A, EP-B and EP-C) between late September until mid-May 2013, when the blanket started to be fragmented.The polygons adjacent to the gully had a snow cover which was, at most, 10 cm deep during the whole winter 2012-2013.In the depressions (gully), accumulations could either be absent or thicker than 1 m deep, depending on the channel sections.Flat surfaces outside the gully acted as a source for snow and the gully channels acted as sinks where snow accumulated.Shallow snow cover, exposing the ground to increased temperature variations during winter, enabled more heat to be extracted from the ground and diminished the insulating role of snow, with a potential for a thinner maximum active layer depth.A slightly greater n f was obtained at site EP-B (n f 2012-2013 = 0.99 and n f 2013-2014 = 0.88) compared to IP-A (n f 2012-2013 = 0.96 and n f 2013-2014 = 0.85), indicating a thinner snow cover and greater heat extraction during winter.This resulted in a higher DDF soil value for EP-B compared to IP-A for both monitored winters (Table 2).

Summer dynamics
Snowmelt enabled the formation of a shallow pond each year in the center of the intact polygon (IP-A, Fig. 3a).The pond presence was reflected by a ground saturation of 100 % in the IP-A center until 20 June 2013, as indicated by sensor no. 3, and remained greater than 90 % until 21 June 2013 for the other sensors (IP-A, Fig. 4).
The water levels in the area occupied by the pond fluctuated until its disappearance, between late June and early July   sponse to evapotranspiration and the lowering of the water table following the propagation of the thaw front in the active layer and the consequent enlargement of the reservoir capacity.Between 17 June and 1 July 2013, the thaw front progression was 7 mm d −1 , followed by a rapid progression of 100 mm in 4 days.A few precipitation events (IP-A, Fig. 4, E1, E2, E3 and E5) restored the nearly saturated state for short periods.Moisture for IP-A (Fig. 4) was quasi-uniform between the five TDR sensors deployed in its center at any given time.While the moisture inside IP-A was quasi-uniform between sensors at any given time, the difference in moisture decreased significantly between the initial saturated state and the drier, late summer state (Figs. 4 and 5).The last readings of the summer indicated saturation varying between 47.8 and 63 % for a thaw depth of 0.5 m.Thus, moisture was consistent across time at the scale of this specific polygon center.Median saturation levels varied between 82.5 and 94.4 % during the measurement period, at levels which are relatively close to the saturation levels (Fig. 5).Overall, moisture levels responded homogeneously to either input (rain) and output (evapotranspiration) at all monitored locations inside this polygon and it stayed moist during summer.

Standard deviation between each individual TDR in IP-A
Near-surface ground moisture conditions evolved differently in eroded polygons EP-A, EP-B (Fig. 3a) and EP-C (Fig. 3b) compared to the intact polygon IP-A (Figs. 4 and  5), when considering whether the moisture balance of individual eroded polygons (intra-polygon) and between polygons (inter-polygons).
Eroded polygons EP-A, EP-B and EP-C (Fig. 4) were characterized simultaneously by (a) a strong variability in moisture saturation between sensors in their respective centers (intra-polygon); and (b) a variability in overall moisture saturation between each eroded polygon (inter-polygons).Among the 15 sensors deployed in the 3 eroded polygons, only 2 sensors (EP-A no. 2 and EP-C no. 5) recorded curves which were similar (but nevertheless drier) to the wet polygon IP-A.Those two emplacements inside the eroded polygons were near saturation (∼ 90 % saturation) at the beginning of the logging interval, had moisture peak recorded when precipitation events E1, E3 and E5 occurred, and experienced a progressive decrease in moisture levels in response to evapotranspiration and thaw front deepening (EP-A no. 2 median = 74.4%, min = 63.7 % max = 94.5 % and EP-C no. 5 median = 76.7 %, min = 62.9 %, max = 91.3%).On the other hand, EP-C no. 1 (located ∼ 1 m from EP-C no. 2 mentioned above) had a median moisture saturation of 33.5 %, underlining significant differences in moisture inside the same polygon (Fig. 5).For each time step, EP-C moisture intra-polygon readings were the most variable: standard deviation for the duration of the record was at least 14.2 and at most 18.1 % saturation.
More stable ground moisture levels were also observed in eroded polygons.These moisture readings were characterized by a weak response to (a) precipitation inputs and (b) evapotranspiration and thaw front deepening, when compared to the intact IP-A or single emplacement in eroded polygons (EP-A no. 2 and EP-C no. 5).Data recorded by all the sensors in EP-B underlined the low variability of the readings quite clearly, especially as displayed by nearly linear percent saturation moisture levels (Fig. 4, EP-B) and by whisker boxes for EP-B no. 2, EP-B no. 5 (Fig. 5).Thaw depth was very shallow (∼ 20 cm) in some parts of the polygon EP-B (Fig. 4, EP-B, red dashed line); mean thaw depth for EP-B on 1 July 2013 was 19 cm (SD = 4).Sensors located 1.5 m on each side of EP-B no. 2 (EP-B no. 1 and EP-B no. 3) provided non-overlapping ranges for moisture, underlining a great variability under short distance in this disturbed polygon.Other sensors reported similar low moisture variability in EP-A and EP-C.For instance, in EP-B no. 2, percent moisture median was 92.8, varying during the summer between 90.7 and 96.3 %; thus there was a wet emplacement in this polygon, where an adjacent sensor (EP-B no. 1, located 1.5 m nearby) recorded drier conditions with a median of 69.9 %.

Ground temperature and active layer thickness
Ground thermal regime monitoring, obtained in an intact low-center polygon (BYLOTPD) between the winter of 2010 and the summer of 2013, provided maximum active layer depth of 56, 48.5, 52 and 40 cm for 2010, 2011, 2012and 2013, respectively (BYLOTPD, Figs. 6, 7).During these 4 years, the 0 • C isoline reached a depth of 10 cm between 22 and 30 June.Temperatures recorded at BYLOTPD surface between 2002 and 2007 indicated an average delay of 8 days for the thaw front to lower from 0 to 10 cm.Between 2010 and 2013, the thaw front usually progressed rapidly down to 29 cm deep before slowing.Initial ground thaw progression rates varied 1.5-8 cm d −1 and then slowed to less than 0.5 cm d −1 until reaching the near-maximum thaw depth.While the maximum attained thaw depth varied from year to year (BYLOTPD, Fig. 6), the maximum depth ±2 cm persisted for 40 ± 1 d.In 2012, the active layer refroze, slowly at first and then very rapidly, when all the latent heat of water was extracted from the ground.Freezeback of the active layer was completed during the first week of October.During the winter of 2012-2013, the temperature at 66 cm depth remained below −16 • C for 121 consecutive days (violet, BY-LOTPD, in Fig. 7).The 0 • C isoline contour was akin to an irregular-shaped curve, with a rapid thaw and freeze back at the beginning and the end of the summer, respectively, and a relatively stable, long-standing maximum thaw depth in between.
The intact polygon IP-A located near the gully shared some characteristics with BYLOTPD.In 2012, IP-A and BYwww.biogeosciences.net/13/1439/2016/Biogeosciences, 13, 1439-1452, 2016 LOTPD had a similar maximum thaw depth (52 cm) following a similar input of thawing degree days, 450 and 474 DDT, respectively (Fig. 6).While there was some interannual variability, dispersion was limited and similar to readings from other sites (particularly 2010-2012, IP-A, EP-A, EP-C).All sites responded similarly for each monitored year when considering the direction of change, implying that every site had a shallower maximum thaw depth in 2013 compared with 2012, resulting from a lower DDT that year.At the beginning of the record (7 July) the thaw depth was at 42 cm and reached its maximum depth 52 days later (29 August) at 52 cm deep (0.19 cm d −1 during this interval).A complete thawing season was recorded in 2013, showing a rapid initial progression of the thaw front during the first 9 days at 1.67 cm d −1 (beginning 10 June) and then a slower rate until the maximum thaw depth was attained at 50 cm.Once the deepest point was reached, it froze back in 33 days.Overall, the thaw season lasted 94 days for IP-A during 2013.Temperature at a depth of 63 cm was below −16 • C for 165 continuous days during the winter of 2012-2013 and 143 days during 2013-2014 (violet, IP-A, in Fig. 7).Short peaks of warm temperatures were observed in IP-A's center, exceeding 8 • C; in comparison, warm peaks were not present on BYLOTPD since there were no sensors near the surface (shallowest at 10 cm).The shape of the 0 • C isoline contour of IP-A (Fig. 7) was parabola-like in 2012 and 2013.
When comparing ground thermal regime of intact polygons against sites adjacent to the gully, maximum active layer depths were within the same range for 2012 (EP-A = 47 cm, EP-C = 44.5 cm against IP-A = 52 cm and BY-LOTPD = 52 cm) (Fig. 6).Active-layer depths were similar to those sites (except BYLOTPD) in 2013, with thaw initiation varying between 8 and 12 June.In 2013, BYLOTPD maximum thaw depth was shallower than the other sites, yet followed a similar trend than most of the other sites during this colder year.EP-A and EP-C had a parabola-shaped evolution of the 0 • C isoline, with rapid thaw initiation and refreezing, similarly to IP-A.There was no temporal persistence of the thaw front for IP-A: once the maximum depth was reached, the active layer started thinning on the day following the maxima (for 2012 and 2013).EP-A and EP-C ground temperature at 63 cm was below −16 • C continuously during 2012-2013 for 176 and 170 days, respectively; during the winter of 2013-2014 for 168 and 163 days, respectively (violet., EP-A and EP-C, in Fig. 7).
On the other hand, EP-B had a very shallow maximum active layer depth (ALD) with 21 and 20.5 cm recorded for 2012 and 2013, less than half of the other monitored polygons, either intact or eroded (Fig. 6).Thaw depth progressed at a rate of approximately 1 cm d −1 for 9 days, followed by a clear slowdown toward stabilization of the thaw progression.Maximum thaw depth persisted for 59 days in 2012 and 51 days in 2013.The 0 • C isoline dropped sharply following the initial thaw, followed by a long, stable maximum depth and rapid refreezing.Surface temperature peaked at 11.4 • C dur- During the summer of 2013, the n factor n t was closer to 1 for IP-A than EP-B (1.02 and 1.06, respectively).Thermal gradient was steeper during summer for EP-B than for IP-A (Table 4).
The n factor n f computed for the sites indicated values closer to 1 during winter for EP-B compared to IP-A.A n f of 0.99 for EP-B during 2012-2013 (Table 2) clearly suggested a thin snow cover at this location.Thermal gradient for EP-B was very steep at shallow depths with −26 ± 9 • C m −1 (Table 4) during winter, reflecting the absence of a substantial snow cover and considerably larger than other intact sites, either at Bylot Island or at other undisturbed sites.Thermal gradient in the literature during winter in the active layer or the near surface was similar for all sites except EP-B (Table 4).
Therefore, the proximity of the gully and the consequent shallow snow cover in polygons near the eroded channels could impact heat extraction during winter, compared with an intact polygon (BYLOTPD) where the microtopography of polygon ridges and absence of the gully enabled thicker cover.Near-surface averaged maximum temperatures were generally cooler in the intact polygon (BYLOTPD, 3 • C) in 2012, compared with the other polygons, as shown by the reddish colors delineating the 2 • C isolines in Fig. 7 The year 2012 had a warmer summer than 2013 and the winter of 2012-2013 was warmer than the winter of 2011-2012 (Table 3).The polygons located within a 1 km radius and those near the gully were exposed to similar DDT air .Inter-polygonal differences in the active layer dynamics (e.g., moment of maximum ALD, maximum depth of the active layer, averaged values for near-surface temperature) were due to polygon-specific surface characteristics (absence of snow cover during winter and vegetation, moisture during summer) impacting ground thermal dynamics at each respective site.

Vegetation in intact and eroded polygons
In 2010 and 2014, the center of the intact wet polygon (IP-A) was uniformly vegetated with typical wetland vegetation with low vascular plant diversity (see wetland vegetation in Perreault et al., 2015 in this issue).Perreault (2012) measured a strong cover of living mosses Drepanocladus sp.53.3 % and Polytrichum sp.1.3 %) and Carex aquatilis (27.5 %), a sparse cover of Dupontia fisheri (0.5 %) and traces of Pedicularis sudetica, Arctagrostis latifolia and Salix arctica (Supplement Table S1).
The three other polygons had a higher vascular plant diversity and less uniform vegetation with a mixture of wetlands and mesic species typical of disturbed polygons (Perreault et al., 2015, this issue).Polygon EP-A had 20 vascular plant species with the higher covers for wet habitat species with dried Drepanocladus sp.mosses (77.5 %), Carex aquatilis (6 %), Dupontia fisheri (3 %), Eriophorum angustifolium (2.5 %) and Eriophorum scheuchzeri (1.3 %), as well as traces for a number of typical mesic habitat species such as Arctagrostis latifolia, Cerastium alpinum, Luzula confusa, L. arctica and Stellaria longipes).Polygon EP-B had a fraction of its surface as bare ground (2.5 %), a sign of disturbance; 42.7 % of dried mosses of wet habitat Drepanocladus sp.; and 16 vascular plant species.Wet habitat species and mesic species shared the dominance with 21 % Eriophorum angustifolium and 10.8 % of arctic willow (Salix arctica).Lastly, polygon EP-C with 11 vascular plant species was dominated by mesic habitat mosses Aulacomnium sp. ( 41 typical wet species were not present in this polygon anymore (Perreault et al., 2015).

Intact polygons
The comparison of shallow ground thermal regime between IP-A and BYLOTPD showed that similar terrain located only a few kilometers apart in a similar terrain can behave slightly differently in response to local conditions.For those two intact sites, the active layer depth for most years (Fig. 6) was within similar ranges, but the shape of the 0 • C isoline and the timing of ground thaw initiation differed slightly (Fig. 7), implying site-distinct edaphic factors.Thermal gradient was greater for IP-A than for BYLOTPD (Table 4).This is interpreted as the result of differential lateral subsurface water flow between these two sites.Indeed, the presence of a gully near IP-A favored drainage which reduced subsurface flow, whereas BYLOTPD had no such gully effect during the thaw season.Also, snow cover dynamics may differ in a similar manner between two such polygons, with an increased possibility of snow being blown from a polygon into the gully.This is reflected by colder temperatures at 60 cm depth during winter for BYLOTPD compared to the other sites (violet, in BYLOTPD Fig. 7).Polygon IP-A was self contained (intact rims) which enabled intra-polygonal moisture homogeneity in its center along the monitoring transect (Fig. 5).
Carex aquatilis, the dominant species in this polygon (Supplement Table S1), was an efficient competitor in such stable and wet environments (Billings and Peterson, 1980).

Eroded polygons
Neighboring (IP-A, EP-A and EP-B) and nearby (EP-C) polygons had less environmental variability among each other because of their proximity.The eroded polygons were each located close to a gully, their topographic gradient was alike and their proximity to a topographic low (gully bottom) was similar.What varied was the length of their contour adjacent to the gully, the integrity of their ridge and their intrinsic potential to retain moisture in their center, including the capacity to retain snow during winter.
Prior to 1999 and prior to gully inception, these eroded polygons were low-center polygons (Godin and Fortier, 2012b), with similar surface conditions.When erosion was initiated in 1999, not every site was affected in the same way: EP-A had a secondary ice wedge that melted in the center of EP-A, creating a linear trench and reducing ground moisture at all monitoring sites.However, it was observed that some parts of the center were able to maintain wet conditions even up to 14 years after degradation (EP-A nos. 1, 2 and 3) (Figs. 4 and 5).This result stresses the heterogeneous character of moisture within eroded polygon centers.This situation translated into more microhabitats and higher plant diversity (20 vascular species, see Supplement Table S1) and the presence of plants typical of stable wet environments (Carex aquatilis) and a species efficient in colonizing wet environments (Dupontia fisheri) (Billings and Peterson, 1980) or when water supply was uncertain (Bliss and Gold, 1994).
EP-C's moisture signature in this well-drained polygon shared characteristics with EP-A.Differential thaw subsidence in the center changed the polygon microtopography and gave rise to some wetter and drier areas, although to a lesser extent than EP-A (Fig. 5).One sensor in particular (EP-C no. 1) had a low moisture saturation through the recording interval; such an intra-polygon heterogeneity and consequent moisture conditions favored the growth of plants that were either rare or absent from IP-A (Supplement Table S1: Salix arctica, S. reticulate and Arctagrostis latifolia).
EP-B's thermal regime and the moisture signature were strikingly different from observations in the other polygons.The active layer was thin and with a persistent maximum thaw depth at the borehole emplacement during summer.A thin active layer resulted in extreme thermal gradients, very high during winter and summer (Table 4).Two sides exposed to the gully contributed to a shallow snow cover during winter, enabling heat to be further extracted from its center.Thinning of the active layer resulted in the aggradation of ground ice as the permafrost table rose up.Therefore, thaw depth progression was closely connected with the sum of the DDT air for intact polygons, which is uncertain with an eroded site.The thinning of the active layer further implied local ground ice, carbon and nutrient fixation due to freezing following upward permafrost aggradation, but on a local scale.Moisture conditions varied considerably over very short distances.Thinning of the active layer also contributed to reduce the potential volume to retain water and contributed to keep the water table close to the surface.Intrapolygonal heterogeneity in moisture conditions resulted in a plant distribution including wet and mesic species along with dead mosses and vascular plants, remnants of the previous wetter state of the polygon.These results indicate that degraded low-center polygons switch from relatively homogeneous wet conditions to heterogeneous conditions in moisture, thermal conditions and plant distribution.
When impacted with erosion, polygon centers' physical characteristics do not necessarily tend to all change toward the same trend.Three nearby eroded polygons in the same gully network, eroded at approximately the same time, shared a (more or less severe) diminution of the ground moisture, and not in a uniform way inside each respective polygon.Thermal dynamics during summer were more complex than those during winter; individual sites' albedo, nature of the cover, moisture and ground temperature need to be taken into consideration to precisely identify differences between sites.

The question of scale: time and space
The polygon terraces characterizing the floor of the valley of glacier C-79 is the product of thousands of years of eolian sedimentation, organic and ground-ice sedimentation along with syngenetic growth and cracking of ice wedges under cold conditions.The development of ridges created lowcenter polygons and, under certain terrain conditions, groups of polygons evolved into high-center polygons as their ridges were gradually eroded.
In low-center polygon terraces, when the proper conditions are met, large gullies may form very rapidly (Fortier et al., 2007) with yearly rates of ice-wedge erosion greater than several tens of meters per year.The erosional processes will be very dynamic and occur rapidly, immediately following the triggering events as the components of the landscape in disequilibrium are in transition and unstable (Sidorchuk, 1999).Gully channels ranged in dimensions as small as 1-2 m wide and deep and smaller trenches were observed to further connect polygons to the gully (Fig. 3), to large channels measuring more than 10 m wide and 4 m deep, networked in kilometer-long channels (Godin et al., 2014), thus fragmenting the terrace.Once initiated, not only it could take years (and decades) for erosion to stabilize (Godin and Fortier, 2012a) but the gully may capture and facilitate the water fluxes downstream, enabling enhanced drainage and soil displacement outside the watershed (Fortier et al., 2015;Godin et al., 2014).At medium spatial scale, more than a thousand polygons were directly eroded following gullying in this valley, with various levels of severity.As a result, some polygons were completely eroded (baydzherakhi), and others severely reduced in area as two sides or more were eroded.At the fine spatial scale, perturbed polygons such as EP-A and EP-B, which have a ridge or two eroded and most of their internal area intact, experienced physical surface changes as well.
As the gully evolved over the years, the drainage network found alternative paths downstream.Water entry points were abandoned as the gully evolved upstream and retrogressive thaw slumps were stabilizing near an axis end when a stream was not accelerating the erosion process (Godin and Fortier, 2012b).The capacity of eroded polygon to retain water and moisture remained diminished compared to intact polygons through time, even though they were in the course of stabilization toward their new equilibrium state.Polygons in the current study were in transition toward the new equilibrium -changes in surface conditions caused nonlinear fluctuations in the active layer depth, moisture content and plant species distribution and occurrence.In any case, the eroded polygons definitely will not return to their pre-erosion state.The gully network across the polygons constitutes a fragmentation of an otherwise flat terrain unit with an improved drainage network and drier polygons in a wider area than only the one directly linked to the gully, potentially diminishing the flux of the lateral flow.Instead of a flat, continuous terrace, gullies fragment the terrace and affect stream circulation and water recharge potential.The terrace should be a drier overall place than before gullies formed, with a more heterogeneous (and simultaneously wet) environment, as the new steady state.
Medium-scale studies monitoring gully channels and other erosion morphologies do not provide quantitative information on how polygons adjacent to the gully could be affected.Terrain heterogeneity following permafrost erosion and how the hydrologic system evolves thereafter is essential when modeling methane and carbon emissions: a model by Cresto Aleina et al. (2013) took into account the local heterogeneity of polygonal terrain and satisfactorily compared this against GHG emissions measured on the field.Eroded polygons by gullying, as presented in this paper, underlined how a study at fine scale could result in identifying different impacts of adjacent sites eroded in a similar manner.Finescale changes, such as variable moisture in a single polygon, changing ground thermal regime and consequent heterogeneity of plant cover in an eroded polygon, are major impacts as thousands of polygons are subject to erosion.When the gully and the eroded polygons stabilize and a new equilibrium state is attained, heterogeneous conditions characterize those stabilized surfaces compared to the more uniform intact polygons.

Limitations
Once triggered, as already demonstrated in previous papers (Godin and Fortier, 2012b;Fortier et al., 2007), a gully can change quite fast by very rapid impact from nearby polygons.Yet, gully inception (triggering) observation in the field was exceptional: once in 1999 -others were observed but were not further studied.Historical air imagery was quite useful to assess the period of initiation and evolution of several large gullies (Godin and Fortier, 2012a).However, precise or finerscale observations of active layer depth, ground temperature and ground moisture were only obtained during the last few years.

Conclusions
In the valley of the Glacier C-79 on Bylot Island, the rapid thermo-erosion of permafrost leading to gullies enabled the disturbance of ice-wedge polygons and consequently caused mass transfers such as water, sediment, peat transport, permafrost and massive ice thaw.The gully changed fluxes of heat and matter in the geosystem; once the gully stabilized it began to exist as a new steady state, permanently distinct from the initial conditions.Polygons adjacent to the gully can be partially impacted at various degrees of severity.Gully effects on polygons can be partial erosion and subsidence of the ground, inducing differential microtopography, ground moisture retention capacity (rain and snow), thermal conditions (via water and snow) and thus, active layer thickness.Such E. Godin et al.: Nonlinear thermal and moisture response of ice-wedge polygons changes in a polygon have impacts on ground ice conditions, and the carbon and nutrient cycle.These differential conditions create heterogeneity at a fine scale in eroded polygons, where colder, warmer, wet and mesic conditions can coexist in the proximity of the same polygon center and adjacent polygons.These differential conditions affect plants' distribution, diversity and abundance.Vegetation will have an effect on physical and biogeochemical properties of the active layer dynamics in the long term, but more research and monitoring is needed to understand these trajectories.
The Supplement related to this article is available online at doi:10.5194/bg-13-1439-2016-supplement.

Figure 1 .
Figure 1.The study site is on Bylot Island in the Canadian Arctic archipelago (73 • 09 N, 79 • 57 W), north of Baffin Island (a).The valley of Qalikturvik is located in the southwestern section of the island (b) (background: NRCan Landsat-7 orthorectified mosaic, 3 August 2010).

Figure 2 .
Figure 2. The location of the stations at study site; map background is a GEO-Eye false color satellite image (Near Infra-Red, Red, Green) obtained 2 September 2010.BYLOSIL and BYLCAMP are meteorological stations that are part of the SILA network(CEN, 2014).BY-LOTPD is a reference ground temperature monitoring site installed in a low-center polygon(Allard et al., 2014).The gully R08 is located on the right side of the image, NE of the lake.Stations EP-A, EP-B, EP-C and IP-A are located near the gully margin.

Figure 3 .
Figure 3. Digital elevation model of sites EP-A, EP-B, EP-C and IP-A.Altitudinal isoline contours were digitized on the figure at each 0.1 m.A gully channel (deep blue) was initiated in 1999 near the polygons; rims (yellow and orange) delineate polygon contours; light blue indicates a breach by erosion in EP-A, EP-B and EP-C.The IP-A ridge contouring the polygon is intact.The gully channel floor (dark blue) is approximately 2 m lower than nearby polygons' rims at this position.Boreholes equipped with a string of thermistors (BH in the figure) are identified in red in each polygon and moisture sensors are identified in blue (TDR in the figure).

Figure 4 .
Figure 4. Moisture readings for the near surface of sites EP-A, EP-B, EP-C and IP-A during the summer of 2013 (Fig.3).Daily precipitation readings indicated six rain events through the summer, identified as E1-E6.Propagation of the thaw front in the active layer evolution is identified as the 0 • C isoline for each site.Moisture levels in percent saturation decreased following active-layer depth for all sites except EP-B.

Figure 5 .
Figure 5. Variability of moisture conditions during the summer of 2013 (percent saturation) per polygon, for each sensor.The line located in the center of each box indicates the median, lower box line indicates the first quartile and the higher line indicates the third quartile; points are outliers.Underlying data refer to the variability of Fig. 4.

Figure 6 .
Figure 6.Maximum thaw depth (y axis, cm) related to the square sum of the degree days of thawing (x axis, √ DDT air ) for sites between 2010 and 2013.The relation between those two variables for the sites was similar for all locations except EP-B.Sites measured in 2012 were exposed to more DDT than in 2013, but this tendency was not reflected clearly in the maximum active-layer depth, possibly implying varying edaphic factors between both years.EP-B's depth was virtually the same between 2012 and 2013.

Figure 7 .
Figure 7. Air Temp.(top) represents the mean daily air temperatures from BYLCAMP (yellow) and BYLOSIL (red) recorded between July 2012 to July 2014.The 0 • C limit is highlighted in red as a reference for mean daily air temperature.BYLOTPD, EP-A, EP-B, EP-C and IP-A's top 80 cm ground thermal regime is as follows: purple color indicating colder temperature and reddish color indicating warmer temperature.The 0 • C isotherm of the ground temperature is indicated as a red line.

Table 2 .
Degree days of freezing (DDF air ) as recorded by the on-site meteorological stations, between 2010 and 2014(CEN, 2014), and DDF soil for EP-B and IP-A (the winter of2012-2013 and 2013-2014).The length of the thawing season is indicated for each measured unit under its respective n day column.Wintern day DDF air n day DDF soil IP-A n day DDF soil EP-B

Table 3 .
Degree days of thawing (DDT air ) as recorded by the on-site meteorological stations, between 2010 and 2014(CEN, 2014), and DDT soil for EP-B and IP-A (summer 2013).The length of the thawing season is indicated for each measured unit under their respective n-day column.Summer n day DDT air n day DDT soil IP-A n day DDT soil EP-B

Table 4 .
Thermal gradients reported from the literature (as reference) and from sites in this study (EP-B, IP-A).The gradient i • C m −1 , the depth, the season, the location and the source are mentioned.Winter thermal gradients were common while spring/summer gradients were sparser.