Effects of thermo-erosion gullying on hydrologic flow networks, discharge and soil loss

Thermo-erosion gullies in continuous permafrost regions where ice-wedge polygons are widespread contribute and change the drainage of periglacial landscapes. Gullying processes are causing long-term impacts to the Arctic landscape such as drainage network restructuring, permafrost erosion, sediment transport. Between 2009 and 2013, 35 gullies were mapped in a polygon terrace in the valley of the Glacier C-79 on Bylot Island, Nunavut (Canada), one of which was monitored for its hydrology. A gully (R08p) initiated in 1999 in a low-center polygon terrace. Between 1999 and 2013, 202 polygons over a surface of 28 891 m2 were breached by gullying. Overall, 1401 polygons were similarly breached on the terrace in the valley before 2013. R08p is fed by a 1.74 km2 watershed and the hydrological regime is characterized by peak flows of 0.69 m3 s−1 and a cumulative volume of 229 662 m3 for 2013. Historic aerial photography from 1972 and recent field surveys showed a change in the paths of water tracks and an increase in channelized flow in the gully area from none to 35% of the overall flow path of the section. The overall eroded area for the studied gullies in the valley up to 2013 was estimated at 158 000 m2 and a potential volume close to 200 000 m3. Gullying processes increased drainage of wetlands and the hydrological connectivity in the valley, while lowering residence time of water near gullied areas.


Introduction
Thermo-erosion gullying is a type of permafrost degradation process generally observed in the continuous permafrost landscape. Specifically, thermo-erosion is a thermo-mechanical process in which ground ice is eroded and residual sediments and carbon displaced by hydraulic transport (van Everdingen 1998). Thermo-erosion is a permafrost degradation process contributing significantly to permafrost landscape evolution along littorals (Hoque and Pollard 2009), river (Costard et al 2003) and inland environments (Fortier et al 2007). Thermo-erosion gullies are an integral part of continuous permafrost landscape evolution and as such were observed in Alaska, Yukon, North Western Territories, the Canadian Arctic Archipelago (Bylot Island) (Grosse et al 2011), the Antarctic Dry Valleys (Levy et al 2008) and Siberia (Czudek andDemek 1970, Sidorchuk 1999).

Environmental Research Letters
Environ. Res. Lett. 9 (2014) 105010 (10pp) doi:10.1088/1748-9326/9/10/105010 Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Sudden, repeated or massive input of run-off water over inland polygons can trigger the formation of sinkholes and tunneling of ice wedges, leading to the formation of gully networks (Fortier et al 2007). Following inception, the rate of gully erosion can reach several 10's m per year and subsequent positive feedbacks can keep the gully active and enlarging for more than a decade (Godin and Fortier 2012b). Climatic models (Easterling et al 2000) and recent climate observations indicate that the likeliness of occurrence of extreme events (such as air temperature variability, amount of precipitations) is increasing, which favors thermokarst development (Hinzman et al 2005) and thermo-erosion gullying.
Sheet flow in lowlands polygons, particularly during freshet is the main contribution to the hydrologic recharge through the year. Permafrost gullying, as an incision in the terrace, changes how surface water flows overland and potentially diminishes the extent of the recharge during freshet and afterward because of an enhanced drainage (Woo and Young 2006). Once gullying occurs, not only does the hydrologic connectivity increase but also the capacity of flow network transport to sediments (Poesen et al 2003) and nutrients (Harms et al 2013) in the watershed.
In the valley of glacier C-79 on Bylot Island, 35 thermoerosion gullies were identified and mapped since 1999 (figure 1). The gullies covered in this paper consist of a contiguous network of thawed polygon troughs and edges resulting from the partial erosion of numerous low-center polygons. One of these gullies was closely monitored since its inception in 1999 (Fortier et al 2007). Water tracks existing on the terrace before the formation of the gully were visible at the current gully location on aerial photography from 1972. Water tracks and stream flow were captured by headward gully thermo-erosion, increasing connectivity (Pringle 2003) and further enabling transport of sediments and nutrients (Bowden et al 2008).
The gullying of the terrace in the valley has fragmented the terrace, changed the local hydrology and enabled permafrost thawing. In this paper we assessed the impacts of the permafrost gullying on the integrity of the terrace by first presenting the gully geometry and their rate of erosion in the valley based on field measurements and historic aerial imagery. Second, changes in the hydrologic network associated with a single gully were identified by comparing historical aerial photography with recent observations on the field and excluding changes in the distribution of two flow types (channelized versus overland and sheet flow). Further, the hydrograph of the gully enabled the quantification of the water transiting through the outlet and its flashiness following a rain event could be assessed. Finally, the tri-dimensional geometry of the gully was surveyed to determine the volume of permafrost eroded.

Site
The study site (73°09′ N-79°57′ W) is located in the valley of glacier C-79 on Bylot Island in Nunavut, Canada. The valley's dimension is ∼17 km long by ∼4 km wide from the sea to the current position of the glacier (figure 1). The valley is bordered on each side by shale and mudrock plateaus about 500 m a.s.l. A braided river flowing on the floor of the valley drains meltwaters of C-79 and C-93 glaciers (Inland Water Branch 1969) towards the Navy Board Inlet. A few ephemeral streams originate from sub-perpendicular valleys and alluvial fans and are connected to the braided river. The river is bordered on each side by low gradient terraces with lowcenter polygons associated with ice-wedge formation. The terraces consist of an accumulation of over 4 m of peat mixed with aeolian silts and sands, and in some sectors alluvial/ colluvial deposits, overlying glacio-fluvial sand and gravels and glacial drift accumulated on the valley floor following the retreat of the glacier during Holocene (Fortier and Allard 2004). The depth of the active layer varies between 40-60 cm in peaty-silts and up to 1 m in sands and gravels. The regional depth of the permafrost body was estimated to 400 m (Smith and Burgess 2000).
Climate normal  were recorded at Pond Inlet Airport located 85 km southeast from the study site (Environment Canada 2010). Mean daily air temperatures between Pond Inlet Airport and a meteorological station located on the study site were quasi-identical for the 1994-2002 period (Fortier and Allard 2005). Mean air temperatures for 1981-2010 at Pond Inlet were −14.6°C (61 m a.s.l.) while it was −14.5°C at 24 m a.s.l. in the valley of the glacier C-79 (CEN 2013) for the 1994-2013 period.
Normal yearly average precipitation at Pond Inlet was 189 mm, with an observed 91 mm of rain, and the remaining as snow water equivalent (Environment Canada 2010). Hydrological regime of the stream flowing in the gully is nival, which is characterized by a seasonal high flow following snowmelt and diminishing drastically after, which is common for small streams in the High Arctic (Woo 2012). Streams observed on the valley slopes flow towards the terraces and were either positively contributing to the water balance of the wetlands or drained by the gullies.
The orientation of the main channel of the gullies was generally normal to the river. The recent (2007-2011) erosion rates of the largest and most active gullies in the valley varied between 17 m y −1 and 23 m y −1 for the main axis length and 306 m 2 y −1 and 872 m 2 y −1 in eroded area. (Godin and Fortier 2012a).

Erosion rates
Historical aerial photography, satellite images and field surveys over a decade were used to build a geospatial database of the gully geometric and geomorphological layout, precisely track spatial changes in the landscape and approximate the date of initiation if available (Godin and Fortier 2012a). Thirty five gullies identified on the photo were visited at least once between 2009 and 2013 to record if thermo-erosion recently occurred at each emplacements. Gully margins were mapped using a Differential GPS (or DGPS) model Trimble Pathfinder Pro XRS with a TSC1 data collector (2009)(2010)(2011)(2012) and a Trimble R8 GNSS system (2013). Differential correction was applied to the recorded spatial data obtained with the Pro XRS system by using the Thule (Greenland, DK) reference station located 496 km to the east of the Glacier C-79 Valley. The report of the correction indicated a spatial precision of 0.5 m-1 m for the x,y,z, axis for 99% of the records. Centimeter-level precision was obtained while mapping with the real-time R8 GNSS System. Gullies were manually digitized in a geographic information system (ArcGIS v10) for 1972, 2007, and 2010 based on air photography and satellites images: aerial photography from 1972 (A23052- [61,62,77,209,211,213], A23098-181, 1:15 000) from the National Air Photo Library (NAPL) were used to estimate years of gully development. An IKONOS pan-chromatic image was used as a reference for 2007 (late July 2007, pixel = 1 m) and a true color satellite (GeoEYE) image was used as a reference for 2010 (2 September 2010, pixel = 0.5 m). Gully margins were superposed in the GIS and differences in geometry related to time were obtained. The length of the main axis and the area of each gully were calculated in the GIS. It was therefore possible to determine the gullies as being present or not during the defined intervals.
Superposing the digitized gully boundaries for each lapse on the available images provided the count of breached polygons per interval (1972, 2007, 2010 and 2013). Polygon ridges and troughs contrasted well on the images and those contiguous to a gulling axis were manually interpreted. There was no differentiation between various generations of polygons (when applicable) and the smallest unit identified counted as one.

Hydrology
The hydrological impacts of gully development focused on changes of water tracks and drainage patterns following the formation of a gully during the late 90's. Water tracks and stream paths near R08 were already visible on the aerial photography from 1972 (A23052-61, A23052-77, 1:15000) in the area where a gully later initiated in 1999 were digitized in the GIS software. Recent surveys of the stream network up to 2013 were compared against photography pre-dating the gully using a GIS. The contributing area of the watershed was defined using as reference a topographical digital elevation model [038C04, 048D01, 1:50000] (NRCan 2004) and the GEOEye 2010 satellite image. This contributing area was divided by terrain unit: sub-unit 1 included exclusively the terrace, while sub-unit 2 is the colluvial hillslope connecting with the plateau and the southern valley wall (figure 1). Changes in the network were calculated using the drainage density relation (Dingman 2002) for 1972 and 2013: where D d is the drainage density of the contributing area, ΣL the sum of the length of the streams in the area and A d the area of the contributing area. Lengths were separately calculated for incised channels (Bledsoe et al 2002), water tracks or streams flowing on the terrace and streams originating from a connecting secondary valley. Precipitation readings were manually obtained daily on site from a Hellmann Rain Gauge, compact version (CEN 2013).
Streamflow was measured at the gully outlet morning and evening daily between June and August in 2012 and 2013. A velocity-area method (Dingman 2002) was used, where Q is the discharge, V i the velocity and A i the area of the measured section. The cross-section micro-topography at outlet was measured at an interval of 20 cm. Therefore, a precise measurement of the wet area of the cross-section A i could be obtained for any level of water in the channel by using: where W i represent an interval of 20 cm of the cross-section width and D i the average depth of the channel at W i . Average velocity V i was obtained at 0.5D at center of the stream using a Global Water FP-101, and a tape was used to measure the cross section width W i .
A cumulative discharge graph related to time with a separation of rain events display how the streamflow respond to those events. Bi-daily discharge Q records were used to model the volume of water flowing thought the outlet per 12 h by transforming the sampled discharge per second to a related volume per 12 h Q 12 . The estimation of the rain volume contribution to the watershed is the product of the watershed area with the amount of precipitation.
Gully R08p was monitored with precision for linear erosion rates of gully branches, and total eroded volumes of permafrost. Gully margins were mapped with a GNSS in 2013. Gully topographical cross-sections were measured with a Trimble VX Station at 29 emplacements distributed across the whole gully, including shallow retrogressive thaw slumps, stabilized sections and recent, deeper sections. 20 to 40 tridimensional measurements were obtained for each cross sections. Each cross-section profile was simplified from 3d to 2d by summing vectors from point to point. Each profile line was fitted (using a 3rd to 10th degree equation depending on the complexity of the profile). Depth values were obtained using the resulting equation at intervals of 20 cm, and thereafter averaged. The average depth of a given cross section was multiplied with a section's area to obtain an eroded volume (figure 2). The volumes of all the sections of the gully were summed to obtain the total eroded volume. Only the channels were taken into account for volume calculation, which excluded non-collapsed polygon centers for instance.

Erosion rates
Thirty five gullies developed, enlarged or stabilized since at least the 1972's in the deposits forming the terrace found on the valley floor. Of the 35 gullies located in the valley, 19 were actively eroding and 16 were inactive during the 2007-2013 period (figure 1). Of the 19 active gullies, 14 were showing very limited erosion rates (cm y −1 ) and 5 larger erosion rates (m y −1 ) (figure 3). Most gullies were directly connected to the river, either directly or indirectly via a lake.
Yearly rates of erosion for the main axis of gullies (excluding R08p) varied between less than one and up to 27.67 m y −1 in length and from two to 1080 m 2 in area The rate of erosion in main axis length of observed gullies correspond well with the area eroded (R 2 = 0.87 on a power regression line) only for gullies without a secondary axis. Gullies with a secondary axis were excluded from the relation, because secondary channeling does not count in increasing the main axis length.
R08p is ranked second in eroded area (28 991 m 2 ) over the 35 observed gullies totalling 158 500 m 2 (2013), the first being R05 with 29 485 m 2 . R08p is ranked first in average rate of eroded area at 2063 m 2 y −1 over the 14 years since its formation in 1999. 18 gullies were active between 1972 and 2007, five between 2007 and 2010 and four between 2010 and 2013. Three were active through the whole interval .
As gullies evolved through time, polygons ridges, troughs and even polygon centers were eroded or breached. Sizes, shapes and rates of erosion of the gullies were highly variable. Most gullies (63%) in this study were  existing prior to the images from 1972. Coarse images (1:60 000) existed for 1958 and already 16 gullies were identifiable on the air photography but the potential error of interpretation at such a scale (±15-30 m) was already larger than the typical gully channel. To summarize: Table 1 show a steady increase in the amount of gullies in the valley until 2007. Every new gullies formed between 1958 and 1972 were quasi linear, where it was exceptional that more than one or two side of a polygon breached. Except for two gullies that were longer than 1 km, all others were a few hundred meter in length. Between 1972 and 2007 not only new gullies continued to initiate at the same rate, but among those new gullies a few prove to be larger and wider than those already formed (figure 4, yellow bars-newly eroded polygons as a reference). There were no new gullies between 2007 and 2013 (figure 4, orange and red bars), and four were significantly active during that interval: of those four, one formed before 1958, two between 1972 and 2007 and R08p specifically in 1999. Gullies have a limited period of development. 11 gullies were active only between 1972 and 2007. Only one gully was active from 1972 to 2013.

Hydrology
Changes 1972-2010-13. The formation of gully R08p in 1999 (Fortier et al 2007) was the result of a process where streams were already flowing over the surface for a few decades (at least 1972). There was two flow sources feeding the gully, one located 100 m to the south flowing through a lateral moraine from a nearby secondary valley. The second one was the result of water track concentration from the upper slope to the terrace, south-east of the gully head. Both are ephemeral and were often observed to not be flowing during late July, depending on the antecedent moisture conditions. A mesh of water tracks were identifiable on figure 5(a) (1972).
A 1.74 km 2 watershed encompassed the gully R08p (figure 1). About 81% (or 1.41 km 2 ) of the watershed occupied by a sequence of plateau, slope and lateral moraines generally characterizing the southern wall of the valley, south of the gully. Few changes were observed between the 1972 and 2010 observations in this contributing area of the watershed. Since the focus of this research is the gullying of polygons on the terrace, the watershed was divided in two contributing area: sub-unit 1 (0.33 km 2 ) on the terrace and sub-unit 2 on the slope and plateau (figures 5(a), (b) in yellow).
Following formation of the gully in the sub-unit 1, an important fraction of the water track length changed in this section. Of the 4.42 km of water track length observed in 1972, only 2.62 km remained in 2013 as the length of the channelized flow directly resulting from the gully formation and further evolution raised to 1.53 km in length (or 35% of the total flow path length in 2013). There was no channelized flow in this section in 1972.
All watershed metrics, overall and per sub-unit, for 1972 and 2013 (flow path length, drainage density and separation of flow type such as channelized, stream and water track) are available in the supplemental data, available at stacks.iop.org/ ERL/9/105010/mmedia. Entry points and discharge. During mid-June 2010 and 2012, 29 water entry points (sinkholes and gully heads, figure 5(b)) position and discharge were obtained for the gully. Each blue dot represents a water track or stream entry point in the gully. Discharge at any of those points varied between less than 0.01 m 3 s −1 up to 0.1 m 3 s −1 . Four gully head and sinkholes had discharge of 0.1 m 3 s −1 .
Snow accumulating within the 1.58 km 2 watershed and in the gully channel prevents free flowing of the freshet during late spring when snow melts. Pools of water that accumulated behind snow banks eventually breached, freeing restrained water until the next snow obstruction downstream. This cascade of restrained water eventually reached the gully head area and the process continued inside the gully network. The 835 m long multi-channelled gully had about half a dozen snow banks, which broke one after the other as water flowed inside the network towards the outlet.
Discharge and precipitation were recorded from June to August 2013 (figure 6) at gully outlet (red dot, figure 5(b)). Stream flow in the gully following the final snow bank collapse near the red dot (figure 5(b) initiated on 10 June 2013. The yearly peak flow followed the dam collapse by a few hours and lasted until16 June 2013 (E1) with a cumulative flow of 134 600 m 3 in 6 days (59% of the total flow for the year, figure 6). On site, 45 mm of rain were  To assess the efficiency of the gully as a drainage vector at watershed scale, input volume of rain were related to the overall volume output measured at gully outlet, per event (figure 7).
The input was calculated by multiplying the watershed area (1.74 km 2 ) with the amount of rainfall (in mm) which resulted in a total volume for each events. Output volume was strongly influenced by antecedent conditions and amount of precipitation. Apart from E1 and E4, the relation between input and output corresponded well and is thoroughly interpreted in the discussion.

Soil losses
The gully R08p ( figure 1, figure 5) was measured during summer 2013 for its eroded area and volume. The eroded volume of the gully R08p was estimated at 32 093 m 3 for 22 574 m 2 in channeled area. When taking into account the whole gully network, including partially eroded polygons, baydzherakhii and local ground subsidence caused by gullying the area in 2013 was 28 891 m 2 . The mean depth for the gully channel was 1.42 m. Nevertheless the typical crosssections geometry of R08p varied tremendously along the  Carex aquatilis was spreading in the previously wet polygon and following its breaching by gullying, the polygon was drained, therefore inducing a transition in plants locally. A 1 m wide channelized stream flowed at the floor level of the gully ( figure 8, cross-section 1, c.). The other side of the gully (figure 8, cross-section 1, e.) indicate a transition toward stabilization, where pioneer species such as Tephroseris palustris were colonizing the bare ground, while plants typical of the wetland environment were present over the edge of the cross-section as a remnant of pre-gullying (figure 8, cross-section 1, f.).
Cross-section 2 (figure 8) initiated erosion by gullying in 1999. In 2013, the section was mostly stabilized. Dimensions were 17.2 m wide and 2.19 m at maximum depth with an average depth of 1.21 m. On the left side of the cross-section 2 (a. b.) Salix sp. and Poaceae were colonizing directly on the slope inside the gully. Alluvial bars were common on the floor of the gully near this section.
Cross-section 1 was much more incised than crosssection 2: recent gully channels were not stabilized and therefore may not have slopes in equilibrium. Another factor explaining the difference in incision for the two cross-sections is the base level of the large glacial river located in the center of the valley, which is estimated at 7-8 m above sea level (NRCan 2004). More than 400 m separate the cross-section from one another; there is also a difference of 5 m in altitude of the terrace between the two emplacements, which could further explain the differences in incision.
The total soil (active layer plus permafrost) eroded volume by gullying in the valley, taking into consideration the total incised area and the average depth of R08p (1.42 m) was estimated to be 200 000 m 3 .

Increased hydrological connectivity
Apart from polygon troughs and ridges forming in the 1% slope of the terrace and occasional pond located at troughs intersect, this terrain unit was relatively uniform. Between 1972 and 2013, the contributing area 1) (0.329 km 2 ) (figure 1) saw drastic changes. In 1972 there was 4.4 km of water tracks and distinct drainage channels were absent, but by 2013 the channelized flow was 1.53 km in length. While channels shapes varied widely, from less than 1 m to more than 4 m deep and up to 20 m wide, the gully formation was morphing the flat terrace toward a dissected, incised surface. Gully heads were evolving at the rate of the gully erosion, generally toward water sources to the south and east ( figure 5(b)), while simultaneously eroding polygons ridges and breaching lowcenter polygons along the way. Thus, the partially breached polygons had their water conservation capacity reduced and the progression of the gullyhead further prevents any eventual recharge of the polygon as the feeding source is captured headward or rapidly drained away during snowmelt. Several breached polygons were observed and resulted in concentrations of dried out Eriophorum angustifolium and Carex aquatilis as a consequence of a rapid drop in moisture conditions. Since 1999 around gully R08p, 202 polygons at least were either breached or partially to completely eroded ( figure 4).
The year 2013, considered particularly dry, was an excellent period to assess the hydrologic response of the gully following the few, relatively short, rain events (figure 6). By graphing the quantity of rain with the estimated volume of water at the gully outlet for an event (figure 7), it was possible to display how the hydrologic system react for each events.
It is generally assumed that the water volume measured at the outlet should increase when it rains, which was the case here, the figure indicating a flashy discharge at each event. The event E1, while featuring the highest estimated flow, was induced by a small rain event (5.5 mm in 72 h) and a smaller input in precipitation. The reason is there was still a high flow from snowmelt at mid-June (figure 6-E1) and a very shallow active layer (less than 5 cm). We infer that by eliminating the snowmelt partition of the flow, the point would be much nearer the trendline (figure 7-E1). Similarly for E4, a very small rain event (3.5 mm in 48 h) occurred at mid-July when the active layer was generally deeper (20 cm deep) in the polygons surrounding the gully. The event barely showed on the cumulative discharge graph (figure 6-E4). The active layer as a reservoir was possibly depleted or had an important potential for recharge following a dozen days without precipitation, between E3 and E4 (figure 6) as displayed by the difference between precipitation input against actual discharge at outlet (figure 7).
The remaining four events (E2, E3, E5 and E6) fitted the drawn trendline with a R 2 of 0.96 (figure 7) and clearly a relation could be drawn between the amount of rain for a given event, the input at the watershed level and the output at the outlet. The slope of the predicted input on the watershed should be corrected by obtaining more precise rainfall data, as there was only one station deployed at the site. Further, onethird of the watershed had a high plateau (up to 300 m), which could potentially increase rainfall measurements variability. Nevertheless, it was possible to assume that the flashiness increased when precipitation increased, suggesting a saturation of the storage with a larger input. The depth of the active layer which should deepen progressively, did not seem to be a factor for those 4 events.
Usually during snowmelt, the shallow active layer was saturated and ponds were commonly observed in low center polygons due to the relative impermeability of the frozen ground, which is common in this type of environment (Woo et al 2008). The snowmelt contribution to the cumulative flow volume at the outlet (figure 6) is pretty clear: the first 6 days recorded account for 59% of the bulk of all the water that would transit at the outlet during the summer 2013, water that will not be available to recharge the wetland water storage. This new state of hydrologic deficit for the breached polygons near the gully and the other polygons in the area of influence of the gully may be lasting, as a significant part of further input of water will be used to recharge the previous deficit (Bowling et al 2003).
In summary, during the relatively dry year of 2013, the two main sources of water input in this system, both snowmelt and rainfall, transited rapidly toward outlet and out of the polygon terrace, suggesting a low residence time for water on the terrace, and a shift from water tracks and sheet flow to channelized flow.

Erosion rates and feedbacks
It was well demonstrated in the literature that the peak of erosion rate and enlargement of a gully occurs during the first fraction of its entire lifetime (Sidorchuk 1999). This model of evolution concord with the rates of erosion observed on R08p, with a observed erosion length of 390 m of the main axis a few month following its initiation (Godin and Fortier 2012b), down to a few dozen meters and less in the last few years (figure 3). The two other gullies which initiated during the 1972-2007 lapse are enlarging at this time but the rate of erosion trend is slowing to a few dozen meters to a few meters per year. Gullies initiated before 1972 are either stabilized, or eroding a few meters to centimeters per year and on their way to become stabilized.
Every active gullies observed were characterized by a localized erosion area, where erosion processes occurred. Sinkholes, tunnels and recent slumps were usually found near an active gully branch, but as time passed and the gullies lengthened and enlarged and the branches eventually stabilized. A stabilized state was usually reached between 5 and 10 years after the moment of erosion. Therefore only a limited fraction of any given gully was actively eroding after the initial peak of erosion. The size of the section varied essentially depending on the amount of secondary branches and stream entry points for instance.
Further, hydrological networks within a basin in a continuous permafrost environment hardly reach full maturity because of limited mechanical erosion due to the physical properties of the frozen ground and therefore seldom incise the streams toward channelized flow (McNamara et al 1999). This typical lack of maturity in an intact section of the terrace may further contribute to the strong positive feedbacks in the time lapse past initial incision.
The group of four yellow dots (figure 3) representing active gullies during 2007-2010 were all within the same range for development rates, suggesting a relatively recent initiation (1990 to early 2000)-following the aforementioned model of evolution.
These currently active gully, while still not in equilibrium yet, are progressing toward that eventual state, similarly to those with a low rate of erosion like most observations between 1972 and 2007 (figures 3 and 4). Stabilization establish a new steady state compared to the pre-gully state. No negative feedback following stabilization can restore the initial state.

Soil losses
Retrogressive thaw slumps headwall can reach annual erosion rates varying between a few meters up to a few dozen meters (Kokelj and Jorgenson 2013). This figure is similar to a developing gully but possibly less than an initiating gully where collapsing tunnels and sinkholes contribute to a very rapid initial erosion due to an important local disequilibrium. While thermo-erosion gullies more likely require a water source to initiate and further evolve, retrogressive thaw slumps do not. Much of the water originating from a thaw slump will be from thawing ground ice where the thermoerosion gully will be a mix of thawing ground ice and runoff. Over an area of 28 891 m 2 (a bit less than 3 ha), 32 093 m 3 were eroded and transported outside the network between 1999 and 2013 from a single gully, and an estimated cumulative volume of 200 000 m 3 for all the gullies in the valley. When taking into account the fractions of eroded permafrost resulting from gullying (Fortier and Allard 2004), melted ground ice was easily put into circulation, mixed with drained runoff water from streams and snowmelt, flowing poorly decomposed peat and aeolian deposits (or colluvial in certain cases) and ultimately either displaced downstream or out of the gully network. Clean water flowing from water tracks and hillslopes transiting through gullies became turbid while transporting fine sediments resulting from a disequilibrium state. This disequilibrium state is only temporary, as the gullied network will eventually stabilize and reach a new equilibrium state. This state occurs when the collapsed blocks, slopes, internal retrogressive thaw slumps and the active layer become stabilized, and no new stream source are captured by the network.
A gully may partially reactivate if a new source of water reach the gully network or a secondary branch, if an exceptionally rapid snowmelt or warm summer occurs or more broadly, if any disturbing factors induced a reactivation.
In any case water may keep transiting through the gullies as one long term impact, several decades up to centuries past its initiation.

Conclusion
Permafrost gullying near the outlet of a small High Arctic watershed enabled not only a change in drainage pattern but also a more efficient drainage of the adjacent wetlands as a direct result of gullying. This increased connectivity enabled the freshet to flow rapidly outside the watershed toward the sea lessening its potential for moisture recharge in the active layer. A flat and topographically continuous wetland is in the process of being punctured by 35 gullies, more than half of those active to this day, therefore establishing a new and increasingly drained equilibrium state for this landscape.