BAWLD-CH₄: a comprehensive dataset of methane fluxes from boreal and arctic ecosystems

Methane (CH₄) emissions from the boreal and arctic region are globally significant and highly sensitive to climate change. There is currently a wide range in estimates of high-latitude annual CH₄ fluxes, where estimates based on land cover inventories and empirical CH₄ flux data or process models (bottom-up approaches) generally are greater than atmospheric inversions (top-down approaches). A limitation of bottom-up approaches has been the lack of harmonization between inventories of site-level CH₄ flux data and the land cover classes present in high-latitude spatial datasets. Here we present a comprehensive dataset of small-scale, surface CH₄ flux data from 540 terrestrial sites (wetland and non-wetland) and 1247 aquatic sites (lakes and ponds), compiled from 189 studies. The Boreal–Arctic Wetland and Lake Methane Dataset (BAWLD-CH₄) was constructed in parallel with a compatible land cover dataset, sharing the same land cover classes to enable refined bottom-up assessments. BAWLD-CH₄ includes information on site-level CH₄ fluxes but also on study design (measurement method, timing, and frequency) and site characteristics (vegetation, climate, hydrology, soil, and sediment types, permafrost conditions, lake size and depth, and our determination of land cover class). The different land cover classes had distinct CH₄ fluxes, resulting from definitions that were either based on or co-varied with key environmental controls. Fluxes of CH₄ from terrestrial ecosystems were primarily influenced by water table position, soil temperature, and vegetation composition, while CH₄ fluxes from aquatic ecosystems were primarily influenced by water temperature, lake size, and lake genesis. Models could explain more of the between-site variability in CH₄ fluxes for terrestrial than aquatic ecosystems, likely due to both less precise assessments of lake CH₄ fluxes and fewer consistently reported lake site characteristics. Analysis of BAWLD-CH₄ identified both land cover classes and regions within the boreal and arctic domain, where future studies should be focused, alongside methodological approaches. Overall, BAWLD-CH₄ provides a comprehensive dataset of CH₄ emissions from high-latitude ecosystems that are useful for identifying research opportunities, for comparison against new field data, and model parameterization or validation. BAWLD-CH₄ can be downloaded from https://doi.org/10.18739/A2DN3ZX1R (Kuhn et al., 2021).

1 Introduction

Methane (CH₄) is a strong climate forcing trace gas that is naturally produced and emitted from wetlands and lakes, which are abundant in northern regions (Matthews and Fung, 1987; Lehner and Döll, 2004; Messager et al., 2016). Current estimates of CH₄ fluxes from the northern boreal and arctic region (∼ > 50^∘) range between 9 and 53 Tg CH₄ yr⁻¹ from wetlands (Spahni et al., 2011; McGuire et al., 2012; Zhu et al., 2013; Bruhwiler et al., 2014; Treat et al., 2018; Watts et al., 2014; Thompson et al., 2018; Peltola et al., 2019; Saunois et al., 2020) and between 12 and 24 Tg CH₄ yr⁻¹ from lakes (Bastviken et al., 2011; Wik et al., 2016a; Tan et al., 2016; Walter Anthony et al., 2016; Matthews et al., 2020; Saunois et al., 2020). Combined, CH₄ emissions from northern ecosystems make up a significant but uncertain portion of fluxes from natural sources (232 to 367 Tg CH₄ yr⁻¹ for averaged bottom-up and top-down global estimates, respectively; Saunois et al., 2020). One reason for the large range of high-latitude CH₄ emissions estimates is the consistently lower estimates based on top-down approaches compared to bottom-up approaches. Top-down approaches use atmospheric observations of CH₄ concentrations with atmospheric inverse modelling frameworks to estimate regional CH₄ budgets (e.g., Bruhwiler et al., 2014; Thompson et al., 2018), while bottom-up approaches merge land cover datasets and empirical CH₄ flux inventories or process-based models to scale emissions across regional scales (e.g., Wik et al., 2016a; Treat et al., 2018; Peltola et al., 2019). A key issue for bottom-up approaches is the lack of differentiation among different wetland and lake types despite clear evidence indicating differences in both the magnitude and drivers of CH₄ fluxes among wetland and lake types (Olefeldt et al., 2013; Turetsky et al., 2014; Wik et al., 2016a; Treat et al., 2018).

Net CH₄ flux to the atmosphere depends on a suite of physical and biological controls linked to microbial production, oxidation, and transport via diffusion, ebullition, and plant-mediated processes (Bastviken et al., 2004; Whalen, 2005). While the basic underlying CH₄ processes are the same across all ecosystems, the dominance of different production, oxidation, and transport pathways varies within and among terrestrial (wetlands and non-wetlands) and lentic open-water aquatic ecosystems (lakes and ponds), leading to a wide range of reported CH₄ fluxes at the site level with differences of up to 4 orders of magnitude (Olefeldt et al., 2013; Wik et al., 2016a; Treat et al., 2018). Furthermore, drier terrestrial sites may drawdown, or uptake, CH₄ out of the atmosphere (Treat et al., 2018). Despite the wide range in reported CH₄ fluxes, key overarching controls on emissions from wetlands and aquatic ecosystems have been identified through the work of syntheses (Olefeldt et al., 2013; Wik et al., 2016a; Treat et al., 2018), suggesting that different ecosystems can be partitioned based on a handful of key CH₄-emitting characteristics.

For terrestrial ecosystems, CH₄ fluxes across the boreal–arctic region are primarily linked to permafrost conditions and hydrology (Olefeldt et al., 2013; Treat et al., 2018), which encompass other important controls on CH₄ emissions. For example, permafrost condition and hydrology can be directly linked to water table position and redox conditions (Moore et al., 1994; von Fischer et al., 2010; Olefeldt et al., 2017), which in turn influence plant composition (i.e., plant function types including graminoids, Sphagnum mosses, shrubs, and trees; Olefeldt et al., 2013; Bridgham et al., 2013), microbial community composition (McCalley et al., 2014), productivity (Christensen et al., 2003), and organic matter availability (Wagner et al., 2003; Christensen et al., 2003). Both permafrost condition and hydrology can further be used as an indication of soil temperature with typically colder conditions in drier soils and permafrost-dominated landscapes (Olefeldt et al., 2017). Methane fluxes are typically highest from graminoid-dominant wetlands such as marshes and fens, which are frequently inundated. Inundation, in turn, enhances primary productivity (Ström et al., 2012), creates a soil habitat conducive to CH₄-producing microbes (Woodcroft et al., 2018), and facilitates transport CH₄ through aerenchymatous roots and stems (Chanton et al., 1993; Ström and Christensen, 2007). Conversely, CH₄ fluxes are typically low from permafrost bogs and bogs which tend to have colder (in the case of permafrost bogs) and drier soil conditions (Beylea and Baird, 2006; Andersen et al., 2011), which are less conducive to the presence of graminoid species and promote the consumption of CH₄ through oxidation (Bartlett et al., 1992; Moosavi and Crill, 1997).

Methane fluxes from aquatic ecosystems (lakes and ponds) are highly influenced by lake morphology (Rasilo et al., 2015; Holgerson and Raymond, 2016) and lake genesis (Wik et al., 2016a), including underlying permafrost conditions (Walter et al., 2006), which are associated with other key controls and CH₄ fluxes. Lake morphology influences sediment temperature, macrophyte presence (Marinho et al., 2015; Wik et al., 2018), and turbulent transfer (MacIntyre et al., 2018). Lake morphology, permafrost condition, and lake genesis all determine organic substrate availability in sediments (Walter et al., 2006; Wik et al., 2016a) and trophic status (Bastviken et al., 2004; DelSontro et al., 2016). For example, peatland lakes and ponds, which form through degradation and permafrost thaw processes in peatlands, are relatively high CH₄ emitters (Matveev et al., 2016; Kuhn et al., 2018; Burke et al., 2019). These waterbodies are underlain by organic-rich sediments and are typically small and shallow and less likely to be seasonally stratified, allowing for rapid sediment warming and carbon mineralization (Matveev et al., 2016). Glacial and post-glacial waterbodies, on the other hand, have relatively low CH₄ fluxes due to deeper water columns, which limit ebullition by creating cooler sediment temperatures and greater hydrostatic pressures for bubbles to overcome (Bastviken et al., 2004; DelSontro et al., 2016). These waterbodies also tend to have mineral-rich sediments with typically less labile organic substrates (Schnurrenberger et al., 2003; DelSontro et al., 2016; Wik et al., 2016a). Therefore, while there are many physical and biogeochemical controls on aquatic CH₄ fluxes, size and lake genesis can be useful proxies for many of these underlying factors.

There are various methodologies used to measure surface CH₄ fluxes from terrestrial and aquatic ecosystems. Two approaches used in both terrestrial and aquatic ecosystems include micrometeorological eddy covariance (EC) techniques and chamber measurement techniques. Eddy covariance measurements are collected at high temporal frequencies from towers and typically cover a footprint of 100–10 000 m². The near-continuous nature of EC measurements provides valuable insight into the temporal patterns and drivers of CH₄ fluxes; however, towers are geographically limited across the boreal–arctic region, and it can be difficult to attribute flux transport pathways and specific source areas at fine spatial scales (Knox et al., 2019; Delwiche et al., 2021). Conversely, static chamber measurements cover small spatial areas that allow for detailed assessments of environmental controls on fluxes (Bäckstrand et al., 2008; Olefeldt et al., 2013). Chamber-based methods quantify fluxes by calculating the change in chamber headspace concentration over a set time, which varies based on extraction methods (i.e., syringe, automated chamber, or portable gas analyzer). While chamber-based techniques have drawbacks, including surface disturbance, typically low sampling frequency and high labour intensity, they are easily installed, can capture environmental controls of CH₄ fluxes at a sub-metre scale, and are cheaper options compared to installing and maintaining EC towers. Thus, we focus mostly on chamber-based flux measurements in this synthesis because they have been performed at many sites across the boreal and arctic region and represent more of the geographic variation across the region.

In aquatic ecosystems, turbulence-driven modelling approaches, inverted funnels (i.e., bubbles traps), and ice bubble surveys (IBSs) are additionally used to quantify fluxes. Modelling approaches calculate net hydrodynamic flux (herein referred to as diffusion) to the atmosphere by determining the concentration of dissolved CH₄ in the water column and an estimate of the gas transfer velocity k (see Sect. 2.4 for more information). Bubble traps capture the volume of bubble gas released from sediments; ebullitive flux can be estimated by using the concentration of CH₄ found in the bubble (Wik et al., 2013). Finally, IBSs are used to quantify the spatial abundance and types of bubble formations trapped within lake ice over the winter (Walter Anthony et al., 2010). Importantly, these surface-based methods can be used to assess controls of CH₄ exchange at scales of individual ponds, lakes, and portions of open-water wetlands, providing key insights into the environmental processes controlling CH₄ flux to the atmosphere (Olefeldt et al., 2013; Wik et al., 2016a).

Here we expanded, updated, and merged previous CH₄ flux syntheses for northern wetlands (Olefeldt et al., 2013) and lakes (Wik et al., 2016a) to create a small-scale (sub-metre), surface-based dataset for CH₄ fluxes collected from 189 studies across the boreal–arctic region. The dataset was built in parallel with a CH₄-specific land cover dataset for the circumpolar north – the Boreal–Arctic Wetland and Lake Dataset (BAWLD; Olefeldt et al., 2021a, b) – allowing for flux observations and spatial distribution of land cover features to be classified under the same criteria at a pan-arctic scale. This dataset includes surface-based fluxes and associated environmental drivers from aquatic, wetland, and upland (i.e., non-wetland) ecosystems and can be utilized by both field researchers and the modelling community. The boreal–arctic region represents a potentially globally significant, but still highly unknown source of CH₄. This dataset can be used to help constrain boreal–arctic flux estimates, compare field results, identify new research opportunities, or build and test models. This dataset includes and uniformly classifies lake, wetland, and upland (non-wetland) surface CH₄ flux data for the circumpolar north (Kuhn et al., 2021). We show CH₄ flux distributions and environmental drivers from various terrestrial (wetland and upland) and aquatic ecosystems across the north, compare the results to previous CH₄ flux syntheses, highlight key gaps in the data, and suggest future research directions.

2 Dataset description and BAWLD land cover classification

The dataset is composed of two parts, including (1) terrestrial ecosystems (vegetated wetland and non-wetland ecosystems) and (2) lentic open-water aquatic ecosystems (lakes, ponds, and open-water pools, hereafter referred to as “aquatic ecosystems”). This synthesis does not include lotic systems (streams and rivers), which are already synthesized in Stanley et al. (2016). The datasets for terrestrial and aquatic ecosystems are reported as separate components due to differences between both the drivers of CH₄ fluxes and data collection methods. The terrestrial dataset extends the work by Olefeldt et al. (2013), who compiled CH₄ flux estimates for wetlands in the permafrost zones designated by Brown et al. (2002). Our dataset expands on this initial work to include flux data from non-permafrost and non-wetland sites throughout the arctic and boreal region (Olson et al., 2001) and flux data from studies between 2012 and February 2020. We updated the initial dataset to include separate entries for individual sites that reported flux and water table data for multiple years. We expanded the number of site year flux estimates in the original terrestrial dataset by 83 % and expanded the number of independent studies by 86 %, leading to a total of 555 warm-season (∼ May through October, depending on the location) flux estimates and 121 studies (Fig. 1a). The aquatic dataset extends the work by Wik et al. (2016a), which is a compilation of CH₄ flux data for lakes and ponds north of 50^∘ N. We expand on this initial work to include studies between 2016 and February 2020. Additionally, we updated the original aquatic dataset to include the within-lake location for ebullition measurements and the equation used to model the gas velocity coefficient k. We expanded the number of lakes in the dataset by 71 % and the number of studies by 66 %, summing to a total of 1251 lakes and 68 independent studies (Fig. 1b). Finally, each terrestrial and aquatic site was reclassified into a new land cover classification, further explained below.

Figure 1 Maps of the individual sites (orange circles) incorporated in BAWLD-CH₄. (a) Sites included in the terrestrial flux dataset. (b) Sites included in the aquatic flux dataset. The number of “sites” in the terrestrial dataset represents site years, which in some cases represent multiple years of data from one site or data from the same site reported by different studies. “Sites” in the aquatic dataset represent the reported average fluxes for one or multiple lakes. In some cases, studies reported one mean value for multiple lakes; therefore the number of lakes and the number of sites is not the same. ^* Boreal–arctic region boundary from Olson et al. (2001). Permafrost zones are from Brown et al. (2002). Continental shoreline base layers are from Wessel and Smith (1996).

2.1 The Boreal–Arctic Wetland and Lake Dataset

2.1.1 Land cover classification

Land cover classes in the Boreal–Arctic Wetland and Lake Dataset (BAWLD; Olefeldt et al., 2021a, b) were chosen and defined to enable upscaling of CH₄ fluxes at large spatial scales. As such, we aimed to include as few classes as possible to facilitate large-scale mapping while still having sufficient classes to allow separation among groups of ecosystems with similarities in hydrology, ecology, and biogeochemistry and therefore net CH₄ fluxes. The BAWLD land cover classification is hierarchical, with four upland classes, five wetland classes, seven lentic aquatic classes, and three lotic aquatic classes. As mentioned previously, fluxes from lotic ecosystems (streams and rivers) are not included in this dataset but are covered by Stanley et al. (2016).

2.1.2 Wetland classes

Wetlands are defined by having a water table near or above the land surface for sufficient time to cause the development of wetland soils (either mineral soils with redoximorphic features or organic soils with > 40 cm peat) and the presence of plant species with adaptations to wet environments (Canada Committee on Ecological (Biophysical) Land Classification et al., 1997; Jorgenson et al., 2001; Hugelius et al., 2020). Wetland classifications for boreal and arctic biomes can focus on either small-scale wetland classes that have distinct hydrological regimes, vegetation composition, and biogeochemistry or on larger-scale wetland complexes that are comprised of distinct patterns of smaller wetland and open-water classes (Glaser et al., 2004; Masing et al., 2010; Gunnarsson and Löfroth, 2014; Terentieva et al., 2016). While larger-scale wetland complexes are easier to identify through remote sensing techniques (e.g., patterned fens comprised of higher-elevation ridges and inundated hollows), our classification focuses on wetland classes due to greater homogeneity of hydrological, ecological, and biogeochemical characteristics that regulate CH₄ fluxes (Heiskanen et al., 2021).

Several boreal countries identify four main wetland classes, differentiated primarily based on hydrodynamic characterization: bogs, fens, marshes, and swamps (Canada Committee on Ecological (Biophysical) Land Classification et al., 1997; Masing et al., 2010; Gunnarsson and Löfroth, 2014). The BAWLD classification follows this general framework but further uses the presence or absence of permafrost as a primary characteristic for classification and excludes a distinct swamp class, yielding five classes: Bogs, Fens, Marshes, Permafrost Bogs, and Tundra Wetlands (see Figs. 2 and 3). The swamp class was omitted due to the wide range of moisture and nutrient conditions of swamps as well as the limited number of studies of swamp CH₄ fluxes. We instead included swamp ecosystems in expanded descriptions of Bogs, Fens, and Marshes. The presence or absence of near-surface permafrost was used as a primary characteristic to distinguish between Permafrost Bogs and Bogs and to distinguish Tundra Wetlands from Marshes and Fens. The presence or absence of near-surface permafrost is considered key for controlling CH₄ emissions given its influence on hydrology and for the potential of permafrost thaw and thermokarst collapse to cause rapid non-linear shifts to CH₄ emissions (Bubier et al., 1995; Turetsky et al., 2002; Malhotra and Roulet, 2015; Fig. 3). Finally, while some classifications include shallow (e.g., 2 m depth), open-water ecosystems within the definition of wetlands (Canada Committee on Ecological (Biophysical) Land Classification et al., 1997; Gunnarsson and Löfroth, 2014), we have included all open-water ecosystems without emergent vegetation within the lake classes (see below) due to the strong influence of emergent vegetation in controlling CH₄ emissions (Juutinen et al., 2003).

Figure 2 Conceptual diagram of the terrestrial land cover classes and their CH₄-emitting characteristics including permafrost conditions, hydrology, organic layer depth, and associated nutrient and vegetation characteristics. Numbers within the brackets represent the interquartile (IQR) flux ranges. Arrows are scaled based on mean flux values. See Sect. 3.2 for a detailed breakdown of terrestrial fluxes.

Figure 3 Definitions of the five wetland classes in BAWLD along axes of moisture regime and nutrient regime.

Bogs are described as ombrotrophic peatland ecosystems, i.e., only dependent on precipitation and snowmelt for water inputs. Peat thickness is at least 40 cm, with maximum thickness > 10 m. The peat profile is not affected by permafrost, although in some climatically colder settings there may be permafrost below the peat profile. Bogs are wet to saturated ecosystems, often with small-scale (< 10 m) microtopographic variability, with stagnant water and a water table that rarely is above the surface or more than 50 cm below the surface (Fig. 3). Bogs have low pH (< 5), low concentrations of dissolved ions, and low nutrient availability, resulting from a lack of hydrological connectivity to surrounding mineral soils. Vegetation is commonly dominated by Sphagnum mosses, lichens, and woody shrubs and can be either treed or treeless (Beaulne et al., 2021). Our description of Bogs also includes what is commonly classified as treed swamps, which generally represent ecotonal transitions between peatlands and upland forests (Canada Committee on Ecological (Biophysical) Land Classification et al., 1997).

Fens are described as minerotrophic peatland ecosystems, i.e., hydrologically connected to surrounding mineral soils through surface water or groundwater inputs. A Fen peat profile is at least 40 cm thick (Gorham, 1991), although maximum peat thickness is generally less than for bogs. The peat profile is not affected by permafrost. Fens are wet to saturated ecosystems, with generally slow-moving water (Fig. 3). Fens have widely ranging nutrient regimes and levels of dissolved ions depending on the degree and type of hydrological connectivity to their surroundings, ranging from poor fens to rich fens. Vegetation largely depends on wetness and nutrient availability, where more nutrient-poor fens can have Sphagnum mosses, shrubs, and trees, while rich fens are dominated by brown mosses, graminoids (sedges, rushes), herbaceous plants, and sometimes coniferous or deciduous trees (e.g., willows, birch, larch). Our description of Fens also includes what is commonly classified as shrubby swamps, which often are associated with riparian ecotones and lake shorelines.

Marshes are minerotrophic wetlands with dynamic hydrology and often high nutrient availability (Fig. 3). Vegetation is dominated by emergent macrophytes, including tall graminoids such as rushes, reeds, grasses, and sedges – some of which can persist in settings with > 1.5 m of standing water. Marshes are saturated to inundated wetlands, often with highly fluctuating water levels as they generally are located along shorelines of lakes or coasts, along streams and rivers, or on floodplains and deltas. It is common for marshes to exhibit both flooded and dry periods. Dry periods facilitate the decomposition of organic matter and can prevent the build-up of peat. As such, Marshes generally have mineral soils, although some settings allow for the accumulation of highly humified organic layers – sometimes indicating ongoing succession towards a peatland ecosystem. Salinity can vary depending on water sources, with brackish to saline conditions in some areas of groundwater discharge or in coastal settings.

Permafrost Bogs are peatland ecosystems, although the peat thickness in cold climates is often relatively shallow. Permafrost Bogs have a seasonally thawed active layer that is 30 to 70 cm thick, with the remainder of the peat profile perennially frozen (i.e., permafrost). Excess ground ice and ice expansion often elevate Permafrost Bogs up to a few metres above their surroundings, and as such, they are ombrotrophic and relatively well drained (Fig. 3). Permafrost Bogs have moist to wet soil conditions, often with a water table that follows the base of the seasonally developing thawed soil layer. Ombrotrophic conditions cause nutrient-poor conditions, and the vegetation is dominated by lichens, Sphagnum mosses, woody shrubs, and sometimes stunted coniferous trees. Permafrost Bogs are often interspersed in a fine-scale mosaic (10 to 100 m) with other wetland classes, e.g., Bogs and Fens. Common Permafrost Bog landforms include palsas, peat plateaus, and the elevated portions of high- and low-centre polygonal peatlands.

Tundra Wetlands are treeless ecosystems with saturated to inundated conditions, most commonly with near-surface permafrost (Fig. 3). Tundra Wetlands can have either mineral soils or shallow organic soils and generally receive surface or near-surface waters from their surroundings as permafrost conditions preclude connectivity to deeper groundwater sources. Vegetation is dominated by short emergent vegetation, including sedges and grasses, with mosses and shrubs in slightly drier sites. Tundra Wetlands have a lower maximum depth of standing water than Marshes due to the shorter vegetation. Tundra Wetlands can be found in basin depressions; in low-centre polygonal wetlands; and along rivers, deltas, lake shorelines, and on floodplains in regions of continuous permafrost. Despite the name, limited wetlands with these characteristics (hydrology, permafrost conditions, and vegetation) can also be found within the continuous permafrost zone in boreal and sub-arctic regions (Virtanen et al., 2016).

2.1.3 Upland and other classes

Upland and other classes in BAWLD – Glaciers, Rocklands, Dry Tundra, and Boreal Forests – have in common that they are neither wetlands nor aquatic ecosystems. Glaciers are assumed to have neutral CH₄ fluxes; however, to our knowledge there are no published studies with field data from the glacier surface. There are a handful of studies that highlight lateral CH₄ export and emission from glacial outflows and termini (Christiansen and Jørgensen, 2018; Burns et al., 2018; Lamarche-Gagnon et al., 2019); however due to both limited atmospheric flux measurements and information on the spatial distributions of termini features and difficulties in mapping their areas at the circumpolar scale, we did not include these fluxes. Fluxes from glacial outflows and streams are considered to be riverine fluxes, and our flux synthesis does not include riverine fluxes. Rocklands are also expected to have very low CH₄ fluxes (Oh et al., 2020), potentially with more frequent CH₄ uptake than release. No sites included in the database were described as Rocklands. There are five sites described as high polar desert or desert tundra, which were included as Dry Tundra sites (Emmerton et al., 2014).

The Dry Tundra class includes both lowland arctic tundra and alpine tundra, both treeless ecosystems dominated by graminoid or shrub vegetation. Dry Tundra ecosystems generally have near-surface permafrost, with seasonally thawed active layers between 20 and 150 cm depending on climate, soil texture, and landscape position (Heikkinen et al., 2004; van der Molen et al., 2007). Near-surface permafrost in Dry Tundra prevents vertical drainage, but lateral drainage ensures predominately oxic soil conditions. A water table is either absent or close to the base of the seasonally thawing active layer. Dry Tundra is differentiated from Permafrost Bogs by having thinner organic soil (< 40 cm) and from Tundra Wetlands by their drained soils (average water table position > 5 cm below soil surface).

Boreal Forests are treed ecosystems with non-wetland soils. Coniferous trees are dominant, but the class also includes deciduous trees in warmer climates and landscape positions. Boreal Forests may have permafrost or non-permafrost ground, where the absence of permafrost often allows for better drainage. Overall, it is rare for anoxic conditions to occur in Boreal Forest soils, and CH₄ uptake is prevalent, although low CH₄ emissions have been observed during brief periods during snowmelt or following summer storms (Matson et al., 2009) or conveyed through tree stems and shoots (Machacova et al., 2016). The Boreal Forest class also includes the few agricultural and pasture ecosystems within the boreal biome.

2.1.4 Aquatic classes

Lakes in BAWLD include all lentic open-water ecosystems (herein referred to as aquatic ecosystems), regardless of surface area and depth of standing water. It is common in ice-rich permafrost lowlands and peatlands for open waterbodies to have shallow depths, often less than 2 m, even when surface areas are up to hundreds of square kilometres in size (Grosse et al., 2013). While small, shallow open-water pools often are included in definitions of wetlands (Canada Committee on Ecological (Biophysical) Land Classification et al., 1997; Gunnarsson and Löfroth, 2014; Treat et al., 2018), we include them here within the lake classes as controls on net CH₄ emissions depend strongly on the presence or absence of emergent macrophytes (Juutinen et al., 2003). Further classification of lakes in BAWLD is based on lake size and lake genesis, where lake genesis influences lake bathymetry and sediment characteristics (Fig. 4). Previous global spatial inventories of lakes include detailed information on size and location of individual larger lakes (Downing et al., 2012; Messager et al., 2016) but do not include open-water ecosystems < 0.1 km² in size and do not differentiate between lakes of different genesis (e.g., tectonic, glacial, organic, and yedoma lakes). Small waterbodies are disproportionately abundant in some high-latitude environments (Muster et al., 2019), have high emissions of CH₄ (Holgerson and Raymond, 2016), and therefore require explicit classification apart from larger waterbodies. Furthermore, lake genesis and sediment type have been shown to influence net CH₄ flux from lakes (Wik et al., 2016a). In BAWLD we thus differentiate between large (> 10 km²), midsize (0.1 to 10 km²), and small (< 0.1 km²) lake classes and further differentiate between three lake types for midsize and small lakes: peatland, yedoma, and glacial lakes (Fig. 4).

Figure 4 Conceptual diagram of the aquatic land cover classes. Key differences between the three overarching lake genesis “types” and their CH₄-emitting characteristics are shown, including sediment type, permafrost conditions, and water column depth. Fluxes (interquartile ranges, IQRs) for each class size within the overarching types are shown above the lakes for both diffusive and ebullitive transport pathways. Arrows are scaled based on mean flux values. See Sect. 3.3 for a detailed breakdown of aquatic fluxes. Large Lakes are not shown.

Small and Midsize Peatland Lakes are described as lakes with thick organic sediments that are mainly found adjacent to or surrounded by peatlands or in lowland tundra regions with organic-rich soils. Small Peatland Lakes include the numerous small pools often found in extensive peatlands and lowland tundra regions, e.g., including the open-water parts of string fens and polygonal peatlands. Peatland Lakes generally form as a result of interactions between local hydrology and the accumulation of peat, which can create open-water pools and lakes (Garneau et al., 2018; Harris et al., 2020) but can also form in peatlands as a result of permafrost dynamics (Sannel and Kuhry, 2011; Liljedahl et al., 2016). As such, these lakes with thick organic sediments are often shallow and have a relatively low shoreline development index. Peatland Lakes typically have dark waters with high concentrations of dissolved organic carbon.

Small and Midsize Yedoma Lakes are exclusive to non-glaciated regions of eastern Siberia, Alaska, and the Yukon, where yedoma deposits accumulated during the Pleistocene (Strauss et al., 2017). Yedoma permafrost soils are ice-rich and contain fine-grained, organic-rich loess that was deposited by wind and accumulated upwards in parallel with permafrost aggradation, thus limiting decomposition and facilitating organic matter burial (Schirrmeister et al., 2013). Notable thermokarst features, including lakes, often develop when yedoma permafrost thaws, causing labile organic matter to become available for microbial mineralization (Walter Anthony et al., 2016). Small Yedoma Lakes are thus more likely to have actively thawing and expanding lake edges where CH₄ emissions can be extremely high, largely driven by hot spot ebullition emissions (Walter Anthony et al., 2016; Fig. 4). Century-scale development of yedoma lakes can shift the main source of CH₄ production from yedoma deposits to new organic-rich sediment that accumulated from allochthonous and autochthonous sources, resulting in such lakes here being considered as Peatland Lakes.

Small and Midsize Glacial Lakes include all lakes with organic-poor sediments – predominately those formed through glacial or post-glacial processes, e.g., kettle lakes and bedrock depressions. However, due to similarities in CH₄ emissions and controls thereof, we also include all other lakes with organic-poor sediments within these classes. Glacial Lakes typically have rocky bottoms or mineral sediments with limited organic content. Lakes in this class are abundant on the Canadian Shield and in Scandinavia but can be found throughout the boreal and tundra biomes. Many Glacial Lakes have a high shoreline development index, with irregular, elongated shapes. Generally, Glacial Lakes are deeper than lakes in the other classes when comparing lakes with similar lake areas and are more likely to stratify seasonally than peatland lakes (Fig. 4).

Large Lakes are greater than 10 km² in surface area. Most Large Lakes are glacial or structural and tectonic in origin. Lake genesis is not considered for further differentiation within this land cover class.

2.2 Terrestrial Methane Flux Dataset

The Terrestrial Methane Flux Dataset includes warm-season (∼ May–October depending on the location) fluxes and was compiled using data from studies published before February 2020. We identified relevant studies using (1) JSTOR™, Google Scholar™, and Web of Science™ searches with the terms (peatland OR wetland OR bog OR fen OR marsh OR upland) AND (north^* OR boreal OR arctic OR sub-arctic) AND (methane OR CH₄ OR greenhouse gas^*); (2) references from published studies; and (3) contributions of unpublished data (n=2). If multiple, yearly CH₄ flux and water table measurements were reported from one site or if multiple studies reported fluxes from the same site, the data were entered as separate individual lines and were considered each their own “site”. Sites that underwent manipulations (soil temperature, water table, nutrients, etc.) were not included in the dataset; however, any control or undisturbed sites included within manipulation studies were included. Sites that had recently experienced disturbance from thermokarst processes were included. Winter flux measurements from terrestrial sites were excluded from this dataset (winter and ice-out emissions from aquatic ecosystems are included; see Sect. 2.3). A comprehensive synthesis of seasonal winter estimates of CH₄ emissions from northern terrestrial ecosystems is presented in Treat et al. (2018).

The terrestrial dataset includes predominantly chamber measurements (n=519) at the sub-metre scale, which allows for a detailed representation of specific land cover classes (i.e., one land cover class per chamber measurement). However, a handful of eddy covariance measurements were included if the authors could clearly partition fluxes based on specific land cover classes (n=36). For more information on EC-based CH₄ synthesis, we direct the readers to the FLUXNET-CH₄ Community Product (Knox et al., 2019; Delwiche et al., 2021) and additional FLUXNET resources (https://fluxnet.org, last access: 15 April 2021). We grouped chamber measurements from specific studies by “site”, which we defined as surfaces with similar vegetation composition (dominant, present, absent) and physical characteristics (including water table position, permafrost conditions, organic layer depth, soil moisture, and pH) within proximity to each other (typically 1–100 m radius). In most cases, chambers and sites were already classified by these standards; however, sometimes it was necessary to combine or split chamber measurements presented by the authors into our site and classifications. By combining and splitting sites this way, we were able to classify sites into BAWLD land cover classes. Average daily warm-season fluxes were then calculated from the average CH₄ flux from each site over the study's measurement period.

In addition to CH₄ flux data, we extracted various site descriptors and categorical and continuous environmental variables (See Table 1 for detailed attribute information and additional variables not discussed here). For all sites, we included information on the site name (Site), location (LatDec/LongDec, Country), the months measurements were taken (SampMonths), the flux measurement method (Meth), the author's description of the site (SiteDescrip), and vegetation composition. Most studies did not classify land cover types with similar BAWLD criteria; therefore we assigned BAWLD land cover classifications. Permafrost zone was assigned according to Brown et al. (2002). When reported by the authors, we also extracted continuous variables including mean annual air temperature (MAAT), mean annual precipitation (MAP), growing season length, net ecosystem productivity (NEP), ecosystem respiration (ER), gross ecosystem photosynthesis (GPP_Per), air temperature (T_Per), soil temperature at 0–5 cm (TSoilA) and at 5–25 cm (TSoilB), water table depth (WT_Av), organic layer depth (Org), active layer depth (AL), pH, and soil moisture (SoilMoist), all averaged over the same period as the flux measurements. The categorical variables collected include absence or presence of permafrost within the top 2 m (PfConA), permafrost thaw (PfTh), and vegetation composition (absent, present, dominant) for graminoid (Carex spp. and Eriophorum spp.; referred to as “Sedge” in the dataset), Sphagnum moss (Sphag), non-Sphagnum moss (Moss), tree, and shrub species. Vegetation composition of the functional plant type was considered dominant if that type made up greater than 50 % of the reported biomass or areal coverage or was one of only two species present at the site. Trees were assigned as the dominant vegetation type if the canopy was described as closed. Gridded (0.5 by 0.5^∘) climate variables including mean annual temperature (referred to as GRID_T) and mean annual precipitation (CD_Pcp_An) were extracted from WorldClim2 (http://www.worldclim.com/version2, last access: 15 April 2021).

Table 1 Attribute information for the terrestrial flux dataset.

2.3 Aquatic Methane Flux Dataset

The Aquatic Methane Flux Dataset includes ice-free season (∼ May–October depending on the location) and winter and ice-out fluxes and was compiled using data from studies published before February 2020. We identified new studies using (1) JSTOR™, Google Scholar™, and Web of Science™ searches with the terms (lake^* OR pond^*) AND (north^* OR boreal OR arctic OR sub-arctic) AND (methane OR CH₄ OR greenhouse gas^*); (2) references from published studies; and (3) contributions of unpublished data (n=1). If multiple, yearly measurements were given for one site by the same study, we averaged the flux values (following the initial protocol taken by Wik et al., 2016a). If different studies reported fluxes from the same lake, then these data were reported as separate entries. In instances where ice-free season fluxes and storage and ice-out fluxes were reported for the same lake, those data were entered on separate lines, but the number of lakes was designated as NA for the winter measurement so as to not add to the total lake count. We defined sites based on reported average CH₄ fluxes. For example, some studies reported one average flux value for a group of lakes, and this was considered one “site”; however, the number of lakes was noted. Studies that only reported CH₄ concentrations and not a flux estimate were not included.

Like the terrestrial dataset, the aquatic dataset focuses on small-scale measurement techniques that allow for flux estimates to be attributed to one specific land cover class. Therefore, ice-free season diffusive fluxes included in this dataset were measured using dissolved CH₄ concentrations and modelling approaches (n=254) or floating chambers (n=181), while ebullitive fluxes were measured by bubble trap (n=187) or floating chamber (n=34). Diffusive modelling approaches include an estimate of the gas transfer coefficient, k. Gas transfer velocity estimates are commonly calculated using equations (e.g., Cole and Caraco, 1998). However, more recent efforts with EC systems, chambers, and either calculation or measurement of the near-surface turbulence that enables flux across the air–water interface indicated that fluxes calculated using the wind-based model of Cole and Caraco (1998) of gas transfer velocities underestimate fluxes from non-sheltered and sheltered waterbodies by a factor of 2 to 4 (Heiskanen et al., 2014; Mammarella et al., 2015; MacIntyre et al., 2020). Highly sheltered waterbodies, such as small lakes surrounded by trees, may be an exception and can have reduced mean lake k values (Markfort et al., 2010). While we do not recalculate fluxes in this synthesis, we indicate which k calculations were used so that future studies can easily identify and recalculate fluxes when required. Only a handful of eddy covariance (EC) measurements (n=5) were included in the dataset. We included a limited number of EC measurements due to difficulties that most studies had in attributing the fluxes to lakes specifically. We classified all EC fluxes as diffusive fluxes as it is hard to separate between ebullition and diffusion within this measurement technique; however, for this reason, EC measurements were excluded from statistical analysis for ice-free season fluxes.

We further delineated aquatic fluxes by transport pathway including ebullition (bubbles), diffusion (hydrodynamic flux), and winter storage and ice-out flux. Ebullition and diffusion measurements were averaged over the ice-free season to represent a mean daily flux estimate across a lake. In some cases, if measurements were only taken from one zone of the lake (i.e., just lake edge or just lake centre) we averaged the fluxes and assumed whole-lake fluxes. Some studies only reported a seasonal ice-free flux estimate. If they also reported the number of days in the ice-free season, we then calculated the average daily flux rate. Storage and ice-out flux includes the annual release of CH₄ that accumulates within and under the ice over the winter and is released during spring turnover and includes estimates from ice bubble surveys (IBSs). Our storage flux estimate does not include estimates of fall circulation fluxes, wherein CH₄ that is stored in the deep portion of the water column is released upon seasonal turnover of the water column (Karlsson et al., 2013; Sepulveda-Jauregui et al., 2015). We also include an estimate of the ice-free season ebullition and diffusive fluxes if provided by the authors or if the authors provided the number of ice-free days. Note that flux measurements that include the transport of CH₄ through littoral vascular plants were not included as aquatic fluxes, but as Marsh or Tundra Wetland fluxes within the terrestrial dataset.

In addition to aquatic CH₄ flux data, we also collected various site descriptors and categorical and continuous environmental variables (see Table 2 for detailed attribute information and additional variables not discussed here). For all sites we extracted information about the site name and location (latitude and longitude and country), the number of lakes for a reported flux estimate, sampling season (SEASON) and within-lake sampling location (E.LOCATION), sampling pathway (PATHWAY), the general sampling dates (YEAR/MONTH), and the number of times sampled (D.DAYS/E.DAYS). When available, we added a column for the equation used to estimate the gas transfer velocity constant (k) using modelling approaches (K600_EQ). Categorical variables included lake sediment type (BOTTOM), permafrost zone (PERMA.ZONE), presence of talik (TALIK), ecoregion (ECOREGION), and the original lake types outlined by Wik et al. (2016a) (LAKE.TYPE). BAWLD-specific categorical variables include the overarching lake genesis type (TYPE), binned waterbody size (SIZE), and BAWLD land cover class (CLASS). BAWLD land cover classes were assigned based on author descriptions of the waterbodies. If the authors did not provide information indicating the lake type, we used the coordinates provided to find the waterbody on Google Earth™ and used yedoma permafrost (Strauss et al., 2017) and organic soil maps (Hugelius et al., 2014) to determine the land cover class. In a handful of cases, the land cover class could not be determined, and we left the “class” field blank. When reported, we extracted the following continuous variables: surface area (SA), waterbody depth (DEPTH), water temperature (TEMP), dissolved organic carbon concentration (DOC), and pH. Gridded (0.5 by 0.5^∘) climate variables including mean annual temperature (GRID_T) and mean annual precipitation (CD_Pcp_An) were extracted from WorldClim2 (http://www.worldclim.com/version2, last access: 15 April 2021).

Table 2 Attribute information for the aquatic flux dataset.

2.4 Statistics

All statistical analyses were performed in R statistical software (Version 1.1.383; https://www.r-project.org, last access: 10 April 2021). We tested for significant relationships between log-transformed warm-season (terrestrial sites) or ice-free season (aquatic sites) average CH₄ fluxes and several covariates using a combination of linear regression and linear mixed-effects models when necessary (R Package 3.3.3; Lme4 Package; Bates et al., 2014). To include sites with CH₄ uptake or near-zero fluxes we added a constant of 10 (terrestrial fluxes) or 1 (aquatic fluxes) before log transformation. Mixed-effects modelling was used when a given model included sites with multiple yearly measurements or if multiple studies reported fluxes from the same site (R “nmle” package; Pinheiro et al., 2017). In these cases, site ID was included as a random effect in the analysis to help account for lack of independence across repeated measurements and to weight potential biases (as in Treat et al., 2018). Almost no studies in the terrestrial or aquatic datasets provided information on all the variables; therefore, individual statistical analyses have different sample sizes. However, the same subset of data was used to select the best-performing mixed models (n=206 and n=149 for the terrestrial and diffusive aquatic mixed models, respectively). The significance of individual predictor variables in the mixed models was evaluated using forward model selection. Model performance was conducted using the size-corrected Akaike information criterion (AICc; “AICcmodavg” package; Mazerolle and Mazerolle, 2015), wherein a decrease in AICc by 2 or more is an indication of a superior model (as in Olefeldt et al., 2013; Dieleman et al., 2020). All models were tested against each other and the null model. The null model only included the random effects. Non-parametric Tukey's honestly significant difference (HSD) post hoc tests were performed to assess differences in median fluxes among sub-categories if the overall model was determined significant. All aquatic diffusive and ebullitive fluxes were analyzed separately. Eddy covariance CH₄ flux estimates for aquatic ecosystems (n=5) were not included in the statistical analysis since ebullitive and diffusive fluxes could not be partitioned. We modelled the temperature dependence (Q₁₀) of CH₄ fluxes following Rasilo et al. (2015).

2.5 Limitations

Due to limitations of the studies where we extracted data from, some parts of the annual period are not considered in our dataset. Thus, this dataset focuses on small-scale, surface-based spatial patterns in CH₄ fluxes associated with specific land cover classes and does not represent temporal patterns in fluxes. For both terrestrial and aquatic datasets, we extracted data on the average CH₄ fluxes over warm periods or ice-free periods. While we do include an estimate of ice-out and winter fluxes from aquatic ecosystems, our dataset does not include autumnal turnover fluxes from aquatic ecosystems, which may represent a substantial portion of annual emissions (Fernández et al., 2014; Klaus et al., 2018). Nor do we include shoulder season or winter fluxes from terrestrial ecosystems, which can represent substantial components of the annual flux (Zona et al., 2016; Treat et al., 2018). Furthermore, our data extraction methods were not designed to assess inter-annual changes in fluxes as this dataset compiles the data of multiple studies over a large range of years (1986–2020). Despite data limitations, the datasets represent an important step forward regarding the spatial variability in fluxes among different land cover types.

3 Results

3.1 Summary statistics

In total, we extracted 555 site year CH₄ estimates from terrestrial (wetland and non-wetland) ecosystems. Most reported fluxes (site years) were from Canada and Greenland (34 %), followed by Russia (27 %), Alaska (25 %), and Scandinavia (14 %) (Fig. 1a). Terrestrial fluxes followed a bimodal distribution, split by net positive fluxes (82 % of all reported fluxes) and net uptake or zero-emission (18 % of all reported fluxes; Fig. 5a). The median number of measurement days per site year flux for chamber measurements was 10, and the median number of collars per site measurement was 5 (Fig. 6a). Of the site year fluxes reported from aquatic ecosystems, there were 441 diffusive estimates and 175 ebullitive ice-free season estimates as well as 125 estimates of winter and ice-out fluxes (including storage, winter ebullition, ice bubble surveys, or a combination of the three). Aquatic sites were distributed throughout the boreal–arctic region with a greater density of sites in Alaska and eastern Canada (Fig. 1b). Diffusive fluxes showed a unimodal distribution, while ebullition showed bimodal peaks near 100 and 0 mg CH₄ m⁻² d⁻¹ (Fig. 5b, c). The median number of measurement days per site year flux was 3 and 15 for diffusion and ebullition, respectively (Fig. 6b, c). Winter and ice-out fluxes were reported as annual estimates and are shown in Table 5.

Figure 5 Histograms of site-specific average CH₄ fluxes. (a) Terrestrial fluxes. (b) Aquatic diffusive fluxes. (c) Aquatic ebullitive fluxes. Grey bars represent net zero or net uptake fluxes.

Figure 6 Histograms for the number of sampling days contributing to the average warm-season or ice-free season flux value. (a) Terrestrial flux sampling days. (b) Aquatic diffusion flux sampling days. (c) Aquatic ebullition sampling days. The dotted orange lines in (b) and (c) represent the number of recommended sampling days needed to arrive at a flux estimate within 20 % accuracy (11 d for diffusion and 39 d for ebullition; Wik et al., 2016b). The dotted red lines represent an updated estimate of the number of sampling days needed including 14 and 135 d for diffusion and ebullition, respectively (Jansen et al., 2020).

3.2 Correlations with terrestrial fluxes

Of the continuous variables, water table (WT_Av) and soil temperature (TSoilA at 5–25 cm) were significantly and linearly correlated with CH₄ (WT_Av: χ²=121, P<0.0001, R²m=0.28, df = 380; TSoilA: χ²=54.6, P<0.0001, R²=0.21, df = 283), and gross primary productivity (GPP) was logarithmically correlated with CH₄ (χ²=5.8, P=0.016, R²m=0.15; df = 56; Fig. 7). However, given the relatively low sample size for GPP (n=57), we do not include GPP in mixed model analyses. The temperature sensitivity (Q₁₀) for all terrestrial emissions was 2.8 (Table S1 in the Supplement). Of the categorical variables, there was no difference between the different permafrost zones (χ²=0.88, P=0.83, df = 539), but CH₄ fluxes were higher from sites without permafrost present in the top 2 m (χ²=16.37, P<0.0001, df = 482; Fig. 8). For vegetation composition, sites dominated by shrubs had lower fluxes than those sites with shrubs present or absent (χ²=34.66, P<0.001, df = 2; Fig. 8). The strongest relationship between vegetation composition and CH₄ flux was emergent graminoid cover. Sites with dominant graminoid composition had higher fluxes than sites where graminoids were present or absent (χ²=148.95, P<0.0001, df = 2; Fig. 8). The best explanatory model for terrestrial CH₄ emissions was an additive model that included site-level predictors of water table, soil temperature, and graminoid cover alongside the broader classification of land cover class (R²m=0.69; P<0.0001, df = 224; Table S2). There was no effect on model performance using interactive effects (DeltaAICc = 0.84); however, the R²m did increase to 0.73 (Table S2). Notably, on their own, individual models with just the site-level predictors or with just land cover type explained close to the same amount of variation in CH₄ fluxes (R²m=0.55 and 0.54, respectively). Methane uptake fluxes, when analyzed separately, were positively correlated with thaw depth (i.e., more uptake with greater thaw depths; R²m=0.55, χ²=19.61, P<0.0001, df = 22; Fig. S1 in the Supplement). No other continuous variables were correlated with CH₄ uptake; however, sites where shrubs were present had significantly higher uptake than sites where shrubs were absent or dominant (Tukey's post hoc, P<0.001 for both, df = 2; Fig. S2).

Figure 7 Relationships between site-averaged warm-season CH₄ flux and environmental variables. Environmental variables include water table, soil temperature at 2–25 cm depth, active layer depth, latitude, air temperature, organic layer thickness, ecosystem respiration (ER), gross primary productivity (GPP), and soil pH. Regression lines and R² values are shown for significant relationships. Note the log scale. CH₄ flux was linearly related to water table and soil temperature and was logarithmically related to GPP (dotted line). Points below the dotted red line represent net uptake fluxes. ^* P<0.05, ${}^{**}$ P<0.01, ${}^{***}$ P<0.001.

Figure 8 Warm-season CH₄ fluxes classified by categorical variables. Orange circles represent mean flux values. The number of sites for each category is represented in the column to the right (n), and statistical differences among the categories are indicated by the letters (Sig), wherein bars with the same letters are not significantly different. Permafrost zones are from Brown et al. (2002). Permafrost condition represents the presence of permafrost in the top 2 m as reported by the authors. See text for definitions used to classify vegetation cover. Outlier fluxes greater than 380 are not shown.

There were significant differences in fluxes among the BAWLD terrestrial land cover classes (χ²=253.69, P<0.001, df = 6; Fig. 9a, Table 3). Median fluxes were highest from Marshes, Tundra Wetlands, and Fens (mean water table = +2, −0.4, and −6 cm, respectively). Median fluxes from Bogs were lower than the Marshes, Tundra Wetlands, and Fens but higher than Permafrost Bogs, Dry Tundra, and Boreal Forests. Permafrost Bogs were the only wetland class that fell into the lowest-emitting group of classes. However, the frozen and elevated nature of Permafrost Bogs typically leads to lower-water-table conditions more like Dry Tundra and Boreal Forests (mean water table = −22, −15, and −40 cm, respectively). However, it must be noted that in most Boreal Forest sites the water table is not in the top 2 m; therefore water table is not commonly measured or reported. The mean water table depth presented here is likely an overestimate that represents wetter Boreal Forest sites that had measurable water tables in the top 2 m. Boreal Forest ecosystems were the only class to have negative median CH₄ flux for the entire class (net uptake). Permafrost Bogs and Dry Tundra classes also included net uptake site year CH₄ estimates (n=17 and 31, respectively). One Wetland Tundra site in the Canadian High Arctic had net CH₄ uptake for 1 of the 3 years it was measured (Emmerton et al., 2014). Notably, the apparent temperature sensitivity from the drier terrestrial sites (Boreal Forest, Dry Tundra, and Permafrost Bogs: Q₁₀=3.7) was higher than from the wet terrestrial sites (Marshes, Tundra Wetlands, Bogs, and Fens: Q₁₀=2.8).

Figure 9 Relationship between methane flux and BAWLD land cover classes. (a) Terrestrial fluxes per class. (b) Aquatic fluxes including diffusion and ebullition per class. Orange dots represent the arithmetic mean flux values, and black lines represent median flux values. Boxes represent 25th and 75th percentiles. Outlier fluxes over 380 are not shown. The letters represent significant differences in fluxes among classes, wherein bars with the same letters are not significantly different.

Table 3 Characteristics of BAWLD terrestrial classes based on environmental variables. The number of sites (site years) and contributing studies are shown for each class. Also shown are the mean, median, and quartiles for site-average CH₄ flux, water table, soil temperature between 5 and 25 cm (TSoilB), sedge cover, pH, ecosystem respiration (ER), and gross primary productivity (GPP). ^* In some cases, one study contributed flux data for multiple classes.

3.3 Correlations with aquatic fluxes

Diffusive CH₄ fluxes from aquatic ecosystems were negatively correlated with the continuous variables logged surface area (χ²=73.0, P<0.0001, R²m=0.20, df = 235; Fig. 10a), logged waterbody depth (χ²=23.5, P<0.0001, R²m=0.09, df = 275; Fig. 10b), and latitude (F=54.6, P<0.0001, R²=0.13, df = 361; Fig. 10c) and positively correlated with DOC (F=71.7, P<0.0001, R²=0.21, df = 261; Fig. 10d) and water temperature (F=57.2, P<0.001, R²=0.19, df = 236; Fig. 10e). The apparent Q₁₀ for diffusive emissions was 4.3 (Table S1). Diffusive CH₄ fluxes were highest from the sporadic permafrost zone (χ²= 17.2, P=0.002, df = 3; Fig. 11). Furthermore, diffusive fluxes were significantly higher from small lakes compared to midsize and large lakes (χ²=30.5, P<0.0001, df = 2; Fig. 11) and from lakes with peaty or organic-rich sediments compared to lakes with Yedoma and Glacial sediment types (χ²=103.9, P<0.0001, df = 2; Fig. 11). The best explanatory model for diffusive CH₄ fluxes was an additive model including an interaction between lake surface area (continuous) and type (i.e., overarching lake genesis) alongside water temperature as predictor variables (F=14.9, P<0.0001, adj.R²=0.41, df = 149; Table S3). Land cover class on its own explained 25 % of the flux variation (F=22.8, P<0.0001, df = 149).

Figure 10 Relationships between site-averaged ice-free diffusive and ebullitive CH₄ fluxes (note the log scale) and environmental variables. Environmental variables include surface area, waterbody depth, latitude, dissolved organic carbon (DOC) concentration, water temperature, and pH. Regression lines and R² values are shown for significant relationships. Log diffusive CH₄ flux was linearly related to surface area, depth, latitude, water temperature, and DOC. Log ebullitive fluxes were linearly related to surface area, latitude, DOC, and water temperature. ^* P<0.05, ${}^{**}$ P<0.01, ${}^{***}$ P<0.001.

Ebullitive CH₄ fluxes from aquatic ecosystems were positively correlated with logged DOC (F=12.25, P=0.0008, adj.R²=0.14, df = 71; Fig. 10d), negatively correlated with surface area (F=13.88, P=0.0003, adj.R²=0.08, df = 164; Fig. 10a) and latitude (F=5.38, P=0.02, adj.R²=0.03, df = 160; Fig. 10c), and weakly correlated with water temperature (F=5.55, P=0.02, adj.R²=0.06, df = 67; Fig. 10e). The apparent Q₁₀ for ebullitive emissions was 2.4 (Table S1). There was no apparent relationship with lake depth and ebullitive fluxes (F=0.02, P=0.91, df = 151; Fig. 10b). There were no differences in ebullitive emissions between the permafrost zones with the exception of lower ebullitive emissions from the continuous zone compared to the sporadic zone (Tukey's HSD, P < 0.001; Fig. 11). Like diffusive fluxes, ebullitive fluxes were higher from the small lake classes compared to midsize lakes (Wilcoxon rank sum, P=0.0006; note that Large Lakes did not have a large enough sample size (n=1) to be included in the post hoc analysis). Finally, ebullitive fluxes were similarly higher from Peatland and Yedoma lakes compared to Glacial Lakes (Tukey's HSD, P=0.006 and 0.001, respectively). The best explanatory model for ebullitive fluxes using a subset of the data with complete information for predictor variables of interest (i.e., SA, log.CH4.E.FLUX.plus1, SITE, CLASS, SIZE, DOC, TYPE, LAT, GRID_T) included just waterbody surface area (continuous) as a predictor variable (F=19.85, P=0.0001, adj.R²=0.21, df = 68).

Figure 11 Ice-free season diffusion (left) and ebullitive (right) CH₄ fluxes as described by categorical variables. Orange circles represent mean flux values. The number of sites for each category is represented in the column to the right (n) in the representative colours for diffusion (light blue) and ebullition (dark blue). The letters (sig) indicate statistical differences among the categories, wherein bars with the same letters are not significantly different. Lake size represents binned surface areas for < 0.1 km² (small), 0.1–10 km² (midsize), and > 10 km² (large). Lake type represents the BAWLD classification of waterbody types including Peatland, Yedoma, and Glacial Lakes. Fluxes higher than 380 are not shown.

There were clear differences in diffusive CH₄ fluxes among the aquatic class types, but few differences were observed for ebullitive fluxes. Diffusive fluxes were higher from the Peatland and Yedoma Lake classes (both small and midsize), associated with organic-rich sediments, compared to mineral-rich glacial and large lakes (χ²=119.8, P<0.001, df = 6; Fig. 9b, Table 4). While ebullition fluxes appear to follow a similar trend to diffusive fluxes, the only significant difference was between Small Yedoma Lakes and Midsize Glacial Lakes (Tukey's HSD, P<0.001; Fig. 9b). However, the lack of statistical differences found for ebullition between lake classes may in part be due to fewer and more variable ebullition measurements compared to diffusion (Table 4). Reported winter ice-out emission estimates (including storage flux and ice bubble survey (IBS) flux) were scarce in comparison to reported ice-free season emissions. Small Glacial Lakes and Midsize Glacial Lakes had the most reported winter ice-out emission estimates (n=20 and 31, respectively). Average winter emissions (storage flux + IBS) generally were lower than annual estimates of ice-free diffusive and ebullitive emissions (Table 5); however, statistical tests were not performed across all the classes due to low sample sizes from some of the classes. Winter ebullition estimates (i.e., direct ebullition emission to the atmosphere from seeps during the ice-cover winter season) were not included in winter emission sums because of the non-uniform spatial nature of these emission types (Sepulveda-Jauregui et al., 2015; Wik et al., 2016a) but are shown in Table 5. In the future, more estimates of winter emissions from aquatic systems are needed to more accurately estimate total annual emissions.

Table 4 Characteristics of the BAWLD aquatic classes based on CH₄ and environmental variables. The number of sites and contributing studies are shown for each class and flux pathway. Also shown are the mean, median, and quartiles for site-average diffusive and ebullitive CH₄ flux, waterbody surface area, waterbody depth, and dissolved organic carbon concentrations (DOC). ^* In some cases, one study contributed flux data for multiple classes and pathway types. One ebullition outlier point (flux = 1815 mg CH₄ m² d⁻¹) was excluded from the Midsize Glacial class as it was influenced by beaver activity (Sepulveda-Jauregui et al., 2015).

Table 5 Winter fluxes, including storage, ice bubble storage (IBS), and winter ebullition for each class type. Annual estimates of ice-free diffusion and ebullition are included for comparison. ^* Winter ebullition from constant seeps not included in sum of winter and ice-out emissions.

3.4 Joint analysis of terrestrial and aquatic fluxes

We performed joint analysis of fluxes from both the aquatic and terrestrial datasets with regional predictor variables (Class, MAAT, MAP, Permafrost Zone, and Biome) using mixed models to assess the potential for universal drivers across all boreal–arctic ecosystems. The best model included class and MAAT (χ²=345.6, P<0.0001, R²m=0.47, df = 18: Table S4). However, class alone explained 44 % of the variation in fluxes (compared to 47 % in the best model; Table S4), suggesting that ecosystem classification based on CH₄-emitting characteristics, alongside corresponding spatial extent, is one of the most important variables to consider when scaling CH₄ fluxes across the boreal–arctic region.

4 Discussion

4.1 Flux variation largely explained by land cover classes

In this review, we assessed the controls on CH₄ emissions from 189 studies across terrestrial and aquatic ecosystems in the boreal–arctic region. A central component to this study was the inclusion of new land cover classes split by CH₄-emitting characteristics common across terrestrial and aquatic ecosystems, respectively. Terrestrial classes were split by permafrost conditions and hydrology (and vegetation and nutrient conditions therein), whereas aquatic classes were split by size and lake genesis (i.e., type). We found that much of the observed CH₄ flux variability from terrestrial and aquatic ecosystems could be explained by this land cover classification system (Fig. 9). When modelling fluxes for both aquatic and terrestrial ecosystems together with regional-level predictors (variables assigned to sites based on the gridded product including biome, permafrost zone, MAAT, and MAP), land cover class explained most of the variation (44 %), with significant but small contributions in explained variation from gridded MAAT (3 % of 47 % total variation explained; Table S2). This suggests that spatial differences in land cover classes are the most important consideration for estimating CH₄ flux at this scale, with some influence of MAAT.

For terrestrial fluxes alone, land cover class as a predictor variable explained 55 % of the flux variation. Site-level predictors, including water table, temperature, and vegetation conditions, explained 54 % of the variation in the fluxes when analysed separately. The best model for terrestrial fluxes included these site-level variables and land cover class and explained 69 %–73 % of the variation (depending on additive or interactive effects; Table S2). This model likely performed better than land cover class on its own because the extra information added from the continuous soil temperature and water table variables captured the variation in these conditions within each class. While permafrost presence came out as a non-significant term in our best model (Table S2), the effects of permafrost presence and absence, including confounding temperature effects, were already intertwined into the land cover classes.

For aquatic ecosystems, the best models for diffusive and ebullitive fluxes contained different predictor variables. The best model for diffusive fluxes explained 41 % of CH₄ flux variability and included an interactive effect between surface area and lake type (Peatland, Yedoma, and Glacial) and water temperature. Land cover classes (i.e., lake types split by small and midsize categorical sizes) did not come out as significant in this model because the continuous variable of surface area captures the size variation within each lake type. However, land cover class modelled on its own explained 25 % of the flux variation. The significant effect of surface area is consistent with previous global synthesis efforts that found that small waterbodies tend to have higher CH₄ fluxes, likely due to the compounding effects of higher substrate availability and warmer temperatures compared to larger waterbodies (Holgerson and Raymond, 2016; DelSontro et al., 2018). Notably, previous synthesis efforts also found that waterbody depth was a significant predictor variable of diffusive fluxes (Wik et al., 2016a; Li et al., 2020). While depth did not come out as significant in our model, the effect of waterbody depth is taken into account with the lake types. For example, we found that diffusive fluxes are typically higher from Peatland Lake types compared to Glacial Lakes, which have average depths of 1.6 and 6.7 m, respectively. Waterbody depth is also an important factor contributing to waterbody temperature (i.e., warmer waters in shallower waterbodies); thus the effect of waterbody depth may also be confounded with that of the temperature variable.

The best model for ebullition contained waterbody surface area as a predictor and explained 21 % of the variation in the fluxes. Previous synthesis efforts have linked ebullition fluxes to both temperature (Aben et al., 2017) and waterbody depth (Wik et al., 2016a). There are a few potential explanations as to why we did not find similar relationships between ebullition and temperature or waterbody depth. First, Aben et al. (2017) include global data that encompass sites across broad temperature ranges from the north to the tropics. It is possible that the range of temperatures represented by our dataset is not wide enough to capture this relationship. It is also possible that the summary data collected, including average temperature and average flux over the ice-free season, are too coarse to show a relationship. It is likely that temperature and depth influence is clearer over time and space in each respective waterbody and that a higher resolution of data would show these relationships. Regarding waterbody depth, it is also possible that in the absence of detailed surveys, estimated mean and max depths may be less reliable. The effects of depth may also be confounded with surface area as the two metrics are highly correlated (Fig. S5). While this dataset represents one of the largest collections of ebullitive emissions from northern lakes so far, this emission pathway is still largely underrepresented, and waterbody depth and temperature are not always reported with the flux estimates. Furthermore, we collected information on surface water temperature for this dataset because it was the most widely available temperature metric. Sediment temperature is a better metric to collect in hand with ebullition due to production and transport directly from the sediments (Wik et al., 2013; Aben et al., 2017). Future studies should work to report sediment temperature and water column temperature alongside their flux measurements.

Figure 12 Geographical distribution and flux frequencies and for each terrestrial class. Relative land cover for each type is represented in green on the map. Site locations are represented by orange circles. Note the log scale for CH₄ flux. Land cover distributions from Olefeldt et al. (2021a, b). Histograms of non-transformed flux data can be found in Fig. S3.

4.2 Directions for future research

While our small-scale, surface CH₄ flux datasets for northern ecosystems are the most extensive datasets compiled to date for the boreal–arctic region, we identified key gaps in the data and areas of improvement that future studies should focus on. While the geographical gaps represented in Fig. 1a suggest widespread geographic under-representation of terrestrial ecosystems, especially across central Russia and the Canadian territories of Nunavut and Northwest Territories, these regions are comprised primarily of Boreal Forest and Dry Tundra ecosystems, respectively (Fig. 12e, g). Study sites for many of the other land cover types, for example, Bogs and Fens, were relatively well distributed across the boreal and arctic region (Fig. 12a, b). However, to assess how well or poorly represented a land cover class is, class area and flux magnitude must also be considered (Fig. 13a). For example, Fens are a high-emitting land cover class and are spatially abundant, leading to a high total flux contribution across the study region (∼ 41 %, Fig. 13a); however, the relative number of Fen sites represented in the available literature is not proportional to the total flux contribution (∼ 26 %). This, alongside the large spread of reported flux magnitudes (Fig. 9a), suggests that future studies should focus on Fens to better constrain the flux magnitude. Conversely, Permafrost Bogs are low contributors to the total wetland flux (∼ 4 %), and sites are well represented throughout the literature (∼ 19 %), suggesting that fewer direct flux measurements are needed from these ecosystems.

Figure 13 Relative total flux contribution (mean flux ⋅ total class area) for each land cover shown with the relative contribution of flux measurements for each class. (a) Wetland classes. (b) Aquatic classes. The bars represent the percent of total flux contribution and percent of reported flux sites for each class. Aquatic flux contributions represent average ebullition + average diffusion fluxes.

For aquatic ecosystems, there are key data gaps in geography and flux pathway representation with relatively few measurements of ebullition and ice-out fluxes compared to diffusive fluxes. Geographically there are very few flux measurements from Midsize Glacial Lakes and Large Lakes in the western Canadian Shield (Fig. 14a, d), despite this region containing the most lakes per unit area throughout the north (Messager et al., 2016). Notably, Large Lakes are the least represented of all the aquatic classes (∼ 2.4 % of measurements) but could potentially contribute ∼ 17 % of the total flux, mostly from diffusive emissions. Interestingly, while Small Peatland Lakes are well represented (∼ 42 % of measurements and 37 % of potential total flux contribution), Midsize Peatland Lakes are under-represented (∼ 5 % of measurements) compared to their estimated flux contribution (∼ 28 %). Thus, Large Lakes and Midsize Peatland Lakes may be important focal points for future research; however, more empirical scaling-based uncertainty analyses should be explored.

There are fewer ebullition measurements compared to diffusive flux measurements from aquatic ecosystems (21 % and 79 % of ice-free fluxes, respectively). Average ebullitive fluxes were greater than diffusive estimates for all the land cover classes except Large Lakes (Fig. 7b) and thus represent an important component of total CH₄ fluxes from these systems; however, none of the models performed in this study could explain a large amount of the variation in ebullitive fluxes. More ebullition measurements, across all the land cover classes, will help to constrain our understanding of CH₄ transport mechanisms and drivers. However, it is important to note that more representative ice-free season flux estimates are needed for both ebullition and diffusion. Wik et al. (2016b) suggest that ∼ 11 diffusive day flux measurements and ∼ 39 ebullition day flux measurements are required to calculate a mean ice-free flux estimate within 20% of the true value. A total of 86 % of diffusive estimates were under the recommended 11 d mark, and 58 % of ebullition estimates were below the recommended 39 d mark (Fig. 5b, c). Jansen et al. (2020) posit that an even higher frequency of sampling is required (14–22 and 135 d for diffusion and ebullition, respectively). Further, Wik et al. (2016b) recommend that in addition to the number of sampling days, measurements should be distributed spatially across the waterbody using a depth-stratified approach including ∼ 3 and ∼ 11 locations for diffusion and ebullition, respectively. While we did not collect data on the number of sampling locations across each waterbody, it is likely that many of the average fluxes included the dataset also represent spatially undersampled measurements. Undersampling potentially reduces the accuracy of mean CH₄ flux estimates, leading to the relatively poor fitness and explanatory power of the aquatic regression analysis in this study compared to the terrestrial models' performances. This is especially true for ebullitive emissions, which were poorly explained by the reported predictor variables available for this dataset. Calculation of mean ice-free fluxes from a greater number of flux measurements is an important step forward towards better constraining CH₄ fluxes from aquatic ecosystems. Finally, there are very few flux estimates for lakes over the shoulder seasons and winter and ice-out compared to the ice-free season (Table 5). While shoulder season flux estimates, including autumnal turnover, were not included in this dataset, only 7 % of the lakes in this synthesis had winter and ice-out measurements collected. Winter and ice-out emissions could potentially contribute a significant portion of annual fluxes from aquatic ecosystems (Karlsson et al., 2013; Sepulveda-Jauregui et al., 2015) and therefore represent an important gap in CH₄ flux data.

5 Data availability

The BAWLD-CH₄ flux dataset is available for download at the Arctic Data Center at https://doi.org/10.18739/A2DN3ZX1R (Kuhn et al., 2021). The companion land cover spatial dataset is also available at the Arctic Data Center at https://doi.org/10.18739/A2C824F9X (Olefeldt et al., 2021a).

6 Conclusions

Methane fluxes from northern ecosystems represent an important component of the global CH₄ cycle (Saunois et al., 2020). BAWLD-CH₄ is a comprehensive flux dataset that uniquely represents flux data from both terrestrial and aquatic ecosystems across the boreal–arctic region. BAWLD-CH₄ has many potential applications including benchmarking for process-based models, use in empirical scaling models, and informing future research directions. Importantly, we show that land cover class, split CH₄-emitting ecosystem characteristics, is a significant flux predictor variable across terrestrial and aquatic ecosystems, and we suggest that future studies should scale CH₄ emissions based on CH₄-emitting land cover characteristics. We show that while land cover class explains most of the flux variation for wetland and aquatic ecosystems when analyzed jointly, MAAT significantly explains ∼ 3 % of the variation, which has important implications for future scaling efforts. Finally, we found that a higher percentage of terrestrial CH₄ fluxes could be explained by land cover class and site-level variables than for diffusive and ebullitive fluxes from aquatic ecosystems (73 % vs. 41 % and 21 %, respectively). Undersampling of aquatic ecosystems is likely responsible for the lower explained variation observed in our models compared to terrestrial ecosystems. Future studies should increase the number of sampling days for both diffusive and ebullitive fluxes from aquatic ecosystems to arrive at more representative ice-free flux estimates and total CH₄ emissions from the boreal–arctic region.

Figure 14 Flux frequencies and geographical distribution for each aquatic class. Relative land cover for each class type is represented in blue on the map. Site locations are represented by orange circles. Note the log scale for CH₄ flux. Land cover distributions from Olefeldt et al. (2021a, b). Histograms of non-transformed flux data are shown in Fig. S4.

Appendix A

Table A1 Table of references in the dataset.

NA – not available.