Seasonality buffers carbon budget variability across heterogeneous landscapes in Alaskan Arctic tundra

EC study sites are located on continuous permafrost tundra on the North Slope of Alaska, near Utqiaġvik (figure 1(a)). The sites include the Barrow Environmental Observatory (US-Beo) (figure 1(b)), Biocomplexity Experiment South (US-Bes) (figure 1(c)) and a site near the NOAA Earth System Research Laboratory (US-Brw) (figure 1(d)). These sites were chosen for long-term continuous data acquisition as they capture dominant landscape variability of the region. US-Bes is in a drained lake basin containing the wettest soils (table 1), with the water table above the surface for most of the growing season and is dominated by wet sedges and sphagnum moss (Davidson et al 2016a). US-Brw is a moist upland tundra containing the driest soils (table 1) and is dominated by graminoids and lichens (Kwon et al 2006). US-Beo is characterized by ice wedge polygon formations that arise from the freeze thaw cycle (Webber 1978). Due to these polygon formations, US-Beo is a mixed landscape, exhibiting both inundated and drained areas and consists of wet sedge/sphagnum moss dominated areas as well as drier graminoid/lichen dominated areas (Davidson et al 2016a).

Zoom In Zoom Out Reset image size

CH₄ and CO₂ fluxes were estimated at half-hourly intervals from year-round data collected at 10 Hz following the procedures outlined from LI-COR® EddyPro®. A double rotation was applied to the axis rotations of three-dimensional wind speeds according to Wilczak et al (2001) and a block averaging interval was used to define turbulent fluctuations. An in situ/analytic correction, according to Ibrom et al (2007), was applied to the gas analyzer data as the greenhouse gas analyzer (GGA) has a closed path. Quality flags were output within datasets according to Mauder and Foken (2011) and data that did not pass the quality requirements were removed. An internal chamber pressure of ⩾20.67 kPa (155 torr) in the GGA indicates line blockage or instrument failure, and these data were eliminated. Additionally, a turbulence threshold was applied, identifying conditions with insufficient turbulence (indicated by low friction velocity (u* ⩽ 0.1 m s⁻¹)), and those data were removed in accordance with Reichstein et al (2005). A moving window of 2 weeks was applied and fluxes that were three standard deviations away from the mean were removed as outliers for CH₄ and CO₂ fluxes. EC tower site and instrumentation information can be found in table 1 (Goodrich et al 2016, Arndt et al 2020).

Meteorological data were obtained at 30 s intervals and averaged into half-hourly means. Each of the three sites had independent meteorological instrumentation and measurements used in defining the peak thaw and zero-curtain period. Meteorological instrumentation used in analysis included soil water content (SWC) (Campbell Scientific® CS616 Water Content Reflectometer), soil temperature (Omega Engineering™, type T thermocouples), air temperature and humidity (Vaisala, HMP 45), and photosynthetically active radiation (PAR; LI-COR® LI-190R quantum sensor). At each site, soil temperature was measured at 0, −5, −15, and −30 cm from the soil surface and soil moisture probes were inserted from the soil surface to a depth of 20 cm, providing the average soil moisture in the top 20 cm of the soil column. Measurements of soil moisture were taken at one location at US-Brw and US-Bes, and two locations at US-Beo (a high center polygon and polygon trough) for better relief representation then averaged over the growing seasons within the study period. Data from both EC and meteorological instrumentation were collected using datalogger/multiplexer arrays from Campbell Scientific® (i.e. CR-3000, CR-23X and AM-1632). Thaw depth was measured weekly during the growing season by probing the land surface to the extent of the active layer along a transect with a small diameter metal rod.

EC tower footprints were estimated with the analytical footprint model of Korman and Meixner (KM) (2001) using the R package 'FReddyPro' v1.0 (Xenakis 2016). The KM model calculates the density function of the footprint contribution for a two-dimensional area surrounding the EC tower and was used to estimate the landscape area in which 80% of fluxes originated by averaging half-hourly single flux footprints during 2016. WorldView-3 (Maxar Technologies) imagery (1.24 m multispectral resolution) of Utqiaġvik, AK was used to show variability in the WorldView normalized difference water index (NDWI) at each of the three sites. As SWC is measured at one to two locations within each EC tower footprint and site hydrology can be variable, NDWI is used in conjunction with SWC to characterize differences in site moisture regime. The WorldView NDWI is calculated as the normalized difference between the coastal band (R_c, 400–450 nm) reflectance and the second near infrared band (R_NIR2, 860–1040 nm) reflectance (equation (1)). This is because it has been shown that the R_c and R_NIR2 bands show a better soil-water separation than the typical green (510–580 nm) and the first near infrared band (770–895 nm) combinations and is indicative of surface water moisture levels (Maglione et al 2014). The NDWI results in values between −1 and 1 where more positive values represent wetter landscapes. Pixels outside of the footprint of the EC towers as well as those representing structures within the tower footprints were masked for statistical analysis. Before calculating statistics, pixels were aggregated into 3 × 3 pixel grids using the nearest neighbor approach to avoid bias by over sampling and to reduce high variability. A pairwise Wilcoxon rank sum tests was used to compare NDWI values among sites. Additional imagery for NDWI analysis across study period can be found in supplementary information (figure S1 (available online at stacks.iop.org/ERL/16/035008/mmedia))

$$\begin{equation}{\text{NDWI}} = { }\frac{{{R_{\text{C}}} - {R_{{\text{NIR}}2}}}}{{{R_{\text{C}}} + {R_{{\text{NIR}}2}}}}.\end{equation} \tag{ 1 }$$

Daily average fluxes were calculated in R V 3.6.2 (R Core Team 2019) and R Studio software using the 'data.table' package and were calculated with a minimum of 30 half-hourly samples per day to ensure proper representation of the diurnal patterns of CH₄ and CO₂ fluxes. Daily averages of CH₄ and CO₂ fluxes for purposes of the examination of site variability and C budgets are beneficial as they show an accurate representation of the systems while qualifying data complexity and size. Data gaps are unavoidable in the harsh conditions characteristic of Arctic environments. These gaps can be a result of power or network outages as well as instrument failure. Total data coverage is between 61% and 71% depending on the site with the best coverage during the summer and fall periods. Detailed data coverage information by season and annual totals can be found in the supporting information (table S1). Data gaps were filled using random forest machine learning (R package, 'missForest') utilizing a 300-decision tree design. Model validation can be found in supporting information (figure S2). Comparisons of model validation between the default half-hourly data output and daily averages show models perform much better when using a daily average. This method reduces the 'noise' and is therefore better equipped to inform machine learning processes.

The beginning of the growing season was defined as the period where the top five cm of the soil are above 0 °C, ending at the onset of the zero-curtain. The zero-curtain was defined as the period during the fall shoulder beginning when soil temperature of the top 5 cm of the soil are less than 0 °C for 3 or more days, ending when the temperature at −15 cm (roughly the middle of the active layer, table 1) dropped below −0.75 °C for 3 or more days. Non-growing season, as defined here, includes both winter and spring, beginning at the end of the zero-curtain period and ending at the beginning of the growing season.

ER and GPP were partitioned from net ecosystem exchange (NEE) according to Lasslop et al (2010) using the 'REddyproc' package in R (Wutzler et al 2018), as nighttime data is unobtainable during Arctic summer. Temperature response relationships were calculated using a weekly mean to reduce noise and to better represent annual trends. Models showing temperature response curves for CH₄ were calculated using soil temperature, rather than air temperature, as this has been shown to act as a better predictor for CH₄ fluxes (Arndt et al 2019a) (figure S3). All temperature response curves were linearized by log transformation and compared using analyses of covariance (ANCOVA) to test for homogeneity among the regressions with the 'car' R package (Fox and Weisberg 2019) and Q₁₀ values were calculated from temperature response regressions with the 'respirometry' R package (Birk 2020).

NDWI as well as the KM model footprint of each EC tower were used to establish differences in surface water content within each EC tower footprint (figure 2). Using a pairwise Wilcoxon rank sum tests, NDWI was found to be significantly different (p < 0.001) among each of the three sites. The results of the NDWI, showing levels of surface moisture, agreed with soil moisture data at the sites (table 1) with US-Bes showing the wettest conditions (NDWI (mean ± standard error) = −0.079 ± 0.005), US-Beo showing intermediate NDWI levels (NDWI = −0.144 ± 0.003) and US-Brw with the lowest NDWI supporting its position as the driest site in the study (NDWI = −0.204 ± 0.003). Further analyses of NDWI show that imagery acquired 24 July 2016 was representative of site differences and while some variability occurs, these positions are maintained (figure S1) as changes in vegetation community or hydrology at the landscape scale happens over longer periods of time (Liljedahl et al 2016, Arndt et al 2019b).

Zoom In Zoom Out Reset image size

Peak growing season CH₄ emissions were higher at US-Bes (1.43 ± 0.11 mg C–CH₄ m⁻² h⁻¹) and US-Beo (1.35 ± 0.31 mg C–CH₄ m⁻² h⁻¹) in comparison with US-Brw (0.76 ± 0.10 mg C–CH₄ m⁻² h⁻¹), however, CH₄ fluxes during the zero-curtain period showed lower variability across the three study sites (figure 3(a)). The annual average of peak uptake in NEE was greatest at the driest site, US-Brw (75.26 ± 8.9 mg C–CO₂ m⁻² h⁻¹), followed by US-Beo (60.5 ± 7.9 mg C–CO₂ m⁻² h⁻¹) and least pronounced at wettest site, US-Bes (39 ± 6.4 mg C–CO₂ m⁻² h⁻¹). The annual average of peak emission in NEE followed the same order of US-Brw (44.25 ± 4.9 mg C–CO₂ m⁻² h⁻¹), US-Beo (24.74 ± 5.3 mg C–CO₂ m⁻² h⁻¹) and US-Bes (13.75 ± 1.7 mg C–CO₂ m⁻² h⁻¹) (figure 3(b)).

Zoom In Zoom Out Reset image size

Seasonal fluctuations (growing season to zero-curtain) of CO₂ budgets are most pronounced at the driest site, US-Brw, and least pronounced at wettest site, US-Bes (figure 4(a)). Due to the higher zero-curtain CO₂ emissions dampening growing season CO₂ uptake, US-Brw is the weakest CO₂ sink on average. This trend holds with drier sites generally emitting larger amounts of CO₂ during the zero-curtain offsetting much of the uptake during the growing season (US-Brw: 75%; US-Beo: 62%; US-Bes: 25%). Growing season CH₄ emissions are highest at US-Beo and US-Bes (figure 4(b), p < 0.05). However, zero-curtain CH₄ emissions are roughly equal across all landscapes. This indicates that the percent contribution to local CH₄ budget varies by site during the zero-curtain (US-Brw: 45%; US-Beo: 34%; US-Bes: 32%) and total non-growing season (US-Brw: 56%; US-Beo: 48%; US-Bes: 43%).

Zoom In Zoom Out Reset image size

Site variability in combined C budgets are most pronounced during the growing season, yet still occur, driven by CO₂ emission, during the zero-curtain period (figure 4(c)). As US-Beo is characterized by low centered polygons, it exhibits characteristics of both US-Brw and US-Bes. This accounts for the larger values observed in CH₄ emission relative to US-Brw and the larger values in net CO₂ uptake relative to US-Bes during the growing season. The effect is that the mixed landscape, US-Beo, is the largest mean C contributor (table 2). More detailed information regarding C budget by year/season can be found in supplementary information (table S2).

Yearly cumulative emissions indicate that interannual variability in CH₄ fluxes is limited in inundated sites (figure 5(a)). US-Bes exhibited the lowest interannual variability (standard deviation—CO₂-eq (σ) = 3.96), followed by US-Brw (σ = 15.2) then US-Beo (σ = 23.8). Interannual variability in CO₂ fluxes is highest at US-Brw (σ = 34.98), followed by US-Bes (σ = 18.11) then US-Beo (σ = 10.65) (figure 5(b)). Most of the variability in combined cumulative C emissions is therefore controlled by variability in CO₂ at US-Brw and US-Bes, and in CH₄ at US-Beo (figure 5(c)).

Zoom In Zoom Out Reset image size

ER temperature response relationships show that air temperature increases in predictive strength as site wetness decreases (figures 6(a)–(c)). In the inter-site comparison of linearized ER regressions, US-Beo and US-Brw show significant similarity (p < 0.05) while US-Bes is significantly different than both US-Brw (p = 0.12) and US-Beo (p = 0.076). Q₁₀ is highest at US-Beo (3.5), followed by US-Brw (2.5) and US-Bes (2.2). Methane temperature response relationships show similar predictive strength across all sites (figures 6(d)–(f)). Q₁₀ for CH₄ is again highest at US-Beo (4.6), followed by US-Bes (4.2) and US-Brw (3.1). Contrary to ER, in the inter-site comparison of linearized CH₄ regressions, US-Beo and US-Bes show significant similarity (p < 0.05) while US-Brw is significantly different than both US-Beo (p = 0.11) and US-Bes (p = 0.14).

Zoom In Zoom Out Reset image size

Seasonality appears to be the dominant factor in annual variation in NEE; however, we find the magnitude of that effect is dependent on local variation in site hydrology (figure 4). US-Brw exhibited the highest GPP yet was the weakest mean annual CO2 sink in comparison with the other wetter sites (table 2). Conversely, US-Bes exhibited the lowest GPP yet was the strongest mean annual CO2 sink (table 2). This is because sites with higher summer GPP exhibited larger zero-curtain CO2 emissions. The lower SWC at US-Brw likely leads to a substantially larger portion of the soil column under oxic conditions, therefore supporting aerobic respiration, increasing CO₂ emissions. While it is probable that this constitutes much of the budget disparity, the contribution of CH₄ oxidation to ER in primarily methanogenic areas has been found to be up to 35% (Nielsen et al 2019). As US-Brw contained the deepest active layer (table 1) and largest GPP (table 2), soils may contain higher amounts of photosynthates and labile C. This can increase methanogenesis rates deeper in the soil column (Dorodnikov et al 2011) and fuel CO2 producing methanotrophs closer to the surface that would be more active under the oxic conditions in drier areas relative to waterlogged areas (Megonigal and Schlesinger 2002). Alternatively, as US-Bes is inundated, methanotrophy is likely substantially lower during the zero-curtain period when the highest amount of CO2 loss is observed. Moreover, the lower GPP at US-Bes may result in more recalcitrant C and thus lower ER rates. The US-Beo site exhibited intermediate growing season GPP and zero-curtain CO2 emission relative to the other two sites likely due to US-Beo exhibiting a mixed landscape regarding the prevalence of drained and inundated areas (figure 4(a)).

ER temperature response relationships show that US-Beo had the strongest temperature dependence (Q₁₀ = 3.5) relative to US-Brw and US-Bes (Q₁₀ = 2.5 and Q₁₀ = 2.2) (figure 6). This may be linked to the intermediate soil moisture of US-Beo. Higher SWC can limit the temperature sensitivity of soil respiration in wetland regions due to the restriction of oxygen and thus, aerobic respiration (Chen et al 2018). Alternatively, respiration can be limited by lower SWC through a reduction in microbial mobility and substrate diffusion (Grant and Rochette 1994). The temperature sensitivity of belowground respiration can also be dependent on productivity by providing photosynthates as substrates (Hartley et al 2006). This may contribute to US-Bes having the lowest temperature sensitivity, as this site exhibited the lowest productivity (table 2).

Similar studies of annual CO₂ budgets show that the largest annual net CO₂ loss is seen during the non-growing season, particularly associated with early winter respiration (Oechel et al 2014, Commane et al 2017, Euskirchen et al 2017), agreeing with data presented here. As early winter respiration comprises a large part of the CO₂ budget, further increases in zero-curtain duration will likely result in winter CO₂ emissions that exceed growing season uptake (Arndt et al 2019b). Each of the sites in this study were found to act as a weak sink for CO₂, yet others have reported relatively strong annual source signals from similar systems (Commane et al 2017, Eurskirchen et al 2017). This highlights the need for the monitoring and greater representation of mesoscale (⩽1 km) processes in climate projections, particularly at sub-grid scales.

Growing season CH₄ emissions are lowest at US-Brw (figure 4(b)). This is likely due to lower soil moisture increasing the volume of soil experiencing aerobic conditions not conducive to methanogenesis (Garcia et al 2000). The polygon tundra site, US-Beo, had similar CH₄ emissions relative to the inundated US-Bes site despite US-Bes being the site with the higher SWC (figure 4(b)). This may be related to the vegetation community composition. As US-Bes contains a lower percent cover of sedges (Davidson et al 2016b), US-Beo may produce more photosynthates, like acetate, that could leach into surrounding waterlogged soil in polygonal environments, further fueling methanogenesis via acetoclastic methanogenic pathway (King et al 2002, Dorodnikov et al 2011). Further, sedge density is positively correlated to CH₄ emissions as sedges provide a pathway for CH₄ through the vegetation to the atmosphere (Lai 2009, Andresen et al 2017). Interannual variability in CH₄ emissions is low at US-Bes compared with the other sites (figure 5). As US-Bes is consistently inundated, interannual differences in snow melt and rainfall would have a larger impact on SWC and by proxy, oxygen availability at US-Brw, and US-Beo, possibly explaining this variability.

Contrary to zero-curtain CO₂ emission having site dependent variability, zero-curtain CH₄ emissions are roughly equal across all sites. This may be due to the frozen surface soils, creating an ice 'cap', and limiting oxygen diffusion into the soil column thereby equalizing oxygen availability and by extension, methanogenesis, across sites. This shows that variability in the growing season, rather than the zero-curtain, may have a stronger impact on annual CH₄ variability across different landscapes. However, zero-curtain emissions may increase as the zero-curtain extends longer into the winter with a warming climate (Arndt et al 2019a). Zero-curtain CH₄ contributions were found to be higher than in previous works due to the length of the study period capturing the range of annual variability (Zona et al 2016). Methane temperature response relationships also indicate that temperature dependence was strongest at US-Beo (Q₁₀ = 4.6), followed by US-Bes (Q₁₀ = 4.2), and US-Brw (Q₁₀ = 3.1) (figure 6). As US-Beo, like US-Bes, contains large amounts of anaerobic soil, temperature sensitivity of methane production would be stronger than US-Brw where soil moisture is a limiting factor.

The largest mean combined C emissions were from the mixed landscape US-Beo, exhibiting both inundated and drained areas. These polygonized landscapes comprise close to 65% of the Alaskan coastal plain (Lara et al 2018) and contain both anaerobic areas that produce large quantities of CH₄ as well as drained areas where aerobic respiration can readily occur. US-Beo exhibited the strongest temperature response, for both CO₂ and CH₄. On this basis, it is possible that further climate change may disproportionately increase C emissions from polygonized landscapes as rising temperatures will support increased production and emission of CO₂ and CH₄. However, rising temperatures will likely coincide with polygon succession and hydrologic transitions (Liljedahl et al 2016). These hydrologic transitions can significantly alter annual carbon budgets (Kittler et al 2017), stressing the importance of monitoring landscape heterogeneity for in these regions for carbon budget estimation.

The interplay of CO₂ and CH₄ dynamics are affected strongly by both seasonality and by mesoscale landscape variability. Though the variability in summer C emissions is significant among the landscapes studied, zero curtain releases of CO₂ and CH₄ tend to offset site differences. This acts as a buffer to variability and leads to similar annual combined C budgets across the sites studied. However, the differences in timing and magnitude of CO₂ and CH₄ fluxes elevate the importance of mesoscale processes for restricting uncertainty in Arctic model projections. Arctic regions make up the largest portion of uncertainty in climate global climate models (IPCC 2014). Pan-Arctic models are typically run at coarse scales that describe landscape heterogeneity by the dominant landscape. It has been stressed that a higher degree of spatial and temporal coverage is needed (Natali et al 2019) and that representation of wet and dry tundra at a finer scale (⩽4 km²) can result in a threefold reduction in model error (Lara et al 2020). Data presented here demonstrate the need for this improvement, particularly for models that can represent mesoscale landscape heterogeneity and subsequent differences in seasonal carbon emission patterns.