Environmental controls on ecosystem-scale cold season methane and carbon dioxide ﬂuxes in an Arctic tundra ecosystem.

. Understanding the processes that inﬂuence and control carbon cycling in Arctic tundra ecosystems is essential for making accurate predictions about what role these ecosystems will play in potential future climate change scenarios. Particularly, air–surface ﬂuxes of methane and carbon dioxide are of interest as recent observations suggest that the vast stores of soil carbon found in the Arctic tundra are becoming more available to release to the atmosphere in the form of these greenhouse gases. Further, harsh wintertime conditions and complex logistics have limited the number of year-round and cold 5 season studies and hence too our understanding of carbon cycle processes during these periods. We present here a two-year micrometeorological data set of methane and carbon dioxide ﬂuxes, along with supporting soil pore gas proﬁles, that provide near-continuous data throughout the active summer and cold winter seasons. Net emission of methane and carbon dioxide in one of the study years totalled 3.7 and 89 g C m − 2 a − 1 respectively, with cold season methane emission representing 54% of the annual total. In the other year, net emission totals of methane and carbon dioxide were 4.9 and 485 g C m − 2 a − 1 10 respectively, with cold season methane emission here representing 82% of the annual total – a larger proportion than has been previously reported in the Arctic tundra. Regression tree analysis suggests that, due to relatively warmer air temperatures and deeper snow depths, deeper soil horizons – where most microbial methanogenic activity takes place – remained warm enough to maintain efﬁcient methane production whilst surface soil temperatures were simultaneously cold enough to limit microbial methanotrophic activity. These results provide


Introduction
Active-layer soils and permafrost soils in the Arctic permafrost region contain significant stores of terrestrial organic carbon.
These ecosystems account for an estimated 1307 (1140-1476) Pg of organic carbon, with ∼1035 Pg found within soils between 20 0 and 3 m depth (Hugelius et al., 2014). Recently-observed increases in surface air temperature within these regions (Polyakov et al., 2002) have sparked interest in the biogeochemical cycling of this carbon store, as substrate metabolic activity -shown to be positively correlated to temperature -can break down organic compounds in the soil, releasing soil organic carbon to 1 the atmosphere in the form of carbon dioxide and methane (Lai, 2009). Furthermore, the "active layer" horizon, within which most soil carbon decomposition takes place, has been observed in places to be expanding as the underlying permafrost thaws 25 under the influence of a warming atmosphere, thus exposing larger quantities of organic carbon to decomposition (Schuur et al., 2009;Romanovsky et al., 2010;Schuur et al., 2015;Vonk and Gustafsson, 2013).
Production of methane in the carbon-rich soils of the Arctic tundra takes place as a result of microbial metabolic activity (Lai, 2009). The break-down of organic carbon to form methane is a complex process that requires contributions from various taxa of microorganisms (Whalen, 2005). Methanogens form the last step in this process of producing methane from organic 30 polymers (Conrad, 1999); whilst these methanogens encounter oxygen in the breaking down of acetate and carbon dioxide, they cannot survive in oxic environments (Whalen, 2005;Kamal and Varma, 2008;Lai, 2009). As such, methanogenesis in peatlands is an obligate anaerobic process that takes place largely within deeper, anoxic layers of the soil, generally below the water table level (Le Mer and Roger, 2001). Methane production at these depths creates concentration gradients that lead to upward diffusion of methane through the soil to the surface (Preuss et al., 2013). Furthermore, methane can be transported to 35 the surface via ebullition and through aerenchymatous tissues within some vascular plants (Joabsson et al., 1999;Lai, 2009).
In the upper soil layers, another subset of microorganisms known as methanotrophs consume a portion of this produced methane in the presence of oxygen, eventually oxidising it to carbon dioxide (Lai, 2009). These methanotrophs are largely exposed to methane diffusing through the soil pore space, as ebullition is too quick to allow exposure to methanotrophs (Boone, 2000) and plant vascular transport shields methane from methanotrophic activity (Schimel, 1995;King et al., 1998;Verville 40 et al., 1998). The rate of methane consumption is usually highest immediately above the water table level, where high concentrations of methane formed from the underlying methanogens meet with sufficient oxygen levels from the overlying atmosphere (Dedysh et al., 2002). Rates of both methanogenesis and methanotrophy are highly dependent on temperature, with optimal metabolic rates (as determined in the laboratory) occurring at temperatures of around 25 • C (Dunfield et al., 1993). Of the two competing processes, methanogenesis has shown to be more temperature-dependent with higher reported rate changes per unit 45 warming (e.g. per 10 • C (Q 10 ): 5.3-16) compared to methanotrophy (Q 10 : 1.4-2.1) (Dunfield et al., 1993).
Estimating the methane exchange budget in Arctic tundra ecosystems and how it relates to temperature are challenging objectives, currently subject to considerable uncertainties. Ambient observations in northern Alaska over 29 years showed no clear increase in ambient atmospheric methane concentration enhancements during this period, despite noticeably warmer air temperatures . Direct observations of methane exchange, however, during the Carbon in Permafrost 50 Experimental Heating Research (CiPEHR) project, showed significant increases in methane emission under warming soil conditions . Multi-year carbon exchange data sets are rare, and challenging winter conditions at Arctic sites has led to many studies focussing largely on summer and early autumn periods (Euskirchen et al., 2012). The first yearround micrometeorological methane exchange measurements reported showed that cold season methane emission dominated the annual exchange budget, suggesting a predominant role of cold-season processes (Zona et al., 2016). At five Alaskan sites 55 (four on coastal plains and tundra dominated by sedges, grasses, and mosses within the northern coastal region surrounding Utqiaġvik and one over tussock-sedge dwarf-shrub, moss tundra at Ivotuk on the North Slope of the Brooks Range, Walker et al., 2005), substantial methane emission was reported to have occurred throughout the cold season, and although emission rates were lower than those measured during the warmer summer period, prolonged wintertime activity amounted to 50% ± 9% (mean ± 95% confidence interval) of annual emission. The authors suggest that this cold season emission may become more 60 important under forecasted climate conditions that include higher air temperatures and deeper snowpacks (Hay and McCabe, 2010;Zona et al., 2016).
We present here a two-year micrometeorological methane and carbon dioxide exchange data set, undertaken over an acidic tussock tundra site near the Toolik Field Station, Alaska, USA, on the north slope of the Brooks Range. Complimentary to airsurface exchange measurements, we report soil pore space methane, carbon dioxide and oxygen concentrations and soil water 65 content in the upper 40 cm, as well as soil temperature profiles at and near the site to a depth of 150 cm. All measurement systems were deployed year-round, providing near-continuous data coverage throughout both the summer growing seasons and the cold winter seasons in both years. The goal of this study was to investigate environmental controls that significantly impact the magnitude and direction of methane fluxes in this environment -particularly over the colder winter months -as environmental control-methane flux relationships during these periods are relatively poorly understood. These results expand 70 our coverage of year-round methane and carbon dioxide exchange data sets across different bioclimates and landscapes, as well as add to our growing understanding of carbon exchange dynamics in the Arctic tundra soils, and provide insights into how these dynamics may evolve under forecasted changing conditions in this region.  (Shaver and Chapin III, 1991). Vegetation in other areas of the measurement footprint were characterised as wet graminoid tundra (sedge and moss tundra). Both mineral and organic soil profiles were closely present together with different horizon depths. Soil organic carbon content, based on eight soil sampling pits to a depth of 90 cm around the flux tower, was highly variable, with A-horizon organic carbon concentrations averaging 10.3% (range of 7 to 14%) and B-horizon organic carbon concentrations averaging 2.4% (range of 1 to 4%) (Olson et al.,85 2018). Estimates of soil organic carbon stores in the flux footprint range from 5 to 25 kg C m −2 (Fig. S2).

Instrumentation
An aerodynamic gradient approach was utilised for observing air-surface methane and carbon dioxide fluxes. The aerodynamic gradient method was chosen for compatibility with a concurrent atmospheric mercury flux study (Obrist et al., 2017). This method has been shown to have greater variability compared to the more widely-used eddy covariance method over diel 90 time scales (Muller et al., 2009) though with reasonable agreement over longer time scales, with the caveat that the concentration gradients are precisely quantified using high-precision gas analysers (Zhao et al., 2019;Karlsson, 2017;Fritsche et al., 2008). Turbulent characteristics were measured using a Metek USA-1 sonic anemometer (Metek GmbH, Elmshorn, Germany), positioned 2.36 m above the tundra soil. Atmospheric sampling of trace gas concentrations was performed at heights z 1 = 0.61 m and z 2 = 3.63 m. In addition to atmospheric trace gas sampling, a soil trace gas system consisting of four to six 95 soil inlets in two vertical profiles installed between 15 and 20 m to the north of the flux tower (herein "flux tower profiles") allowed for determination of near-surface trace gas concentration gradients in soil pore air (Fig. S3). These profile locations were chosen to represent an organic matter-rich soil profile (herein "organic") and a less organic matter-rich profile (herein "mineral"). Data from these two soil inlet profiles were collected at two depths in Year 1 (10 and 40 cm), and three depths in Year 2 (10, 20 and 40 cm). Perfluoroalkoxy Teflon tubing from both the atmosphere and soil inlets were directed in a heated 100 conduit to an onsite field laboratory, and a solenoid valve system allowed sequential sampling between all inlets with switching interval of 10 minutes (Obrist et al., 2017). Methane and carbon dioxide concentrations were quantified using a Los Gatos 915-0011 ultra-portable greenhouse gas analyser (Los Gatos Research, Mountain View, CA, USA), factory-calibrated prior to installation and operating at 1 Hz. The analyser was zeroed using methane-and carbon dioxide-free zero air approximately every 6 weeks. Span calibrations were achieved using the internal calibration routine, as recommended by the manufacturer. Line 105 intercomparison tests were also performed with the same frequency by moving both inlets to the same height and sampling for between 12 and 24 hours (average 17 hours). Concentration differences between the sample lines during the intercomparisons were <0.001 µmol mol −1 (or 0.025% of the mean concentration gradient) for methane (mean methane gradients were -4 ± 7 µmol mol −1 ) and <0.1 µmol mol −1 (0.1% of the mean concentration gradient) for carbon dioxide (mean carbon dioxide gradients were 100 ± 200 µmol mol −1 ). Although the emphasis here is on methane flux magnitudes and dynamics, carbon 110 dioxide fluxes are discussed at length in order to understand corresponding respiration processes that help us constrain the influence of microbial activity on observed methane fluxes.
Within the flux tower soil profiles, temperature and volumetric water content (VWC) were also measured at depths of 10, 20 and 40 cm. Temperatures were measured using soil temperature probes (Model 107, Campbell Scientific Inc., Logan, UT, USA) and VWC was measured at the same depths using time-domain reflectometry (Model CS615-L Soil Volumetric Water 115 Reflectometers, Campbell Scientific Inc., Logan, UT, USA). Further from the flux tower (430 m to the north-east), Toolik Field Station operates two profiles of soil temperature (thermocouple) measurements to a depth of 150 cm (0, 5, 10, 20, 50, 100 and 150 cm). Despite the increased distance, these profiles were included in the current analysis as they provide a longer time series (measurements have been taken continuously since 1988) and information at deeper depths than the flux tower profiles.
A snow tower (Seok et al., 2009;Faïn et al., 2013) was installed prior to the first snowfall and recorded temperatures at 0, 10, 120 20, 30, 40 and 110 cm above the soil surface, thus measuring temperatures within the snowpack as it developed above each measurement height. The average snowpack depth over the site was measured daily using a camera set to automatically record images of reference snow stakes (Agnan et al., 2018). These depth measurements began in November 2014, and so the first snowfalls in that year were not recorded. Fluxes of methane and carbon dioxide were calculated using the aerodynamic gradient approach described by Edwards et al. (2005): where F represents the flux of either methane or carbon dioxide, k the von Kármán constant, u * the friction velocity, C i the concentration of atmospheric trace gas species in question at height i = [1, 2], z i the sampling height, d the displacement 130 height and Ψ i the stability-dependent integrated similarity functions for heat, as given by Businger et al. (1971), that are wellrepresented for scalars within a range of Obukhov lengths between -2.5 and 2. Values outside of this range (highly stable or highly unstable) were filtered from analysis, resulting in a loss of 1.1% of available data. Additionally, periods for which the friction velocity was below 0.08 m s −1 (Muller et al., 2009) were excluded from analysis. This resulted in a loss of 7.1% of all available data, with a slight seasonal bias (7.1% of winter (DJF) data compared to 2.9% of summer (JJA) data). Herein we 135 follow the convention of positive flux values representing emission, whilst negative values represent deposition.
Atmospheric turbulent characteristics (friction velocity and Obukhov stability) were calculated using the flux processing software EddyPro v.6.2.0 (Li-COR, Lincoln, NE, USA) using 30-minute averaging periods. Rotation of sonic data into mean wind vectors was accomplished using the double rotation technique and quality control tests for steady state and developed turbulent conditions were implemented according to Foken et al. (2004). Apparent sonic anemometer sampling height was 140 altered according to daily observed snow depth in increments of 5 cm, as were gradient intake sampling heights. Displacement height d was set as 0.7h c (canopy height, ∼0.2 m) during snow-free periods and at 0 m during snow-covered periods (Oke, 1978). As a single instrument was used for trace gas sampling at both intake heights, leading to a loss of temporal coverage within each 30-minute period (Woodruff, 1986), gaps in the concentration time series were estimated for each averaging period using a 4th order polynomial fit to the observed concentration time series. Average concentrations at each height were then 145 calculated from a truncated mean (10th-90th percentile) in order to reduce effects of outliers (Fig. S4).
Two-dimensional footprint analyses were undertaken for each 30-minute period using the method of Kljun et al. (2015) and fluxes for which the footprint intensity over the adjacent Toolik Lake was shown to be >20% of the total were removed from analysis. Analysis of energy balance closure showed that calculated turbulent and soil heat fluxes for snow-free periods, excluding fetches in the direction of Toolik Lake, accounted for approximately 88% of net radiative fluxes (linear least 150 squares, p <0.001). Gaps in both methane and carbon dioxide fluxes, resulting from quality control and instrument downtime/maintenance, were filled using the MDSGapFill function (Reichstein et al., 2005) within the R package REddyProc (Wutzler et al., 2018). The efficacy of this gap filling was tested against a randomly-selected validation set of size equal to 10% of available flux values (Fig. S5). Ecosystem respiration was approximated using half-hourly, gap-filled carbon dioxide fluxes, filtered to exclude times when incoming photosynthetically active radiation (PAR) was above 5 µmol m −2 s −1 (Natali et al., 155 2015). Periods during which no data fitting this criteria are available (i.e. polar day) were not gap-filled, resulting in incomplete temporal coverage for ecosystem respiration. Regression tree analysis (Sachs et al., 2008;De'ath and Fabricius, 2000;Breiman et al., 1984) was undertaken on 80% of observed flux data (20% validation fraction) using the TreeBagger function in Matlab 2016a (MathWorks, Natick, MA, USA). 500 cross validations were ran with a minimum leaf size of 1% of the training set size, with the tree with lowest mean squared error chosen as our predictive model. Deeper snowpacks are able to provide an increased thermoinsulation effect from cold air temperatures, particularly in the 180 early cold season (Maksimova et al., 1977;Sokratov and Barry, 2002), thus leading to warmer and less variable surface soil temperatures. This effect can be seen in the temperature pulses shown in the snow tower thermocouple data ( Arctic tundra ecosystems are highly heterogeneous within the scale of micrometeorological flux footprints (typically 10s 185 to 100s metres, Kljun et al., 2015;Fox et al., 2008), and during the winter, the combined effects of wind and topography lead to even greater spatial heterogeneity in snow depths and snow physical properties (Agnan et al., 2018). Sub-surface soil temperatures, which are further influenced by air temperature and downwelling radiation; overlying vegetation and snow; and soil properties and moisture, are likely highly spatially variable within the footprint as well. As a result, the limited number of soil temperature profiles within the flux measurement footprint may not be fully representative of the average temperature 190 6 within the flux footprint. The four soil temperature profiles we have available are separated both in space and in the soil properties in which they were installed. Fig. S9 shows the time series for each of these profile measurements at two common depths (10 and 20 cm). This time series shows that, whilst the absolute range of temperatures between measurements can be pronounced, the correlation between these soil temperature measurements is reasonable enough throughout most of the study period (mean R 2 = 0.64) to warrant use of Toolik Field Station data for further investigation of temporal trends. The decision 195 to use Toolik Field Station data was driven primarily by the deeper profiles measured here, as this information is vital to the primary outcomes of this study (see Section 3.4).

Annual and seasonal flux patterns
Half-hourly methane and carbon dioxide fluxes are shown in Figure 1, along with total cumulative values over the ∼two-year (727-day) study. Overall, based on combined raw and gap-filled data, the mean half-hourly methane flux showed an emission of 200 0.5 ± 0.5 mg C m −2 h −1 (herein, uncertainty is expressed as one standard deviation of the measured values). The distribution of methane fluxes (Fig. S10) showed a positive skew, as well a secondary peak in values close to zero. Cumulative diel sums gave a mean net daily flux of 11 ± 8 mg C m −2 d −1 . Over the study period as a whole the site acted as a net source of methane, with a cumulative methane emission of 8.6 g C m −2 (4.9 g C m −2 and 3.7 g C m −2 in the first and second years, respectively).
The overall mean carbon dioxide flux across the two years of measurements was 0.0 ± 0.2 g C m −2 h −1 (mean net daily flux of , with the site acting as a net source of carbon dioxide during both years of measurements. The distribution of carbon dioxide fluxes across the study (Fig. S10) did not show skewness as seen in the methane flux distribution, though it did show a higher level of kurtosis. During the 24-month measurement period, the site emitted a net carbon dioxide flux equivalent to 583 g C m −2 (485 g C m −2 and 89 g C m −2 in the first and second years).
Seasonality can be defined in a number of different ways depending on the processes of interest (Mastepanov et al., 2013); 210 initially, we followed similar definitions as those described by Zona et al. (2016) who investigated changes in methane fluxes based on surface (10 cm) soil temperatures. Periods where surface soil temperatures were above 0 • C were defined as "active" seasons (yellow shading in Fig. 1) and those where soil temperatures were below 0 • C were defined as "frozen" seasons (blue shading in Fig. 1). Zero curtain periods, where surface soil temperature remains close to 0 • C (± 0.5 • C) for prolonged time periods due to latent heat released or absorbed from soil water, were separated into "freezing" or "thawing" seasons.  Tables 1 and 2 give summary methane and carbon dioxide flux data, respectively, for these seasons and years so defined.
220 Table 1 shows marked differences in the magnitude and seasonality of methane fluxes between the two years. Cumulative methane emission in Year 1 was 1.3-fold that of Year 2, across a slightly shorter period (347 days compared to 380 days). All seasons showed net methane emission across the study period, with statistically significant differences (Student's two sample t-test, p <0.05) in the net daily flux between Year 1 and Year 2 for all seasons. The largest seasonal differences in net daily methane fluxes between years were for the active and frozen seasons (6 and 10 mg C m −2 d −1 , respectively). For Year 2, the 225 active season showed significantly higher methane emission compared to Year 1, releasing 1.7 g C m −2 (15 mg C m −2 d −1 ), or 46% of the annual total, compared to 0.9 g C m −2 (9 mg C m −2 d −1 ), or 18% of the annual total for Year 1. Conversely, the frozen season showed higher emission in Year 1, releasing 2.5 g C m −2 (16 mg C m −2 d −1 ), or 51% of the annual total, compared to 1.1 g C m −2 (6 mg C m −2 d −1 ), or 30% of the annual total in Year 2. Year 1 also showed higher methane emission in the freezing and thawing seasons, though these represented similar percentages of the annual total across both years (∼25% 230 of annual total for freezing and 3-4% of annual total for thawing). Zona et al. (2016) reported average cold season methane emission from five Alaskan Arctic sites of 1.7 ± 0.2 g C m −2 , accounting for between 37 and 64% of the total annual methane budget at these sites. The authors note that these contributions are higher than those estimated from previous models and periodic chamber observations. In our study, observations in Year    Carbon dioxide fluxes (Table 2) showed significant net emission throughout the entire cold season in Year 1 (471 g C m −2 or 1.9 g C m −2 d −1 ), followed by minor net emission in the subsequent active season (14 g C m −2 or 0.1 g C m −2 d −1 ). Year 2 freezing and frozen seasons showed lower carbon dioxide emission than Year 1, and the thawing season showed net carbon dioxide deposition, resulting in a combined cold season emission of 294 g C m −2 or 1.1 g C m −2 d −1 (38% lower than Year 1). Net active season carbon dioxide fluxes in Year 2 showed significant deposition, at -211 g C m −2 or -1.9 g C m −2 d −1 .

260
Annual and multi-season micrometeorological flux studies are rare for the Alaskan Arctic (Commane et al., 2017); however, the net annual carbon dioxide flux for Year 2 is within the range of values reported for wet sedge tundra (2 to 147 g C m −2 a −1 , Euskirchen et al., 2012Euskirchen et al., , 2017, and larger than for heath tundra (21 to 61 g C m −2 a −1 , Euskirchen et al., 2012Euskirchen et al., , 2017 or tussock tundra (13 to 15 g C m −2 a −1 , Euskirchen et al., 2012;Oechel et al., 2014). for Year 1 was significantly above the range previously reported for wet sedge tundra. Arctic tundra ecosystems are highly heterogeneous both physically and biogeochemically (Fox et al., 2008) and the area examined here is no exception. Seasonal two-dimensional footprint analyses (Fig. S11) showed a prdeominantly southerly footprint during all seasons, where fens and moist tundra are more abundant (Fig. S1). Importantly, the homogeneity in all seasonal footprints shown in Fig. S11 excludes the possibility that observed differences in seasonal flux magnitudes (i.e. higher cold season contributions in Year 1 relative to 275 Year 2) are due to flux footprint differences over heterogeneous surfaces. latitudinal comparison from three Alaskan tundra sites, Grogan and Chapin III (1999) reported significantly higher wintertime carbon dioxide efflux from Toolik, relative to Fairbanks (south of Toolik) and Sagwon (north of Toolik). This wintertime carbon dioxide efflux was correlated with warmer surface soil temperatures (quantified as 5 cm soil temperatures greater than

9
-5 • C), which at Toolik were relatively high due to thermal insulation by a substantial early snowfall. The extended period of 295 increased carbon dioxide emission in the Year 1 freezing season is likely also associated with the insulating effects of an early substantial snowfall and the associated warmer surface soil temperatures (explored in greater detail in the following section).
Importantly, in both years, MA28 carbon dioxide fluxes do not completely cease during the cold season and always maintain a small emission throughout winter (up to 0.1 g C m −2 h −1 ).
In summary, the two years showed substantial temporal differences in methane fluxes, with Year 1 showing higher methane 300 emission throughout most of the cold season (100% greater), contributing a high fraction (82%) of annual net methane emission. This is in contrast to Year 2, which experienced a continued decline in methane emission that began early in the freezing season, resulting in a relatively low contribution (54%) to the annual total. Annual carbon dioxide flux magnitudes were most similar to other wet sedge tundra measurements; show strong seasonal trends with relatively high respiration in the freezing season and prolonged but low carbon dioxide emission in the frozen season; and carbon dioxide deposition during the thawing 305 and/or active season, largely in agreement with carbon dioxide flux patterns reported for northern tundra ecosystems before.
Inter-annual comparison, however, showed cold-season carbon dioxide fluxes that were 38% higher in Year 1, also largely driven by slower and later declines in carbon dioxide emission fluxes during the freezing and frozen periods.

Soil temperature relationships
Continuous cold season methane flux data at the ecosystem level are rare for Arctic tundra ecosystems. Given the strong dom- with substantial variability and a strong hysteresis between freezing and thawing periods. Zona et al. (2016) suggested that 330 temperature-dependent decreases in the near-surface methane oxidative capacity were largely responsible for the slow attenuation of methane fluxes in the early frozen period (here F-G), noting that sites with the largest and warmest active layers displayed the slowest decrease in methane fluxes. Soil pore gas concentration measured in our study show that oxygen levels within the upper 40 cm were sufficient (>17%) to ensure methane oxidation in this zone across the entire study period and hence that the top 40 cm soils were continuously oxic (Fig. S3). We also note that soil pore gas measurements took place in an 335 elevated, drier tussock region and that the thickness of the upper oxidative region is expected to be smaller in lower depression areas and more water-saturated wet sedge regions (Gebauer et al., 1996).  whilst in Year 1 these values were more variable, declining from 0.1 g C m −2 h −1 to near-zero before increasing again to 0.05 g C m −2 h −1 . Respiration fluxes increased in both years prior to and during the thawing periods as surface soils warmed.

365
It is noteworthy that the beginning of these increases coincided with turning points in the MA28 methane flux-surface temperature relationship (points B and G). In the active season, as discussed previously, Year 2 showed higher net carbon dioxide emission than Year 1. Similarly, since surface soil temperatures were much warmer in the Year 1 cold season, heterotrophic respiration in this upper oxic soil region remained high relative to Year 2 leading to higher cold season cumulative carbon dioxide losses ( Table 2). The strong relationships between cold season ecosystem respiration fluxes and surface soil temperature, 370 and the relative similarity between the two years, is consistent with patterns reported in previous research (Lüers et al., 2014;Euskirchen et al., 2012;Björkman et al., 2010), and largely explain differences in the temporal trends of ecosystem respiration flux between the two years. This result suggests that changing heterotrophic microbial respiration in the upper soil region is not a suitable explanation for the differences in methane flux-surface soil temperature relationships observed between Years 1 and 2. including temperature, water table depth, oxygen availability and Eh, soil organic matter content, soil pH, soil texture and soil 380 mineralogy, though soil temperature and water table depth are often identified as the major of these controlling factors (Le Mer and Roger, 2001;Yvon-Durocher et al., 2014;Gulledge et al., 1997). We employed a regression tree approach (Sachs et al., 2008) to explore non-linear relationships between observed methane fluxes and variables identified in the literature as influential to methanotrophic/methanogenic activity. Of the known variables, soil organic carbon content, pH, texture and mineralogy cannot explain changes in fluxes over short time periods and hence were not included. Additionally, pore-space 385 oxygen concentrations were not included in the analysis since it was measured only in the upper 40 cm and remained oxic throughout the entire study period. For the regression tree, we hence used soil temperature data from the surface to 150 cm depth, as well as surface VWC. Daily values were chosen in order to reduce the influence of diel variability, with net daily sums (in the case of methane fluxes) and mean daily values (for temperature and water content) as inputs into the model. The outcome of this analysis is shown in Fig. 3a, along with the time series of net daily methane fluxes used to build the model 390 (Fig. 3b). Horizontal lines and coloured shading in Fig. 3b show mean ± one standard deviation of the input methane flux data, as grouped by the predictive model. The predictive capability of this model, tested against a randomly-selected validation set representing 20% of the available input data, shows an R 2 value of 0.69 (p <0.001).
The two variables that most effectively cluster methane flux values within the hyper-dimensional data space are soil temperatures measured at 100 cm and at 10 cm depths. Critically, these two temperature variables separate, and likely explain, 395 substantial methane flux differences observed during the frozen seasons between Year 1 and Year 2. Specifically, the frozen season in Year 2 was largely represented by regression tree outcomes when temperatures at 100 cm soil depth were below -2.4 • C, marked in Fig. 3b in blue and red. This is the only season during the study period when 100 cm soil temperature fell below this threshold (see also Fig. 4). Measurements throughout the active layer by Gebauer et al. (1996) at nearby Imnavait Creek suggest that it is highly likely that soils at these depths are anoxic and thus methanogenesis is the dominant relevant 400 microbial process taking place here. In contrast to Year 2, the Year 1 frozen season was largely separated from the remainder of the study period (with 100 cm soil temperatures above -2.4 • C) when 10 cm temperatures were simultaneously below -0.6 • C (purple and green shading). During this season, methane emission values were amongst the highest observed throughout the entire study period. As discussed previously, our own soil pore gas measurements show that these surface soils are oxic and thus methanotrophy is here likely the dominant relevant microbial process. It must be noted that these temperature thresholds 405 do not represent mechanistic limits but instead the most effective clustering of observed data, based on the chosen environmental parameters. Even so, the reasonably good predictive capability of this model provides strong evidence that a critical reason for strong flux differences between frozen season methane fluxes was differences in deep soil temperature between the two years. More precisely, these first two results from the regression tree analysis show a threshold temperature value at a depth of approximately 100 cm that is linked to low methane emission in Year 2, and an additional threshold temperature value at a 410 depth of approximately 10 cm that is linked with high frozen season emission in Year 1.
A decrease in methanogenesis below a temperature threshold around -2.4 • C as suggested by the model is in reasonable agreement with several experimental laboratory studies. Incubation studies investigating the temperature dependency of methanogenesis in Arctic soils have shown that it can take place at sub-zero temperatures, though at greatly reduced rates. Rivkina et al. (2004) reported substantial methane production at temperatures of -1.8 • C in Siberian permafrost soils, as well as 415 methane production in these soils at temperatures as low as -16.5 • C, though at a rate 100 times lower than at -1.8 • C. Similarly, Panikov and Dedysh (2000) observed minor methane emission (∼0.1 mg C dm −3 h −1 ) from Siberian peat bog soils at -20 • C that increased by an order of magnitude after thawing. Chowdhury et al. (2015), using soils from Barrow, Alaska, showed evidence of substantial methanogenesis in organic and mineral active layer soils kept at 4 and 8 • C, yet this was not observed in permafrost soils or in active layer soils kept at -2 • C. Similarly, for methanotrophic activity, temperature dependencies have 420 also been observed, again with lower microbial activity reported at lower soil temperatures. Jørgensen et al. (2015) reported an exponential relationship between temperature and methane uptake in unsaturated Arctic tundra soils. At 18 • C they observed a deposition flux of 192 µg C m −2 h −1 -this decreased to 24 µg C m −2 h −1 in soils kept at -4 • C. Richter (2019) also observed a temperature-related decrease in methane oxidation in A-and B-horizon soils sampled near Toolik Lake, to below detection levels at temperatures below -2 • C. Based on this evidence and our regression tree analysis, we suggest that the separation of 425 methane fluxes in the Year 2 frozen season (the period with the lowest methane fluxes across the entire study period) is linked to an inhibition of methane production due to low soil temperatures in deep, anoxic soil horizons. Further, we suggest that the separation of methane fluxes in the Year 1 frozen season (the period with some of the highest methane fluxes across the entire study period) is linked to an inhibition of methane oxidation due to low soil temperatures in oxic, surface soil horizons.
The full time series of interpolated soil temperature profile measurements at Toolik Field Station is shown in Fig. 4, along   430 with the -2.4 • C isotherm identified in the first grouping of the regression tree analysis. These profile data show clearly that the deeper soil horizons never reached the cold temperatures in the Year 1 frozen season that they did in the Year 2 frozen season -13 in fact, the -2.4 • C isotherm did not move below 70 cm depth in Year 1. Taking as a starting point the only methane production rate in the aforementioned laboratory studies given per soil volume (0.1 mg C dm −3 h −1 at -20 • C, Panikov and Dedysh, 2000), the 80 cm of soil below this depth known to be above -2.4 • C could presumably sustain a methane production rate on 435 the order of 80 mg C m −2 h −1 . The MA28 methane flux (black line, lettering corresponds in time to that given in Fig. 2) shows cold-season decreases that accompany cold temperature pulses into the deeper soil horizons. This is most readily seen in Year 2 (F-G) with the greater contrast in soil temperatures, though it is also noted that in Year 1 the onset of a decrease in methane emission (B-C) corresponds in time to a lowering of the -2.4 • C isotherm and a still-perceptible cold temperature pulse to lower soil horizons. This highlights that the limits identified in the regression tree analysis (i.e. -2.4 • at 100 cm depth) are not claimed 440 to be mechanistic, yet they still provide valuable insight into competing methanogenic/methanotrophic processes within the soil profile. Further, whilst soil methane concentrations were not measured deep enough to confirm methane production at 100 cm depth, the sustained methane gradient in the upper 40 cm of the organic profile during the Year 2 cold season, and especially the Year 1 cold season, provide supporting evidence for production of methane below this surface layer. Frozen season differences in the MA28 respiration flux (green line) and, particularly, the ratio of methane flux to respiration (yellow line) further reiterate 445 the disconnect between methane production and respiration that was highlighted in the discussion surrounding Fig. 2.
The shading in Fig. 3 shows that, whilst cold season methane fluxes could not be grouped together within the hyperdimensional data space for Years 1 and 2, freezing and thawing periods largely were. This suggests that the remaining seasonal methane flux dynamics can be related to the balance of continued methanogenesis at depth and methanotrophic activity near the surface. Specifically, predicted methane fluxes outside the frozen seasons (orange, yellow, and brown shading) are separated 450 according to temperatures at 150 cm and 20 cm, yet follow a similar pattern whereby colder temperatures at depth (suggesting inhibited methanogenesis) and warmer surface temperatures (suggesting enhanced methanotrophy) jointly lead to smaller predicted methane fluxes (see also Fig. 4). One exception to these general patterns is during active seasons at times when surface VWC is greater than 0.65. During these periods the largest mean methane fluxes (as well as the largest methane flux variability) of any grouping were observed (pink shading). This observation is consistent with studies reporting increased methanogenic 455 activity relative to methanotrophy under water-saturated soil conditions due to reduced oxygen diffusivity and highly reducing conditions in otherwise oxic surface soils (Le Mer and Roger, 2001). Evidence of this reduced oxygen diffusivity, as well as inhibition of gas diffusion through the soil profile, can be seen in the soil pore gas measurements in Fig such evidence in the Year 1 thawing season. Decreased gas diffusivity during these periods likely contributed to a suppression of the methane flux, which were amongst the lowest seen throughout both years (leaf group 6 in Fig. 3).

Implications for Arctic methane fluxes
Considerable debate exists over the potential future of methane fluxes in the Arctic tundra under future climates . Hydrologic modelling under IPCC forecasts by Hay and McCabe (2010) predicted warmer air temperatures with greater 465 precipitation, leading to the suggestion that methane fluxes may increase as labile carbon becomes available due to permafrost thaw. Warming experiments undertaken in the field in Alaska have also shown that warmer and wetter soils resulting from increased snow cover emit considerably more methane during the active period . Zona et al. (2016) stated that cold season fluxes made up to 64% of their reported annual methane emission, due to relatively low but consistent emission over a large portion of the year. Our Year 2 cold season methane emission accounts for only 54% of the annual total, yet for 470 this study site is still greater than the cold season total given by Zona et al. (2016) for their study sites (2.0 g C m −2 here compared to 1.7 ± 0.2 g C m −2 ). It is significant that Ivotuk, the upland tundra site in this study, exhibited the largest cold season, and the largest annual, net methane emissions. Euskirchen et al. (2017), in their four-year study of methane exchange at Imnavait Creek (another upland tundra site), reported lower net cold season methane emission of 1.4 ± 0.3 g C m −2 , yet note that they only had partial zero curtain (freezing season) data for one year. The methane emission dynamics we observed Modelled forecasting of Arctic methane fluxes is typically undertaken using air temperature data, due to its relative ease of measurement and prediction, and the assumption that air temperature is closely linked to soil temperature (Riley et al., 2011;Koven et al., 2013;Zhu et al., 2014). An analysis of a 29-year record at Barrow, Alaska, however, showed no correlation 485 between increasing air temperature and methane concentration anomaly , suggesting that air temperature is an inadequate variable for predicting methane fluxes. Air and soil temperature measurements at Toolik Field Station taken since 1988 (Fig. S7) show that, whilst Year 1 was the shortest winter (defined here once more as the period where MA28 air temperatures <0 • C) on record, both Year 1 and Year 2 were unusually warm and similarly ranked in regards to total FDDs (3rd and 4th warmest on record, respectively). Similarly, minimum cold season soil temperatures in the upper (20 cm) 490 and deep (100 cm) horizons were unusually warm in both years, comparative to the long-term record (2nd and 3rd highest values, respectively). Given the relative similarity in the temperature anomaly of both years compared to the last 28, the large differences in cold season methane emission (2 times larger in Year 1 compared to Year 2) is unlikely related to a simple linear relationship with increasing air, or even soil, temperatures.
Instead, we suggest the presence of a deep soil temperature threshold in anoxic horizons above which cold season methano-495 genesis -and hence net methane emission -remains high. Climatologically, there was little difference between Years 1 and 2 in terms of cold season FDDs, or minimum soil temperatures, relative to the previous 28 years. Yet, as suggested in the regression tree analysis and the temperature profiles in Fig. 4, any such temperature threshold was not crossed in the Year 1 cold season, allowing methanogenesis to continue relatively unabated. Snow profile measurements (Fig. S8) show that, in addition to both winters experiencing relatively warm air temperatures, deeper snow in Year 1 likely insulated the underlying soil such 500 that anoxic soil horizons cooled at a much slower rate. If, as has been predicted (e.g. Hay and McCabe, 2010), the Arctic continues to warm and precipitation increases, high methane emission winters will likely become more prevalent in the future, particularly also if enhanced summertime warming pulses penetrate deeper in the soil profile. Our observations highlight the need for more sophisticated modelling of temperature regimes in the forecasting of methane emission. More importantly, we suggest that the increasing number of year-round ecosystem flux measurement sites operating in Arctic regions should monitor 505 soil temperatures throughout the entire active soil region, rather than limit observations to the upper surface horizons. This is particularly important for large-scale soil monitoring networks such as the Soil Climate Analysis Network (SCAN), the outputs from which are important for enabling gridded modelling products for quantifying regional-scale carbon fluxes. Temperature data throughout the entire active soil profile, preferably in conjunction with estimates of soil redox conditions, would help to further elucidate the competing microbial processes that drive methane fluxes at the surface.

Conclusions
Year-round measurements of ecosystem-scale methane and carbon dioxide fluxes were undertaken at Toolik Field Station in the Alaskan Arctic over two years. Annual carbon dioxide exchange budgets suggest that these observations are representative of wet sedge tundra, with seasonal patterns that are characteristic of the Alaskan North Slope generally. Net methane and carbon dioxide fluxes in the Year 2 cold season (2.0 g C m −2 and 294 g C m −2 respectively, over 269 days) were similar in 515 magnitude to values reported in similar studies, and positive correlations between surface soil temperature and methane were observed as previously reported by Zona et al. (2016). Year 1 cold season net methane and carbon dioxide fluxes, however, were 100% and 38% higher, over a shorter cold season (22 days shorter). Relationships between respiration fluxes and surface soil temperature were similar between years and with those reported in the literature, suggesting that warmer soil temperatures in the oxic surface horizon can largely explain the differences in annual cold season carbon dioxide budgets between the two 520 years. Methane flux and surface soil temperature, by contrast, showed almost reversed relationships between the two years, suggesting that surface soil temperature was not always sufficient to explain methane emission dynamics over the course of this study.
Whilst cold season soil temperatures and FDDs were similar across both years (relative to the 28-year record), we observed that deeper snow pack in Year 1 led to significantly warmer soil temperatures, particularly in the deeper portion of the active 525 soil profile. A regression tree analysis shows that high Year 2 frozen season methane fluxes were clustered from other data along the deep (100 cm) temperature axis at a threshold of -2.4 • C, suggesting inhibited methanogenesis in deeper, anoxic soil horizons. The highest cold season fluxes (among the highest of the two-year study) were observed during Year 1, when deep soil temperatures remained above this threshold whilst surface temperatures were simultaneously below -0.6 • C, suggesting limited methanotrophy in upper soils being unable to offset methane production at depth. From our data we cannot reliably 530 state that these thresholds represent mechanistic limits, only that they highlight a pattern of temperature dynamics between the upper, oxic layer and the deeper, anoxic layer that are key to controlling surface methane fluxes via their limiting influence on the competing processes of methanotrophy and methanogenesis. These temperature dynamic patterns further explain methane flux dynamics outside of the frozen season during both years. Our results suggest that high cold season methane emission may be associated with warmer atmospheric temperatures and deeper snowpacks, and highlight a need for measurement and mod-535 elling of soil temperatures throughout all seasons, and throughout the entire active soil profile. Such expansion in observation capacities will allow more accurate prediction of potential changes in the annual methane exchange budget in Arctic tundra regions.
Data availability. Data of this project, including atmospheric mercury, carbon dioxide and methane concentrations and mercury concentration measurements in vegetation, soils and snow are archived with NSF's Arctic Data Center (https://arcticdata.io/), DOI: 10.18739/A21Z41S5S.

540
Author contributions. Dean Howard led the analysis of data and manuscript preparation. Yannick Agnan, Detlev Helmig, Yu Yang and Daniel Obrist designed the field measurement campaign. Yannick Agnan, Detlev Helmig and Daniel Obrist conducted field measurements.
All authors contributed to writing of the manuscript.
Competing interests. The authors declare that they have no competing interests, financial or otherwise.
Acknowledgements. This study was funded by a by the U.S. Department of Energy (DE-SC0014275) and the U.S. National Science Foun-545 dation Office of Polar Programs (# OPR 1304305 and 1739567). The authors would like to acknowledge the help of Toolik Field Station for measurement support and the Toolik Field Station Environmental Data Center for providing soil temperature data. We also thank Donald A.
(Skip) Walker for allowing the use of the Toolik field maps reproduced in the supplemental material. We wish to also thank the anonymous reviewers whose constructive feedback contributed significantly to the development of the manuscript. 26 Table 1. Overview of methane fluxes split according to season. Values in parentheses in cumulative columns are percentages of the total for that year. Daily differences are the differences in mean daily fluxes for that season between the two years, with p-statistics from Student's two-sample t-test giving the significance with which the null hypothesis (values come from distributions with same mean) can be rejected. Year