Influence of Tundra Polygon Type and Climate Variability on CO2 and CH4 Fluxes Near Utqiagvik, Alaska

Arctic tundra has the potential to generate significant climate feedbacks, but spatial complexity makes it difficult to quantify the impacts of climate on ecosystem‐atmosphere fluxes, particularly in polygonal tundra comprising wetter and drier polygon types on the scale of tens of meters. We measured CO2, CH4, and energy fluxes using eddy covariance for 7 yr (April to November, 2013–2019) in polygonal tundra near Utqiagvik, Alaska. This period saw the earliest snowmelt, latest snow accumulation, and hottest summer on record. To estimate fluxes by polygon type, we combined a polygon classification with a flux‐footprint model. Methane fluxes were highest in the summer months but were also large during freeze‐up and increased with the warming trend in August–November temperatures. While CO2 respiration had a consistent, exponential relationship with temperature, net ecosystem exchange was more variable among years. CO2 and CH4 exchange (June–September) ranged between −0.83 (Standard error [SE] = 0.03) and −1.32 (SE = 0.04) μmol m−2 s−1 and 13.92 (SE = 0.26)—23.42 (SE = 0.45) nmol m−2 s−1, respectively, and varied interannually (p ≤ 0.05). The maximum‐influence method effectively attributed fluxes to polygon types. Areas dominated by low‐centered polygons had higher CO2 fluxes except in 2016–2017. Methane fluxes were highest in low‐centered polygons 2013–2015 and in flat‐centered polygons in subsequent years, possibly due to increasing temperature and precipitation. Sensible and latent heat fluxes also varied significantly among polygon types. Accurate characterization of Arctic fluxes and their climate dependencies requires spatial disaggregation and long term observations.

The Arctic coastal plain of the North Slope of Alaska is characterized by ice wedge polygonal tundra, small ponds or lakes, and drained lake basins. For example, on the Barrow Peninsula near Utqiagvik (formerly Barrow), roughly half the area is covered by low-centered polygons and lakes (24% each), with other areas covered by drier flat-centered (17%) and high-centered (11%) polygons (Lara et al., 2015;Wainwright et al., 2015). Ice wedge polygon formation is related to cycles of freezing and thawing of the active layer (Brown, 1967;Liljedahl et al., 2016), creating units that extend from a few meters to tens of meters across. The different polygon types differ in soil microbial community (Graham et al., 2013;Jansson et al., 2012;Taş et al., 2018), soil organic carbon content (Bockheim et al., 1999;Dennis et al., 1978;Michaelson et al., 1996), hydrology and inundated area (Brown, 1967;Liljedahl et al., 2016;Zona et al., 2011), and vegetation cover (Dennis et al., 1978;Webber, 1978). These diverse properties of tundra landscape combine to spatially diverging CO 2 and CH 4 fluxes across the landscape as shown in many high latitude studies (e.g., Pirk et al., 2017 and references therein).
The Arctic is undergoing a warming trend (Box et al., 2019;Jeong et al., 2018;Moritz et al., 2002) causing later active-layer freeze-up (Box et al., 2019) as well as a deepening of the active layer (Andresen et al., 2020;Gewirtzman et al., 2018). Along the Arctic coastal plain of Alaska, active layer depth (ALD) varies at fine scale due to microtopography, vegetation, snow dynamics, and belowground properties, such as soil moisture for example, but a trend of deepening ALD has been observed at some locations. For example, at the Barrow Environmental Observatory, near Utqiagvik and environs where this study was based, ALD deepened between 2000 and 2015 as found by the recent modeling study by Yi et al. (2018), integrating remote sensing and field observations. This said, a declining trend (ALD shallower) was observed from manual annual probing at the site between 1995 and 2007 by Streletskiy et al. (2008), coinciding with a period of global warming hiatus (Dai & Wang, 2018) with an average maximum depth of 0.39 m (SE = 0.01) in 2013-2018 (D. Nicolsky, University of Alaska, Fairbanks [personal communication]; Romanovsky et al., 2017) usually reached in mid to end August (Shiklomanov et al., 2010;Streletskiy et al., 2008). Precipitation may also be trending, with recent studies showing an increase across Alaska (Wendler et al., 2017) and the greater pan-Arctic region (Box et al., 2019). In Utqiagvik, winter and spring precipitation are increasing (McAfee et al., 2013). Models project that climate changes such as these will alter the functioning and composition of polygonal tundra (Lara et al., 2018), including altering greenhouse gas and energy fluxes (Grant et al., 2019;Grant, Mekonnen, Riley, Wainwright, et al., 2017;Jeong et al., 2018;López-Blanco et al., 2017).
While net CO 2 , CH 4 , and energy exchange at the landscape scale are measured in the Arctic using the well-established eddy covariance (EC) technique (Aubinet et al., 2012;Baldocchi, 2003), partitioning Arctic CO 2 fluxes into ecosystem respiration and gross primary productivity (GPP) is less straightforward than in many other locations due to the lack of true night during the growing season.
One approach to measuring fluxes from small landscape features is the chamber method (Kutzbach et al., 2007;Norman et al., 1997;Pavelka et al., 2018). Chamber measurements of greenhouse gas fluxes in polygonal tundra show high variability in relative CO 2 flux rates among polygon types and saturation. For example, some studies report lower ecosystem respiration in low-centered polygons when the ground was wetter (Olivas et al., 2011) in comparison to drier periods (study only looked at low-centered polygons) while Pirk et al. (2017) measured higher respiration in wetter areas (no polygon type indicated) at their polygonal tundra site. Wilkman et al. (2018) reported consistently higher respiration rates in troughs and low centered polygons in comparison to high-centered polygons. However, Vaughn et al. (2016) found no significant difference among polygon types. These results illustrate the high variability in respiration in this spatially and temporally heterogeneous environment. In contrast, methane chamber measurements in tundra ecosystems appear more predictable, with consistently higher efflux in low-centered polygons and troughs, and efflux in flat-centered polygons only during wetter periods (Lara et al., 2015;Sachs et al., 2010;Vaughn et al., 2016;von Fischer et al., 2010;Wagner et al., 2003). While chamber measurements are useful for characterizing controls of fluxes by plant, moisture, and temperature conditions, their fine scale (<1 m 2 ) and intermittent deployment-continuous chamber measurements especially tend to be poorly replicated-makes scaling up over space and time challenging. In contrast, the EC method captures the turbulent exchange continuously and over a larger area up to 1 km 2 depending on atmospheric stability (Burba, 2001). The incommensurate nature of these measurements complicates scaling chamber fluxes to the 10.1029/2021JG006262 3 of 23 ecosystem level or using flux data to provide mechanistic explanations for observed EC fluxes (Fox et al., 2008;Zhang et al., 2012). Additional means of attributing EC fluxes to landscape features are needed.
Footprint studies (Kljun et al., 2004;Kormann & Meixner, 2001) help in estimating the field of view of EC measurements at larger scale, but have not been applied to quantify fluxes from smaller or distinct landscape features, such as different polygon types. The polygon classification work by Wainwright et al. (2015) opens up a new measurement avenue in which polygon distribution, footprint analysis, and the maximum-effect source location within the footprint (estimated for each half-hourly flux value) and its distance from the EC tower can be combined to attribute flux components to locations within the footprint or to specific landscape types.
Although EC measurements cover a spatially small area (footprint) of up to a few hundred meters with sites across the Arctic still limited in numbers, their data are nonetheless used in estimating regional fluxes of CO 2 (Ueyama et al., 2013), CH 4 (Sayres et al., 2017;Zona et al., 2016), or energy fluxes (Cristóbal et al., 2017). The spatial heterogeneity of the Arctic polygon tundra has large implications for regional fluxes and land surface models and the fidelity and uncertainty of their predictions. A synthesis of 40 terrestrial models' simulations of Alaska (Fisher et al., 2014) showed almost no consistent spatial patterns in net ecosystem exchange (NEE). Different models showed the entire state as a strong carbon sink, while others showed the opposite. Some predicted carbon neutral emissions and some models showed a spatial distribution of sinks and sources, as also did a newer study by Commane et al. (2017).
A model study on energy fluxes of the North Slope of Alaska (Lynch, Bonan, et al., 1999; detailed the difficulty of model validation and performance using measured data from multiple meteorological and/or EC sites and research groups, highlighting the uncertainties attached to measurement methods, instruments, and model parametrization. Uncertainties that accompany such measurements vary from EC tower/ met station location to the exact placement and depth of the individual soil temperature probe or soil heat-flux plates and their respective spatial representativeness. These observational challenges and model studies highlight the difficulties in understanding and predicting the complex and heterogeneous Arctic polygon tundra, and its dynamics under a changing climate. To better understand the temporal and spatial heterogeneity of CO 2 , CH 4 , and energy exchange at the Next Generation Ecosystem Experiments (NGEE Arctic) Barrow research site near Utqiagvik, Alaska, we analyzed seven consecutive years of EC and meteorological measurements covering April to November, 2013-2019. This period represents some meteorological extremes, including the earliest snowmelt dates, latest snow accumulation dates (Cox et al., 2017), and hottest summer on record (National Oceanic and Atmospheric Administration [NOAA], 2019). The objectives of this study were to (a) document the inter-and intra-seasonal variation in CO 2 , CH 4 , and energy fluxes at a polygonal tundra site; (b) explore how variability and trends in meteorological variables affect these fluxes; and (c) investigate the spatial variability of these fluxes, to ultimately (d) assess the possibility of separating these landscape-scale fluxes by polygon type.

Site Description
This study was carried out in the Barrow Environmental Observatory (BEO, 71°16′48.26″N; 156°36′33.12″W) on the coastal plain near Utqiagvik (Barrow), Alaska, as part of the United States Department of Energy Next Generation Ecosystem Experiments (NGEE-Arctic) project ( Figure 1). The mean annual temperature is −12°C and mean annual precipitation is 114.3 mm. The site is characterized by polygonal structures classified into low-, flat-, and high-centered polygons, each representing a different microtopography and hence surface, soil, and hydrological properties (Brown, 1967;Liljedahl et al., 2016;Wainwright et al., 2015;Zona et al., 2011). These features are identifiable in aerial images (Cherry & Crowder, 2016;Svensson, 1963) and in high-resolution elevation models (Wilson & Altmann, 2017). According to Lara et al. (2015) and Wainwright et al. (2015), low-centered polygons occupy the largest area across the Barrow Peninsula and the BEO, followed by flat-centered ones, and high-centered polygons only cover around 11% of the tundra. Areas north of the EC tower site (AmeriFlux ID US-NGB, Figures 1 and 2) are dominated by low-centered polygons while the southern and western sections by flat-and high-centered polygons. The main vegetation types found across the BEO are Carex aquatilis subsp. Stans, Saxifraga cernua communities, Luzula confusa, Dupontia fisheri as well as non-acidic tundra Carex spp. and Eriophorum spp. (

EC Measurements
We measured CO 2 , CH 4 and energy exchange at the landscape level at the NGEE-Arctic eddy flux site  using EC for 7 yr (April to November, 2013-2019) with a 74% overall April-November data coverage (30 min data) and 79% for the June-September data range. The system consisted of a Gill R3-50 ultrasonic anemometer (Gill Instruments, Lymington, UK), and an open path LI-7500A CO 2 / H 2 O gas analyzer (LI-COR Biosciences, Lincoln, NE, USA), along with an LI-7700 CH 4 gas analyzer (LI-COR Biosciences) at 4 m above ground. Digital raw data streams from all three of these instruments were synchronized and recorded by a program (HuskerFlux) developed by D. Billesbach; running on a small, single-board computer, which also provided communications for data retrieval, and served as a platform for the Campbell Scientific LoggerNet program which archived raw data from other (slow response) instruments. The following processing steps were taken in calculating 30-min fluxes from the raw 10 Hz data: spike detection was implemented using the standard Vickers and Mahrt (1997) de-spiking algorithm. We applied the coordinate rotation using the second rotation method described in Baldocchi et al. (1988), sonic virtual temperature correction (Kaimal & Gaynor, 1991), frequency response correction (Massman, 2000), and the buoyancy flux correction defined in Schotanus et al. (1983). The computed co-spectra for vertical wind and temperature, CO 2 , CH 4 , and H 2 O, respectively, for each half hour using the Fast Fourier Transform were used to verify the magnitude of the frequency response correction, reflecting the ability of the system to measure high-frequency gas transport. We applied the Webb-Pearman-Leuning density correction (WPL; Webb et al., 1980) to the data originating from both gas analyzers, as well as the LI-7700-specific corrections for spectroscopic effects (LI-COR, 2010;Mc-Dermitt et al., 2011). Furthermore, data quality control measures were carried out following guidelines given in Foken et al. (2004) by removing flux values with quality flags of seven and higher, in addition to excluding flux data recorded under low signal strength and rainy periods as rain has an impact on the performance of open-path gas analyzers as well as sonic anemometers (Burba, 2013). These quality control procedures, together with fluxes filtered for low friction velocity led to a removal of 38% of measured CO 2 fluxes (27% for June-September time periods). In addition to the friction-velocity filter, one of the common reasons for data removal was low signal strength.
Due to its proximity to the Chukchi Sea, the site experiences prolonged periods when the foggy air is enriched with salty sea air leading to low signal strength values, sometimes within hours of cleaning the optical windows of the open path analyzers. Furthermore, the site experiences extensive periods without power.

Meteorological Measurements
Meteorological measurements were made along with the EC measurements (April to November, 2013-2019). The suite of instrumentation included air temperature and relative humidity (HMP60, Campbell Scientific, Logan, UT, USA) measured at 3 m above ground (overall data coverage 80%, 93% for June-September time periods, and similar coverage for all other met variables), incoming and outgoing short-and longwave radiation (CNR4, Kipp & Zonen, Delft, The Netherlands), and incoming and reflected photosynthetic photon flux density (PPFD) that was measured with two LI-190 quantum sensors (LI-COR Biosciences) at 2 m above ground. Surface temperature was measured using a SI-111 Infrared Radiometer (Apogee Instruments, Logan, UT, USA) with a 22° half-angle field of view attached at 2 m above the ground. These IR sensors use a portion of the longwave thermal radiation emitted by the surface (i.e., tundra) within the field of view. Soil heat flux (G) was measured with four soil heat flux plates (Hukseflux, Delft, The Netherlands) at 7 cm below the surface in the vicinity of the EC tower. No attempt was made to estimate the heat flux above the soil heat flux plates or at surface level. These meteorological data  were quality controlled following methods briefly introduced in Christianson et al. (2017) and Pastorello et al. (2017) by removing observations outside realistic margins, as well as correcting for timestamp drifts, instrument drift, or use of wrong coefficients, for example, soil temperatures at 10 cm depth were measured at five locations along a transect across polygons near the EC tower by Romanovsky et al. (2017). Daily thaw depth values were estimated from soil temperature readings recorded along a vertical profile in the center of a polygon to a depth of 1.5 m .

Flux Partitioning, Gap Filling, Calculation of Cumulative Fluxes, and Energy Balance
Due to the late snowmelt season (mostly June) and early autumn snow accumulation that is often already occurring in September (Figure 3, see incoming and reflected radiation) only June-September data were partitioned into ecosystem respiration and GPP. To do so ecosystem respiration was calculated from relationships between nocturnal NEE and soil temperature measured within the footprint area. Windy (friction velocity > 0.15 m s −1 ) nighttime flux data (PPFD < 50 μmol m −2 s −1 , as also used by Euskirchen et al., 2012;Nordstroem et al., 2001 at their Arctic sites) were plotted against soil temperature and an exponential function (Lloyd & Taylor, 1994) fitted. Currently there is no consent on use of soil versus air temperature in Arctic ecosystem respiration studies (Euskirchen et al., 2012(Euskirchen et al., , 2014Grøndahl & Ministry of the Environment-Denmark, 2006;Kutzbach, 2006;Shaver et al., 2013;Vourlitis et al., 2000;Zamolodchikov et al., 2003) with air temperature generally being used due to its more frequent availability, as stated by Keenan et al. (2019). Ecosystem respiration represents the sum of heterotrophic (soil) respiration that is driven by soil temperature and moisture, and autotrophic respiration that is driven by air temperature (Ruimy et al., 1995) so a combination of both appears most appropriate. Wohlfahrt and Galvagno (2017) pointed out that the use of soil temperature as the driving temperature underestimates daytime respiration while the use of air temperature as the driving temperature to overestimate daytime respiration. In order to obviate such an under-or overestimation of daytime respiration we plotted nighttime fluxes against soil temperature and used air temperature to estimate daytime respiration and GPP. We believe this to be the appropriate method at this site as this method avoids this bias together with the additional uncertainty introduced by thermal inertia/thermal offset and the resulting time lag of Arctic soil temperature in relation to air temperature (Aalto et al., 2018;Biskaborn et al., 2019;Lyu, 2015;Lyu & Zhuang, 2018;Zhang et al., 2005) and its effect on respiration (Shaver et al., 2013).
Following Causton and Dale (1990) and Ruimy et al. (1995), the Michaelis-Menten rectangular hyperbola was used to describe the relationship between photosynthesis (GPP) and PPFD, representing the apparent photochemical efficiency (light use curve): where Q g is the incident PPFD, A max the maximum gross photosynthesis at infinite irradiance, and b is the PPFD value at which half of the light-saturated photosynthesis is occurring. A max /b is the slope of the curve at A max = 0, representing the apparent photochemical efficiency (light use efficiency [LUE]; Causton & Dale, 1990).
In order to calculate the net carbon (g C m −2 ) cumulative uptake of the site, gap-free data are required. Following the same procedure as with partitioning NEE, we only gap-filled NEE for June-September by applying the Amer-iFlux ONEFlux Marginal Distribution Sampling (Pastorello et al., 2020) gap-filling procedure incorporating incoming shortwave radiation, air temperature, and vapor pressure deficit which act as drivers of net fluxes.
Missing meteorological values have implications for gap-filling methods/procedures and associated uncertainties (see Reichstein et al., 2005). The site had no power in September 2016 and to avoid introducing very high uncertainties no attempt was made to gap-fill missing September 2016 values and so this month was not included in the cumulative sums. Gap-filled averages and associated standard errors for each June-September (2013-2019) period are reported in Table 1. A detailed description of the method, procedures, and uncertainties can be found in Reichstein et al. (2005; Appendix A) and Pastorello et al. (2020). CH 4 flux data were not gap-filled.
Energy balance was calculated as the fraction of the total sum of latent, sensible and soil heat flux, and the net incoming and outgoing radiation measured at the site. While latent and sensible heat were measured with the EC suite of instrumentation, net radiation was measured off the EC tower and ground heat flux plates installed near the tower. This common measurement method is well known to be accompanied with uncertainty introduced by the combination of a varying footprint (potentially each half hour) and fixed point measurements.

Polygon Classification and Footprint Determination: Maximum Influence Estimation
Since source and sink locations are dependent on the underlying surface, in our case polygonal features, we include the polygon type classification introduced in Wainwright et al. (2015) in conjunction with a simple application of the Kljun et al. (2004) Lagrangian footprint model.
The ability of this model to predict the area of influence (footprint/field of view of the EC system) and the maximum-effect source location (location with the highest contribution) within the footprint (estimated for each half-hourly flux value) and its distance from the EC tower are key to our analysis. The spatial extent of an individual footprint and distance to the maximum contribution location depends on measurement height, surface roughness, atmospheric stability, wind direction, and friction velocity. Due to the fine-scale heterogeneity of the site no attempt was made in classifying each individual footprint in the current study. Note. In addition, number of datapoints for each month are added in brackets as well as the standard error of the means (SE), as a measure of uncertainty, are included for all fluxes (June-September averages).

) and Sensible Heat (H-W m −2 ), Together With Monthly and June-Septemeber Total Precipitation (Prec-mm) and the June-September Averaged Gap-Filled (FC gf ) CO 2 Fluxes
Furthermore, footprint models assume stationarity of the flux across space and time. In order to combine these spatial and surface characteristics we overlaid the polygon type (LCP, FCP, and HCP) across coordinates on a map, allocating coordinates to each maximum flux contribution location derived from wind direction and distance allowing an attribution of fluxes to polygon type ( Figure 2) while acknowledging the potential uncertainty introduced due to the conversion of microtopography to a plane surface, a common feature at complex sites.

Statistics
We applied the bootstrapping approach (Banjanovic & Osborne, 2016) by resampling the flux data 10,000 times (with replacement) in order to compare the true/population means by looking at their confidence intervals together with the difference in means across all June-September periods as well as by polygon type and the degree of significance of differences between individual comparison pairs. The Mann-Kendall (Kendall, 1975;Mann, 1945) and the modified Mann-Kendall trend analysis (Hamed & Rao, 1998) commonly applied in hydrometeorological time series (Wang et al., 2020 and reference therein) but also EC time series  was applied to the data in order to investigate existing trends in fluxes or meteorological conditions over time and the statistical significance thereof. Prior to applying the modified Mann-Kendall trend analysis we assessed the data for existing autocorrelation and partial autocorrelation, a procedure common in time series analysis, as autocorrelation is known to increase the chances of detecting significant trends where in fact there are none. Results summarize the Kendall rank correlation coefficient (τ) indicating the trend strength (−1 < τ < 1) as well as the magnitude (S) of the trend (Sen's slope [Sen, 1968], change in flux) together with the significance at the 95% level (p ≤ 0.05) of the observed trend. These analyses were applied to daily averages treated separately by the three time periods: snowmelt (April and May), June-Septemeber, and freeze-up (October and November) as well as separated by polygon type.

Temporal Variability in Net CO 2 Exchange, GPP, and Respiration Over 7 yr
NEE averaged over April to November (Table S1) Table 1). There were significant differences in fluxes between some years (bootstrapped, 95% confidence level using the percentile method; Figure 4), that are partially explained by variation in climate over these seven years (see Section 3.4). For example, we can say with 95% confidence that 2013 fluxes were somewhere between −0.30 up to −0.60 μmol m −2 s −1 higher (sign convention, more uptake) than 2014, but 2013 was not different from 2015, having been 0.17 μmol m −2 s −1 lower than 2015 up to −0.14 μmol m −2 s −1 higher than 2015 (see confidence interval values in Table S2 for all CO 2 [by year] comparisons).
The tundra was a net sink of CO 2 in each of the June-September periods (Figure 5) with 65.4-102.1 g C m −2 cumulative net uptake. Late summer/early autumn uptake seen in Figure 5 can possibly be explained by warm temperature spells promoting continuation in biomass growth and carbon uptake (see September averages for air temperature and CO 2 fluxes in Table 1). No clear trend with length of the snow-free season was found.
There was a weak but significant decreasing trend in April-November NEE over the 7 yr (becoming more positive) with τ = +0.066, p = 0.080 with an average daily change rate of S = 0.0003 μmol m −2 s −1 but not for the June-September values.
The weak decreasing trends for April and May (snowmelt) and Octoberand November (freeze-up/post photosynthetic season) periods were not significant (Table S3).
GPP, based on partitioning NEE, was highest in 2013 and lowest in 2014. To investigate ecosystem efficiency in converting light into photosynthate, we estimated the LUE as the relationship between PPFD and GPP. The response is asymptotic (fitted rectangular hyperbola, Figure 6), with parameters depending on leaf photosynthetic properties and the environment: In our case it describes gross photosynthesis, and so the curve passes through the origin. Parameters estimated from the Michaelis-Menten rectangular hyperbola (A max , b, and LUE) for each month separately are summarized in Table S4. June-September GPP and LUE by polygon type were not significantly different from each other.  No relationship with amount of precipitation nor a lag between respiration and precipitation was found. As with NEE, no significant trend was observed in GPP nor respiration over the seven June-September periods.

Methane Fluxes
CH 4 fluxes averaged over each April to November (Table S1) (Table 1) (Table 1). Averaged over all years, June-September CH 4 flux was 18.44 (SE = 0.14) nmol m −2 s −1 . Figure 8 shows the cold-season (April and May/spring snowmelt and October and November/fall freezeup) distribution of CH 4 fluxes outside the typically reported snow-free June-September period as well as the whole April-November data range.
According to Mann-Kendall analysis, there appears to have been an increasing trend in June-September methane fluxes over these seven years (Mann-Kendall analysis, τ = +0.111, p = 0.01) with an average change rate of 5.47 nmol m −2 s −1 per year. Spring (April andMay) fluxes were very low, averaging 1.1 (SE = 0.1) nmol m −2 s −1 with no significant trend across years. In contrast, fall (October and November) CH 4 fluxes, which averaged 13.93 (SE = 0.34) nmol m −2 s −1 , had a significant increasing trend with an average change in rate of 0.125 nmol m −2 s −1 day −1 (p = 0.029; Table S2).
As with CO 2 fluxes, the June-September CH 4 flux averages varied significantly between some but not all years (Figure 4). Methane fluxes in 2013 were higher than 2014, somewhere between 4.75 and 6.68 nmol m −2 s −1 (95%   Table S2).
We saw a positive relationship between CH 4 efflux and surface temperature (Figure 9, r 2 = 0.23, p ≤ 0.0001), as well as very high daily average CH 4 fluxes during periods of frequent thawing and refreezing of the surface in July and August. No relationship with amount of precipitation was found.
Furthermore, CH 4 fluxes correlated with a deepening of the active layer (Figure 10) during the summer month (maximum thaw depth mid-end August) followed by the "zero curtain" (Hinkel et al., 2001) period when the upper surface layer starts to refreeze (data not shown here) but with still substantial CH 4 production in September and the freeze-up season.

Climate Variability 2013-2019
There was significant climate variability over the study period in snow-free season dates, rainfall, and temperature. The length of the snow-free season varied from 76 days in 2014 to 159 days in 2016. The variation was due to   (Table 1), while 2019 had the highest precipitation (176.7 mm, data extracted from the NOAA Barrow Atmospheric Baseline Observatory according to Diamond et al., 2013). In the June-September period each year June was the coldest and July was the warmest month. Although there was high variability among months (see Table 1), the trend across these seven years is that June-September air temperature was slightly increasing (Table S3, τ = 0.08, p = 0.09). The 2014 had the coldest and 2019 by far the warmest June-September temperatures (Table 1)  Incoming spring (April and May) radiation (PPDF) showed a very strong increasing trend over the years (τ = 0.305, p ≤ 0.0001, S = 0.578; 2017-2019 were sunnier then previous years and hence warmer spring temperatures), but no significant trend over the October and November fall months. Regarding albedo, no changes were observed over the years.
Spring temperatures increased more than did June-September temperatures (τ = 0.104, p = 0.045, S = 0.01°C) while the observed increase in October and November temperatures (Table S3) were not significant.

Attributing Variation in EC Fluxes to Different Polygon Types
To attribute fluxes by polygon type, we mapped the maximum flux-contribution location obtained from our footprint analysis (Figure 2) for April to November, 2013-2019, overlaid on a polygon classification map (see map in Figure 2).
Fluxes varied among polygon types ( Figure 11) and the relative flux rates among types changed over time.
Measured NEE and CH 4 fluxes from low-and flat-centered polygons can be compared across all June-September periods in Table S5. Within a polygon type, there was no significant trend in spring nor June-September period NEE (Table S3) while during fall months NEE from FCP showed a weak decrease (sign convention for CO 2 fluxes, more negative = more uptake, therefore a positive trend means fluxes are becoming more positive; p = 0.11), while the τ value of 0.114 for HCP (p = 0. 19) hints at a strong decrease in cold season NEE (i.e., less CO 2 uptake). For methane, both polygon types showed a positive trend in June-September CH 4 fluxes (Table S3). There was no detectable trend in spring CH 4 fluxes. Fall CH 4 fluxes, however, exhibited a significant increasing trend in LCPs and marginally significant increasing trend in FCPs. Based on the bootstrapping approach, fluxes across LCP and FCP were not significantly different (at the 95% confidence level) for all June-September periods ( Figure 11). Mean CO 2 fluxes in low centered polygons were higher in 2013 and 2015, while in 2017 low centered polygons CO 2 flux means were lower somewhere between 0.12 and 0.42 μmol m −2 s −1 (Table S6). In contrast, CH 4 fluxes appeared to be significantly different in four out of 7 yr ( Figure 11) when comparing low centered with flat centered polygons.

Energy Fluxes and Energy Balance by Polygon Type
Energy fluxes were similar across the different polygon types in some years and different in others years (Figure 11 and Table S6 (Table 1).
Considering all June-September data (Figure 12), energy balance ((LE-H + G)/net radiation) averaged 70% with large fluctuations among polygon types, years, and months (range 35%-94%). Such energy balance estimates at EC sites include uncertainties introduced by the lack of representativeness of net radiation and soil heat flux measurements (fixed point measurements) in the wider area near the eddy flux tower (change in footprint, regarding LE and H), as well as lateral energy losses via subsurface water flows.
No significant difference in energy balance closure was found amongst the different polygons in any year. Across years, however, energy balance closure sharply increased for the spring months (τ = 0.24, p ≤ 0.0001) but no significant trend in the June-September nor the fall/freeze-up period was observed. One uncertainty we found in energy balance measurements at our site were the behavior of ground heat fluxes over the years. The two continuously running soil heat flux plates (both in LCP) showed no trend in the spring but opposing trends during June-September (τ = −0.09, p = 0.06 and τ = 0.101, p = 0.02) and clear increasing trends in the fall months (τ = 0.193 and 0.109) with only one being weakly significant (p = 0.07).

Discussion
Fluxes of CO 2 , CH 4 , and energy over 7 yr (2013-2019), including early spring/snowmelt (April and May), June-September and the post photosynthetic period/freeze-up season (October and November) from a polygon tundra, EC flux site near Utqiagvik, Alaska are discussed here. An overview of the meteorological conditions and flux distribution shows inter-and intra-seasonal variation but also notable spatial variability in partitioned flux values.

CO 2 Fluxes
The timing of snowmelt and spring temperatures had a large impact on tundra carbon sequestration. While the cumulative sink by the beginning of August 2018 was only 50 g C m −2 ( Figure 5), in 2015 (snow-free since May and with the warmest spring temperatures of the 7 yr [ Table 1]), approximately twice as much carbon was sequestered. This said, 2016, the year with the longest snow-free season but lowest precipitation values had intermediate carbon uptake (see Figure 5) Our findings show no trend in CO 2 uptake linked to snow-free season length, an observation also made in High Arctic Siberia by Parmentier et al. (2011) from 8 consecutive years of CO 2 EC flux measurements.
The tundra was a net CO 2 sink in every of the seven June-September periods, driven by summer uptake ( Figure 5) at rates similar to those measured by Euskirchen et al. (2012) for Imnavait Creek in the northern foothills of the Brooks Range, Alaska, but slightly lower than 119 g m −2 June-August average in the Lena Delta, Siberia (Kutzbach et al., 2007). Net CO 2 uptake commenced shortly after snowmelt, when air temperatures are still relatively low but insolation levels are highest, thus favoring GPP over soil respiration. The strong influence of PPFD (Figure 6) on Arctic GPP has been documented with EC in high-Arctic polygon tundra in Siberia (Dennis et al., 1978;Kutzbach et al., 2007) and at leaf level at our site (Rogers et al., 2019). Ecosystem respiration increased exponentially with surface temperature (Figure 7), as seen in other studies at the BEO (Harazono et al., 2003;Oechel et al., 1995) and Lena Delta (Kutzbach et al., 2007).

Changes in CH 4 Fluxes Over Time
June-September CH 4 fluxes showed an increasing trend (Mann-Kendall analysis) between 2013 and 2019, in concert with a slight increase in June-September air temperature (Table S3) and a potential deepening of the active layer (measurements near the site are done in up to 5 cm increments up to a depth of 0.4 m and 10-20 cm increments from 0.4 m downward, see Figure 10). The strong link with air and/or soil surface temperature (Figure 9) was also found in other Arctic studies (Arndt et al., 2019;Friborg et al., 2000;Sachs et al., 2010;Wille et al., 2008;Zona, Oechel, Kochendorfer, et al., 2009). Generally Arctic tundra soil temperatures at deeper depths vary too little (thermal inertia) however, to explain the high variability in CH 4 fluxes. Zona et al. (2016) and Taylor et al. (2018) estimated that around half of annual Arctic CH 4 emissions occur outside the snow-free season, with substantial contributions during freeze-up (Mastepanov et al., 2008;Pirk et al., 2015;Tagesson et al., 2012) and spring thaw (Raz-Yaseef et al., 2017;Wille et al., 2008) originating from individual pulse events, a phenomenon also observed in boreal permafrost regions (Hargreaves et al., 2001;Song et al., 2012). We observed large CH 4 fluxes during freeze-up in all years (Figures 3, 8, and 9), and substantial spring-thaw pulses only in 2014 (Raz-Yaseef et al., 2017).

Spatial Heterogeneity of Tundra Fluxes
Inundated areas, ponds, and lakes in high-latitude tundra can act as hotspots for both CO 2 and CH 4 fluxes (Jammet et al., 2015(Jammet et al., , 2017Langer et al., 2015), contributing to the spatial heterogeneity of high-latitude greenhouse gas fluxes. In addition, studies of vegetated drained thaw lake basins (DTLBs) near Utqiagvik found that the age of the basins (i.e., time since lake drainage) affects flux rates too (Lara et al., 2015;Sturtevant & Oechel, 2013), with higher NEE and GPPs in young DTLBs (Lara et al., 2015;Zona, Oechel, Peterson, et al., 2009). Changes in land cover type, age, permafrost thaw, and wetter or drier areas have a potential to impact future polygon tundra properties and fluxes in a warming climate (Lara et al., 2015(Lara et al., , 2018.

Impact of Polygonal Features on Tundra Greenhouse Gases
What impact do polygonal features have on tundra greenhouse gases? NEE fluxes from the two dominant polygon types, low-and flat-centered polygons, were significantly different during only three years ( Figure 11). Regarding CH 4 fluxes, in the first 3 yr (2013-2015) moister LCPs had higher fluxes (Table S5), followed by 4 yr (2016-2019) in which FCPs had higher values, with the increase in FCP CH 4 fluxes partially explained by an increase in seasonal air temperature (p < 0.001) and variation in precipitation (p = 0.067) dynamics, both promoting a deepening of the active layer (Douglas et al., 2020).
Polygon type alone did not explain the differences in fluxes. In addition, snow cover distribution is not uniform in these ecosystems (Gouttevin et al., 2018;Kępski et al., 2017;Wainwright et al., 2017), with snow depth during the first half of the June-September period potentially being more persistent in low-centered polygons than in flat-or high-centered polygons (Zona, Oechel, Kochendorfer, et al., 2009).
The result could be that due to wind and snow drift insulating the soil and exposing some flux drivers much later in the season than in drier high-and flat-centered polygons, thus contributing to differences in fluxes among polygon types. This also applies to spatially changing soil temperatures in the upper most layers of the tundra soil during snow melt, as well as the "zero-curtain" (Hinkel et al., 2001) effect later in the season.

Gas and Energy Fluxes at Polygon Level
This is the first study to attribute eddy fluxes to different polygon types in the footprint. Chamber data can provide an independent check on the relative patterns we found. Studies near our site (Lara et al., 2015;Sturtevant & Oechel, 2013;Vaughn et al., 2016;von Fischer et al., 2010;Wainwright et al., 2015) and in Siberian polygon tundra (Sachs et al., 2010;Wagner et al., 2003) report higher CH 4 fluxes from low-centered polygons than from drier areas such as high-and flat-centered polygons. These studies also found high fluxes from wetter features within the polygon landscape, such as depressions and troughs compared to the drier rims, but these fine-scale variations cannot be resolved by EC data. Studies have also found high CH 4 fluxes from the normally drier FCPs during wetter periods, consistent with our observations of the impact of an increase in precipitation.
The Arctic tundra not only shows a high spatial heterogeneity in greenhouse gas flux but also in surface energy fluxes (Boike et al., 2008;Langer et al., 2011;Lynch, Bonan et al., 1999;Muster et al., 2012;Ohmura, 1981Ohmura, , 1982, an interplay between climate drivers (e.g., available radiation/meteorological conditions) and ecosystem properties (soil, hydrology or wetness, and vegetation cover; Ling & Zhang, 2004;Sedighi et al., 2018). Our estimation of the site energy balance was problematic, as it is for nearly all EC systems, due to a mismatch in spatial scale in different measurements. While outgoing radiation and soil heat flux were estimated as fixed point measurements, turbulent fluxes were averaged over larger areas with a changing footprint, and influenced by microtopographical features, dry versus wet areas, and lateral losses of energy. Nevertheless, our half-hourly energy balance values (mostly above 60%, range 35%-94%) were within the ranges found in other energy balance studies in the High Arctic (Lynch, Bonan, et al., 1999 and authors therein;Grachev et al., 2019;Langer et al., 2011;Lund et al., 2014).

Heterogeneity of Polygon Tundra and Land-Surface Models
The heterogeneity of polygon tundra poses a challenge to land-surface models, predictions, and model fidelity. For example, a comparison of 40 models (Fisher et al., 2014) as well as a newer comparison study of 11 models by Commane et al. (2017) showed that the spatial patterns in mean NEE for Alaska varied widely among models, essentially showing no consistent distribution and some models completely contradicting the outcome from other models. A model study of energy fluxes limited to the North Slope of Alaska (Lynch, Bonan et al., 1999; detailed the difficulty of model validation and performance using data measured at multiple sites across similar terrain by multiple research groups, highlighting the uncertainties attached to measurement methods and instruments used. Model input affected model outputs by under/over estimating individual land cover features (Wainwright et al., 2015), flux components (Grant, Mekonnen, Riley, Wainwright, et al., 2017), and fluxes (Grant et al., 2019;Grant, Mekonnen, Riley, Arora, & Torn, 2017), ultimately illustrating the complexity of small-scale variations of the greenhouse gas and energy balance in ecosystem processes at regional scales. These numerous studies and models highlight the lack of complete understanding of the complex and heterogeneous Arctic polygon tundra, its dynamics, and future development under a changing climate.

Conclusion
Seven years of CO 2 , CH 4 , and energy fluxes (April to November, 2013-2019) revealed an increasing trend in CH 4 fluxes from Arctic tundra, including increasing fall (freeze up) fluxes as well as a consistent CO 2 sink during June to September months in each year, at the Barrow Ecological Observatory near Utqiagvik Alaska. There was high temporal and spatial variability in NEE, CH 4 , LE, and H fluxes linked to climate variability and site heterogeneity. To disentangle the heterogeneity of the complex polygonal tundra, we combined a polygon type classification (Wainwright et al., 2015) and footprint analysis that allowed us to separate the maximum flux-contribution location by polygon type. It revealed which polygon type had higher NEE or CH 4 fluxes that varied by time period investigated and year. Multi-year, multi-season data were necessary to quantify the long term trends and the relationship between spatial patterns (e.g., among polygon types) and climate variability. We hope this approach will pave the way for further such studies, extending understanding of the heterogeneity of Arctic ecosystems and the use of flux data in land model development, parameterization and validation.

Data Availability Statement
All eddy covariance and meteorological data included in this study are freely obtainable from the AmeriFlux database (https://ameriflux.lbl.gov/data/data-availability/) and additional auxiliary data from the US Department of Energy Next-Generation Ecosystem Experiments (NGEE Arctic) Data Portal (https://ngee-arctic.ornl.gov/data/).