Controls on carbon dioxide and methane fluxes from a low-center polygonal peatland in the Mackenzie River Delta, Northwest Territories

The EC system measured fluxes of CO₂ () and CH₄ () over the study period. It consisted of a closed-path infrared CO₂/H₂O gas analyzer (IRGA, model LI-7200, LI-COR Inc., Lincoln, NE, USA; LI-COR), an open-path CH₄ analyzer (model LI-7700, LI-COR), and a CSAT3 sonic anemometer (Campbell Scientific Inc, Logan, UT, USA; CSI) mounted on a tripod at a measurement height (z_m) of 2.87 m (Fig. 1). The EC data and air pressure (P_a) were logged on a LI-7550 Analyzer Interface Unit (LI-COR). The CSAT3 was oriented to the northeast (35°) because northerly and easterly winds are most common for this time of year (NWT Water Resources 2019). Due to an error setting up the LI-7550, fluxes were only logged at 1 Hz until 12:00 DOY 193; the settings were corrected after this and the fluxes were logged at 10 Hz for the remainder of the study period. In principle, the low sampling frequency does not cause a bias in the mean flux, but does increase the variance (Bosveld and Beljaars 2001). This is supported by experimental observations (Rinne et al. 2008; Holl et al. 2019). To confirm this error had limited impact on our results we down sampled the July flux data for the period following the error correction (13–29 July). We artificially generated a 1 Hz data set by selecting every tenth value from the raw data, and processed it using the same steps as detailed below. Student’s t tests indicated no significant difference between the mean values of the 10 Hz data and the down sampled 1 Hz data for either CO₂ or CH₄ fluxes. Boxplots of the results 1 and 10 Hz fluxes are shown in Supplementary Fig. S1¹.

Half-hourly fluxes were calculated with EddyPro V.6.2.0 (LI-COR). The software performed statistical assessments (Vickers and Mahrt 1997), low and high frequency spectral corrections (Moncrieff et al. 1997 and 2004), and a double rotation (Wilczak et al. 2001). For the open path LI-7700, the Webb, Pearman, and Leuning (WPL) density correction was applied (Webb et al. 1980) with spectroscopic correction following McDermitt et al. (2011). For the closed path LI-7200 density corrections were calculated following Ibrom et al. (2007). Quality control flags (0–2) were assigned following Mauder and Foken (2004) and fluxes with a flag of 2 were discarded (4.6% and 10.8% of and data, respectively).

Post processing treatments were conducted as follows. (1) Removing LI-7700 observations during precipitation events (0.2% of data). (2) Removing LI-7700 observations when signal strength (received signal strength indicator, RSSI) from the LI-7700 was below 20% (13.3% of all data). Methane concentrations and were plotted against RSSI to confirm this was an acceptable threshold (Supplementary Fig. S2¹). Low RSSI values can result from condensation or dust on the mirror of the LI-7700. We elected not to use a cleaning pump at this site because of the limited power supply. This resulted in more missing data but extended the lifespan of the power supply for measuring fluxes. (3) Removing observations with mean wind direction from 215° ± 30° (8.3% and 7.1% of and data, respectively) to avoid uncertainties associated with the wake of the sonic anemometer. (4) Storage correcting by calculating net fluxes (NEE and NME) as the sum of the observed scalar flux ( and ) and the rate of change in their scalar concentrations at z_m following Aubinet et al. (2012). (5) Removal of spurious half-hourly measurements (4.1% and 3.3% of NEE and NME data, respectively) using the median absolute deviation about the median method following Papale et al. (2006). (6) Filtering fluxes by friction velocities (u_*) below an 0.1 m s–1 (4.4% and 2.3% of NEE and NME data). We confirmed that this value was suitable by calculating an iteratively determined u_* threshold (0.094 m s⁻¹) following Papale et al. (2006). In total, quality control and post processing removed 22% and 55% of and observations, respectively. Additionally, there was a 14 day gap in flux observations from 00:30 DOY 240 to 16:30 DOY 254 due to insufficient power supply from the solar panels. The CSAT3 and climate station remained operational, but the Li-7200 and Li-7700 were shut down.

The flux footprint was calculated following Kljun et al. (2015). All variables needed for footprint calculations were collected onsite, except planetary boundary layer height. Inverse distance weighting was used to interpolate these values from 3 hour reanalysis data for the 25 nearest grid points spanning the domain 68°30′N to 70°N and 134°W to 135°40′W. Half-hourly values were then approximated from the 3 hour value using linear interpolation. The flux footprint was calculated at a spatial resolution of 2 m² over a 2000 m × 2000 m grid centered on the tripod for all valid half-hourly measurement periods: u_* > 0.1 m s⁻¹ and mean wind direction ≠ 215° ± 30°. The half-hourly footprints were intersected with the landscape classification to determine the relative flux contribution of low-center polygons (F_Cnt), rims (F_Rim), and troughs (F_Tro) to each measurement. The contribution from outside the landscape classification (F_Out) was calculated as 1 – (F_Cnt + F_Rim + F_Tro). To allow for a direct comparison with the source area fractions, the LF values weighted by (1 – F_Out), using the median value of F_Out to calculate the adjusted landscape fraction (ALF).

Neural networks (NN) were used to identify relevant environmental controls over NEE and NME, gap fill the time series, and partition NEE into its component fluxes ER and GPP. Various studies have applied NNs to flux data to identify and analyze flux drivers (Moffat et al. 2010; Briegel et al. 2020; Skeeter et al. 2020), investigate the influence of spatial heterogeneity (Morin et al. 2014; Skeeter et al. 2020), and to gap-fill time series of NEE and NME (Papale and Valentini 2003; Dengel et al. 2013; Knox et al. 2016; Skeeter et al. 2020). The steps for identifying relevant controls, gap filling, and flux partitioning are discussed in the Feature identification and gap-filling section.

Climate data at the EC site were logged on a CR1000 datalogger every 1 s and averages/totals stored at 5 minute intervals. A NRLite net radiometer (Kipp & Zonen, Delft, Netherlands) measured net all-wave radiation (R_n) and a SQ-110 quantum sensor (Apogee Instruments, Logan, UT, USA) measured photosynthetic photon flux density (PPFD) at 3 m a.g.l, and a HMP35 (CSI) measured T_a and humidity (RH) in a shielded, naturally ventilated screen all at 2 m a.g.l. on the same tripod as the EC system. A tipping bucket rain gauge (R.M Young Company, Travers City, MI, USA) was mounted 2 m south of the tripod at a height of 0.4 m.

Soil temperature, moisture, and water table depth data were recorded at 30 minute intervals on three separate, automatic soil sites, each operated by a CR10x datalogger (CSI) at three locations near the tripod (Fig. 2a): one in polygon a center (S1), one on a polygon rim (S2), and one on a border between the rim and the center (S3). At each of the soil sites, soil temperatures (T_s) were recorded with custom made type-T thermocouples at depths of 2.5, 5, and 15 cm, volumetric water content (θ_w) was measured with CS616 water content reflectometers (CSI) at 0.1 m depth, and water table depth (W_td) was sampled with PLS probes (CSI). The climate and soil stations S1 and S2 operated uninterrupted from 23 June (DOY 175) until 13 September 2017 (DOY 257). On 17 July an animal dug up and destroyed the thermocouples and CS616 at S3. The PLS probes were installed at the base of the active layer and adjusted to account for increasing thaw depth on 10 July and 1 August. Thaw depth (TD) was measured during each site visit by inserting a graduated steel probe into the ground to point of refusal. Thaw depths were measured ten times for both polygon rims and centers and the median value calculated for each set of observations. Canopy heights were also measured for polygon rims and centers on each site visit using the median value of 10 replications.

The sun remains above the horizon much of the summer at 69°N. The date of the first sunset and timing subsequent sunrises and sunsets is important for understanding NEE and partitioning GPP and ER. We obtained the timing of sunrises and sunsets from a sunrise/sunset calculator (National Research Council Canada 2012). We used this information to one-hot encode a Daytime variable. Every half-hourly measurement period occurring before the first sunset or between sunrise and sunset was assigned a value of one. Measurements between sunset and sunrise were assigned a value of zero. If a sunrise or sunset occurred during a measurement period, the encoding was set to 1, as the purpose of this coding is to identify periods with no GPP.

In addition to EC measurements, ER was independently measured using the closed chamber method on two site visits (DOY 191, and DOY 233). Closed chamber measurements were sampled using a LI-700 (LI-COR) mounted in a custom made portable automated chamber as described by Christen et al. (2017). ER was measured at polygon centers (ER_Cnt) and rims (ER_Rim) using eight and four replications, respectively. Opaque PVC collars were inserted 10 cm into the ground on DOY 173. Before each measurement, chamber heights were measured with four repetitions per collar. Then the transparent PVC sensor head was placed over the collar and covered with and opaque canvas bag to block out all sunlight. The sampling period lasted 2 min and fluxes were then calculated following Christen et al. (2017). Due to limited time and resources, we were only able to collect a small number of samples and we were unable to collect chamber observations of GPP. The small sample size limits the conclusions we can draw from the chamber data, but this data are still valuable because they provide context to our EC observations and allow for comparison with other sites.

NN are flexible machine learning methods that make no prior assumptions about functional relationships within a data set and are ideally suited to perform non-linear, multivariate regression (Hornik 1991; Melesse and Hanley 2005; Desai et al. 2008). The goal of a NN is to approximate a target function as:

where t(X) is the target (e.g., NEE’s response to environmental drivers), ε(X) the noise, and is the set of input variables; M denotes the number of inputs variables, and x₀ = 1 is a bias term (Khosravi et al. 2011). In the context of EC measurements, this means that if all relevant climate and ecosystem information are available to a network, the remaining variability can be attributed to noise in the measurement (Moffat et al. 2010).

Here, we used feed-forward dense NN with a single hidden layer:

Where g(·) is a non-linear transfer function, here we used the rectified linear activation unit (ReLu) (Anders and Korn 1999). H denotes the number of hidden nodes in the network. The weights are randomly initialized and after each model iteration are updated by backpropagating the error through the network, β are the weights of the output layer and γ are the weights of the hidden layer. Weights are adjusted in the direction that decrease the error and training continues. We used early stopping to terminate training when the mean squared error (MSE) of a test data set failed to improve for 10 consecutive iterations. The test set consisted of 10% of the training data, and was not used if fitting f(X,w). Early stopping prevents f(X,w) from overfitting the training data (Weigend and Lebaron 1994; Sarle 1995; Tetko et al. 1995; Sarle 2014). When using early stopping, the choice of H is somewhat arbitrary, but H must be large enough to ensure the model has sufficient flexibility (Sarle 2014). We set H to 1/30^th the number of training samples.

To account for the random initialization of w and further minimize overfitting, stratified bootstrapping (with replacement) was used to generate a set of 30 bootstrapped training data sets (B). Separate f(X, w)_i were trained on each bootstrapped training set and for each model and validation metrics: root mean squared error (RMSE) and the coefficient of determination (r²) were calculated using the out of bag (OOB) bootstrapped samples. Thus, the final response model response was calculated as:where F(X) is the response of NEE or NME to a set of inputs (environmental drivers). The variance of the model outputs is:

A confidence interval (CI) for F(X) can be calculated as , where is the students t score, 1-α is the desired confidence level, and df are the degrees of freedom that are set to the number of bootstrapped samples B.

Our NN analysis used the Weights method (Gevrey et al. 2003) to quantify the relative influence of various inputs on NEE and NME and to prune the NN to reduce the number of input variables. This method has been used to study the influence of inputs for ecological modelling (Lee et al. 2003; Olden, et al. 2004; Fischer 2015; Liyanaarachchi et al. 2020), but to our knowledge, it has yet to be applied to EC flux data. For both NEE or NME, F(X) were initially trained on 21 initial inputs X including: radiation (R_n, PPFD, Daytime), atmospheric conditions (T_a, vapor pressure deficit (VPD), P_a), wind (U, u_*), flux contribution (F_Cnt, F_Rim, F_Tro), and sub-surface (T_s (all depths, polygon center and rim), VWC (polygon center and rim), water table depth (W_TD), thaw depth (TD)) variables. These X were chosen because they potentially influence NEE or NME.

After training, we calculated the partial derivatives (d_ji) of each input factor x_m in X for every sample (j=1, …, N), where N is the total number of flux observations:

where S_j is the derivative of the output with respect to the input, β_h and γ_mh are the weights of the output and hidden neurons, and I_hj is the response of the of the h^th hidden neuron for the input x_j,m. The max function is the first derivative of the ReLu activation function. d_ji were averaged over the bootstrapped data sets and plotted to visualize the influence of factors. Because the inputs and the target must be standardized (mean = 0, standard deviation = 1) before training a NN, d_ji are not in the units of NEE or NME, rather they can be interpreted in terms of relative magnitude (i.e., compared between inputs). When d_ji < 0 the respective input has a negative influence on the target and when d_ji > 0 it has a positive influence.

Next, the sum the of partial derivatives (SD_m) and sum of the squared partial derivates were calculated (SSD_m) for each input m:

The sign of the S_Dm indicates whether the variable has a net positive or negative influence on the output (over the input domain). The SS_Dm describes the relative magnitude of the variables influence (if d_ji changes sign, S_Dm is reduced but SS_Dm is not). The SS_Dm values for each input were then normalized by the sum of all SSD values to quantify the relative influence (RI) of each input factor. The 21-input model served as a benchmark to gauge the performance of a “pruned” model using only the inputs X with over 2.5% influence. Amiri et al. (2020) suggest removing inputs that have less than a 5% influence on the model; we used twice as many initial inputs as they did, so we chose 2.5%. The SS_Dm of the pruned models were calculated and used to gauge the relative influence of the flux drivers. Additionally, the partial derivatives for some the most important drivers were plotted to visualize their relative influence.

The B = 30 pruned models for both NEE and NME were then run on the time series of drivers over the full study period. The 30 bootstrapped model outputs were averaged and used to estimate mean NEE and NME by gap-filling the NEE and NME time series. Most drivers just had nearly complete records (>99%); brief gaps in the drivers were gap-filled with linear interpolation (up to 60 minutes). Wind speed and u_* had slightly more missing records (99% and 98% complete, respectively) and were also gap-filled with linear interpolation. The footprint drivers had more missing values (95% complete), F_Cnt, F_Rim, and F_Tro were gap-filled using their respective mean values, binned by wind direction in 10° intervals. Confidence intervals for NEE and NME were calculated using the total variance of multiple imputations, which is the sum of variance between and within imputations.

Typically, NEE is gap-filled using flux-partitioning algorithms that model ER using T_s or T_a and GPP using PPFD (e.g., Aubinet et al. 2012; Lee et al. 2017) where NEE = ER − GPP. However, these methods require nighttime observations to estimate ER and, thus, do not perform well for Arctic summertime measurements due to the limited number of samples available during low light conditions (Kutzbach et al. 2007; Runkle et al. 2013). We tested a Q10 T_a response curve (eq. 9.4 in Aubinet 2012) and a logistic T_s response curve (eq. 1 in Lee et al. 2017) using each of the six T_s measurements to gauge their effectiveness when trained on the limited number (n = 186) of ER samples.

A NN trained on CO₂ fluxes will just estimate NEE rather than the component fluxes ER and GPP. To compensate for this, we estimated daytime ER with the pruned NN by using an artificial “dark” input set (Skeeter et al. 2020). The inputs Daytime and PPFD were both set to zero and any R_n > 0 were also set to zero. This approximation projects ER well beyond conditions under which it was actually observed, so caution needs to be taken when interpreting the output. However, the NN performed better than the traditional methods when estimating nighttime ER (see Ecosystem respiration section). Therefore, we feel this is the best option for estimating ER and is useful for comparing to the ER chamber observations. The ER estimate and confidence interval were obtained by “gap-filling” the sparse record of nighttime NEE (n = 186) observations using the same procedure as for NEE and NME.

The model trained on NEE performed very well on the OOB validation data (r² = 0.93, RMSE =0.39 μmol·m⁻² s⁻¹) and performance only dropped slightly compared with the benchmark (r² = 0.95). Pruning selected 8 inputs for NEE, which are listed in Table 2. PPFD is the dominant driver of GPP, and is the primary control over NEE (64% RI). The partial first derivatives of PPFD (Fig. 6a) show a strong negative influence that begins to decrease in magnitude at 100 μmol·m⁻² s⁻¹. The function reaches a minimum around 750 μmol·m⁻² s⁻¹ (dependent upon the other drivers), above which it has a weak positive influence. The PPFD response is modulated by the other drivers, primarily VPD (8% RI). It had a strong negative effect over NEE at low VPD, that reached an “optimal” value around 340 Pa (Fig. 6b). At higher VPD, it was a strong limiting factor over GPP. The response curve of NEE to PPFD shows that all else being equal, high VPD (1000 Pa) can limit CO₂ uptake by more than 0.5 μmol·m⁻² s⁻¹ relative to the optimal value (Fig. 7a).

Thaw depth (7% RI) gave the model a seasonal signal that influenced both GPP and ER. Low thaw depths (early season) had a negative influence, NEE trended more negative as vegetation matured (Fig. 6c). The model response shows a minimum around 0.35 m, below this point, all else being equal, NEE trends positive. Figure 7b shows the model response to TD for “typical” daytime and nighttime conditions. During daytime, maximum CO₂ uptake can be seen around the timing of the minima identified in Fig. 6c. With radiative inputs fixed to “dark” values, ER decreases steadily from a maximum during the early season, then leveled out around 0.4 m (corresponds to DOY 226) at about 1.0 μmol·m⁻² s⁻¹. We had limited nighttime NEE observations, and none before DOY 206, so the confidence interval around the modeled ER was broader.

The model identified soil temperature at 5 cm (T_Cnt5 RI 7%) as the main driver of ER. Above 3 °C, increasing T_Cnt5 consistently increased ER (Fig. 6d). This was modulated by T_Rim5 (5% RI) and T_Rim15 (3% RI) that had moderately negative influences at lower temperatures (Fig. 6e and f). At higher temperatures the negative T_Rim5 influence was negligible and T_Rim15 became a positive driver. Polygon center and rim temperatures at 5 cm were highly correlated (r = 0.81) and moderately to weakly correlated with T_Rim15 (T_Cnt5 r = 0.58, T_Rim5 = 0.19). Rims were consistently warmer than centers, but the magnitude of the temperature gradient (T_Δ5 = T_Rim5 – T_Cnt5) decreased through the study (Fig. 3a). Daily maximum T_Δ5 was in the late afternoon and minimum was in the late night/early morning. The response of NEE to typical late afternoon and late night T_Δ5 is shown in Fig. 7c. There was a strong response to T_Cnt at higher temperatures, but the response was muted or even reversed at lower values. Temperatures at 15 cm also had a diurnal signal, but they were much more muted. In the model, T_Rim15 had more impact on the seasonality of NEE, with higher temperatures at depth indicating increased respiration throughout the soil profile.

Mean estimated ER was 1.54 (CI_95% ± 0.87) μmol·m⁻² s⁻¹, which suggests a mean GPP around 2.14 μmol·m⁻² s⁻¹. This is an extrapolation well beyond conditions the model was trained on (see Fig. 8b), so there is more uncertainty around this prediction. The nighttime ER validation statistics for the model (r² = 0.58, RMSE = 0.23 μmol·m⁻² s⁻¹) show this method outperformed the Q₁₀ response curve (r² = 0.32, RMSE = 0.28 μmol·m⁻² s⁻¹) and the best logistic soil temperature curve fit (T_Rim2.5: r² = 0.45, RMSE = 0.26 μmol·m⁻² s⁻¹). The NN does a better job of approximating ER, especially at higher and lower NEE (Fig. 8a). The chamber data gave us increased confidence in our NN derived ER estimate (Fig. 8b). Due to the limited sample size, here, we report the median and inner quartile range (IQR). Median ER_Rim over the two collection days was 2.77 (IQR 2.33 – 4.05) μmol·m⁻² s⁻¹ and was more than double ER_Cnt 1.09 (IQR 0.90– 1.39) μmol·m⁻² s⁻¹. Median values between sites also decreased from 1.57 μmol·m⁻² s⁻¹ n DOY 191 to 1.08 μmol·m⁻² s⁻¹ n DOY 233. Our NN derived ER estimate is between median ER_Cnt and ER_Rim and also showed a decreasing seasonal pattern.

In low-center polygonal landscapes, polygon centers are depressed relative to rims. Relief was about 10–20 cm at Fish Island, whereas a site on Samoylov Island in Siberia’s Lena Delta had up to 50 cm (Sachs et al. 2010). Anaerobic respiration is favored in the depressed areas, promoting CH₄ production, and inhibiting ER and CH₄ consumption (Lai 2009). Our chamber ER data show higher ER from rims than centers at Fish Island. Similar patterns have been found in chamber ER studies of low-center polygon sites near Barrow (Utqiaġvik), AK (Olivas et al. 2011) and Samoylov Island, Siberia (Eckhardt et al. 2019). At Samoylov Island, NME differed by an order of magnitude between centers and rims as well (Sachs et al. 2010). However, rim GPP was higher and lower than center GPP at the Alaskan and Siberian sites, respectively (Olivas et al. 2011; Eckhardt et al. 2019). The rims at Fish Island had more shrub cover than either of the other sites, so we cannot assume the GPP response here would be similar to either the Alaskan or Siberian sites.

At the plot scale, there was significant spatial heterogeneity at Fish Island. The NN analysis indicated polygon rim and center fraction were both relevant for NME, with a combined influence of 6%. Polygon rims had lower NME than centers and both had lower NME than troughs. The source area fractions are not drivers of landscape scale CH₄ per se. Rather they highlight substantial spatial heterogeneity in CH₄ fluxes at Fish Island and indicate location bias had a minor impact on observed NME (Schmid and Lloyd 1999). Following a three-tier site-footprint representative index based on the dominant landcover fraction and chi-square analysis, our site’s footprint would be classified as Medium (Chu et al. 2021). Projecting to the ALF shown in Table 1 estimates that landscape scale NME is about 1 nmol·m⁻² s⁻¹ lower. This projection can also be used to infer that if permafrost degradation leads to the development of more troughs at Fish Island in the future, landscape scale NME will go up.

For NEE the temperature gradient between rims and centers was relevant, which indicates that microtopography has an effect on landscape scale CO₂ fluxes at Fish Island. Footprint source area inputs were not identified as important drivers (RI > 2.5%) of NEE by the benchmark model, so they were pruned. In the benchmark NEE model, polygon rim (RI = 1.9%) and center (RI = 1.7%) ranked 12^th and 14^th among the initial 21 inputs, indicating they had little influence on footprint scale NEE compared with the selected drivers. At very low values (e.g., F_Rim < 10%), the benchmark model predicts more negative NEE, which supports the chamber observations of reduced ER in polygon centers, but conditions like that accounted for less than 4% of our flux observations. For further comparison, trough fraction ranked 21^st (RI = 0.5%) and 9^th (RI = 2.4%) in the benchmark NEE and NME models, respectively. The low importance assigned to the source area fractions suggests the placement of the tripod relative to microtopographic features did not have an appreciable impact on the NEE observed.

The resolution of the landscape classification and footprint function (Kljun et al. 2015) may have been limiting factors. However, the average polygon was only 270 m² (∼16.5 m across), whereas the 50% and 90% contours for the typical flux footprint were 2300 m² and 6500 m², respectively. Given the relatively small size of the polygons and the regularity with which they repeat across space, the spatial heterogeneity at Fish Island mostly averages out at the footprint scale.

The NN analysis identified PPFD as the primary control over GPP, and highlighted VPD as the main limiting factor of GPP, which was expected (Aubinet et al. 2012). This relationship has been identified by NN analysis in other studies (Moffat et al. 2010; Skeeter et al. 2020) and has been noted at other wet tundra sites across the Arctic (Kwon et al. 2006; Fox et al. 2008). The model selected polygon center temperatures at 5 cm as the dominant driver, which makes sense given that polygon centers cover the majority of the land area (66%) at the site. It also indicated that the center–rim temperature gradient played an important role in modulating ER. Other studies have found surface temperature rather than soil temperatures to be the primary control over ER in tundra environments (Kutzbach et al. 2007; Runkle et al. 2013). We lacked an observation of surface temperature and would consider adding it in a future study. However, temperatures throughout the upper layers of peat were shown to be strong predictors of ER at a polygonal peatland site in Salluit, Quebec (Gangnon et al. 2017). Further, given the flexibility of NN to account for interactions between variables we are confident this provides a better estimate than soil surface temperature alone could.

When modeling NEE, the network does not inherently differentiate between ER and GPP. Rather it conflates the signals of the two responses into one output. Some studies have trained the NN separately on nighttime and daytime conditions to resolve ER (Papale and Valentini 2003). This approach was not well suited for our site. Instead, we used the Daytime variable to give the model the ability to differentiate between night and day without splitting our training data. Skeeter et al. (2020) approximated ER using a NN by setting radiative inputs to “nighttime” values, but they did not use a Daytime variable, and their NN underperformed relative to a Q₁₀ curve. We found adding Daytime drastically increased the model’s ability to resolve ER, which we support with comparisons to the Q₁₀ and logistic temperature response curves. It also allowed us to estimate ER over the full study period, despite the limited number of nighttime observations, using an artificial “dark” data set.

Our NN derived ER (1.54 μmol·m⁻² s⁻¹) was between the ER_Cnt and ER_Rim, which gives us confidence in our estimate. Using the half-hourly footprint source area fractions, we can spatially weight our chamber observations. Assuming polygon center/rim distribution outside the LCM mirrors the distribution within it, and ER from troughs is between ER_Cnt and ER_Rim we get an estimate of 1.62 (IQR1.34 to 2.22) μmol·m⁻² s⁻¹. This is a rough, back of the envelope comparison, but it does further support using this method to estimate ER.

Mean gap-filled NEE, estimated ER and GPP at Fish Island all exceeded ranges observed at a low-center polygon sites in Alaska (Olivas et al. 2011) and Siberia (Eckhardt et al. 2019) using flux chambers. We attribute the differences in GPP to (1) The lower latitude of Fish Island (69°) compared with the Alaskan (71°) and Siberian (72°) sites. (2) More productive vegetation at Fish Island, e.g., Salix spp. were absent at the Alaskan and Siberian sites. Differences in ER are attributable to colder soil temperatures and shallower thaw depths (25 cm) and (35 cm) at the Alaskan and Siberian sites respectively (Olivas et al. 2011; Eckhardt et al. 2019).

There was a significant diurnal cycle in NME throughout the season, and this was identified by the NN selecting R_n as the dominant driver. Many studies in Arctic, subarctic, and peatland sites have not found a distinct diurnal cycles in CH₄ fluxes (Rinne et al. 2007; Sachs et al. 2008; Nadeau et al. 2013; Lee et al. 2017; Rößger et al. 2019b). However, mid-day maxima in NME have been observed at subarctic fens in the Hudson Bay lowlands (Chanton et al. 1992) and Siberia (Veretennikova and Dyukarev 2017). In the Siberian fen, mid-day NME was three times the nighttime values (Veretennikova and Dyukarev 2017), which is even greater than observed at Fish Island. They proposed diurnal temperature variations explain the cycle. We postulate the same is occurring at Fish Island, given the model both polygon center temperatures at 15 cm (T_Cnt15) as an auxiliar drivers and the diurnal cycle in T_Cnt15. However, we cannot overlook the potential role of plant-mediated transport as another driving factor (Wang and Han 2005). Given the strong positive correlation between NEE and NME, we think this merits further investigation as well.

W_TD is known to be a primary determinant of CH₄ production and consumption in peatlands (Kutzbach et al. 2007; Lai 2009). However, Arctic sites with little variation in W_TD have found its influence to be negligible (Sachs et al. 2008; Rößger et al. 2019b). At Fish Island, decreasing water table depth had a weak negative influence, which was unexpected. We propose three explanations for this. (1) The range of water table depth was fairly small over the study period. This reasoning was given to explain a similar pattern at a bog in the James Bay lowlands (Nadeau et al. 2013). (2) Data during/after precipitation events when the water table rose rapidly were generally filtered out due to low signal strength (RSSI). Using the unfiltered data, we see strong peaks associated with rising water table, but these signals are unreliable so more studies would be needed to verify this. Gagnon et al. (2017) found precipitation events created conditions that suppressed ER in a polygonal peatland. It is probable that the same thing occurs at Fish Island and leads to elevated NME. (3) Significant loss of data in the second half of the study period means we lack NME observations for many specific combinations of thaw depths and water table levels, limiting the model’s ability to resolve the relationship between these drivers. To adequately resolve the influence of water table depth and its interaction with thaw depth, a multi-year data set would be needed.

Perhaps the most interesting finding of our NN analysis is the complex relationship of NME with wind speed (U) and friction velocity (u_*). Spikes in NME associated with high friction velocity have been observed at a low-center polygon site in Siberia (Sachs et al. 2008) and a peatland in the James Bay lowlands (Nadeau et al. 2013). One explanation is that gas transfer between open water and the atmosphere increases proportionally to the third power of the windspeed (Sachs et al. 2010; Wanninkhof and McGillis 1999). This likely had some effect because of the ponded troughs within the footprint. However, most of the land area was consistently above the water table. Instead, we frame the hypothesis that high winds and turbulence may enhance pressure pumping (Laemmel et al. 2017; Mohr et al. 2016) that ventilates the aerobic peat layer in the polygon center. We hypothesize that pressure pumping reduces the time CH₄ is subject to methanotrophy within the aerobic layer of peat and, therefore, enhances net CH₄ emissions by transporting it more quickly from the anaerobic layer to the atmosphere. Given that high wind events ventilate CH₄ that would otherwise be consumed, a change in storm frequency paired with an increase in the duration of the snow free season could lead to an increase in emissions from sites like Fish Island, a hypothesis that requires further exploration.