Cold season emissions dominate the Arctic tundra methane budget

Fig. 2 shows continuous eddy flux data for five tundra sites in Alaska: three in Barrow (CMDL, BEO, and BES), one in ATQ, and one in IVO (Materials and Methods). Methane emission rates from the cold seasons (September to May) were comparable to (e.g., BEO and ATQ; Fig. 1 C and D) or higher than (e.g., CMDL; Fig. 1B) emissions in summer over a prolonged period. Cumulative emissions for the cold season averaged 1.7 ± 0.2 [mean ± confidence interval (CI)] g C-CH₄ m⁻² at our five sites, accounting on average for 50 ± 9% (mean ± CI) of the annual budget (BES, 37%; BEO, 43%; CMDL, 64%; ATQ, 47%; IVO, 59%). Cold-season emissions dominated the annual CH₄ budget in the driest sites (CMDL, ATQ, IVO), representing a notably higher contribution than previously modeled (6) in other continuous permafrost sites (35%) and also higher than observed year round in boreal Alaska [40%, using periodic sampling of static chambers (20)]. The boreal systems are underlain by discontinuous or sporadic permafrost and are therefore subject to different soil processes than Arctic sites underlain by continuous permafrost (which prevents drainage for extended areas for example).

The highest fall and winter CH₄ fluxes were observed at IVO, an upland tundra site (with a water table below the surface for most of the summer), which had the longest zero curtain period (101 d; Table S1), the warmest soil temperatures during the cold season (Fig. 3 and Fig. S1), the deepest snow depth (SI Materials and Methods), and the deepest active layer (Fig. S2 A and B). Soil temperatures were also poised near 0 °C for more than 90 d at much wetter sites near Barrow (BES). In both cases, the zero curtain period lasted as long as, or longer than, the summer season (Fig. S1 and Table S1). Based on direct measurement of the active layer depth and on soil temperature data, the maximum thaw depth did not begin to decrease appreciably until November or later in all of the sites measured (Fig. S2 A and B), even though the surface froze in September. During the zero curtain period, we observed strong CH₄ emissions from all five sites, 0.3–2.4 g C-CH₄ m⁻² (Fig. 2), albeit somewhat lower than the peak summer season CH₄ fluxes observed. The overall contribution of these zero curtain periods to annual emissions was important because of their extended duration (Fig. 2, Fig. S1, and Table S1): emissions of CH₄ during the zero curtain period alone contributed ∼20% of the annual budget (BES, 18%; BEO, 20%; CMDL, 20%; ATQ, 16%; IVO, 32%).

A few previous studies reported measurements of Arctic CH₄ fluxes during the fall (6, 7, 9, 10), but the measurements did not extend to winter and spring. We found that sites with similar summertime CH₄ fluxes had different zero curtain emissions because of different durations and depths of unfrozen soil (Fig. 2 and Fig. S2). For example, summertime cumulative emissions in IVO were 1.9 g C-CH₄ m⁻² in 2013 and 2.7 g C-CH₄ m⁻² in 2014, similar to the 2.3 g C-CH₄ m⁻² (in both years) at BES. However, cumulative CH₄ emissions during the zero curtain were much higher in IVO (2.4 and 2.1 g C-CH₄ m⁻² in 2013 and 2014, respectively) than BES (0.9 and 0.7 g C-CH₄ m⁻² in 2013 and 2014, respectively) probably because of interacting effects of greater CH₄ production at IVO, the inhibition of surface oxidation in the fall (Fig. 1), and the deeper thaw depth delaying the complete soil freezing in IVO (Figs. S1 and S2). The emissions of CH₄ produced deeper in the soil continued during the cold season, presumably through cracks and pathways in the near-surface frozen soils (7).

Linear mixed effects modeling (SI Materials and Methods) suggested that the depth of the active layer was a critical control on CH₄ fluxes during the summer. The presence of this unfrozen soil layer in the fall and early winter was also a major control on cold season CH₄ emissions; warmer soils resulted in greater CH₄ emission over the entire year. The importance of warm soil temperatures and deep active layer is consistent with the observed higher winter emissions in IVO, where soil temperature at 15 and 30 cm below the surface never dropped below approximately −8 °C compared with at or below −15 °C at the northern sites (e.g., BES and ATQ). The observed CH₄ emissions during fall and winter are consistent with data showing significant microbial populations and metabolic activity at and below 0 °C in the Arctic (16, 21), reflecting the availability of unfrozen water films (22) under these conditions (16). Measurable metabolism has been observed down to −40 °C (23), and CH₄ production has been observed down to −16 °C (21, 24). Soil particles maintain liquid water films until a temperature of at least −10 °C (25), and this unfrozen water can sustain microbial metabolism and greenhouse gas production (26), even as the soil bulk water freezes (25). The direct effect of higher temperature on metabolic activity and the indirect effect of temperature through greater liquid water volume should result in a larger population size and more activity in the methanogenic (i.e., methane-producing) community in the winter at IVO compared with the other, colder, sites. Unfortunately, IVO is the only tower collecting CH₄ fluxes and environmental variables continuously year round over upland tundra at this latitude in Alaska. Therefore, we encourage the establishment of similar upland sites in the Arctic to confirm these observations.

Across all our sites, areas of lower inundation (i.e., less surface area with water table above the surface for most or all of the growing season) had the greatest percentage of total emissions from the cold season, with the highest emissions from IVO with <5% inundation (Fig. 2). In contrast, most modeling studies limit CH₄ emissions to areas with inundated or saturated soils (27). The observed CH₄ emissions that persisted, even when temperatures were well below 0 °C (Fig. 2), present a remarkably uniform temperature response with a decrease in emission rates as soil temperatures drop (Fig. 3). The fall fluxes show clear relationships with declining soil temperature in the active layer, with little discontinuity in the flux relationship with soil temperature as the soils freeze (Fig. 3). It is likely that freezing of the surface soils decreases near-surface CH₄ oxidation (Fig. 1), maintaining net soil CH₄ emissions even as decreasing soil temperatures results in decreasing CH₄ production rates. At IVO, warmer soil and deeper thaw depth (and therefore greater metabolically active soil volume) resulted in the highest cold season emission rates. This seasonal pattern is very different from that reported by Mastepanov et al. (7, 10), who showed a drop in emissions in late summer/early fall from Greenland tundra, followed by large late-fall CH₄ emissions peaking during complete freezing of the active layer. We instead found fall emissions were persistent until the soil temperatures were well below 0 °C (Fig. 2), with a few instances of sporadic, exceptionally high emissions, e.g., in IVO (Fig. 2) contributing just ∼15% of the zero curtain emissions and ∼5% of the total annual CH₄ emissions. The underlying sensitivity of CH₄ fluxes to temperature at our sites was, on average, a factor of 2.7 (Fig. 2) for a temperature rise from 0°, to 5 °C, slightly more sensitive than the global mean described by Yvon-Durocher et al. (2).

Spring CH₄ fluxes also increased with increasing active layer temperatures (Fig. 3). The northern sites (e.g., BES and BEO; Fig. 3) showed prompt, steep increases in CH₄ emissions coincident with increasing soil temperatures. The southernmost site (IVO) showed a very different pattern, with apparently much lower temperature sensitivity of net fluxes in the spring vs. fall (Fig. 3). Unlike the wet tundra sites, there is substantial seasonal hysteresis at IVO, likely reflecting a combination of CH₄ oxidation in the spring and summer in the warmer, dry surface layers and CH₄ storage in the deepening, porous active layer. Also, methanogenesis may be stimulated by reduced oxygen in the unfrozen active layer, because the frozen surface (Fig. 1) slows diffusion of oxygen into the soil column (28).

Microbial consumption of CH₄ in the near-surface soil layer (methanotrophy) can be very active in summer (28) but is inhibited by near-surface soil freezing (28, 29). Thus, the fraction of CH₄ escaping to the atmosphere likely increases as the soil surface freezes in the fall. The wettest sites, such as Barrow-BES, where the water table was on average above the surface for the entire measuring period (Fig. S2 C and D), presumably had low levels of surface oxidation of CH₄. Therefore, this site showed the greatest relative decrease of cold season CH₄ fluxes compared with summer (Fig. 2) because decreasing temperatures reduced CH₄ production, but because oxidation rates were low, there was little benefit from suppression of oxidation in the surface layer in fall.

Our measurements of CH₄ emissions from Arctic tundra are more extensive in both time and space than what have been used to develop and test existing models. Annual CH₄ emissions rates from noninundated Arctic tundra (<20% surface water; Fig. 2) are comparable to those of inundated environments. Most models map CH₄ fluxes to the Arctic landscape using inundation (27), thus dramatically underestimating the emitting area in the Arctic, including during the cold season. The zero curtain interval in fall and winter, and even the period of frozen soils in winter, produce significant, previously underestimated, CH₄ emissions (27). Our work provides the basis for parametric representation of these fluxes and highlights the critical importance of driving models with subsurface soil temperature, and not air temperature.

Regional CH₄ fluxes calculated from aircraft observations (30) show a strikingly consistent pattern to our eddy flux data (Fig. 4), notably including the persistence of CH₄ emissions into the cold season. The regional aircraft fluxes derived from the CARVE (Materials and Methods, SI Materials and Methods, and Fig. S3) flights were at times lower than the mean of the EC tower fluxes, as has been observed previously in point-scale and regional-scale flux comparisons (SI Materials and Methods). Global-scale measurements (HIPPO; Materials and Methods) detected a large enhancement of CH₄ in the Arctic in early November, peaking in the boundary layer of the northern high latitudes (Fig. 5). Because of the flight plans of the HIPPO flights conducted in 2009 to 2011, fluxes could not be calculated from the HIPPO data. However, the HIPPO data are important to understanding whether the CH₄ fluxes calculated at the flux towers and during CARVE are relevant to CH₄-mixing ratios on the global scale. In the North Slope vicinity (71° N > latitude >65° N), CH₄ is enhanced compared with the global mean, but there is no corresponding elevation of CO, indicating that the CH₄ sources are not associated with transported pollution or fossil fuel burning (Fig. 5B; we have only considered CH₄ data between 65° N and 71° N to remove the influence of CH₄ enhancements observed over open leads in sea ice (32)]. By contrast, in January, there were air parcels with high CH₄ consistently associated with CO enhancement, indicating a dominant anthropogenic source of CH₄ compared with the global mean. During this time CH₄ was likely transported from lower latitudes (31). Overall, the HIPPO data are consistent with a substantial biogenic CH₄ source over northern Alaska in fall and with our finding of strong late season biogenic emissions on both a local and regional spatial scale.

Recent estimates using inverse modeling of atmospheric concentration data give CH₄ emissions from Arctic tundra wetlands in the range from 16 ± 5 Tg CH₄ y⁻¹ [from CarbonTracker (32)] to 27 (−15 to 68) Tg CH₄ y⁻¹ (8). Extrapolating our average CH₄ emissions rates to the Circumpolar Arctic tundra (SI Materials and Methods) yields an estimate of 23 ± 8 Tg CH₄ y⁻¹ from Arctic tundra, similar to these previous estimates (8, 32). Our estimated CH₄ cold-season emissions as well as those from inverse analysis (27, 32) are significantly higher than that estimated by land-surface models (27, 32). This difference was thought to be linked to anthropogenic emissions, because biogenic emissions were assumed to be negligible during the cold season (27, 32). Overall, the seasonal patterns estimated by models (27) are very different from ours and generally do not include the substantial cold season CH₄ emissions found here. Our finding of large cold-season biogenic emissions from tundra reconciles the atmospheric observations and inverse model estimates without the need to invoke a large pollution influence.

Data for this study were collected from five EC towers across a 300-km latitudinal transect in the North Slope of Alaska. The three northernmost towers (BES, BEO, CMDL) are in the vicinity of Barrow, AK, where mean annual temperature is −11.3 °C and summer precipitation is 72 mm for the 1948 to 2013 period (37). The fourth site, ATQ, is about 100 km south from Barrow. Mean annual temperature and summer precipitation in ATQ are −10.8 °C and 100 mm, respectively, for the 1999 to 2006 period. The most southerly site (IVO) is located near the IVO Airstrip at the foothills of the Brooks Range Mountains, about 300 km south of Barrow, with a mean annual temperature and summer precipitation of −8.9 °C and 210 mm, respectively, for 2003 to 2008. The average snow depth was about 0.3 ± 0.1 m (mean ± SE) in BEO/BES, 0.2 ± 0.2 in ATQ, and 0.4 ± 0.1 in IVO. The vegetation is classified as W1 in Barrow [wet coastal plain dominated by sedges, grasses, and mosses (38)], as W2 in ATQ (tundra dominated by sedges, grasses, mosses, and some dwarf shrubs <40 cm tall), and as G4 in IVO [tussock-sedge dwarf-shrub, moss tundra (38)]. The land-cover types in Barrow and ATQ together are representative of about 60% of all arctic wetlands (38), whereas IVO represents the dominant vegetation type in Alaska. The zero curtain period (19) was defined as the fall period, when the soil temperature at 15- to 20-cm depth (the last soil layer to freeze in our system) was between 0.75 °C and −0.75 °C (Fig. S1 and Table S1). Among the Barrow sites, CMDL presented the least ice-wedge polygon development and is the driest. The presence of low center polygons results in the presence of wet waterlogged ponds interspersed with drier microtopographic features (high center polygons and/or polygons’ rims) in the BEO site; BES is located in a vegetated drained lake with restricted drainage, the low topographic results in which being the wettest site.

A wide range of meteorological variables were measured at each of the five EC towers, including photosynthetic active radiation (PAR), which was measured with quantum sensors (LI-190; Li-COR) in all sites; net radiation was recorded using a net radiometer [REBS Q7 (Radiation & Energy Balance Systems, Inc.) in BES, BEO, CMDL, and ATQ; and an NR Lite (Kipp & Zonen) in IVO]; incoming solar radiation was measured using pyranometers (CMP3; Kipp & Zonen) in all sites; air temperature and relative humidity (RH) were measured with a Vaisala HMP 45 (in CMDL, BES, BEO, and ATQ), and a Vaisala 155a in IVO (and in BEO after 2013); soil heat flux at −2- to −5-cm depth was measured in four to six locations in all sites with REBS HFT3 (Radiation & Energy Balance); soil temperatures at different depths (at surface −1-, −5, −10-, −20-, and −30-cm depth in BES; at −5, −15, and −30 cm in four profiles in ATQ; and at the surface −5, −15, −30, and −40 cm in four profiles in IVO) were measured with thermocouples (either type-T or type-E; Omega Engineering); soil moisture was measured with a Water Content Reflectometer CS616 (Campbell Scientific) inserted perpendicularly (0–30 cm) or diagonally in different soil layers (0–10 cm and −20/−30 cm) in BES, or horizontally at different depths in the soil (−5, −15, and −30 cm in two different profiles in ATQ; and in three profiles in IVO). More details on these measurements can be found elsewhere (15, 18, 34). Snow depth was measured with Sonic Ranging Sensor in BEO/BES, ATQ, and IVO with an SR50A-L snow senor (Campbell Scientific).

Thaw depth and water table were measured manually during the summer and autumn on a weekly basis in the most accessible sites (BES, BEO, and CMDL) and once every 1–2 months in the more remote sites (ATQ and IVO). A graduated steel rod was used for the thaw-depth measurements, and PVC pipes with holes drilled every centimeter on their sides allowed for the water table measurements (15). Thaw-depth measurements were performed about every 5 m at 45 points in CMDL (three transects of 15 m each in the footprint of the EC tower); at 20 points in BEO, ATQ, and IVO (one transect in the footprint of each EC tower); and at 50 points every 4 m in the footprint of the EC tower at BES. Water table measurements were performed in these same plots for all sites with the exception of CMDL (where the installation of PVC was not allowed).

Half-hourly fluxes were calculated using the EddyPro software (LI-COR), applying the following procedures and corrections: a despiking procedure of fast raw signals (39); the time lag between vertical wind velocity and scalar concentrations was computed based on the maximization of the covariance; turbulent departures from the means were calculated using linear detrending (40); a double-axis rotation and tilt correction was applied according to ref. 41; a correction of the high-pass filtering effect was applied (42); low pass-filtering effect resulting from instrumental attenuations was corrected using different procedures depending on the setup: the analytical method (43) was adopted for the open-path systems (LI-7700); the in situ spectral correction method (44) was used for the closed path (FGGA-24EP) analyzers, more suitable to describe site dependent spectral attenuation along the tubing system. For the closed path and enclosed analyzers (45), the compensation for density fluctuations was not needed as mixing ratio data were measured and used for flux computations (44). For the open-path LI-7700, a spectroscopic correction was computed (46). Spectral correction for instrument separation was applied (47); data quality control (QA/QC) flagging was computed based on stationarity and integral turbulence tests (48, 49), resulting in three flags (0, good; 1, intermediate; 2, poor); random uncertainty attributable to sampling errors was finally estimated following (50). One-point CH₄ storage term was computed based on the concentration measurements of the gas analyzer (51).

To prevent possible biases from the heating system of the METEK in IVO, we modified the activation of the heating of the sonic anemometer to only activate when the data quality was low (as indicated by a quality flag), instead of when temperature was below a set temperature threshold, as commercially available. From September 2014, the CSAT-3D were also externally heated using Freezstop Regular heating cables (Heat Trace), operated at 12 V direct current (DC). These insulated heating lines were cut to length to cover the support arms of the anemometers and yet were far enough removed from the transducer mounting arms to minimize flow distortion and other contamination of wind data. Control of the heating elements was done using the CR3000 data logger and a normally closed solid-state DC relay. The data logger program activated a relay when the sonic transducers were blocked by snow and/or ice, as reported by the diagnostic output by the CSAT-3D. All data when the heating was active were removed for the heated anemometer. Additional data cleaning was performed accordingly to the following criteria: data were removed when the quality flags (48, 49) of H, LE, and CH₄ fluxes were 2; when the internal pressure of the LGR was ≤132 torr (corresponding to malfunctioning of the instrument); or when identified as an outlier (when exceeding a moving-window weekly 1 SD of the individual fluxes).

We used Los Gatos Research FGGA gas analyzers at all sites, except in IVO, where low power availability restricted use to the open-path LI-7700 analyzer. A second LI-7700 was also implemented in CMDL alongside the FGGA to ensure comparability of the results using these two instruments. To this end, in CMDL, we cross-compared the CH₄ fluxes estimated using these two instruments from October 1, 2013 to October 1, 2014 (Fig. S5). The difference between the LICOR and LGR-FGGA-24EP CH₄ fluxes was calculated directly and the resulting distribution used to elucidate features of their uncertainty (Fig. S5). Half-hourly fluxes were used to provide the largest possible sample size. The mean difference (LICOR F_CH4 – LGR-FGGA-24EP F_CH4; ΔF_CH4) was −0.0226 mg m⁻² h⁻¹. Considering that the overall mean values of the LICOR and LGR-FGGA-24EP fluxes were 0.216 and 0.239 mg m⁻² h⁻¹, respectively, this amounts to about a 10% difference between the two sensors, with the average LICOR flux slightly lower than the LGR-FGGA-24EP. The data were heteroscedastic, whereby the uncertainty estimated from Laplace distributions (Fig. S5) of data binned by flux magnitude increased with the value of the fluxes (52, 53). The average annual data coverage of the CH₄ fluxes (after removal of the data as indicated above) for the entire measuring period was 52% (CMDL LGR), 29% (CMDL LI-7700), 54% (BEO LGR), 58% (BES LGR), 49% (ATQ), and 38% (IVO LI-7700). The seasonal distribution of the data coverage is shown in Fig. S4.

Missing CH₄ flux data were gap-filled by applying artificial neural networks as described in refs. 54 and 55. The following meteorological variables, that most commonly act as drivers for CH₄ emissions, were used following the principle of parsimony and avoiding unnecessary input variables of a cross-correlative and cross-dependent nature: air temperature, air pressure, solar radiation, vapor pressure deficit, soil temperature, soil moisture, and the decomposed wind speed and direction, together with the fuzzy sets representing the annual seasons, as discussed by Dengel et al. (55). A training and a testing dataset for the neural networks were produced, including both daytime and nighttime periods in an equal manner covering all meteorological conditions from both years.

Several combinations of input variables and neurons were tested before we decided to include the same input variables and eight neurons for all datasets to keep the method uniform across all five sites. Each site analysis included 500 repetitions of which the 25 best [correlation coefficient values of the predicted (testing stage) values] runs were included in the gap filling itself. The application of error analysis, the agreement or disagreement between measured and modeled data helped to estimate the overall performance of the neural network models. The overall performance of all five models show no under- or overestimation of the true fluxes. Fig. S6 presents the gap-filled data (red) and the original measurements (gray), and Fig. S7 presents the influence of the gap filling on the daily averaged CH₄ fluxes. Final error analysis results can be found in Table S2, including the correlation of determination (R²), together with the root mean square error (RMSE) converted to true physical units of mg C-CH₄ m⁻² h⁻¹ indicating the uncertainties of the models.

For their central estimate in the CH₄ emission of Arctic tundra, McGuire et al. (8) used estimated wetland areas of 772,076; 7,540; 18,139; and 812,969 km² for North America, North Atlantic, Northern Europe, and Eurasia subregions, respectively, for a total of 1,610,724 km²; whereas we used the Circumpolar Arctic Vegetation map and included the land cover types B3, B4, G1, G2, G3, G4, P1, P2, S1, S2, W1, W2, and W3, so encompassing all tundra types excluding glaciers and lakes (38), which resulted in a total area of 5,070,000 km² across the entire Arctic. Our choice was justified by the importance of upland tundra (e.g., ATQ, CMDL, IVO). The land area we used was greater, but the measured fluxes used were lower than those used by McGuire et al. (8), who included data from boreal wetlands in Alaska (20, 56) in estimating rates of CH₄ emissions. Fortuitously, the differences in land area and fluxes compensated, resulting in similar annual estimates of CH₄ flux to those reported here.

The aircraft fluxes were at times lower than the mean of the EC tower fluxes, as has been observed previously in point-scale and regional-scale flux comparisons. The influence of the Brooks Mountain Range has been excluded (Fig. S3 C and F, gray shading). Additionally, we note that WRF estimates of planetary boundary layer (PBL) ventilation rates are difficult to assess quantitatively and might be subject to particular bias in the fall, when heat fluxes are low. A ∼28% difference of CO₂ eddy fluxes from aircraft compared with towers has been reported (57) and attributed in part to differences in the aircraft and tower footprints. As seen in Fig. S3, the footprint of the regional flux is much larger than that of the flux towers and includes areas assumed to be less productive particularly in autumn and winter, such as frozen lakes. Nevertheless, the regional CH₄ fluxes strongly support the view that our EC fluxes capture relevant cold season CH₄ dynamics and the response of CH₄ emissions to soil climate across the wider North Slope area.

Regional fluxes of CH₄ were estimated with aircraft data from the CARVE 2012 to 2014 (30). CH₄-, CO₂-, and CO-mixing ratios were measured using two independent cavity ring-down spectrometers: one operated wet (58) (G1301-m; Picarro) and one dry (30) (G2401-m; Picarro). Each analyzer was calibrated throughout the flight, ensuring a continuous 5-s time series. Ozone-mixing ratios were measured using a 2B Technologies model 205. The aircraft data were aggregated horizontally every 5 km and vertically in 50-m intervals below 1 km and 100-m intervals above 1 km. Each aggregated position was treated as a receptor in a Lagrangian transport model (WRF-STILT), which calculated the back trajectory of 500 particles from each receptor location. WRF-STILT represents the Stochastic Time-Inverted Lagrangian Transport (STILT) model coupled with meteorology fields from the polar variant of the Weather and Research Forecasting model (WRF) [v3.5.1 (59)]. The WRF-STILT calculation allowed for the quantification of the space and time where upstream surface fluxes influenced the measured mixing ratios. A total 24 h 2D surface influence field (i.e., footprint) was calculated for each flight (e.g., Fig. S3 C and F), representing the response of the receptor to a unit surface emission (ppb/mg C-CH₄ m⁻² h⁻¹) of CH₄ in each grid square (0.5° × 0.5° grid). The systematic uncertainty of the calculated surface influence is estimated at 10–20% (59). For comparison with the flux towers, aircraft data were carefully selected and a number of assumptions made to calculate a regional flux of CH₄. Only data collected north of 68° N, west of −153° W, with CO <150 ppb (to remove impacts of anthropogenic influence in the Deadhorse/Prudhoe Bay area), below 1,500 m above ground level and with over 60% surface influence from the North Slope were selected. CH₄ emissions from the higher altitudes of the Brooks Mountain Range were assumed to be negligible (gray area in Fig. S3 C and F). Assuming a uniform land surface emission, CH₄-mixing ratios should scale linearly with the total land surface influence observed at that receptor point. The flux of CH₄ for each flight day was calculated from the correlation of CH₄ with the STILT land surface influence [ordinary least-square (OLS) regression; Fig. S3 A and D], where the slope represents the regional flux and the intercept is the regional background CH₄-mixing ratio, which was assumed not to vary greatly during the flight. The ozone deposition velocity was also calculated in a similar manner (Fig. S3 B and E), and only flights with an ozone deposition velocity consistent with the expected seasonal cycle were used as a valid CH₄ flux.

To understand the environmental control on CH₄ fluxes over the year, weekly median CH₄ fluxes were modeled as a function of the weekly averaged environmental data. Environmental variables included were soil temperature at 0–5, 10–15, and 20–30 cm; soil water content (SWC) at 0–10 and 20–30 cm; and thaw depth. Only weeks with more than 20 flux and environmental measurements were included in the analysis. Because thaw-depth measurements were collected once a week in the Barrow sites, the weekly average of all available flux and environmental measurements were used in the statistical analyses. The remoteness of the IVO and ATQ sites limited the frequency of visits and therefore of thaw depth measurements. As a result, values were linearly interpolated between measurements during the summer. Because we measured a consistent thaw depth from the end of the summer until October 2014, we assumed a stable depth of the active layer from the end of summer throughout the zero curtain period, and we used the soil temperature profiles and the duration of the zero curtain to extrapolate the thaw depth until soil freezing. After soil freezing, we set the thaw depth to zero until the beginning of the following spring (again defined by the soil temperature profiles, as defined in Table S1). Because CH₄ fluxes presented a skewed distribution, they were log-transformed for all of the statistical analysis. Because of the data loss in the soil environmental data in 2013 and 2014 (e.g., the soil moisture data from IVO and soil temperature data in ATQ), we extended the dataset used for the statistical analysis to include data until May 2015, covering a full year for each site. Only the three sites with the most complete environmental datasets were used in the statistical analysis (ATQ, IVO, and BES).

Linear and nonlinear mixed effects models were used for this analysis (using the lme4 and the nlme in R; R Developing Core Team). The mixed-effects models included the “week” of measurement and “site” as continuous and categorical random effects, respectively, to account for the pseudoreplication and the different sites measured. The nonlinear mixed-effects model was used to test whether an exponential or power fit of soil temperature, and SWC were the best predictor of CH₄ fluxes. However, the increase in the complexity of the model did not justify the use of nonlinear mixed effects modeling, as shown by Akaike information criterion (AIC) values and the partial F test of the two models. Therefore, the results of the linear mixed effects models are reported here. The model performance was evaluated based on the AIC values, on the significance of the partial F test (used to compare two models), and on the marginal coefficient of determination (similar to the explanatory power of the linear models) for generalized mixed-effects models as output by the r.squaredGLMM function within the MuMIn package in R (60, 61).

The best univariate model explaining the variability in the CH₄ fluxes during the entire year was soil temperature (T) at 20–30 cm, with an AIC of 5.5 and a marginal coefficient of determination of 0.85. During the cold period (from September to May), soil water content at 20- to 30-cm depth was the most important variable explaining CH₄ fluxes, with an AIC of −27 and a marginal coefficient of determination of 0.89. During the summer period, instead, the best univariate model included thaw depth, presenting an AIC of 1.7 and a marginal coefficient of determination of 0.89. The best multivariate model for the CH₄ fluxes during the entire year included soil water content at 20–30 cm in addition to soil T at 20–30 cm, with an AIC of −23 and a marginal coefficient of determination of 0.89. This multivariate model was significantly different in its explanatory power from the univariate model that only included soil T at 20–30 cm (as shown by the significant partial F test of the difference in the two models). Similarly, during the cold season, the best multivariate model included soil T at 20–30 cm and soil water content at 20–30 cm, with an AIC of −55 and a coefficient of determination of 0.89 (which was significantly different from the univariate model as shown by the partial F test). No multivariate model was statistically different from the univariate one (which only included thaw depth) in the summer.

Seasonal patterns of aerial proportion (%) of surface water inundation within 25 × 25 km footprints extending over the greater Alaska domain is derived from K-band passive microwave satellite remote-sensing (62, 63). Surface inundation at the tower sites persists from late May following surface thaw through November, when colder air temperatures minimize the presence of liquid water above the soil surface. Wet subsurface soil conditions at ATQ and IVO contribute to peak summer CH₄ emissions of 20–50 mg C-CH₄ m⁻² d⁻¹, despite lower surface water inundation relative to the CMDL, BES, and BEO tower sites, indicating that substantial CH₄ emissions are not confined to wet, inundated tundra.