Methane in Zackenberg Valley, NE Greenland: Multidecadal growing season fluxes of a high Arctic tundra

Abstract. The carbon balance of high-latitude terrestrial ecosystems plays an essential role in the atmospheric concentration of trace gases, including carbon dioxide (CO2) and methane (CH4). Increasing levels of atmospheric methane have contributed to ~20 % of the observed global warming since the pre-industrial era. Rising temperatures in the Arctic are expected to promote the release of methane from Arctic ecosystems. Still, existing methane flux data collection efforts are sparse and highly scattered, and further attempts to assess the landscape fluxes over multiple years are needed.Here we use multiyear monitoring from automated flux chambers located on the fringe of a fen area in the center of Zackenberg Valley, northeast Greenland, from July and August (2006–2019). Direct measurements of methane fluxes showed high variability, with mean July–August fluxes ranging from 0.26 to 3.41 mg CH4 m−2 h−1. Methane fluxes based on manual chamber measurements are available from campaigns in 1997, 1999–2000, and in shorter periods from 2007–2013 and have been summarized in several published studies. Fluxes from the multiyear monitoring were combined with fluxes from the most common vegetation types, measured in 2007, and a detailed vegetation cover map to assess the methane flux on a landscape-scale and its variability over time.July–August landscape fluxes, estimated in the current study for the 2006–2019 period, were low compared to previous estimations. For the full study area covering the valley floor, the net methane source during these months was estimated as 0.06 to 0.83 mg CH4 m−2 h−1 and as 0.26 to 3.45 mg CH4 m−2 h−1 for the central fen-rich areas.A 2017–2018 erosion event indicates that some fen and grassland areas along the river in the center of the valley are becoming unstable following pronounced fluvial erosion and a prolonged period of permafrost warming. Although such physical disturbance in the landscape can disrupt the current ecosystem–atmosphere flux patterns, even pronounced future erosion along the river is unlikely to impact methane fluxes at a landscape-scale significantly. Instead, projected changes in future climate in the valley play a more critical role. The results show that multiyear landscape methane fluxes are highly variable at a landscape-scale and stress the need for long-term spatially distributed measurements in the Arctic.


The results show that multiyearmulti-year landscape methane fluxes are highly variable aton a landscape -scale and stress the need for long-term spatially distributed measurements in the Arctic.

Introduction 40
For decades, tundra ecosystems have been subject to attention concerning changingIt was early in climate. The change research recognized that the Arctic is known to be particularly prone to increasing temperatures, and considerable changes to ecosystems at high latitudes were to be expected (IPCC, 1990). Further, high latitude wetlands have, for a long time, been recognized as a significant contributor to the global atmospheric budget of methane (CH4). In the first global budget, the wetland source was estimated at 130 to 260 Tg y −1 of a total 45 atmospheric burden of 529 to 825 Tg y −1 . (IPCC, 1990). For decades, hence, tundra ecosystems have been subject to attention concerning changing climate. Further, high latitude wetlands have, for a long time, been recognized as a contributor to the global atmospheric budget of methane (CH4). In the first global budget, the tundra source was an estimated 1.3 to 13 Tg y −1 of a total atmospheric burden of 529 to 825 Tg y −1 . This term and the overall atmospheric budget have changed remarkably little during nearly 50 years of 50 research, with a more recent estimate ranging from 140 to 280 Tg y −1 . It still forms the background for the concern that with Arctic ecosystems warming, including the wetland source regions, increasing emissions can start a positive feedback in the climate system.emission estimates for all global wetlands ranging from 140 to 280 Tg y −1 in the 1970s as well as in the 2010s . This still also forms the background for the concern that with Arctic ecosystems warming, the so far modest wet tundra emissions may 55 increase and start a positive feedback in the climate system (Knoblauch et al., 2018). In addition to Arctic warming, landscape changes such as permafrost thaw affecting tundra ecosystems and lake formations have been identified as possible hotspots for increased emissions (Schuur et al., 2015;Walter Anthony et al., 2018). (Schuur et al., 2015;Walter Anthony et al., 2018). Also, coastal and further offshore marine sources are subject to possible changed emissions in the Arctic  60 Thornton et al., 2020), but these are not dealt with here. Overall, these concerns led to early and still continuing efforts at quantifying more closely the Arctic natural emissions and their sensitivity and dynamics in relation to climate change.
Studies of tundra methane emissions in several parts of the Arctic were initiated between the 1970s and 1990s Christensen, 65 1993). The . Zackenberg Valley in northeast Greenland was one of the first high Arctic sites to be added to the circumpolar map of methane flux studies in 1997  after the start of the Zackenberg Ecological Research Operations (ZERO) in 1995 .  after the start of the Zackenberg Ecological Research Operations (ZERO) in 1995 . Since then, the valley 70 has seen several different methane flux studies, and methane monitoring became part of the Greenland Ecosystem Monitoring (GEM) program in 2007 . . With numerous short-term research projects, the monitoring has led to the availability of a unique, large number of years with observations of fluxes compared to any other location in the Arctic.
The overarching background for the sizeable emissions in the Arctic is that waterlogged undisturbed soil 75 environments host stable anaerobic environments with optimal conditions for methanogenic activity. Compared with tropical wetlands influenced heavily by the seasonality of flooding (Nisbet et al., 2019), the arctic source areas tend to be more stable geographically. Their emissions subject to the balance between the production at depth and the microbial oxidation in the aerobic surface layer. However, through recent decades, many factors such as nutrients, plant species composition, topography/hydrology have been found to influence will modulate 80 the size of the emissions (AMAP, 2015). From a landscape perspective, the constantly emitting wet soil environments are surrounded by and intermixed with uplands, glaciers, lakes, and rivers, all with their distinct and, in some cases, very different methane flux characteristics. It is rare that a comprehensive analysis of both these small-scale controls and spatial heterogeneity is possible simply due to a lack of locations where enough studies of different kinds have been conducted. Zackenberg Valley is here unique, with such a wide range of 85 studies available.
Here we compile all methane flux studies conducted so far in theZackenberg valley with the objectives to 1) review the combined information from these studies on temporal and spatial variability of methane fluxes in a composite high Arctic landscape and 2) assess the sensitivity of the measured fluxes as they respond to climate warming or local changes. TheA large number of studies and multiple years of observations, along with the 90 addition of flux measurements from a recent gully, provide a unique opportunity to disentangle the effect of different processes and quantify their relative influence on the fluxes at theon a landscape -scale. These processes can broadly be grouped as 1) climate variability and the projections for a gradual warming, 2) increased erosion of vegetated surfaces.
The challenge is to quantify the sensitivity of the landscape fluxes to these factors individually to allow for a 95 quantitative analysis of how the factors combine and compare in terms of sensitivity to established climate warming scenarios. To conduct this study, we use a combination of published data and new measurements of methane fluxes from a recently eroded gully near the Zackenberg Research Station and its immediate surroundings.
A subdivision of the valley study area is Rylekaerene, partly covers the valley floor study area and consists of a 1.3 km 2 patterned wet tundra ecosystem dominated by fen areas and divided by drier patches of heaths and grasslands in the central part of the valley . . Rylekaerene has been a subject of several methane flux studies since the mid-1990s, and several sites within the area have been in use (Fig. 2). 1). 135 The first methane flux measurements in Rylekaerene were carried out in 1997 in studies by Christensen et al. (2000) 145 andFriborg et al. (2000) using manual flux chamber measurements and the eddy covariance (EC) method in the center of the fen area (Fig. 2a). In the south end of Rylekaerene, less than 600 m south of the EC site,  measured methane fluxes in 1999). Close to this site, the AC setup has measured methane and CO2 fluxes along a topographic gradient from a wet fen area and a drier grassland area starting in 2005. The AC site is acentral part of the GeoBasis subprogram of GEM (Mastepanov et al., 2008;Mastepanov et 150 al., 2013).  valley (Fig. 1). The measured methane and CO2 fluxes by gradient EC methods at a site in the center of Rylekaerene in 2008 and 2009. This site (Fig. 2c) is located ~250 m north of the AC site and was later promoted to a permanent EC CO2 flux installation for measurements in the GeoBasis subprogram under the MM2 (MicroMet 2) name (Lund et al., 2008(Lund et al., -2011. 155 In 2007, methane fluxes of the most common vegetation types were measured at 55 plots in a ~600 m 2 area in the center of Rylekaerene, about 300 m north of the AC site , and 50 to 150 m east of the current location of MM2. More recently, several studies focused on areas outside the permanent AC and MM2 sites, including   (Fig. 2d) and   (Fig. 2e) towards the eastern border of Rylekaerene, and Ström et al. 160 (2015) in between the AC site and MM2 (Fig. 2f).
A large dendritically shaped gully formed rapidly during July 2018, most likely triggered by substantial lateral erosion of up to 4.7 m in an outer bend of the river after a glacier lake outburst flood (GLOF) in August 2017 , , right where the gully ends. This gully, previously referred to as 'thermokarst', extends from Gadekaeret, a small fen area close to Zackenberg Research Station (Christensen et al., 165 2020b). The gully is up to 1 to 4 m deep, up to 8 m wide, and ~50 m long in the longest flow direction. This process is different from thermokarst and is similar to the active development of a gully in the northern part of Zackenberg Valley in 1999 (Christiansen et al., 2008). (Christiansen et al., 2008). The retrogressive, branching pattern of the gully most likely occurs along with ice -wedges in the area. Meltwater from a large snowpack located in Gadekaeret in early 2018 has most likely saturated the thawed sediments and thermally eroded the ice -wedges, causing 170 sediment transport from the site and into the river.
The gully was monitored closely during the 2019 field season, but its development had come to a halt.ceased.
Standing water covered the bottom of the gully during August 2018, and elevated methane concentrations were detected in the area (Christensen et al., 2020b). (Christensen et al., 2020b). In 2019, the bottom of the gully had dried out, and methane concentrations were no longer elevated. 175

Fluxes from mobile and automated chambers
Most of the previous studies in the research area were based on mobile flux chambers and stationary automatic chambers, which utilized measured changes in methane concentration over time inside a chamber to estimate surface fluxes. The chamber was sealed off to all sides but the bottom during each measurement, isolating the flux 180 estimates to a specific time and area, analogous to the approach described by, e.g., Hutchinson (1995). Methods vary between the studies, with differences in gas analyzers, chamber designs,   1 shows the coverage of each surface class compared to the valley floor and Rylekaerene areas in both absolute and relative terms, based on hyperspectral aerial data from , which improved on a 195 manual land cover classification by . Zackenberg River has changed its course since 2000 after multiple GLOF events , and we adjusted the land cover classification map to fit the extent of the river in a 2014 orthomosaic of the entire study area (COWI, 2015).

Climatology
The study area is as high Arctic , with mean annual air temperature and precipitation of −8.6 °C and 253 mm (2008-2018) (López-Blanco et al., 2020). Continuous permafrost underlies the valley (Christiansen et al., 2008), and the active layer (AL) reaches a maximum depth of 0.58 to 0.85 m, increasing at a 205 rate of 0.74 cm year −1 (1996-2019) at the ZEROCALM-1 site (Christensen et al., 2020a) near the climate station at the valley floor.

Measurements in Rylekaerene
The first methane flux measurements in Rylekaerene were carried out in 1997 in  and 210 Friborg et al. (2000). They used manual flux chamber measurements and the eddy covariance (EC) method in the center of the fen area (Fig. 1, site a). In the south end of Rylekaerene,  measured methane fluxes in 1999-2000 ( Fig. 1, site b). An automated chamber setup was added nearby in 2005. The automated chambers (AC) measure methane and CO2 fluxes along a topographic gradient from a wet fen to a heath area. The AC site is a part of the GeoBasis subprogram of GEM (Mastepanov et al., 2008;. 215  measured methane and CO2 fluxes by gradient EC methods at a site in the center of Rylekaerene in 2008 and 2009. This site (Fig. 1, site c) is located ~250 m north of the AC site and was later promoted to a permanent EC CO2 flux installation for measurements in the GeoBasis subprogram under the MM2 (MicroMet 2) name (Lund et al., 2008(Lund et al., -2011.
In 2007,  measured methane fluxes of the most common vegetation types in a ~600 m 2 area 220 in the center of Rylekaerene, 50 to 150 m east of the current location of MM2.
More recently, several studies focused on areas outside the permanent AC and MM2 sites, including   (Fig. 1, site d),   (Fig. 1, site e) toward the eastern border of Rylekaerene. Ström et al. (2015) measured fluxes in between the AC site and MM2 (Fig. 1, site f). The length of these campaigns, their onset compared to the beginning of the growing season, the sampling area, and strategy vary between studies. 225 Figure 2 shows a timeline of methane flux campaigns from published studies and data from the GEM database.

230
Most of the previous studies in the research area were based on mobile flux chambers, and stationary automatic chambers, which utilized changes in methane concentration measured over time inside a closed chamber to estimate surface fluxes. The chamber was sealed off to all sides but the bottom during each measurement, isolating the flux estimates to a specific time and area, analogous to the approach described by, e.g.,  and Livingston and Hutchinson (1995). Methods vary between the studies, with differences in gas analyzers, chamber designs, 235 measurement time, replication numbers, sampling frequency, and the length of the study periods. Details on the methods are summarized in Table 2.  LGR Greenhouse Gas Analyzer For the aims of the current study, we analyzed existing data of July-August CH4 fluxes from manual mobile chamber measurements from 2007, published in , and GeoBasis Zackenberg automatic chamber measurements (2006( published in Mastepanov et al. (2013, and 2011-2019 in Mastepanov et al. (in prep.). Details on the methods are described in these publications. Only parts of the AC flux time series were 245 used in this study; the July-August flux includes the peak of the growing season in the valley, and it matches the timing of most of the previous studies.
For the aims of the current study, we combined existing data of growing season methane from the studies shown in Fig. 2 from 2006 and forward. Only parts of the AC flux time series were used in this study; the July-August flux includes the peak of the growing season in the valley, and it matches the timing of most of the previous studies. 250 In addition to this data, we carried out 2.

Measurements in the gully area
We conducted 113 measurements from 23 August 2019 to 1 September 2019 in a grid covering the recent gully (see Fig. 2) and extending toward a small lake in the Gadekaeret fen area, close to Zackenberg Research Station.
The measurements took place between 9 a.m. and 7 p.m. on 43 plots located 4 to 5 m from each other in the three 255 main surface classes in the area: fen, grassland/barren area, and the recently eroded area. Several flux measurements at each surface class were performed every day.
A metal collar was carefully installed on the ground on each plot; a dark acrylic chamber was placed over the collar for a minimum of five minutes. The footprint of the collar was 0.07 m 2 ; the height from the soil surface to the top of the chamber was recorded for each measurementThe measurements took place between 9 a.m. and 7 p.m. on 43 260 plots located 4 to 5 m from each other in the three main surface classes in the area: fens, grasslands, and in the gully. Several flux measurements at each surface class were performed every day.
A metal collar was carefully installed on the ground on each plot; a dark acrylic chamber was placed over the collar for a minimum of five minutes. The footprint of the collar was 0.07 m 2 ; the height from the soil surface to the top of the chamber was recorded for each measurement (. The height ranged from 0.26 m to 0.48 m, depending on the 265 surface topography). The chamber was equipped with a fan inside and had a 3 mm vent on the side. A gas analyzer (Ultraportable Greenhouse Gas Analyzer, Los Gatos Research, USA) was connected to the chamber with a pair of 15 m long HDPE tubes. The chamber was equipped with a fan inside and had a 3 mm vent on the side. A gas analyzer (Ultraportable Greenhouse Gas Analyzer, Los Gatos Research, USA) was connected to the chamber with a pair of 15 m long HDPE tubes and. The gas analyzer was continuously measuring CH4methane concentration in 270 the chamber headspace at 1 Hz frequency. After each sample, the chamber was ventilated for at least two minutes until the methane concentration was down to ambient concentration.
After each sample, the chamber was ventilated for at least two minutes until the methane concentration was down to ambient concentration.

Hyperspectral classification
The hyperspectral remote sensing data were collected on 8 August 2000 by an airborne HyMap campaign (Palmtag et al., 2015). The image data were processed into a 5 m × 5 m resolution land cover map of Zackenberg Valley , improving on a manual land cover classification by . Zackenberg River has changed its course since 2000 after multiple GLOF events Tomczyk and Ewertowski, 280 2020), and we adjusted the land cover classification map to fit the extent of the river in a 2014 orthomosaic of the entire study area (COWI, 2015).

Environmental changes
Monitoring of air (2 m above ground) and soil temperature (0.20 m below ground) were summarized from 1997-285 2019 data (60 min resolution) from the nearby climate station (

Measurements from mobile and automated and manual chambers
Data from several sources were included in the calculation of a timeline of landscape fluxes in Zackenberg Valley.
Previously published and unpublished data were compiled to estimate fluxes and SE on the six surface types (Table   3), combining mobile flux measurements and flux measurements from AC. The measurement methods are described in detail in their respective publication, but for the AC, additional steps were added after applying the 300 same approach as . The AC flux time series was separated into two datasets: one representing a long time series from a 10 m wide transition zone at the fen fringe (chambers 1 to 6). In this area, chambers with higher numbers are generally drier. Another time series represents the changes since 2012 in the four outer chambers, located further out into the fen (chambers 7 to 10). The flux data do not show a diurnal pattern for July and August, so all available data were used in those two months. The flux measurements were first averaged 305 for each chamber. Mean chamber fluxes were then further averaged into a single mean methane flux for each year for all the six inner chambers. The same was done for the four outer chambers from 2012-2019.
In 2006, only four chambers were operating, and no data were available for chambers 4 and 6. Potential differences in methane flux are corrected by multiplying the mean of the available 2006 data by a coefficient based on 2007 data. This coefficient was found by dividing the mean in chambers 1, 2, 3, and 5 by the mean of all six chambers. The mean July-August flux for chambers 7 to 10 were then combined with the mean growing season fluxes from the other studies listed in Table 3. In the valley-wide vegetation cover map, hummocky and continuous fen were not separated into different classes, even though mean fluxes differ substantially for these two surface types (seein, e.g., Table 3 in ). ATagesson et al. (2013) and Table 1 in .
Each mean flux measured in fen areas was paired with the mean flux measured at the fen fringe. Using R (R Core 320 Team, 2021), we applied unweighted Deming linear regression on the data (n = 18, jackknife method, error ratio = 0.44, Pearson's r = 0.64, p-value threshold = 0.05). The approach accounts for uncertainties in both the fen fringe data and in the fen data to estimate a time series for the fen surface types for 2006-2019. The measurements used in the regression for the fen areas are summarized in Table 3.

Fen fringe flux 325
The mean flux from the fen fringes was estimated from the mean July-August flux measured every year in the 2006-2019 period using all the available data from chamber 1 to 6 ( Table 3). The timeline of mean fluxes at the fen fringe represents the outer 10 m of all fen surfaces in the valley, shown in Fig. 1 without further processing.

Grassland flux
The grassland fluxes are held constant over the time series, with data from  as input to the 330 calculation (Table 3), while grassland fluxes from  are omitted due to high spatial variability and a higher average flux.  found a significant methane uptake on unvegetated fell and barren surfaces in Zackenberg Valley in 2012. These measurements were grouped into a broader 'dry tundra' class that includes Dryas heath. 335

Fell and barren fluxes
Here, we use the mean methane flux from the 'dry tundra' class but only applying it for the fell and barren areas.
The mean flux from these surfaces is held constant in our landscape flux time series.

Heaths and Salix fluxes
The heath class in this study includes both Cassiope, Dryas and Vaccinium areas, and data from , , and  are combined to calculate an average estimate for fluxes 340 in these areas, which are held constant over time.  used different groups in their study, where different types of heath areas fall into 'dry' and 'moist' tundra (Cassiope and Salix heath). We calculated the weighted mean flux, taking the different sample sizes for the separate land cover classes into account, was calculated and used for this study. for each study for the heath class, and then we calculated a mean and pooled SE.
Likewise, the Salix snowbed class was calculated in the same way as with heath, with data from the same three 345 studies, but with only the flux from 'moist' tundra from .
The fluxes at the AC site were calculated from each automatic chamber measurement, using the same approach as . The flux time series was separated into two datasets: one representing a long time series from the fen fringe (chambers 1 to 6). Another time series represents the changes since 2012 in the four outer chambers, located further out in the fen (chambers 7 to 10). The flux data do not show a diurnal pattern for July 350 and August, so all available data for each chamber were used in those two months. Data from chambers 1 to 6 were first summarized for each chamber as temporal means, highlighting the differences in fluxes from each of the six inner chambers. July-August data were then further summarized into a single mean methane flux for each year. In 2006, only four chambers were operating, and no data were available for chambers 4 and 6. The difference is corrected by multiplying the mean of the available 2006 data by a coefficient based on 2007 data. This coefficient 355 was found by dividing the mean in chambers 1, 2, 3, and 5 by the mean of all six chambers.  (1)

Land-cover classes
The fourlandscape flux time series 370 The six combined surface classes representoccupy most of the study areas (Rylekaerene and the valley floor, Table 1). In Eq. (2), theThis static areal coverage (Areaclass) of the surface classes was calculated using QGIS v.
The area-weighted flux was calculated for the two study areas, the valley floor, and Rylekaerene (Areatotal)., with 375 total areas of 15,905,003 m 2 and 1,271,303 m 2 . This approach assumes no fluxes of methane flux in the remaining parts18.1 % and 3.7 % of the study areas, including the river, delta, riverabrasion plateaus, and primarily barren areas (Fig. 2).boulder fields in the 'other' category in Table 1. The area-weighted landscape flux (Fluxarea weighted, 2007) calculation is shown in Eq. (2),is calculated for each year, with fluxes in the fens and classes include fen, grassland, Salix snowbeds, and heathfen fringes being the only to change over time. (2)

Combining data into a landscape flux time series
An area-weighted landscape methane flux for each measurement year (Fluxyear) is estimated by combining Eq. (1) and Eq. (2)

Fluxes in the gully area 390
In the recently eroded gully area, the fluxes in 2019 were calculated using the ordinary least square linear (OLS) regression described in Pirk et al. (2016a). Of the 113 measurements, 102 had a significant (p < 0.05) regression slope. The remaining 11 measurements were found on both the grassland areas and on recently eroded surfaces.
They showed an increase in concentration over time close to zero, and they are interpreted as areas with no flux.
The 11 measurements are included in the calculation of the mean flux of the recently eroded surfaces in the gully. 395 The mean flux and the standard error (SE) for all measurements were calculated by calculating the mean flux for repeated measurements for each plot. In the recently eroded gully area, the fluxes in 2019 were calculated using the linear flux model (Pirk et al., 2016a). Of the 113 measurements, 102 had a significant (p < 0.05) regression slope. The remaining 11 measurements were found on both the grassland/barren areas and on recently eroded surfaces. They showed an increase in concentration over time close to zero, and they are interpreted as areas with 400 a zero flux. The 11 measurements are included in the calculation of the mean flux of the recently eroded surfaces in the gully.
The mean flux and the standard error (SE) for all measurements were calculated by calculating the mean flux for repeated measurements for each plot. Further averaging for their respective surface class (recent erosion, barren/grassland, fengully, grasslands, fens) was done afterward, showing how fluxes can change in an area after 405 an erosion event. Flux measurements from the gully area are limited to 10 days in the late growing season. The impacts on methane are estimated from a linear development in eroded surfaces from 0 m to 25 m and from 0 m to 100 m on each side of the Zackenberg River. The development change areas that act as sinks or sources of methane into eroded areas, using the mean flux measured on the recently eroded surfaces in the gully in 2019. The Zackenberg River was digitized to its extent in August 2014, when the valley was mapped in high detail (COWI, 2015). Eighty-five buffer zones were created on each side of the river, each representing one year. These buffer 425 zones had a total width of 25 m and 100 m. The land cover types in each buffer zone were progressively eroded for the sensitivity study, and the associated change in methane flux was subtracted from the mean valley methane flux.

Changes in methane flux from increasing temperatures
In their study,  compared the potential methane flux in the Zackenberg area modeled for the late 21 st century with present-day conditions. They found an exponential growing season temperature-methane flux 430 relationship based on Zackenberg and Kobbefjord data from 2008-2015. Present-day temperatures were modeled for the 1991-2010 period, and late 21 st century temperatures were modeled for 2081-2100 under the Representative Concentration Pathway 8.5 (RCP8.5). Both present-day and future mean temperatures were modeled using the ECHAM5 general circulation model to limit the cold bias in the model, and the relative increase in methane flux was +141 % (Geng, personal communication, 24 April 2020). The relative increase ranges from +114 % to +171 435 % when considering the lower and upper 95 % confidence bounds, see Fig. 3 in .
We assume a linear increase in valley methane flux of +141 % for this sensitivity study, ranging from +114 % to +171 % from 2016-2100, even though this period does not fully match the modeled periods, i.e., from 1991-2010 to 2081-2100.

Sensitivity to increasing erosion activity 440
We establish three pathways to quantify how erosion could affect the landscape flux in Zackenberg Valley. In the first pathway, we calculate the impact on the mean valley flux if the eroded areas are growing at an annual rate equivalent to the size of the recent gully (720 m 2 y -1 ). In the second and third pathway, the eroded area starts at 720 m 2 every year and grows to 5 and 10 times that area per year, respectively, i.e., a linear increase over 85 years from 720 m 2 y -1 to 3,600 m 2 y -1 , and from 720 m 2 y -1 to 7,200 m 2 y -1 . In this calculation, the erosion can only occur in 445 zones with excessive ice-rich permafrost near rivers and streams. To identify these zones, we use GIS data available from Cable et al. (2018). Inside our study area (Fig. 1), we identify parts of the landscape that likely have 'excess ice-rich top 1 m permafrost' and are located near 'rivers', which include both rivers and tributaries, see Fig. 12 in Cable et al. (2018). Based on observations from the recent gully, erosion can occur up to 50 m from the rivers.
Inside the identified zones, we model erosion of the size specified above and the fractional cover of surface classes. 450 The reduction in mean flux for the eroded areas is subtracted from the landscape flux. As time goes, the eroded areas increase, causing larger impacts on the landscape flux, while all fluxes increase at the same rate as when the mean temperatures increase.

Revegetation of eroded areas
A similar gully in the northern part of the valley developed in 1999 (Christiansen et al., 2008), and it shows signs 455 of revegetation on ~40 % of the eroded areas to grassland after 20 years. This estimate is based on a visual interpretation of 100 points randomly scattered over eroded parts within the gully in an orthophoto from August 2019. From the estimate, we set the regrowth rate equal to 2 % year -1, and this rate of regrowth is included in the overall calculation. Figure 3a shows the development of July-August soil and air temperatures since 1997 at 0.2 m depth and 2 m height measured at the climate station. The mean July-August temperature is 5.9 °C for air temperatures and 4.7 °C for the available soil temperatures. Air temperatures increased by 0.07 °C year -1 (1997-2019, n = 23, Pearson's r = 0.43, p = 0.04), while no significant linear trend is observed in the soil temperature data during the same period. 465

Temporal variability in fluxEnvironmental changes
The timing of snowmelt in the study area (Fig. 3b) shows substantial variations between years, from 30 (DOY 150) to 27 July (DOY 208). During July and August, there is so far no clear trend toward drier or wetter conditions for both the heath areas (Fig. 3c) and the fens at the AC site (Fig. 3d) from the mid-2000s and forward. Measurements from chamber 6 show generally drier conditions than from chamber 1, which is located further out in the fen.  .
Treating all six chambers as replicates reveal high temporal variability for the area (red points and line in Fig. 3).  These chambers, named chambers 7 to 10, generally show higher fluxes than those measured at chamber 1 in 485  and Mastepanov et al. (in prep.).

Surface cover
The area of four surface classes of the valley floor and the fen-dominated Rylekaerene was calculated based on 495 HyMap data from August 2000. The surface classes used in this study are a subset of the HyMap dataset, combining heath classes and excluding 'Other' surface classes. This class includes barren areas, lakes, and rivers. The four surface classes cover most of the valley (Fig. 2), with the notable exemptions of the Zackenberg River system and the landing strip. The four surface classes cover 83.1 % of the valley study area and 96.4 % of the Rylekaerene study area, where most methane studies in the valley were carried out (Fig. 2). 500 Table 2 shows the coverage of each surface class compared to the valley floor and Rylekaerene areas in both absolute and relative terms. The fen surface class takes up 12.1 % of the valley floor while accounting for 50.4 % of the Rylekaerene study area. Grassland and heath areas are common (28.5 % and 21.8 %) in the valley floor area while being relatively less important in Rylekaerene. Salix snowbeds take up 20.7 % and 21.1 % in the two study areas. 505 The relative composition of the two study areas is used for estimating the area-weighted methane landscape flux for the time series.

Gully methane fluxes
The methane flux at the exposed, eroded surfaces of the gully was different from the flux on the nearby, undisturbed surfaces (Fig. 4). from a small methane sink to a small source for this area. The source followed an initial substantial episodic release of methane stored in exposed ice in the year when the gully appeared (Christensen et al., 2020b).    the fen fringe, with mean fluxes ranging from 1.75 ± 0.27 to 8.3 ± 0.31 mg m −2 h −1 in 1997  23 and 2000 . Methane fluxes in the fen fringe ranged from 0.26 mg ± 0.07 m −2 h −1 in 2011 to 3.41 ± 0.57 m −2 h −1 in 2007 during July-August Mastepanov et al., in prep.).

Gully methane fluxes
Methane fluxes at exposed, eroded surfaces of the gully were different from the fluxes on the nearby, undisturbed surfaces (Fig. 5). The late growing season mean methane flux of the recently eroded surfaces in the gully in 2019 600 was 0.05 ± 0.02 mg m −2 h −1 , including positive and negative fluxes. The grassland surface cover in areas not disturbed by erosion shows a negative methane flux of −0.06 ± 0.01 mg m −2 h −1 . The mean methane flux in the fen was 3.83 ± 0.76 mg m −2 h −1 , more than 75 times higher than the mean flux in the gully. For several plots in the fen, the flux was highly variable over time, reaching 20 mg m −2 h −1 . Generally, the methane flux decreased with air temperature in the fen. In the gully, plots with patches of live vegetation showed mostly negative flux, while barren 605 plots exhibited a positive flux. In the undisturbed grassland areas, vegetated and barren plot fluxes were negative.
Methane emissions generally increased closer to the open water body of the nearby fen. The emergence of the gully changed the area from a small methane sink to a small source, following an initial substantial episodic release of methane stored in exposed ice in 2018 when the gully appeared (Christensen et al., 2020b).

Estimation of an integrated flux of methane in Zackenberg Valley
The July-August landscape flux of methane exposed large interannual variability over the valley floor and Areas with a high positive flux are, in relative terms, more dominant in Rylekaerene than the valley floor, which instead has a more considerable negative flux from heaths, Salix snowbeds, and fell and barren areas, which decrease the mean landscape flux. 630

Methane flux measurements
Several methane flux methods have been in use since the two first methane flux studies were carried out in 1997 ( Fig. 4) in the Rylekaerene area Friborg et al., 2000), including both manual chambers (Table 2) and EC measurements. Flux measurements over 14 years at AC show large interannual variability in methane fluxes even at a single site. Substantial spatial variability in fluxes is seen in Fig. 4 when studies are 660 available in the same years. Differences in methane fluxes measured from these studies can be explained by site characteristics, equipment and sampling strategies, and both length and timing of the measurement period. Two studies are omitted in Fig. 4: Pirk et al. (2016b) studied the methane fluxes in the fen areas near the AC site under snow-covered conditions before the growing season in 2012 and 2014.  added flux measurements from 2012 covering only 'dry tundra' (Dryas heath, abrasion plateau, and fell field) and 'moist 665 tundra' (Salix snowbeds and Cassiope heath) sites, both of them methane sinks. These findings exceeded the methane uptake on similar tundra surface classes in the area in . Higher soil and air temperatures in 2012 compared to 1997 could explain the higher fluxes (Fig. 3a), which also fits with relatively high methane uptake in 2007 when  did their measurements on corresponding vegetation types. Spatial differences and the inclusion of unvegetated surfaces (abrasion plateau and fell field) and more 670 advanced measurement equipment (Table 2)  and fell field) and moist tundra (Salix snowbeds and Cassiope heath). These findings confirmed the atmospheric sink of the dry tundra heath, first reported by . , which provide data for the landscape fluxes in this study, used a modeling approach for calculating a flux time series in Rylekaerene covering 1997, 2000, 2002, and 2007-2010. Their study was based on both in situ flux measurements and modeling using environmental parameters and satellite imagery, meaning their results are not directly comparable to the in 690 situ measurements shown in Fig. 8. TheChristensen et al. (2000) measured fluxes for an intensive study area, a part of Rylekaerene, in 1997. In Fig. 4, only the combined measurements from hummocky and continuous fen plots are used. Alongside, Friborg et al. (2000) applied the EC method and found a much lower methane flux compared to  in 695 approximately the same area. The fetch of the EC tower included surface classes with lower flux, which is a likely reason for lower mean flux Rylekaerene reported in Friborg et al. (2000). The length of the measurement period was longer in Friborg et al. (2000), starting 1 June (Fig. 2), more than three weeks earlier than , which further explains some of the differences as methane fluxes were low at the beginning of the measurement period.  measured methane fluxes two times per week per plot, while 700 methane fluxes from EC were available every day for several weeks in July and August. The SE for the fen measurements in  in Fig. 4 is high due to the high flux variability within and between the two combined groups, hummocky and continuous fen.
Control plot flux measurements from  are illustrated in Fig. 4 for a study site in 705 the southern part of Rylekaerene in 1999 and 2000. As shown in Table 2, the methods are comparable to , and the mean flux in both years is similar. The SE are smaller, likely due to both fewer and less diverse plots, i.e., six continuous plots compared to a total of 18 hummocky and continuous plots in .

710
Two studies provide a basis for comparison in 2007, the growing season with the highest estimated landscape flux in the 14-year study period.  and  applied similar methods (Table 2), using plots near MM2 in the same year with comparable methane flux as a result. For , the SE reported in Fig. 4 is comparatively high. This difference can be explained by the high SE associated with measurements of plots with high flux and high Eriophorum scheuchzeri coverage, combined with relatively low methane flux plots 715 with low Eriophorum scheuchzeri coverage. In , measurements were done at 25 randomlyplaced plots in hummocky and continuous fen areas, which resulted in a lower SE for the combined fen classes.
The key difference is the range of both high and low Eriophorum scheuchzeri plots in  against the randomly selected placement in two fen classes in . reported mean methane fluxes larger than those measured at the fen fringe at the AC. Measurement methods are 730 comparable to  and . They use an identical gas analyzer and similar measurement frequency. Different chamber sizes are used, but this is not expected to impact the measurements significantly, as chamber volumes are considered in the calculation of fluxes. Differences in mean flux between the three years were attributed to variable measurement periods in the three years and interannual variations in water level depth . 735 Ström et al. (2015) measured methane fluxes with the same approach as . The main differences between the two studies are the timing and the site. The two sites are situated ~600 meters apart, which explains possible differences in fluxes due to water levels and substrate availability. Measurements in 2011 in  began 14 days earlier, which may have caused some of the difference in mean fluxes. In 2012, fluxes were 740 similar between the two studies. The two campaigns started one day apart in that year, indicating that conditions are similar despite the distance.
The outer four chambers of the AC installation show interannual variability in fluxes of a size similar to the manual chamber measurements of , Ström et al. (2015), and EC measurements in . 745 The number of observations for chambers 7 to 10 is high (Table 3), and flux measurements from each chamber are available every 90 minutes during July and August. In some years, measurements are not available for the entire period. Differences in chamber fluxes may be explained by variability in water levels, but data not available from the outermost chambers. Standard errors are higher due to relatively high variability in fluxes during the measurement period. 750 The previous studies show large differences in methane flux within the same fen ecosystem. Local spatial variability, interannual differences may explain these differences, and possibly also by differences in methods.
Spatial variability in methane flux can be pronounced in tundra ecosystems with complex microtopography  and differ significantly even at a meter-scale (Fig. 3). Methane fluxes vary with local 755 differences in substrate, water-table depth, grazing, and vegetation composition and productivity, each of which can either increase or limit the methane flux. Several of these interactions have been studied in Zackenberg Valley: increased substrate availability, mainly acetate, contributes to higher methane fluxes as shown by, e.g., , while a low water table limits methane fluxes . Grazing can either increase  or decrease  the methane flux, depending on the vegetation cover. Vegetation 760 composition and primary productivity are strong drivers of methane fluxes Ström et al., 2015).
Spatial variability in methane flux can be pronounced in tundra ecosystems with complex microtopography  and differ significantly even at a meter-scale (Fig. 4). Methane fluxes vary with local differences in substrate, water-table depth, grazing, vegetation composition and productivity, each of which can 765 either increase or limit the methane flux. Several of these interactions have been studied in Zackenberg Valley: increased substrate availability, mainly acetate, contributes to higher methane fluxes as shown by, e.g., , while a low water table limits methane fluxes . Grazing can either decrease  or increase  the methane flux, depending on the vegetation cover. Vegetation composition and primary productivity are strong drivers of methane fluxes Ström 770 et al., 2012;Ström et al., 2015).
Temporal variability in methane fluxes can be caused by differences in environmental parameters, both within a growing season and from one year to another. Soil temperature was found to explain less variability than species composition and primary productivity , , while soil temperatures showed a high correlation with methane flux within most individual years 775 . . Late-lying snow delays the beginning of the growing season , , which controls several of the above-mentioned parameters mentioned above.. Only a few studies in Zackenberg Valley were conducted outside the growing season. However, the fall season could substantially impact the annual methane budget, mainly through emissions associated with the onset of soil freezing (Mastepanov et al., 2008). (Mastepanov et al., 2008). In contrast, wintertime (November-780 May) emissionemissions may only have a limited impact on the annual methane budget (Pirk et al., 2016b). (Pirk et al., 2016b).
Our analysis focuses on the mean of the July-August fluxes in the valley, but the automated chambers are running from snowmelt to the end of the field season (Fig. 72). The fixed period matches the timing of the previous studies, and the period showed a good representation of the mean flux of the entire measurement dataset. The first 30 to 40 785 days after snowmelt have been shown to express the main differences between years , , which is covered by the July-August meansmean fluxes to a large extent.

LandscapeCombining flux measurements from multiple sources
Comparable mean methane fluxes across several existing studies indicate that the differences in methods and placement have an impact on the flux in fen areas. Differences arise mainly when measurement periods vary by 790 several weeks and when measurement plots are distributed over several different surface classes, i.e., hummocky and continuous fen. However, the fen fringe time series appears to express some of the variability in fluxes as the separate measurement campaigns when these are combined into one. Regression analysis shows a significant correlation between methane fluxes at chambers 1 to 6 and the grouped studies from the fen areas of Zackenberg Valley when pairing data for the same years. Deming regression is a reasonable choice of analysis when both the 795 mean values of chamber 1 to 6 and the fen methane fluxes are associated with measurement errors and uncertainties.
The composite fen dataset, made from six separate studies, increases the temporal coverage and the spatial coverage of the data used in the regression analysis. The SE of the calculated fluxes is estimated with jackknife resampling, providing a reasonable SE estimate for Deming regression (Linnet, 1990). When comparing flux data to OLS regression, the Deming regression slope is steeper. The steeper slope means that the estimated fen flux in 2007, 800 when the flux was high, is ~20 % higher when compared to OLS, which means that the difference in regression methods has a substantial impact on landscape fluxes.
The AC site is located by the outlet of the Rylekaerene fen. The substantial flow of water through the area affects the water level, particularly at the innermost six chambers, which are located along the slight topographic gradient at the fen fringe. The changing water level may control the methane flux and explain its high variability 2006-805 2019, but neither the water levels, the air, nor soil temperatures correlate significantly (p-value threshold of 0.05) with methane when analyzing the interannual variability for July-August mean values. Over the 14 years, p-values range from 0.29 to 0.56. The lack of correlation over the time series illustrates a complex interaction between methane and environmental conditions when analyzed on a decadal scale. Figure 8 shows the methane flux of previous studies, which have mainly focused on wetter parts of the fen where the flux is high. The flux appears to 810 vary relatively less over time further out in the fen, as the soil moisture conditions change less between years than the chambers 1 to 6 at the AC site. One example is the smaller relative variability at chambers 7 to 10. The sizeable interannual variability shown in Fig. 3 could be more common for the smaller, discontinuous fen areas in the valley seen in Fig. 2. These areas are characterized by a less stable inflow of water than Rylekaerene, which may therefore cause high interannual variability in the transition zones between different vegetation types. 815 Methane fluxes measured in chambers 1 to 6 are consistently lower than the fluxes in chambers 7 to 10. Therefore, data from chambers 1 to 6 are split into a separate surface class, the fen fringe. The addition of these lower flux fen areas impacts the landscape flux, but the width of the zone does most likely fluctuate across the valley. Here, 10 m 820 is a simple estimate based on the situation at AC, but the width of the zone ultimately depends on the topographic gradient along the edge of the fens.
The remaining surface classes contribute with the same negative fluxes throughout the period, with a mean value representing each of the surface classes. Their net effect is offsetting some of the positive landscape methane flux, particularly for the valley floor. Several surface classes were studied in detail over a single growing season (Table  825 3), and the specific sites differ between studies. The fluxes from these areas are held constant in lack of more detailed data, and fluxes from different studies are averaged into a single mean flux, which represents a range of sites and environmental conditions. We omitted surface flux data from the grasslands, measured by . The measured mean flux was significantly higher at 2.1 ± 1.6 mg m −2 h −1 , and the measurement site was located near the border of the fen, see Fig. 2 in , making the data more comparable to the 830 fen fringe than the grasslands surface class.
Manual measurements in the gully area were limited to 10 days in the late growing season. Still, the size of the mean fluxes is in good agreement with the concurrent fluxes measured at AC 7 to 10 (Mastepanov et al., in prep.) and the grassland flux reported in . Negative fluxes were measured on plots with patchy 835 vegetation in the gully, while mineral-rich plots had low positive fluxes of methane. The highest methane emissions in the gully were found on two unvegetated plots, located where the gully had recently eroded (up to 0.85 mg m −2 h −1 ). These sparse observations suggest a considerable reduction in flux on eroded gully surfaces over time when limited organic soils remain in the gully.

Methane flux upscaling 840
Vegetation plot measurement on the dominant vegetation classes in 2007 combined with long-term AC site data from 2006-2019 enables the calculation of landscape flux time series for Rylekaerene and the valley floor. This approach is more direct than modeling based on physical parameters that have otherwise been shown to correlate well with methane flux for Rylekaerene . . However, a limitation of the approach is its lack of a dynamic component, which could take spatial differences across the valley into account, 845 e.g., snow cover or changing soil moisture conditions from one year to another, affecting the valley floor and the AC site differently. )Both Christensen et al. (2000 and  found a significant difference in the methane flux in hummocky and continuous fen areas in Rylekaerene and treated the two groups separately.
These groups were combined in this study to match the existing, valley-wide surface cover classification. 850 Comparing the mean flux infound for Rylekaerene over approximately the same period shows good agreement between the flux estimates in this study and the flux in .fluxes reported in . Four growing seasons (2007)(2008)(2009)(2010) are common for the two studies, and in three yearsof the seasons, the flux estimates of this study lie within the model uncertainty in . . The one notable exemption is the 2007 mean flux (3.45 ± 0.3126 ± 1.15 mg m −2 h −1 ) in this study, and it differs a lot 855 from the 1.6 ± 1.0 mg m −2 h −1 modeling result. in . The difference can be explained with a relatively high flux measured at the AC site in a year without extreme temperature or moisture conditions in the valley (Fig. 12), which are central parameters in the modeling of . .
The AC site is located by the outlet of the Rylekaerene fen. The flow of water through the area affects the water 860 level at the site, which is shown in Fig. 3d. Water level data was measured manually, once per day in the growing  Figure 4 shows the methane flux of previous studies, which mainly focused on wetter parts of the fen where the flux is high. The flux appears to vary relatively less over time further out in the fen, as the soil moisture conditions change less between years than the chambers 1 to 6 at the AC site. One example is the smaller relative variability at chambers 7 to 10. Similar, large variability could common for the smaller, discontinuous fen areas in the valley seen in Fig. 2. These areas are characterized by a 875 less stable inflow of water than Rylekaerene, which may therefore cause high interannual variability in the transition zones between different vegetation types.  based their study of Zackenberg Valley methane flux on chamber measurements from a c. ~0.1 km 2 area of the northern part of Rylekaerene (Fig. 21). The methane flux was measured 880 with chambers on the five dominant vegetation types, and a mean flux for the valley was calculated by scaling fluxes to match the land cover classification of . . The upscaled methane flux for the entire valley floor was 1.9 ± 0.7 mg m −2 h −1 , which is higher than in any of the years in this study (, which range from 0.0618 ± 0.0305 mg m −2 h −1 and 0.8467 ± 0.2124 mg m −2 h −1 ).. The distribution of land cover classes differs slightly between   and the HyMap dataset, which explains only some of the differences. However. 885 Still, a primary cause for the much higher valley estimate is the higher measured fluxes in the widespread Grasslandgrassland class with a flux of 2.9 ± 1.6 mg m −2 h −1 in , , compared to 0.1 ± 0.04 mg m −2 h −1 (SE converted from SD) in .   (Table 2). As a result, the magnitude of the fluxes in Fig. 5 may be underestimated, as both positive and negative fluxes from  are lower than found in both  and . Additionally, the year 2007 was chosen for the comparison because it exists in both datasets. However, the methane flux was exceptionally high in that year (Fig.  900 3), which means that all other years in the time series were attributed to a lower flux. The valley methane fluxes are highly variable between years and do not show an increasing trend from the available data. Zackenberg Valley shows the potential for increased methane fluxes in the 21 st century, as methane shows a positive correlation with temperatures which are expected to increase under the RCP8.5 climate scenario , similar to the trend of the rest of Greenland and the Arctic (AMAP, 2017). 905  found a relatively high methane uptake on dry tundra in 2012, e.g., Salix snowbeds and Heath surfaces, compared to what has been found earlier. When including the uptake from , the valley flux is reduced, but the effect could vary between years with larger sinks in dry, warm years, as these surface classes combined account for more than 40 % of the valley floor (Table 1). Soil temperatures were lower 910 in 2012 than in 2007 in July-August, which does not explain the higher uptake on dry tundra soils in 2012, indicating that some dry tundra surfaces have a higher methane uptake than others. Additionally, flux data from 2007 are used for both heath, grassland, and Salix snowbeds surface classes. While flux measurements in the fen areas were exceptionally high in that year, this was not the case for the drier surface types. Landscape fluxes from both the valley floor and Rylekaerene are highly variable between years and do not show an increasing trend from 915 the available data. Both study areas show the same pattern, but the areal distribution between classes differs, making the Rylekaerene study area almost entirely dependent on the variability of only the fen and fen fringe surface types.
In the larger valley floor study area, grassland, Salix snowbeds, heaths, and fell and barren surface classes dampen the variability between years.
Zackenberg Valley shows the potential for increased methane fluxes in the 21 st century, as methane shows a positive 920 correlation with temperatures , similar to the rest of Greenland and the Arctic (AMAP, 2017).

Methane flux upscaling in Arctic landscapes
Landscape-scale methane flux estimations are available from other subarctic and tundra sites (Table 34) in North America, Scandinavia, and Russia. The upscaled growing season fluxes range from 0.5 km 2 , extending from covering a few different ecosystems in a single site (Christensen et al., 2004) to large-scale estimates of fluxes 925 covering up to 320,000 km 2 for the Hudson Bay Lowlands (Roulet et al., 1994). (Christensen et al., 2004) to largescale estimates of fluxes covering up to 320,000 km 2 for the Hudson Bay Lowlands (Roulet et al., 1994). Mean methane fluxes in these landscapes range from 0.3 to 3 mg m −2 h −1 across different scales and landscape types,  Table 34. All the included studies are either fully or partly based on chamber measurements and upscaling with areal coverage of the surface classes, making them comparable to this study.
Several of the studies Roulet et al., 1994;Hartley et al., 2015)Several 935 studies Roulet et al., 1994;Hartley et al., 2015) estimated the landscape methane fluxes based on observations from a single year. The remaining studies combine observations from multiple years or studies to a flux estimate from each surface cover class (Christensen et al., 2004;. (Christensen et al., 2004;. Flux measurements from several years lead to more robust landscape 940 flux estimates, as the fluxes are highly variable between years, as the present study also shows.
The large differences in study area size and composition are ultimately determiningdetermine the mean methane flux estimates of the landscape, which makesmaking direct comparisons between sites difficult. For instance, the mean landscape flux found in this study is nearly fourfive times greater for the fen-rich Rylekaerene study area, which is entirelyfully contained in the valley floor study area. On an even larger scale, the entire northeast 945 Greenland acts as a net sink of methane, as  found a mean methane flux of ~−0.08 mg m −2 h −1 in their 10,675 km 2 study area.  Ecosystem fluxes from these studies range between ~0.1 to 6.2 mg m −2 h −1, with the lowest growing season fluxes measured at an upland tussock tundra site by Eight Mile Lake in Alaska (Taylor et al., 2018) and the highest fluxes measured in a mire in Stordalen, north Sweden (Jackowicz- . Mean growing season fluxes 960 found in Zackenberg using the EC method are within this interval . (Taylor et al., 2018) and the highest fluxes measured in a mire in Stordalen, north Sweden (Jackowicz- . Mean growing season fluxes found in Zackenberg using the EC method are within this range .

Landscape methane flux in a changing climate 965
A warming trend in both air and soil temperatures has been observed for Zackenberg Valley Christensen et al., 2020b). Christensen et al., 2020b). The increase in temperatures has contributed to the destabilization of permafrost, leading to several active periglacial landforms in recent years Cable et al., 2018;Christensen et al., 2020b). Cable et al., 2018;Christensen et al., 2020b). Modeling results show higher soil temperatures and a deepening of the AL in Zackenberg 970 in the future (Christiansen et al., 2008;. (Christiansen et al., 2008;. Increasing temperatures are expected to impact both positive and negative methane flux and surface erosion in the Arctic Schuur et al., 2015), (Schuur et al., 2015;Oh et al., 2020), which is also likely for Zackenberg Valley. The emergence of several active erosion sites along thein Zackenberg RiverValley in recent years could be an initial step towards increased erosion activity induring the 21 st century in 975 the valley.. In 2019, after the disappearance of methane-rich ice -wedges in the previous year, carbon-rich soils had been washed out from the gully, leaving a silt-organic mix with limited potential for methane emission in the area (Christensen et al., 2020b). This (Christensen et al., 2020b). Our study defines twothree erosion scenariospathways to illustrate the sensitivity of methane flux to land cover changes on a valley -scale. We hypothesize large-scale linear growth in eroded areas along the banks of the Zackenberg River by the late 21 st century (i.e., 2081-2100), 980 characterized by the changes in fluxes likein parts of the valley that are likely to have a high ground ice content and are located near streams and rivers identified by Cable et al. (2018). These erosion areas are assumed to share the characteristics of those observed in the gully area in 2019. Increased gully erosion wouldcould transform large areas with limitedsurfaces with both methane emission and uptake into well-drained, low emission eroded surfaces.
These two erosion scenariospathways are unlikely, but they are valuable as a sensitivity study illustrating the 985 difference in importance of increasing temperatures relative to eroding surfaces. In this example, a backward gully erosion area was studied, but erosion caused by the Zackenberg River, and especially GLOFs, appear to have a larger impact on the erosion of the riverbanks  and would most likely have a similar effect on methane fluxes. Further, erosion of the riverbanks is more pronounced on the outer banks of the river bends . Backward gully erosion may be limited to ice-rich parts of the riverbanks, but the extent 990 of these is largely unknown along the Zackenberg River inside the study area (Cable et al., 2018).
Even large-scale erosion alongon the rivervalley floor would have a limited impact on the mean valley methane flux, reducing the fluxit by less thanup to 1.2 % on average between 2041-2060, while surface erosion gradually progresses toward the first third of 25 m or 100 m, respectively. for the most extreme erosion pathway (Fig. 7, f). 995 The reduction becomes more pronounced (4.9 % reduction) between 2081-2100, as eroded areas would develop further inland and cause disturbances to areas dominated by fen. When fen areas are eroded, the flux is expectedcontinue to decline, causing the area erode at a faster pace. Disturbances of this magnitude are comparable to become a relative sink of methane. The erosion could cause the mean methane fluxes of the the size of the edge trimming during the GLOF in 2017 in the lower river section , located inside the valley floor 1000 study area. The impact on fluxes from one added gully, similar in size to decrease by −0.7 %the recent gully per year, is minimal (less than 1 % reduction) at the end of the 21 st century (Fig. 7, d) when compared to the uncertainty range. In our calculation of landscape methane fluxes after erosion, eroded areas become revegetated at a rate of 2 % y -1 . This rate may change over time and −7.8 % in the late period, with 25 mfrom one location to another, dependent on species composition and 100 msoil conditions in the eroded areas along the river. The change is 1005 assumed to be linear in this case but may accelerate over time, as plant communities established. In the calculation, revegetated areas were considered to have the same methane flux as the grasslands in the gully area. Limited data from the measurements in the gully area shows that the undisturbed vegetated areas take up slightly more methane than the disturbed vegetation patches in the gully, but the 10-day measurement period did not allow for more detailed measurements on temporal dynamics in the area. 1010 The resulting changes in flux following erosion, particularly for the less extreme pathways (Fig. 7, d and e), are minor relative to the considerable uncertainty in the general shift in methane for the valley following an increase in temperature. The mean methane flux is shown in Fig. 7 with a wide uncertainty range, partly caused by a lack of perfect fit between fluxes and temperatures and the uncertainties in the estimated landscape fluxes in Fig. 6. The 1015 variability between years in this study is substantial. The mean valley flux for 2016-2019 exceeds the confidence bounds, which shows that this should be seen as a sensitivity study of changes occurring on decadal scales.
In this case, methane fluxes were measured in a gully, which becomes more drained with a loss of organic soils from the surface after an erosion event (Christensen et al., 2020b). Fluvial erosion of vegetated riverbanks would most likely have a similar impact on methane fluxes when organic soils and vegetation are eroded. (Christensen et 1020(Christensen et al., 2020b.
The gully described in this study has likely formed as a direct consequence of pronounced lateral erosion from the river, which steepened the gully, allowing for increased drainage of water and sediments. The increased drainage exposed more and more of the ice -wedges and frozen soils, which later thawed and flowed out through the steeper gully. Lateral erosion of a similar scale has not been recorded to occur along the smaller rivers and streams in 1025 Zackenberg, but several smaller erosion sites have been described in recent years Cable et al., 2018). The gully shares characteristics with thermokarst gullies, a common type of thermokarst erosion, including extent, depth, and shape . . Abrupt thaw, including both gullies and thermokarst areas, can take different forms and affect the surface methane flux through disturbances in both vegetation and hydrology .  use a definition of thermokarst 1030 landscapes, including .  describe thermokarst landscapes, which include thermo-erosion gullies characterized by lateral movement of sediments, similar to the gully described in our study.   estimate that ~20 % of land areas in the northern permafrost zone are thermokarst landscapes, meaning they are either currently characterized by soil settlement or erosion or prone to developing into thermokarst landforms in the future. Thermokarst landforms have diverse impacts on a landscape, 1035 dependent on their type . They could form abruptly and rapidly responding to increasing temperatures (Farquharson et al., 2019;Lewkowicz and Way, 2019). Wickland et al. (2020) . They could form abruptly and rapidly responding to increasing temperatures (Farquharson et al., 2019;Lewkowicz and Way, 2019). Wickland et al. (2020) found a 42 % increase in growing season methane flux from 1949-2018 in a wet polygonal tundra after an increase in thermokarst erosion of ice -wedges. Thermokarst 1040 lakes cover large areas across the Boreal zone and in the Arctic and have been reported to be substantial emitters of methane Walter Anthony et al., 2018;. These Walter Anthony et al., 2018;. The findings from thermokarst lakes contrast the results from this study, where erosion causes draining and loss of organic material.   (Fig. 5).  found a mean flux, combining diffusion and ebullition, of ~3.7 mg m −2 h −1 from thermokarst lakes in the ice-free season. 1045 Walter Anthony et al. (2018)Walter Anthony et al. (2018) found an accelerating increase in methane fluxes from thermokarst lakes in their modeling of methane in the 21 st for the circumpolar region, which also contrastsshows that the impacts from on methane in the gully found in Zackenberg Valleyevent of erosion are diverse. It should be noted that there is a difference in the 'receiving end' of the gully formation between our study and these thermokarst lake studies as the drainage in the Zackenberg valley goes straight out ininto the major river systems with few or 1050 no stagnant reservoirreservoirs between.
In our calculation of landscape methane fluxes to changes in surface cover, eroded areas stay unvegetated for the remainder of the century. Still, eroded areas may partly be revegetated after stabilizing, which could change the size of the relative sink shown in Fig. 6. Further, methane oxidation could increase the sink in upland areas during the 21 st century Oh et al., 2020), which could offset some of the increase in mean flux in 1055 Zackenberg Valley.

Conclusions
In this study, we have summarized 14 measurement-years of methane fluxes and several short-term campaigns, which provide a unique insight into the large variability in methane fluxes in a high Arctic tundra landscape of This study shows the importance of multiyearmulti-year methane monitoring with wide spatial coverage, as interannual variability is substantial when considering a full composite landscape even in a single valley in the 1080 Arctic.  com/app/uploads/2017com/app/uploads/ /11/eBee_saves_day_mapping_greenlands_zackenberg_research_stat ion.pdf, 2015com/app/uploads/ . 1155