Understanding spatial variability of methane fluxes in Arctic wetlands through footprint modelling

The Arctic is warming at twice the rate of the global mean. This warming could further stimulate methane (CH4) emissions from northern wetlands and enhance the greenhouse impact of this region. Arctic wetlands are extremely heterogeneous in terms of geochemistry, vegetation, microtopography, and hydrology, and therefore CH4 fluxes can differ dramatically within the metre scale. Eddy covariance (EC) is one of the most useful methods for estimating CH4 fluxes in remote areas over long periods of time. However, when the areas sampled by these EC towers (i.e. tower footprints) are by definition very heterogeneous, due to encompassing a variety of environmental conditions and vegetation types, modelling environmental controls of CH4 emissions becomes even more challenging, confounding efforts to reduce uncertainty in baseline CH4 emissions from these landscapes. In this study, we evaluated the effect of footprint variability on CH4 fluxes from two EC towers located in wetlands on the North Slope of Alaska. The local domain of each of these sites contains well developed polygonal tundra as well as a drained thermokarst lake basin. We found that the spatiotemporal variability of the footprint, has a significant influence on the observed CH4 fluxes, contributing between 3% and 33% of the variance, depending on site, time period, and modelling method. Multiple indices were used to define spatial heterogeneity, and their explanatory power varied depending on site and season. Overall, the normalised difference water index had the most consistent explanatory power on CH4 fluxes, though generally only when used in concert with at least one other spatial index. The spatial bias (defined here as the difference between the mean for the 0.36 km2 domain around the tower and the footprint-weighted mean) was between ∣51∣% and ∣18∣% depending on the index. This study highlights the need for footprint modelling to infer the representativeness of the carbon fluxes measured by EC towers in these highly heterogeneous tundra ecosystems, and the need to evaluate spatial variability when upscaling EC site-level data to a larger domain.


Introduction
Methane (CH 4 ) emissions from Arctic permafrost soils are a major source of uncertainty in the region's future global warming potential Abbott 2011, IPCC 2013). The Arctic is warming at twice the rate of the global mean (Blunden and Arndt 2019) and its frozen permafrost soils store 1300-1370 Pg of organic carbon (Hugelius et al 2014), twice the current atmospheric stock (IPCC 2013). By the year 2300, thawing permafrost could release between 381 and 616 Pg of carbon to the atmosphere (Schuur et al 2013); and it is imperative to know how much of this carbon will be released as CH 4 , which on a per-molecule basis has a global warming potential 25-28 times greater than carbon dioxide (CO 2 ) (Forster et al 2007, Etminan et al 2016. One method of measuring trace gas fluxes central to understanding the current and future carbon budget is the eddy covariance (EC) technique, as it can bridge the gap between smaller plot scale chamber measurements and larger regional scale remotely sensed (RS) data from aircraft and satellite (Baldocchi 2003, Chen et al 2009. There are, however, challenges to using this technique in the Arctic's heterogeneous landscapes (Davidson et al 2016) as it assumes a uniform sampling area (Foken and Wichura 1996). One approach used to address this issue is footprint modelling (Vesala et al 2008).
The spatial heterogeneity of Arctic wetlands is often related to the presence of cryogenic processes, where the formation and degradation of ice wedges, leads to patterned ground formations such as polygonal tundra (Brown 1967). The polygons consist of distinct microtopographic features, namely rims, troughs and polygon centres, which each have distinct CH 4 emission rates (Sachs et al 2010, Lara et al 2015, Davidson et al 2016, Vaughn et al 2016. These features have differential water drainage and soil moisture, a well-established driver of CH 4 production and consumption, with water-logged anaerobic areas producing CH 4 and dry areas sometimes acting as slight CH 4 sinks (Valentine et al 1994, Segers 1998, Sachs et al 2010, Von Fischer et al 2010, Lipson et al 2012. Typically, rims are well-drained, while troughs are inundated, with polygon centres being either convex (highcentres) or concave (low-centres), depending on age, and thus dry or wet respectively. Microtopography also influences plant (Joabsson et al 1999, Von Fischer et al 2010 and even microbe (Taş et al 2018) distributions which likewise impact CH 4 production and efflux. Sedges grow in water logged areas and release organic acid root exudate that increases CH 4 production, furthermore, a number of species have aerenchymal tissue which allows CH 4 to avoid oxidation by passing through the plant (McEwing et al 2015, Andresen et al 2017). Additionally, Arctic wetlands are further mosaiced due to the thermokarst thaw lake cycle (Zona et al 2009, Sturtevant andOechel 2013). These lakes periodically drain, leaving depressions and drainage channels that vary in the same factors affecting CH 4 emission in polygonal tundra. In addition to spatial variation, CH 4 fluxes vary over time, as water table, thaw depth, soil temperature, plant productivity, and the available carbon pool, all major controls on CH 4 production, change throughout the year (Zona et al 2009, Lipson et al 2012, Zheng et al 2018. A flux tower footprint defines the area being sampled by the EC tower Thurtell 1990, Horst andWeil 1994). Footprint modelling has been used to enable the interpretation of EC data in variable landscapes (Schmid and Lloyd 1999, Göckede et al 2006, Parmentier et al 2011, Tuovinen et al 2018, by assigning relative flux contribution to specific areas based on tower height, wind speed and direction, and turbulence (Leclerc and Thurtell 1990, Schuepp et al 1990, Burba and Anderson 2010. There are a few different types of footprint models including one-and two-dimensional analytic models, Lagrangian, large eddy simulations, and closure models (Vesala et al 2008). Of these, the most commonly used are the analytical models, specifically the Kormann and Meixner (2001), Kljun et al (2004), andHsieh et al (2000) models, due to having relatively low computational complexity and their applicability in a wide array of experiments (Vesala et al 2008, Leclerc andFoken 2014).
Footprint modelling has been utilised in several ways to interpret flux variability in the Arctic. One way footprint modelling can be used is to relate EC measurements to chamber measurements by upscaling them to the EC footprint. In a sub-Arctic mire site in Finland, the Kormann and Meixner (2001) footprint model was applied to get half-hourly footprint-weighted spatial indices, these were used as in-put for a process-based model and improved the correlation between growing season upscaled chamber CH 4 fluxes and EC flux estimates, with an increase in r 2 from 0.41 to 0.72 (Hartley et al 2015). Budishchev et al (2014) also used a footprint model to upscale the chamber measurements to the EC scale, improving the r 2 correlation coefficient from 0.14 to 0.7, in a Russian polygonal tundra. Davidson et al (2017) show a 20%-30% improvement and an r 2 of 0.88, across several tundra sites in Alaska by using footprint modelling to upscale chamber measurements to the EC tower fluxes. One can also use footprint modelling when upscaling fluxes to estimate sensor location bias (Schmid and Lloyd 1999). This metric is the percent difference of a spatial variable's mean within the footprint and the mean of the user-defined area to which one is upscaling. In a study across the entire Canadian flux network (Chen et al 2011), sensor bias was assessed through the distribution of enhanced vegetation index and normalised difference vegetation index (NDVI) showing that four out of twelve sites presented a difference higher than 5% between the EC tower annual footprints (with 90% of the footprint generally within a 1 km 2 radius) and the surrounding area (measured at 1, 2, and 3 km 2 centred at the tower base). Importantly the sites with the highest bias, ranging from −14% to 9%, were classifies as ombrotrophic bog (Chen et al 2011), a vegetation type similar to sites evaluated in this study. More recently, in Siberian tundra Tuovinen et al (2018) used footprint modelling on a half-hourly time scale to parse CH 4 sources and sinks by vegetation type, they also looked at the sites sensor location bias and found a 14% when looking at leaf area index regarding the wider 6.3 km 2 around the base of the tower.
The aim of this study is to evaluate the effect of spatial heterogeneity on the CH 4 fluxes measured by EC method at two tundra sites located on a large wetland area in the North Slope of Alaska. First, we assessed whether the area sampled by the EC measurement, as characterised by the footprint, was representative of the 0.36 km 2 domain around the tower. This was done to determine if the area around the base of the tower was suitable to link to RS data for the purpose of upscaling. We then used statistical modelling, to separate the variability in the measured flux arising from the (i) temporal variability in the environmental controls on the CH 4 fluxes (temperature, soil moisture etc), and (ii) the variability in the footprint, which sampled different parts of the surrounding landscape at different times. We expect the spatial variability in the landscape to have a substantial influence on the fluxes. We also expect that there may be significant sensor bias between the tower footprint and the domain around the tower, due to the inherent heterogeneity of polygonal tundra, which would need to be considered for upscaling the fluxes.

Study sites
Two EC towers located near Barrow (Utqiaġvik), on the North Slope of Alaska, were used in this study (figure 1). Average annual temperature and precipitation measured at the Barrow weather station, between 1948 and 2013, were −11.3°C and 72 mm respectively . The two tower sites utilised for this study include the and Biocomplexity Experiment South (BES) and barrow environmental observatory (BEO). BES and BEO are less than 1 km apart and have similar environmental conditions. The sites share similar vegetation cover, namely sedges (Eriophorum russeolum, Carex aquatilis), grasses (Poa arctica), and mosses (Dicranum spp.) with heights well below 40 cm (Davidson et al 2016). A portion of the footprint area for each of these sites takes the form of well-developed polygonal tundra. Additionally, a substantial part of the footprint in BES is occupied by a drained lake basin (figures 1, 2). The water table within these sites is highly variable given the very different microtopography. Polygon rims are well-drained, while the polygon centres and the drained lake areas have high variability in water table (∼20 cm seasonal range), with some troughs and particularly low-lying areas remaining inundated through most of the summer (Zona et al 2009).

EC measurements and footprint modelling
EC data were collected at 10 Hz with a Campbell Scientific CSAT-3 sonic anemometer and a closed-path LGR analysers (FGGA, Los Gatos Research, Mountain View, CA, USA) at a height of ∼2-3 m (BES: 2.20 m, BEO: 3.12 m). Fluxes were calculated with the EddyPro software package (LI-COR) based on the methodology described in Zona et al (2016). Data were filtered to remove spikes and data unsuitable for footprint modelling as described in Goodrich et al (2016) and Foken et al (2004). Measurement of CH 4 at these sites between July of 2013 and November 2015 is used in this study. The EC data used in this study can be accessed via the Arctic Data Centre (Zona 2019).
The model of (Kormann and Meixner 2001) was applied on a 600×600 m (0.36 km 2 ) grid around the towers to calculate the flux footprint for each halfhourly period (equation (1)). This domain size was chosen based on the assumption that most of the measured fluxes originates within a radius equal to 100 times the height of the tower (Burba and Anderson 2010). During stable atmospheric conditions footprints are elongated and diffuse, leading to some of footprint falling outside of the 0.36 km 2 around the base of the tower. This was tolerated up to the point where >15% of the footprint was outside the 0.36 km 2 area, after which footprints were excluded from the analysis. The grid spatial resolution was set to 1 m 2 , so that the small-scale variability of the polygon troughs was represented. With this model, the probability of a grid cell with co-ordinates x, y contributing to the measured flux is: where s v is the standard deviation of lateral wind speed, L is the Monin-Obukov length, u * is the friction velocity, h is the height of the sensor, and u(h) is the horizontal wind velocity at the measurement height. Thus, w represents the relative contribution of each grid cell to the measured flux. The raster grid of each half-hourly footprint, w, can then be used as a weighting factor to calculate the footprint-weighted average of each of the spatial variables, such as elevation or NDVI.

Temporal and spatial variables
The spatial metrics used in this study were obtained as follows. High resolution topographic data was derived from an airborne LiDAR survey conducted on 12 July 2013 (Wilson and Altmann 2015). The LiDAR derived digital elevation model (DEM) has a horizontal resolution of 0.25 and a 0.143 m vertical resolution. A WorldView-2 (WV2) multispectral satellite image collected on the 6th of July 2013 was used to calculate NDVI and normalised difference water index (NDWI) This image consists of 8 spectral bands with a 1.84 m horizontal resolution. Only this image was used in this study, as it is the only WV2 image collected during the study period, June 2013 to Dec 2015, free of clouds and snow. Sedge cover was obtained from the wetlands map created by the BAID project (Barrow Area Information Database) (Andresen et al 2017), this product is also derived from WV2 data and classifies the Barrow area according to the US Fish and Wildlife Service National Wetlands Inventory Code. For this study, we reclassified the map to presence/absence data for sedge cover, with classes PEM1A-PEM1F indicating sedge cover (where PEM stands for palustrine area with emergent vegetation and the 1A-1F relates to the duration of annual inundation). From these datasets, four indices were calculated to characterise spatial heterogeneity (figure 2). These included elevation (Elev), derived from the LiDAR DEM, sedge cover, and the NDVI, equation (2) and NDWI, equation (3), (Gao 1996), which are band ratios measuring greenness and wetness respectively calculated from the WV2 image.
For each of the four spatial indices (Elevation, NDWI, NDVI, and Sedge Cover), we weighted the values of the index (Iˆt) on the spatial grid with the footprint probabilities w t (equation (5)), at each halfhourly time step t of EC data This yielded a time series of the spatial indices I , (ˆ) which show the changes in ecosystem properties sampled by the tower footprint as it responded to shifts in wind speed and direction. For the metrics of NDVI  and NDWI,  it should be noted that only a single image was available, so we assume the values reflect a pattern of spatial variability which stays consistent over time.
The footprint-weighted spatial indices I (ˆ) were compared with their mean value in the 0.36 km 2 domain around the base of the tower. This domain was chosen because we aimed to examine whether tower footprints were representative of area from which RS data would typically be extracted, as a general assumption for EC towers is that the majority of flux occurs within a radius of 100 times the height of the tower (Burba and Anderson 2010).Thus, one would normally expect the footprint values and the local domain to be very similar. As a summary statistic, we calculated the sensor location bias (sigma, d) developed by Schmid and Lloyd (1999) and Chen et al (2011) Figure 2. Four remotely sensed indices were used to define spatial heterogeneity. They are shown for the two sites, BES and BEO, over the extent of the local domain, 0.36 km 2 , around the base of the EC towers. The 50% and 80% isolines are shown for the full period and each season, to illustrate how the footprint shifts over time. where Iˆis the footprint-weighted spatial index, and I is the mean value of that index over the 0.36 km 2 domain. Therefore, δ characterises the difference between the spatial properties that are sampled by the flux tower and the local mean.

Statistical analysis
We used linear regression modelling to evaluate the explanatory power of the spatial and temporal variables on CH 4 flux. First, we included all terms (i.e. the temporal variables and the spatial indexes) into a full model of the CH 4 fluxes. The temporal variables used were air temperature (T A ), soil temperature (T S ), SWC, PAR, atmospheric stability (Zol), friction velocity (u * ), air pressure (P A ), and NDVI from the MODIS satellite (NDVI MODIS ). While the spatial variables were Elev,  NDVI,  NDWI,  and Sedge.  Model selection used a stepwise regression based on the Akaike information criterion, a widely used goodness of fit criteria. The adjusted R-squared was used to assess the model's explanatory power. To test collinearity in the explanatory variables, the variance inflation factor (VIF) was also calculated and variables that were collinear (VIF>10) were not included in the same model. The variables retained in the model would thus be significant with regards to influencing CH 4 flux variability. Analysis of variance (ANOVA) was used to ensure that model iterations of the model where statistically different from one another. The same model selection process described here for the full model was applied in all subsequent models. Furthermore, we separated the data by site and season, as well as, modelling CH 4 fluxes, using only either the temporal or spatial variables.
Additionally, to more fully isolate the impact of spatial variability in the footprint (i.e. difference between true temporal variability and spatial variability), a two-step analysis was performed as follows. Firstly, we fitted a linear model of CH 4 flux as function of the temporal variables. This thereby accounted for the temporal variation in the fluxes. We then fitted a model to the residuals of this linear model, including only the footprint-weighted spatial indices to examine whether some of the remaining variation could be accounted for by the spatial heterogeneity.
We also examined seasonal and site influences on flux by including them as factors in the full model. Both season and site were significant in explaining the variability in CH 4 fluxes. Three seasonal periods were defined based on the CH 4 flux rates. The snow melt period began with the first non-zero CH 4 flux measurement of the year and ending when NDVI MODIS reached a value of 0.3, which corresponds to a productive grassland (Didan 2015). From that point the growing season period lasted until NDVI MODIS again went below 0.3. Then came the post-growing season, which consisted mostly of the 'zero curtain' period, when soil remains unfrozen around zero degrees, was characterised by stable and substantial CH 4 emissions rates, though lower than those during summer .

Spatial variability and sensor location bias
The footprint modelling showed that there was substantial variability in the footprint-weighted spatial indices (figure 3). The means of the footprint-weighted spatial indices deviated substantially from the means for the 0.36 km 2 domain; the magnitude of the sensor location bias sigma, d, varied among sites and the indices considered, ranging from |18|% to |51|% (table 1, figure 3). BES was more inundated than BEO, with generally higher sedge cover. BES presented a more defined bi-modal distribution of the spatial variables (figure 3), consistent with two very different ecosystem types measured: a drained lake basin (where the EC tower was located), surrounded by polygonised upland tundra. A wet meadow prevails east of the EC tower within the drained lake basin and is characterised by a lower elevation and higher sedge cover. A drier ecosystem prevails west of the tower, with a higher elevation and lower sedge cover (figure 2). The δ was between |18|% and |32|% for nearly all the spatial variables, except for NDWI  at BES, which showed a δ of more than |50|% (table 1). The lowest δ was in the sedge cover and elevation in BES (table 1).

Temporal and spatial controls on CH 4 fluxes
Using a linear model, we could explain 60% of the CH 4 flux variability for the full methane flux dataset (table 2). This model retained the temporal variables, soil temperature (T S ), SWC, and PAR, as well as, the spatial indexes, NDVI  and NDWI.  The model including only spatial variables, was able to explain 33% of the variability in CH 4 emissions; when the model was limited to temporal variables, 53% of the variability in CH 4 emissions was explained, with soil temperature (T S ) having the largest explanatory power (51%) (table 2). In the residual analysis and when the temporal variation in the CH 4 fluxes was accounted for, the footprint-weighted spatial indices were able to explain 7% of the remaining variability in CH 4 fluxes. In all periods, the single temporal variable with the most explanatory power was soil temperature (T S ), except for the snow melt period at BEO, where NDVI MODIS performed better. Generally, when sites and seasons were modelled separately, the explanatory power of the spatial indices doubled for most time periods to ∼20% (table 2). During the curtain period at BES the spatial indexes seemed to have little relevance (table 2). There was, however, little consistency in terms of which spatial index had explanatory power and they often only achieved moderate explanatory power even when used together.

Discussion
The results of this study support the importance of footprint modelling in heterogeneous environments, as found in earlier studies (Budishchev et al 2014, Tuovinen et al 2018. In previous work, Wang et al (2016) used a threshold for δ of 10%, i.e. %d<|10|, when evaluating flux towers in the Canadian flux network, to determine if tower data would be suitable for upscaling fluxes to the regional scale. The lowest observed value of d in this study was |18|% spatial bias and in the case of NDWI, one of the spatial metrics that helped explain methane flux, had a bias of 51% (table 1). In another recent study, Treat et al (2018) found a 20%-65% underestimation of CH 4 emissions flux from a Siberian wetland site unless a high resolution wetland classification was used. While Tuovinen et al (2018), observed a somewhat contrasting result, in a Siberian shrub tundra site, where despite seeing a significant spatial bias (14%) and formally showing a 13% overestimation of methane flux for a 35.8 km 2 area, the results were not Figure 3. Variability in the spatial metrics sampled by the EC towers, due to temporal variation in the footprint data from 2013 to 2015, are shown by the grey bars. The y axis represents the frequency of footprint-weighted means of the spatial index over all halfhourly observations. The footprint-weighted mean is shown by the solid thin black line and the mean for the 0.36 km 2 around the base of the tower is shown by the dashed thick blue line. statistically significant. However, in this study they attributed fluxes to specific land cover classes and upscaled by using a single methane emission value for each class, that were calculated bases on measurements taken over one growning season. They acknowledged their results are heavily dependent on the landcover maps used and it could be a coincidence that the tower footprint was similar to the area to their area of interest. Additionally, they still concluded that in these heterogeneous sites, detailed footprint analysis is necessary, as they observed significant variability within their shrub tundra site. Our results suggest a similar over or underestimation of the CH 4 emissions could be generated if coarse-resolution data is used to upscale the EC data without accounting for footprint variability.
Our results indicate that variability in measured CH 4 flux can partially be attributed to changing footprint areas, in line with previous studies (Sachs et al 2010, Tagesson et al 2013. Overall, when site and season were accounted for, the explanatory power of the spatial indexes varied between 3% and 33% (table 2). During each season, different spatial metrics were better flux predictors for each of the sites. This highlights how polygonal tundra and the drained thermokarst are somewhat distinct habitats. In terms of seasonal variability, during the curtain period at BES the spatial indexes had relatively lower explanatory power, perhaps due to this area freezing up more uniformly. No one spatial index proved to be the most reliable CH 4 flux predictor, as the best results are often achieved when multiple indexes are used concurrently. Even the best spatial model left at least 40% of the variability in CH 4 fluxes unexplained. This could be linked to the complex influence of these variables on CH 4 release not necessarily being captured by the simple models used in this study (Sebacher et al 1983, Herbst et al 2011, Matthes et al 2014. For example, Göckede et al (2019) found that 3%-4% of methane emissions in their wetland study site in the Russian Arctic were in the form of sporadic bursts; a linear model might not be ideal to capture this type of emission pattern. Furthermore, the snow melt season it is a period of rapid transition, where NDVI and NDWI would vary greatly. Subsequent research might significantly improve the explanatory power of the spatial indices by utilising drones to collect high-resolution time series of these metrics. With appropriate development, these results could further develop upscaling methods and improve the extrapolation of the site-level data to the regional scale.

Conclusions
The overarching result from this study is to highlight the necessity of high-resolution footprint modelling when interpreting EC data from heterogeneous environments. One cannot assume that the tower footprint is representative of the even the immediate domain around the tower. Our analysis also shows that the area sampled by the towers differs from the surrounding local domain by |18|% to |51|%, depending on the spatial index used. The results of this study also show that spatial variability in the EC tower footprint in these Arctic wetlands sites accounts for a significant Table 2. Results from linear modelling. The 'Explanatory Variables' column indicates which variables were included as follows: temporal: air temperature (T A ), soil temperature (T S ), soil water content (SWC), photosynthetically active radiation (PAR), atmospheric stability (Zol), friction velocity (u * ), air pressure (P A ), NDVI from the MODIS satellite (NDVI MODIS ), Spatial: Elev,  NDVI,  NDWI,  Sedge  (half-hourly footprint weighted means), remotely sensed (RS): NDVI MODIS. Y is the dependent variable being modelled, either CH 4 emissions or the residuals from the temporal models (R T ) and the next two columns show the selected explanatory variables and the resulting adjusted r 2 value, for the indicated site, and season. All regressions shown are significant with p-values<0.05. amount of the half-hourly variance in the measured CH 4 flux. Models derived using high-resolution data should not be directly applied to lower resolution RS data, such as MODIS data, unless one accounts for footprint variability. The explanatory power of the spatial metrics varies, between 3% and 33%, depending on the site, time period and the modelling approach. Recognising this potential source of error is particularly valuable given that EC data are used to estimate regional and global CH 4 budgets.

Site
To test the larger scale applicability of the results, this analysis should be expanded to include data from additional years and a variety of tundra sites across the Arctic. One might also evaluate whether using different footprint models, such as Vesala et al (2008) or Kljun et al (2004), would help refine the footprint analysis. These steps would help develop an improved methodology for upscaling CH 4 fluxes in the Arctic.