2.1 Study site and data collection

The study took place at Illisarvik, a DTLB on Richards Island (69^∘28^′47.5^′′ N, 134^∘35^′18.7^′′ W), which was drained experimentally in 1978 (Mackay, 1997). Illisarvik has since served as the focus of studies on permafrost growth, active layer development, and vegetation succession (Ovenden, 1986; Mackay and Burn, 2002; O'Neil and Burn, 2012; Wilson et al., 2019). At the nearby Tuktoyaktuk climate station mean annual air temperature (T_a) is −10.1 ^∘C, July is the warmest month with a mean of 11 ^∘C, and January is the coldest at −27 ^∘C. Mean annual precipitation is 160.7 mm yr⁻¹, the majority falling as rain in the summer and autumn. Snow cover typically lasts from mid-September or early October to late May (Environment Canada, 2016). Tuktoyaktuk is 60 km east of Illisarvik and in similar proximity to the coast so the climatology is expected to be similar at Illisarvik.

Table 1Dominant species or landscape feature within the vegetation/cover classes. Unit codes correspond to the map in Fig. 1a.

Download Print Version | Download XLSX

In the 39 years since drainage, Illisarvik has undergone rapid vegetation succession. After drainage, there were two remnant ponds. In the first 5 years after drainage, vegetation colonized the basin margins and wetter areas (Ovenden, 1986). By 1999, low vegetation had proliferated across most of the basin and taller willows had become established along the basin margins (Mackay and Burn; 2002). By 2010, some of the willows had grown to be 3 m in height (O'Neil and Burn; 2012). Current vegetation at Illisarvik is diverse relative to the dwarf-shrub tundra of the surrounding uplands (Table 1); the basin hosts a mix of woody shrubs (Salix spp., Betula spp., and Alnus spp.), wetland vegetation (Carex aquatilis, Arctophila fulva, etc.), and various grasses (Pocacea spp.) (Wilson et al., 2019). The basin is partly ringed by a terrace of peat that formed after a partial drainage event ∼5000 years BP and supports vegetation similar to the uplands (Michel et al., 1989). An ancient DTLB is located 100 m to the south of the Illisarvik basin, and the Arctic Ocean is to the west of the basin, separated by a ridge of upland tundra about 50 m wide at its narrowest (Fig. 1).

Figure 1(a) Map of the distribution of vegetation classes at Illisarvik, with the footprint climatology (F_Clim) over the study period, the locations of the chambers, and the eddy covariance (EC) system. The alphanumeric labels correspond to the unit codes in Table 1. (b) Legend for the map in (a). (c) Oblique drone image of Illisarvik, taken at the 16:40 23 July 2016 view from E of DTLB towards W. The basin and EC system are shown on the image using the same symbology as (a).

A vegetation survey of species composition and abundance was done on a 50 m grid in and around the basin during the 2016 study period (Wilson et al., 2019). A vegetation map was created with 10 units based on plant functional type and vegetation structure, with sub-units denoting sub-canopy vegetation. The unit boundaries between grid points were estimated visually by traversing the grid lines. Additional survey data on vegetation units and canopy height were collected manually with a GPS in the proximity of the EC station because greater resolution was needed for footprint modelling. Aerial imagery was collected on 23 July over two flights using a Phantom 2 drone (DJI, Shenzhen, China). The GPS points and drone imagery were used to cross-reference and modify the map of Wilson et al. (2019). The 10 units were then aggregated into six broader surface cover classes (listed from largest to smallest areal fraction within the footprint climatology (F_Clim); see Sect. 2.3 for definition): shrub, grass, sedge, upland, sparse, and water (Fig. 1 and Table 1).

2.2 Weather and soil measurements

Weather data were logged on a CR1000 data logger (Campbell Scientific Inc, Logan, UT, USA; CSI) at 5 min intervals. Net all-wave radiation (R_n) and photosynthetic photon flux density (PPFD) were measured with a NR Lite net radiometer (Kipp and Zonen, Delft, Netherlands) and a SQ-110 quantum sensor (Apogee Instruments, Logan, UT, USA), respectively, 3.2 m above the grass surface on the main EC system tripod (Fig. 1). A shielded HMP35 (CSI) recorded T_a and relative humidity (RH) 2 m above the surface. A tipping bucket rain gauge (R.M Young Company, Traverse City, MI, USA) was placed 3 m to the west of the main tripod. Soil temperature and moisture were measured within soil pits in two different vegetation types near the tripod: grass (30 m to the east) and shrub (40 m to the north). Measurements were made of ground heat flux (G) with custom-made heat flux plates, soil temperatures (T_s) with custom type-T thermocouples at depths of 0.08 m, and 0–20 cm integrated volumetric water content (VWC) with CS616 water content reflectometers (CSI). The soil measurements were recorded at 30 min intervals on CR10x data loggers (CSI). The climate and soil stations operated uninterrupted from 10 July (day 192) and 11 July (day 193), respectively, until 7 August 2016 (day 220). On 11 July and 6 August thaw depth was measured at each of the 10 chamber sites (see below). Thaw depth was measured by inserting a graduated steel probe into the ground to the point of refusal. Each site was probed five times: the median value has been used as the thaw depth at each location. On 12 and 15 July, a large herd of reindeer (∼500 animals) visited Illisarvik. They mostly avoided the tripod but did graze near it for about an hour on 12 July, which may have affected greenhouse gas fluxes.

2.3 EC fluxes

An EC system was placed in the southwestern portion of the basin (69^∘28^′47.82^′′, −134^∘35^′18.6^′′) and measured fluxes of CO₂ ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) and CH₄ ($F_{{\mathrm{CH}}_{\mathrm{4}}}$) for the full study period between 10 July and 7 August 2016. The EC system consisted of an open-path infrared CO₂∕H₂O gas analyzer (IRGA) (model LI-7500, LI-COR Inc., Lincoln, NE, USA; LI-COR), an open-path CH₄ analyzer (model LI-7700, LI-COR), and a CSAT3 sonic anemometer (CSI) mounted on a tripod at a measurement height (z_m) of 3 m (Fig. 2). The EC data and air pressure (P_a) were logged at 10 Hz on the LI-7550 analyzer interface unit (LI-COR). The CSAT3 was oriented to the northeast (40^∘) because climatology for Tuktoyaktuk indicated northerly and easterly winds are typical for July and August (Environment Canada, 2016).

Figure 2(a) Half-hourly air and soil temperatures, displayed along with photosynthetic photon flux density (PPFD). (b) Hourly soil volumetric water content and daily total precipitation. (c) Half-hourly $F_{{\mathrm{CO}}_{\mathrm{2}}}$ (green) and NEE_NN (grey) and (d) half-hourly $F_{{\mathrm{CH}}_{\mathrm{4}}}$ (red) and NME_NN (grey).

Download

Half-hourly fluxes were calculated with EddyPro v.6.2.0 (LI-COR). The software performed statistical assessments (Vickers and Mart, 1997), performed low- and high-frequency spectral corrections (Moncrieff et al., 1997 and 2004), performed a double rotation (Wilczak et al., 2001), applied the WPL correction to account for density fluctuations (Webb et al., 1980), and computed quality control (qc) flags (Mauder and Foken, 2004). Post-processing treatments included storage correction (calculating the net flux as the sum of the observed scalar flux and the rate of change in scalar concentration at z_m), filtering fluxes by friction velocities (u_∗) below 0.1 m s⁻¹, removing qc flags = 2 (Mauder and Foken, 2004), and the mean absolute deviation spike removal algorithm (Papale et al., 2006). Additionally, observations with mean winds of 220±30^∘ were removed to avoid uncertainties associated with the wake of the sonic anemometer, and observations were removed during precipitation events and when the open-path analyzers indicated there were any other obstructions within the path (Aubinet et al., 2012). The data were gap-filled using neural networks (NNs) which have been applied to $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ in other studies (Moffat et al., 2010; Dengel et al., 2013). Details of the NN methodology are described in Appendix A.

The flux footprint represents the influence of upwind areas on a measured scalar flux, and the footprint climatology is the average of individual footprints over a time period. Evaluation of the flux footprints and climatology helps evaluate the reliability of the dataset and estimate the source area of each individual half-hourly EC flux measurement. A scalar flux F_c sampled at $(\mathrm{0},\mathrm{0},z_{\mathrm{m}})$, where z_m is the height of the EC instrumentation, can be represented as the integral of the flux footprint function f(x,y) and the distribution of sources–sinks (Q_c) over a domain D (Kljun et al., 2015):

The flux contribution of upwind source areas increases sharply upwind from the measurement location to a peak and then decreases gradually with increasing distance (Schmid, 2002). The empirically derived flux footprint function of Kljun et al. (2015) was used to estimate the source area of each half-hourly flux measurement.

The model requires boundary layer heights which were not measured on site. Half-hourly boundary layer heights were interpolated from 3 h estimates obtained from the Global Data Assimilation System of the US National Oceanic and Atmospheric Administration. The model also requires the aerodynamic roughness length (z₀), which is influenced by the canopy height and spacing. Canopy height (C_h) varied considerably within the basin (from >1 m in the north to ∼0 m in the bare-ground areas). Canopy height variability was lower in the vicinity of the EC tripod but ranged from 0.35 to 0.55 m with a few taller shrubs approaching 1 m. Median z₀ was calculated for 30^∘ wind sectors following Paul-Limoges et al. (2013). This calculation was performed for near-neutral conditions: $-\mathrm{0.05}\le \frac{z_{\mathrm{m}}}{L}\le \mathrm{0.05}$, where L is the Obukhov length. The z₀ for each wind sector was found to be insensitive to zero-plane displacement height, d, as z_m≫d, so the mean value of d around the tripod was used, where $d=\mathrm{2}/\mathrm{3}{C}_{\mathrm{h}}$. Zero-plane displacement did not change significantly over the course of the study so z₀ remained fixed over the study period for each wind sector.

For each half-hourly flux observation, f(x,y)_i was solved at 1 m² resolution over a 1 km² domain centred on the EC tripod. Then, f(x,y)_i values were intersected with the surface classes to determine the relative contribution of each surface type to each flux observation (referred to as F_Shrub, F_Sedge, etc.). The footprint function is technically infinite so a fraction of each f(x,y)_i was not contained within the model domain. The out-of-domain source fraction ranged from 1.8 % to 4.9 % with a mean of 3.2 % and was assumed to have minimal impact on the analysis. The flux footprint climatology (F_Clim) was calculated by averaging the half-hourly flux footprints over the study period and is shown in Fig. 1. Table 2 shows the flux contribution of each vegetation class.

Table 2The surface cover class fractions of the basin, along with the mean source area fractions of the footprint climatology (F_Clim) and the range of source area fractions for individual half-hourly observations shown in brackets.

Download Print Version | Download XLSX

2.4 Closed chamber measurements

In addition to EC measurements, fluxes of CO₂ and CH₄ were sampled using a static non-steady-state chamber flux technique on 11 dates between 12 July and 5 August 2016 (Laforce, 2018). Nineteen chamber collars were located at 10 sites, eight sites within and two outside the basin (Fig. 1). Each surface cover class was represented by at least one chamber site, except for open water. At each vegetated site a pair of collars were installed 20 cm apart, except at the “sparse” site where only one collar was installed. The above-ground biomass was removed from one of the collars at each vegetated site. There were three replicates (six collars) for the shrub class; two for the sedge; grass, and upland tundra; and no replicates for the Sparse class. PVC collars 30 cm long and 24.3 cm in diameter were inserted to a depth of approximately 15 cm. The chambers were 34 cm tall and made out of polycarbonate covered in black opaque tape to maintain dark conditions inside the chamber (for more details, see Martin et al., 2018). The chambers contained a small vent (10 cm coiled 1/8 in. diameter copper pipe) to ensure a constant pressure during measurements. The opaque chamber fluxes of CO₂ provided an independent estimation of ER. This helped characterize ER given the challenges with standard NEE partitioning techniques at high-latitude sites during the Arctic summer as noted in Sect. 2.5.1.

Chamber flux measurements were made between 09:00 and 17:00 starting at a different collar set each day to randomize the sampling order to avoid a bias due to diurnal patterns. During gas flux measurements, the chambers were sealed to the top of the collars within a groove filled with water, and five 24 mL air samples were collected into evacuated 12 mL vials sealed with doubled septa. Each vial contained a small amount of magnesium perchlorate to dry the air sample. Samples were collected at 0, 5, 10, 15, and 20 min after the chambers were set on the collars. Air within the chamber was mixed with a 60 mL syringe attached to a three-way stopcock before each air sample was taken. Samples were stored until analysis 1 month later at Carleton University. The integrity of the vials through shipping, storage, and analysis was confirmed using a subset filled with helium before the field season began.

Concentrations of CO₂, CH₄, and N_2O were determined using a CP 3800 gas chromatograph (Varian Inc., Palo Alto, CA, USA) as described by Wilson and Humphreys (2010). Three replicates of five CO₂∕CH₄ standards varying from 383.1 to 15 212.6 ppm CO₂ and from 1.08 to 22.11 ppm CH₄ were included in every set of measurements to create a linear relationship between gas concentration and chromatogram area. The chamber fluxes of CO₂ and CH₄ (F_C) were calculated as follows:

where (dc∕dt) is the linear rate of change in the mixing ratio of the gas, A is the chamber area (0.0464 m³), V is the chamber volume (between 0.0182 and 0.0242 m³ adjusted for collar depth at each collar location), R is the ideal gas constant, P is pressure in pascals, and T is the air temperature in kelvin. P and T values corresponding to the time of each measurement were obtained from the EC station. Visual inspection of the linear trend of gas concentrations (dc∕dt) was used to identify and remove spurious point measurements associated with analysis errors, leaking chambers (isolated decreases in concentration), and contamination or ebullition events (isolated increases in concentration) (0.3 %, 0.7 %, and 2.0 % of CO₂ samples and 2.1 %, 0.5 %, and 1.1 % of CH₄ samples, respectively). In all flux measurements, at least three or more gas samples remained so that dc∕dt and its coefficient of determination (R²) were determined using least-squares linear regression. We did not use R² as an additional quality control criterion as many of our CH₄ fluxes were near zero and tended to have low R² values due to only small variations in the point sample concentrations (see also Clark et al., 2020). A total of 40 % and 32 % of the 227 CH₄ flux measurements and 97 % and 92 % of the 227 CO₂ flux measurements had R² over 0.80 and 0.90, respectively. No flux measurements were removed from the analysis. Positive fluxes indicate emissions of gases to the atmosphere, and negative fluxes indicate uptake by the surface.

2.4.1 Upscaling

Chamber fluxes of ER were upscaled from the plot scale (individual chamber) to the footprint scale using the footprint-weighted average method and to the basin scale using the area-weighted average method (Budishchev et al., 2014). The chamber ER and air temperature from the EC tripod (T_a) were used to determine R₁₀, the base respiration at 10 ^∘C, and Q₁₀, the temperature sensitivity coefficient, using Eq. (3) for five of the six surface classes (Fig. 1) (Laforce, 2018) (Table 3).

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\mathrm{ER}=R_{\mathrm{10}}{Q}_{\mathrm{10}}^{\frac{(T_{\mathrm{a}}-\mathrm{10})}{\mathrm{10}}}\end{array}$$

Half-hourly footprint-scale estimates (ER_FS) were calculated by multiplying ER derived from Eq. (3) for each surface class by the footprint source area fraction and summing over classes. Basin-scale estimates (ER_BS) were estimated the same way but using the mean source area fractions of the basin (Table 2). As there were no open water class ER estimates, ER from open water was assumed to be zero.

Table 3The ER temperature sensitivity (Q₁₀) and base respiration (R₁₀) estimated by Laforce (2018) and estimated from nighttime EC footprint observations.

Download Print Version | Download XLSX

In contrast to ER, there are no standard empirical functions to estimate temporal variations in NME. Instead, we used ordinary least-squares regression (OLS) to estimate NME. The most important environmental controls over $F_{{\mathrm{CH}}_{\mathrm{4}}}$ were VWC and T_s (discussed below). Continuous observations of these factors at the flux chambers were not available; instead chamber NME values were grouped by vegetation class and fit to VWC and T_s measured in the soil pits near the EC station. Half-hourly footprint-scale (NME_FS) and basin-scale (NME_BS) estimates were then made using the OLS parameters for each surface class using the same procedures for ER_FS and ER_BS.

2.5 Factor selection and gap filling

We used an exploratory approach to identify the smallest set of factors that best predicted half-hourly EC-derived NEE and NME without overfitting the dataset using a series of neural networks (NNs). We started with 10 factors: four meteorological variables (PPFD, T_a, vapour pressure deficit (VPD) computed using the T_a and relative humidity (RH) data, and three-dimensional wind speed (U) measured using the CSAT3 sonic anemometer), two soil variables (VWC and T_s averaged between the two soil pits near the EC tripod), and four source area fractions (shrub, F_Shrub; grass, F_Shrub; sedge F_Sedge; and upland, F_Upland). The four source area variables correspond to surface classes sampled by the chambers. We excluded water (F_Water) and sparse (F_Sparse) fractions because its average contribution to the EC observations was only 0.2 % and 2.2 %, respectively, and there were no chamber measurements for the water class while chamber measurements indicated ER was low and NME was not significantly different from zero for the sparse class. A number of these prediction factors were highly correlated, but it was necessary to include them so the model could account for source area heterogeneity.

The NNs were trained iteratively on bootstrapped datasets. First NNs were trained on each factor individually, and the one with the lowest mean squared error (MSE) was selected. Next, NNs were trained on that factor in combination with one of the remaining nine. The best performing additional factor was again selected, and this process was repeated until MSE failed to improve. The most parsimonious model was identified using the 1 standard error (SE) rule. Dybowski and Roberts (2001) give the standard error of a bootstrap estimate of a given error metric (e.g. θ=MSE) to be

where θ_boot is the mean of the bootstrapped samples. The smallest set of factors where θ_boot was within one SE_boot of the minimum θ_boot for both NEE and NME was selected for further analysis. The outputs from the selected models are referred to as NEE_NN and NME_NN. NN modelling was done using the Keras Python library (Chollet et al., 2015); see Appendix A for a more detailed explanation of the NN analysis.

Multiple imputation (MI) was then used to gap-fill the NEE and NME with the NEE_NN and NME_NN, respectively (Vitale et al., 2018). Of the 1296 half-hourly flux observations 28.9 % of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and 31.3 % of $F_{{\mathrm{CH}}_{\mathrm{4}}}$ were missing or filtered out. There were a few gaps in the source area fractions needed to gap-fill the flux time series because the footprint function is not valid when $u_{\ast }<\mathrm{0.1}$ m s⁻¹. When source area fractions were missing, they were gap-filled by using the mean source area fraction observed for winds within ±5^∘ of the observed wind direction. The meteorological and soil data were continuous and did not need to be gap-filled.

3.1 EC observations

Half-hourly observations of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ along with the NEE_NN and NME_NN used to gap-fill the time series are shown in Fig. 2c and d. Gap-filled daily NEE ranged from −3.7 to −0.2 g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ with a mean of −1.5 [CI_95 %±0.2] g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Day-to-day variability was considerable but there was no notable trend in NEE over the peak growing season. The half-hourly NEE during the study period reached a minimum of −10.4 $\mathrm{\mu }\mathrm{mol}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ just before solar noon and peaked at 4.7 $\mathrm{\mu }\mathrm{mol}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ around midnight (Fig. 2c). NEE_NN was used to gap-fill the flux data because it was in good agreement with $F_{{\mathrm{CO}}_{\mathrm{2}}}$ observation (r²=0.91). Daily ER_NN was estimated to be 2.2 [CI_95 %±0.9] g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ with corresponding GPP of 3.7 g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. ER_NN was in poor agreement (R²=0.35, n=95) with night-time $F_{{\mathrm{CO}}_{\mathrm{2}}}$ observations. For comparison, Eq. (3) provided a better fit (R²=0.47) with night-time EC data, and ER${}_{Q_{\mathrm{10}}}$ was estimated to be 3.0 g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. However, NEE${}_{Q_{\mathrm{10}}}$ did not fit $F_{{\mathrm{CO}}_{\mathrm{2}}}$ as well (r²=0.80) as NEE_NN.

Gap-filled daily NME was modest and decreased over the study period. It ranged from 2.0 to 25.1 mg C−CH₄ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ with a mean of 8.7 [CI_95 %±0.4] mg C−CH₄ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Fig. 2d). NME_NN was used to gap-fill the flux data because it provided a reasonable fit (r²=0.62) to $F_{{\mathrm{CH}}_{\mathrm{4}}}$ observations. NME did not constitute a significant component of the carbon balance and thus the flux footprint area was a carbon sink during the peak growing season with negative global warming potential (GWP) after accounting for the greater GWP of CH₄ (IPCC, 2104).

3.2 Chamber observations

ER was highest in the sedge, upland, and grass classes where fluxes were very similar at 5.5 [CI_95 %±1.2], 5.4 [CI_95 %±1.2], and 4.9 [CI_95 %±0.7] g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Shrub ER was significantly less (3.5 [CI_95 %±0.6] g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) than the other vegetated classes, and sparse ER was the lowest among the classes (2.0 [CI_95 %±0.3] g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) (Fig. 3a). The Q₁₀ and R₁₀ values also differed between vegetation classes: ER in the sedge class was the most sensitive to changes in air temperature, and modelled values provided the best fit (R²=0.82) to observations. Upland and grass had the highest base respiration and fit observations moderately well (Table 3).

Figure 3Box plot of (a) ER, (b) NME, and (c) NME fluxes measured using closed chambers, grouped by vegetation class. The orange lines represent the median, blue stars represent means, the boxes indicate the interquartile range (Q₁–Q₃), the whiskers indicate $Q_{\mathrm{1}}-(\mathrm{1.5}\times \mathrm{IQR})$ and $Q_{\mathrm{3}}+(\mathrm{1.5}\times \mathrm{IQR}$), and the circles represent outliers extending beyond the whiskers. Note the scale for (c) sedge is different.

Download

NME was more variable between vegetation classes than ER (Fig. 3b and c). Sedge was a very strong CH₄ source at 114.7 [CI_95 %±15.3] mg C−CH₄ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Shrub and grass were very weak sources, 0.7 [CI_95 %±0.3] and 0.4 [CI_95 %±0.3] mg C−CH₄ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively. Sparse was neutral. Upland was a net CH₄ sink −1.1 [CI_95 %±0.4] mg C−CH₄ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Sedge and shrub NME values were positively correlated with T_s (r=0.61, p<0.01; r=0.35, p=0.04) and VWC (r=0.58, p<0.01; r=0.5, <0.01). They also had a positive correlation with T_a, while upland NME was negatively correlated with T_a. Grass and sparse did not have any significant correlations.

Footprint-scaled chamber fluxes were 59 % and 47 % higher than ER_NN or gap-filled NME, respectively. Mean ER_FS was 3.5 g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ [CI_95 %±0.1], it fit ER${}_{Q_{\mathrm{10}}}$ very well (R²=0.95) as would be expected and ER_NN moderately well (R²=0.46). Mean NME_FS was 12.8 [CI_95 %±1.3] mg C−CH₄ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; it did not fit NME_NN well (R²=0.30). At the basin scale, ER_BS (3.4 [CI_95 %±0.1] g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) was slightly lower than ER_FS because of the exclusion of upland areas. NME_BS was higher (15.2 [CI_95 %±0.1] g C−CO₂ ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) because of the greater sedge fraction in the basin than the footprint (Table 2).

3.3 NEE response to environmental factors and vegetation type

NEE_NN (r²=0.91) was estimated using four factors: PPFD, VPD, VWC, and F_Shrub. PPFD is the primary control over NEE: a NN trained on PPFD alone provided a reasonable fit (r²=0.83). The three additional factors, VPD, VWC, and F_Shrub, helped NEE_NN fit a wider variety of conditions. Examining the partial first derivative of NEE_NN under different conditions provides interpretation of the modelled light response curves (Fig. 4). The minimum values represent the peak light use efficiency and are analogous to α in Eq. (5) (Fig. 4b). With increasing PPFD, light use becomes less efficient and approaches zeros as the light response nears light saturation (Fig. 4b).

Figure 4(a) Modelled NEE response to PPFD under different VPD conditions and (b) the partial first derivatives of NEE with respect to PPFD. (c) Modelled ER (dashed line) and NEE (solid line) response to VWC at different Shrub_% values and (d) the partial first derivatives of ER (dashed lined) and NEE (solid line) with respect to VWC. NEE in (c) was calculated at PPFD = 600 µmol m² s⁻¹. The shaded areas in (a) and (c) are 95 % confidence intervals and grey circles are the EC observations.

Download

VPD was a secondary control over NEE. Increasing VPD increased peak light use efficiency and net CO₂ uptake until a threshold, above which it had a strong limiting effect (Fig. 4a and b). For example, under dry atmospheric conditions (e.g. VPD = 1.5 kPa), peak light use is less efficient (−12 nmol CO₂ µmol⁻¹ photon) than under more humid conditions (−18 nmol CO₂ µmol⁻¹ photon). The value of this VPD threshold was dependent upon soil moisture: from 1 kPa when VWC was highest to 0.25 Pa when VWC was low. Mapping NEE_NN and ER_NN at F_Shrub=100 %, F_Shrub=0 %, and F_Shrub=36 % (F_Clim) shows that VWC and F_Shrub were the primary controls over ER and thus influenced NEE (Fig. 4c and d). We can see from the partial first derivates of NEE_NN that increasing VWC increases ER from shrub areas. In the absence of shrubs, increasing VWC inhibits ER, although it is important to note that variations in VWC were subtle, ranging from 51.7 % to 59.0 %. The partial first derivative of NEE_NN shows that VWC slightly limits NEE from non-shrub areas and significantly reduces it in shrub areas.

3.4 NME response to environmental factors and vegetation type

NME_NN (r²=0.62) was estimated using five factors: F_Sedge, F_Shrub, VWC, T_S, and U. NME was more variable and less dependent on any one factor than NEE, which is why the NME_NN needed an extra factor and had a lower r² score. The source area had a significant effect on NME, and it was encouraging that the model contained F_Sedge and F_Shrub since sedge and shrub were the strongest CH₄ source and largest footprint component, respectively. These two factors can combine to map NME under three general situations: we can extrapolate to F_Sedge=100 % and F_Shrub=0 % or F_Sedge=0 % and F_Shrub=100 %, or we can represent actual F_Clim, where F_Sedge=11 % and F_Shrub=37 % (Table 2). Some upland tundra was included in theF_Clim estimate, which reduced NME.

Figure 5(a) Modelled NME response to VWC at different source area fractions and (b) the partial first derivatives of NME with respect to VWC. (c) Modelled NME response to T_s at different source area fractions and (d) the partial first derivatives of NME with respect to T_s. The shaded areas in (a) and (c) are 95 % confidence intervals and grey circles are the EC observations.

Download

VWC was the primary climatic driver identified by NME_NN. Wetter soils had a consistent positive effect on NME, which was strongest when F_Sedgewas high (Fig. 5a and b). Between the driest and wettest conditions, estimated NME increased: by an order of magnitude at F_Sedge=100 %, 4-fold at F_Shrub=100 %, and from neutral to a source at F_Clim (Fig. 5a). Higher T_s generally had a negative effect on NME (Fig. 5c and d). The negative correlation between T_s and VWC (r=0.54, <0.01) may have contributed to this result. NME_NN performance improved less with the addition of U, indicating the NME_NN was near saturation and its effects are less relevant. Higher U had a weak limiting effect on NME when VWC was high and increased NME when VWC was low (not shown).

4.1 Carbon balance and controlling factors

Compared to other DTLBs, Illisarvik has drier soils and greater shrub and grass cover (Table 4). Peak growing season CO₂ uptake at Illisarvik was greater than at most wet-sedge-dominated DTLBs (Table 4; Zona et al., 2010, Sturtevant and Oechel, 2013; Lara et al., 2015). These differences may be due to differences in the periods of observation and year-to-year variability but may also be due to the presence of more productive shrubs and slightly warmer climate at Illisarvik. Mean 1980–2010 T_a at Utqiaġvik (formerly Barrow, AK) is −11.2 ^∘C (US National Climate Data Centre, 2020). Tuktoyaktuk, the closest station to Illisarvik, is 1.1^∘ warmer. Shrub cover is expected to have a number of impacts on the microclimate and carbon cycle of Arctic tundra (e.g. Myers-Smith et al., 2011). Typically, greater deciduous shrub cover is expected to increase GPP as a result of greater leaf area and photosynthetic potential compared to graminoid-dominated tundra (Sweet et al., 2015; Street et al., 2018). GPP was greater at Illisarvik compared to the young wet-sedge-dominated DTLBs in Alaska (Zona et al., 2010). It was more similar to Katyk, which has significant dwarf-shrub cover, predominately Betula nana and Salix pulchra (van der Molen et al., 2007).

Table 4Growing season (gs) daily range in eddy-covariance-derived NEE and NME from drained thermokarst lake basins (DTLBs) and other select wetland and coastal tundra sites across the Arctic.

The periods of study measurements for the study observations are as follows.
^a Mid-June–end of July.
^b 12 June–28 August 2007, Fig. 4.
^c 11 June–25 August 2011.
^d Upscaled chamber estimates, exact dates not specified.
^e Mean 15 June–31 August 2003–2006.
^f 5 July–4 August 2009.

Download Print Version | Download XLSX

Differences in ER among tundra environments can be related to substrate availability, soil moisture and temperature, and thaw depth, among other factors (Sturtevant and Oechel, 2013). The “snow-shrub hypothesis” (Sturm et al., 2001) describes the potential for greater snow trapping in shrub communities which insulates soils in winter, leads to increased decomposition and nutrient availability, and promotes further shrub growth. At Illisarvik, snow blowing in off the Arctic Ocean results in large snow drifts within the basin where snow depth correlates with vegetation height (Wilson et al., 2019). Wilson et al. (2019) concluded that the soils within the Illisarvik basin were warmer than those of the surrounding dwarf-shrub tundra in part through these snow–shrub interactions. Although our chamber observations suggested shrub ER is lower than ER from other vegetation classes, this may have been an artifact as the taller shrubs (>40 cm) could not fit inside the chambers. In another study, chamber ER increased with greater shrub cover in upland tundra (Ge et al., 2017). ER at Illisarvik was greater than the ER observed at both the young wet-sedge DTLB in Barrow (Zona et al., 2010) and at the shrub–wet-sedge DTLB at Katyk where thaw depth was much shallower (45 to >100 cm at Illisarvik vs. 25 to 40 cm at Katyk; van der Molen et al., 2007). The importance of F_Shrub in describing temporal variations in half-hourly NEE within the flux footprint at Illisarvik is further evidence of the importance of shrub cover on tundra carbon cycle processes in this environment.

PPFD and VPD were the most important factors for predicting half-hourly NEE. This was to be expected as they are typically the primary controls over GPP (Aubinet et al., 2012). The limiting effects of VPD are consistent with another study using NN to analyze NEE at a deciduous forest site (Moffat et al., 2010) and have been found at other tundra sites (Euskirchen et al., 2012; López-Blanco et al., 2017). VWC was also important at Illisarvik. Zona et al. (2010) found VWC could explain 70 % of the variability in daily peak season ER in a young DTLB. Similarly, Kittler et al. (2016) found drier soils increased ER and decreased NEE after a wet tundra drainage experiment in Siberia, consistent with our results at Illisarvik when F_Shrub was low.

As expected, NME at Illisarvik was about half that observed at the Alaskan DTLB sites where soils were wetter with greater sedge cover (Table 4, Zona et al., 2009; Lara et al., 2015). NME at Katyk was even higher than the Barrow DTLBs and had a significant impact on the greenhouse gas (GHG) balance for this site (van der Molen et al., 2007; Parmentier et al., 2011). In our NN modelling of NME at Illisarvik, F_Sedge was the most important factor for predicting half-hourly $F_{{\mathrm{CH}}_{\mathrm{4}}}$. Sedges are aquatic plant species with aerenchymatous tissues that act as conduits for CH₄ from below the water table to the atmosphere and limit CH₄ oxidation by methanotrophs in aerobic surface soils (Lai et al., 2009). The inclusion of F_Shrub further refined the model, allowing it to better fit the site-specific distribution of vegetation types. Budishchev et al. (2014) found shrub and sedge fraction had a significant influence on $F_{{\mathrm{CH}}_{\mathrm{4}}}$ at Katyk. Vegetation type is the dominant control over NME across multiple tundra landscapes and our results further support that (Davidson et al., 2016).

VWC was the second most important factor, which was expected as CH₄ production occurs in anaerobic environments and has been linked to variability in CH₄ emission in many other studies (e.g. Zona et al., 2009; Nadeau et al., 2013; Olefeldt et al., 2013). Soil temperature (T_s) was the third most important factor. Higher T_s values increase the oxidation potential of methanotrophs (Liu et al., 2016; King and Adamsen, 1992), so this result was expected for the drier portions of the basin and upland tundra. However, this was not expected for the sedge areas because most studies find NME in sedges is positively correlated to T_s (Olefeldt et al., 2013). The negative correlation between T_s and VWC may partly explain this.

4.2 Upscaling

ER_FS and NME_FS were about 59 % and 47 % greater than the gap-filled EC estimates. Discrepancies between EC and chamber observations are common and have been attributed to differences in measurement techniques, the small sample size of chamber observations, and sampling bias since all chamber measurements were taken during the day with fair weather (Katayanagi et al., 2005; Chaichana et al., 2018). Meijide et al. (2011) found that chamber NEE could be up to twice as large as EC observations, and Riederer et al. (2014) also found chamber NME estimates were about 30 % higher than EC estimates. Others have been more successful, yielding upscaled chamber NME fluxes within 10 % of EC observations (Zhang et al., 2012; Budishchev et al., 2014; Davidson et al., 2017). A potential reason for the disagreement with ER_FS may be the lack of direct observations by the EC system under low-light conditions. Another potential source of error for the upscaling is inaccuracies in the vegetation map.

4.3 Future trajectories

Presently, peak growing season carbon uptake at Illisarvik is greater than similarly aged landscape features on the Barrow Peninsula, Alaska, and more similar to levels observed at Katyk, Siberia. NME is well below levels observed at any other DTLB studied, making this site a stronger GHG sink than other DTLBs. However, the basin at Illisarvik will continue to evolve, and the trajectory it takes could significantly alter its carbon balance. Historically, DTLBs on Richards Island and the Tuktoyaktuk Peninsula evolve into sedge wetlands, as do DTLBs on the Barrow Peninsula (Ovendend, 1986; Lara et al., 2015). Active maintenance of the outlet channel at Illisarvik has artificially lowered soil moisture and flooding and potentially limited this transition thus far (C. Burn, personal communication, 2016).

If Illisarvik follows the same trajectory as older DTLBs in the area and becomes dominated by sedge wetlands, NME will increase significantly. Figure 5a shows that with extrapolations to full sedge cover (F_Sedge=100 %), NME would be similar to values on the Barrow Peninsula (Zona et al., 2009). If the basin instead transitions into a shrub-dominated DTLB similar to those of Old Crow Flats, Yukon (Lantz and Turner, 2015), NME_NN would remain similar to current levels, meaning the basin would remain a weak source of CH₄. These are projections well beyond F_Clim fractions observed, so confidence in the specific values predicted is low.

The effects of changing shrub and sedge cover on Illisarvik's growing season NEE are less straightforward than on NME, partly because shrub cover had less overall influence on NEE_NN. Figure 4c shows that the model suggests ER decreases and NEE increases with increasing shrub coverage when soils are slightly drier but has the opposite effect under wetter conditions. To our knowledge, only a few winter season (e.g. Zona et al., 2016) and no year-round studies of DTLB NEE and NME have been published to help evaluate the factors influencing carbon losses through the non-growing-season months. Further observation year-round is needed to better understand the implications of continued vegetation change for the carbon balance of DTLBs such as Illisarvik.