Fine-scale spatial variability of winter CO₂ and CH₄ fluxes in Arctic tundra derived from snowpack gradient measurements

Winter carbon dioxide (CO₂) and methane (CH₄) fluxes from soils under seasonal snowpacks make non-negligible, yet poorly constrained contributions to annual carbon budgets across Arctic regions. Quantifying these fluxes and their spatial variance will better constrain uncertainties in simulations of winter carbon fluxes from terrestrial biosphere models. We address this gap by measuring and identifying patterns and spatial variability in CO₂ and CH₄ soil-atmosphere fluxes through late winter snowpacks at an upland tundra site in the western Canadian Arctic. Instantaneous fluxes were calculated from CO₂ and CH₄ microsite concentration gradients at 10 to 20 cm vertical resolution (n=119) through the snowpack across five homogeneous surface covers, representing dominant vegetation types. We measured consistent soil-to-atmosphere CO₂ fluxes but with significantly different rates across surface covers (0.8 to 100 mgC m⁻² d⁻¹), which were strongly influenced by snow depth and soil surface temperature, exhibiting higher emissions under deeper snowpacks and warmer soil surfaces. CH₄ fluxes were also coupled to soil surface temperature and varied between −0.04 and 0.08 mgC m⁻² d⁻¹. Persistent CH₄ uptake was observed in warmer soils (−6.0 to −0.5 °C) in a sparsely populated black spruce and shrub dominated area with deep snow, indicating active methane oxidation during winter. Microsite scale CO₂ and CH₄ fluxes were statistically independent of vertical snow microstructure, indicating that winter fluxes could be reliably calculated from single gas concentration at the soil-snow interface, negating the need for additional snowpack gas measurements. These results open new avenues for quantifying fine-scale spatial variability of wintertime CO₂ and CH₄ fluxes in Arctic tundra, which can constrain biogeochemical process representations in terrestrial biosphere models and inform spatial upscaling methodologies.

1 Introduction

Arctic environments are experiencing up to four times the rate of climate warming than the global average (Rantanen et al., 2022), with greater rates of warming during the cold season months of September to May (Bekryaev et al., 2010; Smith et al., 2019). Longer growing seasons and earlier springtime snowmelt drive changes in vegetation composition and structure (Piao et al., 2019; Shi et al., 2015), such as deciduous shrubs expanding northwards and increasing in size and density (i.e., Arctic greening), especially in the forest-tundra ecotone (Ju and Masek, 2016; Myers-Smith et al., 2011, 2020). Altered Arctic precipitation patterns (Ernakovich et al., 2014) and taller, more plentiful shrubs offer more resistance to wind, favouring the trapping of blowing snow and impacting processes that control the seasonal evolution of snowpack thickness and microstructure (Shirley et al., 2025). This influences the insulating properties of snow on soils throughout winter, here defined as the snow covered period of the cold season (Busseau et al., 2017). Warmer winter soils promote increased microbial activity and more rapid degradation of soil organic matter through heterotrophic respiration. As soil microbes are especially sensitive to temperature and to the availability of liquid soil moisture near and below freezing, warmer wintertime soils could alter emissions of potent greenhouse gases (GHGs) such as carbon dioxide (CO₂) and methane (CH₄) (Natali et al., 2019), with the direction of the net CH₄ response set by the balance between methanogenesis and methanotrophic oxidation under prevailing soil moisture and oxygen availability. Understanding how surface-atmosphere carbon exchanges in different Arctic ecosystems impact net GHG budgets is key to assessing warming potential at a global scale (Ludwig et al., 2024).

Terrestrial Biosphere Models (TBMs) simulate exchanges of energy, water, and carbon at the Earth's surface, and allow us to simulate the effects of warmer air temperatures and land cover change on regional and global climate. TBMs as used in Earth System Models, such as the Community Land Model (CLM5.0) (Lawrence et al., 2019) are used to integrate the effects of dynamically evolving snowpacks on soil thermal and moisture regimes, coupled with biogeochemical processes. Large uncertainties in simulations of cold season carbon fluxes remain, largely because of the poor representation of the effective conductivity of snow and inadequate parametrisation of soil heterotrophic respiration and limited northern latitude GHG flux measurements for model evaluation (Fisher et al., 2014; Natali et al., 2019; Tao et al., 2021).

Simulations urgently require improvement as measurements of cold season GHG fluxes in Arctic ecosystems have been shown to contribute substantially to annual carbon budgets. Field measurements from northern Alaska demonstrated that CO₂ emissions during the snow-covered season can be large enough to offset or exceed growing-season uptake (Fahnestock et al., 1998, 1999; Jones et al., 1999b; Welker et al., 2000), and recent syntheses reinforce these findings at pan-Arctic scales (Natali et al., 2019; Rafat et al., 2021). Relatively small CH₄ emissions during the lengthy cold season have also been shown to contribute > 50 % of annual total sources in Alaskan tundra (Zona et al., 2016). Collectively, such studies established winter as a critical period for Arctic carbon cycling and motivated subsequent efforts to quantify cold-season fluxes across a wider range of Arctic and sub-Arctic ecosystems.

Snow covered landscape wintertime GHG exchanges between soil and atmosphere have been measured using several techniques, yet published methods sacrificed spatial coverage (single sampling sites), spatial resolution by aggregating fluxes over length scales larger than those over which processes vary, or temporal coverage by aggregating flux measurements over several measurement campaigns. Eddy covariance measurements (Baldocchi, 2003) have increasingly been used to measure surface-atmosphere fluxes in the Arctic over footprints of 100 s to 1000 s m² (Virkkala et al., 2022). However, eddy covariance flux systems aggregate fluxes at scales orders of magnitude larger than localised (< 1 m²) soil microbial processes, masking surface cover scale variability driven by localised conditions (i.e., soil temperature, available liquid soil moisture, labile nutrient availability, root zones, and discrete organic inputs such as faecal deposits) (Soegaard et al., 2000; Sullivan et al., 2020; Vourlitis et al., 2000). Conversely, soil chambers provide process scale microsite vertical fluxes by measuring changes in gas concentration within an enclosed volume at the ground-atmosphere interface (Subke et al., 2021). Where wintertime chamber-based CO₂ snowpack flux measurements have been made, they generally underestimate fluxes as the chamber-method excludes lateral diffusive flux driven by wind within the snow porous structure, or require chambers to be installed before snowfall, inevitably changing soil thermal and moisture regimes during the winter (McDowell et al., 2000; Webb et al., 2016). Soil-atmosphere fluxes derived from direct measurements of snowpack CO₂ and CH₄ concentrations and steady-state diffusion models have been used for several decades to estimate winter greenhouse gas fluxes in Arctic ecosystems. Manually extracting air from snowpacks with syringes allows for cost effective and spatially distributed flux measurements. However, large sample sizes (n>10) are difficult to collect within time scales where meteorological conditions are stationary and wind speeds are low, and require several samples to be extracted and analysed per microsite for flux calculations (Jones et al., 1999a; Mavrovic et al., 2023; Seok et al., 2009; Zhu et al., 2014). Faster, high accuracy in situ snowpack gradient measurements of CO₂ flux (and later CH₄) were more recently made possible using portable Infrared Gas Analysers (IRGA), quantifying winter soil-atmosphere fluxes across a variety of Arctic land surface covers (Fahnestock et al., 1998; Jones et al., 1999b; Pirk et al., 2016).

High spatial variability in winter carbon fluxes is due to complex interactions between microtopography, snow structure, and soil moisture. Microtopographic features, such as hummocks and hollows, influence the distribution of snow depth and insulative properties, decoupling soil thermal regimes from cold winter air (Burn and Kokelj, 2009; Quinton and Marsh, 1998; Schädel et al., 2024; Wilcox et al., 2019). Strong soil-air thermal gradients and wind driven snow redistribution in drifts and vegetation result in locally distinct vertical snowpack layering. Dominant snow layers such as wind crusts and depth hoar further affect the insulative capacity and diffusivity of gasses within a snowpack (Domine et al., 2016; Proksch et al., 2015). Additionally, shrubs protruding through snowpacks act as thermal conduits passively cooling soils at their roots, locally diminishing the insulative capacity of snow (Domine et al., 2022). As winter soil thermal regimes vary spatially with snow structure and distribution, winter soil microbial respiration and nutrient mineralization likewise vary (Aanderud et al., 2013; Schimel et al., 2004; Sullivan et al., 2020; Wilkman et al., 2018), because snowpack characteristics control soil temperature and moisture microclimates that regulate microbial activity (Mikan et al., 2002). Wet, low-lying areas with persistent soil saturation are more likely to be hotspots for methanogenesis, whereas drier upland areas promote aerobic methane oxidation (Virkkala et al., 2024; Voigt et al., 2023). However, the spatial distribution of these processes in TBMs remains poorly constrained due to limited high-resolution measurements across tundra ecosystems, especially within frozen soils during winter (Fisher et al., 2014; Levy et al., 2022). Temporally continuous and spatially distributed measurements of winter carbon fluxes remain a necessity to evaluate wintertime soil microbial heterotrophic respiration simulations as a function of ground thermal and moisture regimes (Dutch et al., 2022, 2024; Schädel et al., 2024).

Here we examine the spatial variability of CO₂ and CH₄ wintertime fluxes across local scales (100 m transects) in an Arctic tundra site and develop methodologies for capturing representative mean emission estimates for wider upscaling. Three objectives are to (1) quantify CO₂ and CH₄ flux variability within homogenous surface covers from a broad range of microtopographic and vegetation microsites, (2) quantify and compare errors on measured fluxes using a new snowpack gas concentration gradient measurement approach to current values found in the literature, and (3) identify the required number of samples needed to adequately capture variability within homogenous surface covers for surface cover based spatial scaling. Meeting these objectives, we provide a dataset for constraining uncertainties in biogeochemical processes simulated in TBMs, and statistical insight informing future sampling protocols aimed at improving spatial upscaling from in situ flux measurements and for TBM validation datasets.

2 Methods

2.1 Theoretical framework for gas flux calculation through the snowpack

Winter microbial respiration in soils drive vertical fluxes in GHGs across the soil-atmosphere interface by the production or consumption of CO₂ and CH₄ within the soil column. Snow acts as a porous intermediary medium between soil and atmosphere for gases to traverse, creating a vertical concentration gradient driven by molecular diffusion. Assuming steady state and concentration gradients across the snowpack are linear, then the rate of diffusion of gas (F_gas, mol m⁻² s⁻¹) can be calculated using Fick's 1st law of diffusion as:

$$\begin{array}{}\text{(1)}& F_{\mathrm{gas}}=-\mathit{\phi }\mathit{\tau }D\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}z}}\right)N_{\mathrm{a}}\end{array}$$

where D (m⁻² s⁻¹) is the standard air diffusion coefficient of each gas species (here, $D_{{\mathrm{CO}}_{\mathrm{2}}}$ and $D_{{\mathrm{CH}}_{\mathrm{4}}}$ for CO₂ and CH₄, respectively) through the snowpack, dC $/$ dz (mol m⁻¹) is the vertical gas concentration gradient across the snowpack of height z, φ (m³ m⁻³) is snowpack porosity across dz, τ (m m⁻¹) is snowpack tortuosity across dz, and N_a (mol m⁻³) is the molar volume used to convert units of concentration (from ppm to mol m⁻³). According to the ideal gas law (PV = nRT, where R=8.314 J K⁻¹ mol⁻¹), N_a ($n/V$) is calculated as:

$$\begin{array}{}\text{(2)}& N_{\mathrm{a}}={\displaystyle \frac{P}{RT}}\end{array}$$

where P (Pa) is ambient air pressure and T (K) is the average temperature of the snowpack. Snow porosity is derived from measured snow density along its vertical profile following Kinar and Pomeroy (2015):

$$\begin{array}{}\text{(3)}& \mathit{\phi }=\mathrm{1}-\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{snow}}}{{\mathit{\rho }}_{\mathrm{ice}}}}\right)\end{array}$$

where no liquid content is assumed to be present in the snowpack, ρ_snow is measured snow density and ρ_ice is pure ice density as a function of measured temperature at standard pressure (${\mathit{\rho }}_{\mathrm{ice}}=-\mathrm{0.0001}$ × T+0.9168 g cm⁻³; Haynes, 2016), where T is temperature in °C. Tortuosity through a porous medium is derived from Du Plessis and Masliyah (1991):

$$\begin{array}{}\text{(4)}& \mathit{\tau }={\displaystyle \frac{\mathrm{1}-{\left(\mathrm{1}-\mathit{\phi }\right)}^{\mathrm{2}/\mathrm{3}}}{\mathit{\phi }}}\end{array}$$

Standard CO₂ and CH₄ diffusion coefficients must be corrected for air temperature and pressure (Massman, 1998):

$$\begin{array}{}\text{(5)}& D=D_{\mathrm{0}}\left({\displaystyle \frac{P_{\mathrm{0}}}{P}}\right){\left({\displaystyle \frac{T}{T_{\mathrm{0}}}}\right)}^{\mathrm{1.81}}\end{array}$$

where P₀ and T₀ are standard atmosphere pressure and temperature respectively (P₀=101.325 kPa and T₀=273.15 K), and D₀ is gas diffusion in air at those standards where: $D_{\mathrm{0},{\mathrm{CO}}_{\mathrm{2}}}=$ 0.1381 cm² s⁻¹ (1.381 × 10⁻⁵ m² s⁻¹); $D_{\mathrm{0},{\mathrm{CH}}_{\mathrm{4}}}=$ 0.1952 cm² s⁻¹ (1.952 × 10⁻⁵ m² s⁻¹).

Non-diffusive transport driven by wind pumping is represented through an effective diffusion coefficient altering Eq. (1) to:

$$\begin{array}{}\text{(6)}& F_{\mathrm{gas}}=-\left(\mathrm{1}+\mathit{ϵ}\right)\mathit{\phi }\mathit{\tau }D\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}z}}\right)N_{\mathrm{a}}\end{array}$$

where ϵ is a dimensionless enhancement factor that is zero for pure molecular diffusion and increases with snow-surface wind speed, snowpack air permeability, and barometric variance (Massman et al., 1997; Bowling and Massman, 2011).

Turbulent and barometric pressure pumping is oscillatory with near-zero mean flow, so the net effect is a dispersive enhancement of diffusion rather than a mean advective flux (Massman et al., 1997; Bowling and Massman, 2011), whereas quasi-static (Venturi) advection depends on the wavelength and height of surface relief, such as sastrugi or drifts, and on the resulting surface pressure perturbation (Colbeck, 1997; Massman et al., 1997). Bowling and Massman (2011) directly measured ϵ with values reaching zero in calm air and up to 0.33 at a wind speed of 10 m s⁻¹. Modelling of a layered snowpack driven by such fluctuations at mean wind speeds of 1 to 4 m s⁻¹ shows that pressure pumping does not appreciably enhance through-snowpack diffusion unless it alters the boundary concentrations (Massman et al., 1997). Furthermore, a lower-permeability surface layer such as a wind crust further attenuates the penetration of pressure-driven airflow into the snowpack (Colbeck, 1997). Advective enhancement is therefore bounded near zero (ϵ≈ 0) under which conditions the flux reduces to its diffusive form as Eq. (1) when wind speed is low, the snow surface lacks pronounced relief, and a wind crust is present. As antecedent wind-generated pressure gradients were not directly measured within snowpacks over periods of hours to days before measurement of concentration gradients, their impact on basal snowpack CO₂ concentrations (e.g. Massman and Frank, 2006) cannot be excluded completely. Nevertheless, the contribution to barometric pressure pumping from small atmospheric pressure changes during the measurement period are expected to be negligible.

2.2 Study site

Measurements were undertaken in a 3 km² area surrounding the Wilfrid Laurier University Trail Valley Creek Research Station, (TVC; 68°74^′ N, 133°50^′ W), approximately 45 km north of Inuvik, in the Inuvialuit Settlement Region of western Inuit Nunangat, Northwest Territories, Canada (Fig. 1). TVC has a mean annual precipitation of 210 ± 9 mm for the 1991–2020 period, with approximately two thirds of annual precipitation falling as snow. Mean annual air temperature for that same period is −8.1 ± 0.7 °C (Environment and Climate Change Canada, 2024). Dominant westerly wintertime winds produce stable snow redistribution patterns across years trapping snow in areas with taller shrubs, and into drifts on eastward facing slopes where large eddies form forcing snow deposition (Wilcox et al., 2019). Typically, snow covers persist from mid-October to mid-May depending on accumulation patterns (Wilcox et al., 2019; Pomeroy et al., 1993), with seasonal maximum snow depths ranging from 45 ± 10 cm in riparian shrub and sparsely wooded areas to 71 ± 15 cm in upland tundra and tussock sites (Tutton et al., 2025).

Topography is characterised by gently rolling hills of glacial till, overlain by silty clay and a thin organic top layer. The area is underlain by 350 to 500 m thick ice rich continuous permafrost (i.e., > 90 % in areal extent; Gruber, 2012), generally with active layer thickness varying between 0.5 and 0.8 m (Burn and Kokelj, 2009) and up to 1.5 m on some southward facing white and black spruce (Picea glauca, Picea mariana) populated slopes and below stream channels. Landscape microtopography is dominated by mineral soil hummocks of up to a metre in diameter and peaty inter-hummock hollows with some periglacial tundra landforms and associated cover types including ice-wedge polygons and barren ground (Quinton and Marsh, 1998).

Vegetation is composed of low stature shrub tundra composed of dwarf birch and willow (Betula nana, Betula glandulosa; S. arctica), taller alders surrounding drainage channels (Alnus sp.), graminoids and cotton grass tussocks (Eriophorum spp.), patches of lichen and mosses, with some areas of sparsely populated black spruce (Graveline et al., 2024). Shrub-dominated areas are seen as greener areas, and low-lying lichen and graminoid-dominated vegetation in lighter tones (Fig. 1).

Figure 1 Overview of sample locations in Trail Valley Creek (study site), NWT Canada. Icon colours indicate snow depth at each sampling microsite across sampling transects (Table 1). Land surface imagery from © Google, © 2025 Airbus Maxar Technologies, elevation above sea level contour lines derived from a LIDAR-based digital elevation model produced for Wilfrid Laurier University (KBM Resources Group, 2024).

2.3 Data Collection and Sampling design

All data were collected late winter between 20–26 March 2024 near maximum annual snow water equivalent, but before any snow surface melt had taken place and air temperatures remained below −10 °C. We identified five sampling transects, at length scales of approximately 100 m (T1 to T5; Fig. 1), within different vegetation types representative of the Trail Valley Creek site (Grünberg and Boike, 2019), and snow depth distribution patterns reported in Walker et al. (2021). Topography and vegetation types were used as a visual proxy for spatial delineation of the combined influence of snow depth and soil properties (temperature, moisture and organic content) on potential for carbon fluxes. To account for any preferential orientation of snow resulting from wind redistribution, sampling transects were split into two sections of 100 m, orthogonal to each other. Snow free images of transect vegetative cover are shown in Fig. A1.

Vertical snowpack structure profiles and CO₂ and CH₄ concentration profiles were sampled along five 200 m transects with 10 m spacing between each microsite along each transects (Fig. 1, described in Table 1). Three 1 m spaced samples also taken at the 50 and 150 m mark of each transect (at the midpoint of each section) to further identify highly localised variance at metre-scales. In total, 119 (N_tot) co-located snowpack structure profiles and CO₂ and CH₄concentration profiles were collected across the study site.

Table 1 Samling transect description.

Vertical profiles of snow density, depth, and specific surface area were measured at sub-millimetre resolution using a Snow Micropenetrometer (SMP; Proksch, Löwe, and Schneebeli 2015). Raw SMP measurements were processed using the snowmicropyn (v. 1.2.1) python package with penetration force to snow density relationship coefficients adapted by Dutch et al. (2022) to the arctic snowpacks found in TVC.

Snowpack CO₂ concentration profiles were measured using a temperature controlled portable IRGA (LGR U-GGA-915 Ultraportable, Los Gatos Inc, Los Gatos, CA) deployed in a toboggan. Snow depths at each profile were calculated as the average of three measurements taken within 100 mm of each microsite profile (coloured icons, Fig. 1). Gases within the snowpack were sampled at increments of 100 or 200 mm depending on total thickness of the snowpack, with a minimum of 6 measurement depths per profile. Snowpack air was drawn to the analyser through a metal tube on a graduated avalanche probe by the instrument's fixed internal pump at approximately 0.8 standard L min⁻¹, and each depth was measured for approximately 60 s until CO₂ and CH₄ concentrations reached a stable plateau. The analyser is a flow-through instrument; sampled air was vented to ambient air rather than returned to the snowpack. The volume drawn per depth (∼ 0.8 L over 60 s) is a few percent of the pore-air volume of the surrounding snow, of the order of tens of litres within a radius of at most 0.3 m from the intake at the measured porosities (φ≈ 0.65 to 0.78). Cumulative depletion was therefore limited to the volume within each vertical layer measured. Ambient air concentrations were also recorded at the start of each profile within the transect drawn from approximately 1.5 m height to avoid contamination.

Gas concentration profiles were recorded as one continuous IRGA time series per transect accompanied by manually recorded total profile depths, creating a timeseries of step changes in concentration for each depth where air was sampled in the snowpack. IRGA time series were then partitioned into individual profiles with defined concentrations per depth in a MATLAB environment, using the Signal Processing Toolbox findchangepts.m function (The MathWorks, 2024), forcing known number of step changes in mean gas concentration between the start and end times of each profile. Step changes were first identified on the CO₂ concentration time series as signal noise to change in concentration per step change was lower and then used to partition the CH₄ time series. Concentrations at each step were then associated to independently recorded depths for each profile. Snowpack CO₂ and CH₄ concentration gradients used to calculate soil-atmosphere CO₂ and CH₄ fluxes through the snowpack ($F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$, g C m⁻² d⁻¹, respectively were then calculated by fitting linear models at each profile and testing for significance of gradient (T-test) and model coefficient of determination (R²).

Flux estimates assume quasi-steady-state diffusive transport and linear concentration gradients through the snowpack, computed from the diffusive form of Eq. (1) only. To ensure the validity of this assumption and satisfy the conditions under which the wind-pumping enhancement is negligible (Sect. 2.1), all measurements were conducted during stable meteorological conditions, characterized by wind speeds below 3.5 m s⁻¹, and were restricted to microsites free of pronounced sastrugi or drift relief and capped wind crust free of ice lenses. Barometric pressure was measured at the Environment and Climate Change Canada weather station (Trail Valley Creek) located near transect T2-Mixed, and varied by at most 10 hPa (1.0 %) within each sampling day and by 37 hPa (3.7 %) across the campaign. Pressure changes exchanged no more than the upper 1 % of snowpack pore air with the atmosphere per day (≤ 13 mm for the deepest snowpack), constituting a negligible advective contribution under low-frequency barometric forcing.

Soil-snow interface temperatures (basal temperatures) and snowpack temperatures sampled at 10 cm vertical intervals were measured using a handheld type K thermocouple probe. Missing basal temperature observations due instrument malfunction were inferred using a linear relationship between snow depth and soil temperature (missing = 53, R²=0.76, RMSE = 1.47 °C).

2.4 Snowpack diffusion flux uncertainty quantification

Measurement uncertainties were propagated to $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ according to the generalised error propagation method for nonlinear systems (Taylor, 2022). Gas concentration profile gradient uncertainty is defined as the linear fit 95 % confidence interval using the Wald method for binomial series, assuming that measurements used to calculate gradients are exact as the IRGA concentration measurement error was negligible compared to the error on profile gradients. Error on average profile snow density was defined as the standard deviation of binned densities in 50 mm vertical increments, where error-per-binned densities were propagated to the total average. Soil temperature, snow temperature, and air pressure measurements errors were also fully propagated to $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ using the generalised error propagation method.

2.5 Sub-sampling gas concentration micro-site profiles and transect-level median estimates

Snowpack concentration gradients derived from all vertically spaced measurements in each microsite profile were compared to sub-sampled vertical gradients calculated from; (a) two concentration measurements: one within the air-snow interface layer of snow and one at the bottom of the snowpack (two measurements per microsite), and (b) one measurement from the bottom of the snowpack and the average concentration from ambient air above snow surface for the whole transect (one measurement per microsite and an average of ambient near-surface air).

To assess gas flux spatial variability across transects (homogenous surface covers), robust median estimates with acceptable variances are required to calibrate and test TBMs. However, the ideal number of samples per surface cover cannot be know in advance. Here, we used a bootstrap method where median flux values (Md_full) of the total number of flux measurements (N) from each surface cover transect were compared to median flux values from incrementally smaller subsamples per transect. Subsample size decreased from a total number of flux measurements from each transect $n=N-\mathrm{1}$ to n=2 (T1 N=19; T2 to T5 N=25). Subsamples with n-values less than the maximum in each transect were randomly selected from the total possible. Random selection was repeated 10 000 times, the median calculated for each selection, and then the median of all 10 000 selections was calculated (Md_n). Differences between the medians of subsamples and all measurements in a transect are expressed as a ratio (Md_n $/$ Md_full). Subsampled Md_n were considered different to the median of total number of flux measurements from each transect (Md_full) when the ratio Md_n $/$ Md_full fell within the 19.4 % relative uncertainty reported for the CO₂ concentration-gradient (d[CO₂]/dz) regression of manually sampled snowpack fluxes at Trail Valley Creek by Mavrovic et al. (2023). This tolerance is a method-based precision band rather than a statistical significance level, adopted so that the subsampling result is directly comparable with the accuracy achievable by the gradient method. Profiles exhibiting non-significant concentration gradients (p>0.05) were excluded from subsampling and scaling analyses to avoid introducing noise from profiles indicating negligible net exchange. Profiles exhibiting non-significant concentration gradients were interpreted as representing near-zero flux conditions for the corresponding sampling period.

3 Results

3.1 Snow physical characteristics

Snow density profiles exhibited typical Arctic tundra characteristics of high-density near-surface wind crust and lower density basal depth hoar layers (Fig. 2, a to f). Depth hoar layers were most pronounced in T1-Shrub and T2-Mixed, and T5-Inlet, with median wind crust densities of 297 ±70 kg m⁻³ and depth hoar densities of 223 ± 46 kg m⁻³. Greater spread in T3-Tussock mid-snowpack densities may indicate several harder snow layers forming during consecutive high wind events as the snowpack built over winter. Snowpacks were dry and free of discrete ice lenses or basal ice layers at all microsites throughout transects.

Pronounced microtopography and shallow snow in T3-Tussock created more variable density profiles, corroborating with observed variable wind crust presence and profile depth. T4 forest snowpack were deeper than for other transects but had similar profiles of high-density surface layers overlying lower density depth hoar (wind crust Md: 314 ± 84 kg m⁻³, depth hoar Md: 219 ± 40 kg m⁻³).

Figure 2 Normalised medians snow density profile with 25th/75th percentiles (dark) and 5th/95th percentile (light), for sampling transects T1 to T5 (a to e). (f) shows snow basal temperature and snow depth distribution per transect.

Snow depths and basal (soil-snow interface) temperatures (Fig. 2f) were strongly correlated (Pearson's Linear Correlation Coefficient, r=0.90, α< 0.01). Transect scale basal temperature and snow depth medians were significantly different (alpha < 1 %) to each other, except for T1-Shrub and T5-Inlet that had similar distributions. T3-Tussock basal temperatures were lowest with lowest snow depths (Md: −9.5 ± 0.8 °C, 440 ± 118 mm), contrasting with T4-Forest with highest basal temperatures and snow depths (Md: −2.7 ± 2.6 °C, 1280 ± 343 mm). Basal temperature and snow depth distributions in other transects are all found between T3-Tussock and T4-Forest ranges, with the T2-Mixed significantly colder (α< 0.05) than T1 and T5. Transect density, basal temperature, and snow depths are summarised in Table 2.

Table 2 Transect scale summary of sampled microsites taken at the TVC site in late March 2024.

Notes: ^* Numbers within [square brackets] show interquartile range. ${}^{**}$ CO₂ median fluxes are shown in mgC instead of gC for table readability. Highest flux magnitudes are in bold, with CH₄ uptake in italic.

3.2 Carbon dioxide and methane concentration profiles

Individual snowpack CO₂ and CH₄ concentration profile analysis revealed that most profiles across transects T1-Shrub, T2-Mixed, T4-Forest, and T5-Inlet were explained with high certainty by simple linear models. Transect T3-Tussock microsites contained weak GHG gradients, where the majority of snowpack CO₂ and CH₄ concentrations were not significantly different to atmospheric concentrations with 80 % (CO₂) and 56 % (CH₄) of gradient p-values > 0.05. Although well explained by linear models (gradient p-value < 0.05), some profiles display nonlinearity, suggesting the influence of snowpack vertical structure for some locations.

Transect scale concentration profile linearity were also examined by aggregating normalised CO₂ and CH₄ concentrations as a function of depth (Fig. 3). Linear models were fitted to each normalised gas species (CO₂ and CH₄) with model 95 % confidence intervals illustrating the level of certainty about the model. Normalised T1-Shrub concentration gradients show positive linear trends for both CO₂ and CH₄, with variance well explained by the model (R² > 0.75), indicating a strongly linear trend across all profiles. T2-Mixed has similar relationship to T1-Shrub but shows more spread (R²> 0.55), suggesting either more variance in flux magnitudes or less linearity in concentration profiles. T3-Tussock shows positive gradients for both gas species, with CO₂ variance being poorly explained (R² < 0.1) while CH₄ shows a tighter albeit poor fit (R²=0.37). T4-Forest concentration profiles are the most linear for both CO₂ and CH₄ with very little spread around the linear model (R²> 0.85), indicating the least dependence to vertical snow structure. Contrarily to all other transects, All T4-Forest CH₄ vs. depth gradients are negative, indicating a net CH₄ uptake in the forested area. T5-Inlet CO₂ profiles are also strongly linear (R²=0.75), while CH₄ show more spread (R²=0.35), with 76 % of profile gradients significant (p-value < 0.05). Similarly to T3-Tussock, several profiles across T5 had CH₄ concentrations near atmospheric levels, with a few profiles with high basal concentrations. T5-Inlet was the most varied in terms of land cover type, i.e., dry shallow snowpack slope, and wet high snowpack drainage channel.

Figure 3 Linear model fit through depth normalised CO₂ (a–e) and CH₄ (f–j) concentration profiles across TVC transects. Profiles are aggregated per transect to illustrate consistency of profile shape; individual microsite gradients were fitted separately and used for all flux calculations. All gradients are significant with p-values < 0.01 (T test).

3.3 Carbon dioxide and methane fluxes

Instantaneous CO₂ and CH₄ fluxes ($F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$, respectively) through the snowpack for each vertical profile were calculated using paired snow gas concentrations and snow density profiles, aggregated by transect (Fig. 4). Distributions of CO₂ fluxes closely correspond with snow depth distributions, which strongly control soil surface temperatures. Distributions of fluxes between transects were all significantly different to one another except for T1-Shrub and T5-Inlet. Very low fluxes at T3 coincide with the shallowest snowpacks and coldest basal temperatures. Median CO₂ and CH₄ fluxes at T4-Forest were an order of magnitude higher than other transects, coinciding with the deepest snowpacks and warmest basal temperatures. Five high CO₂ emitting outliers were found in T2-Mixed, and T3-Tussock, T4-Forest.

Figure 4 Transect scale instantaneous CO₂ (a) and CH₄ (b) fluxes ($F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$, respectively) Median snow depths per transect are shown above the CO₂ box and whisker plots in italic. Note the change of unit scale between left and right panels.

Methane fluxes through the snowpack did not follow the same pattern as CO₂ fluxes Instead, significant net CH₄ uptake ($-F_{{\mathrm{CH}}_{\mathrm{4}}}$) was found throughout T4-Forest where snowpacks are thickest and basal temperatures in frozen soils are warmer than −3 °C. Low positive CH₄ fluxes were found across T2-Mixed with T1-Shrub and T3-Tussock, including some high emitting profiles where net fluxes were one order of magnitude higher than profile medians. T5-Inlet CH₄ fluxes had most variability with weak net uptakes and strongly emitting profiles which show little correlation to snowpack depth and basal temperatures (R²=0.06 and −0.02, respectively), contrary to patterns in CO₂ fluxes (Table 2).

We found a significant relationship between T4-Forest microsite rates of net CH₄ uptake (negative $F_{{\mathrm{CH}}_{\mathrm{4}}}$) and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ (Fig. 5). Although significantly related (p-value < 0.01), T4-Forest $F_{{\mathrm{CO}}_{\mathrm{2}}}$ to $F_{{\mathrm{CH}}_{\mathrm{4}}}$ ratios are still poorly explained by magnitudes alone with an R²=0.14. Spread around the linear model is partially explained by basal temperatures, $F_{{\mathrm{CO}}_{\mathrm{2}}}$ increasing proportionally more with warmer soil surfaces than CH₄ uptake rates.

Figure 5 CH₄ uptake rates (negative $F_{{\mathrm{CH}}_{\mathrm{4}}}$) as a function of respiratory $F_{{\mathrm{CO}}_{\mathrm{2}}}$ along T4-Forest. Flux uncertainties are shown by the grey error bars.

Measured CO₂ fluxes followed an exponential relationship to basal temperatures, with R² of 0.79 and RMSE of 0.020, most situated within the model's 95 % regression parameter confidence interval (Fig. 6). A two-parameter exponential regression curve was fitted on $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and basal temperature, and compared to relationships found elsewhere (Natali et al., 2019; Mavrovic et al., 2023). While within this study's 95 % regression parameter confidence interval, Mavrovic et al. (2023) overestimated fluxes with basal temperatures below −4 °C relative to our study, whereas Natali et al. (2019) overestimated low temperature CO₂ fluxes (as noted in Mavrovic et al., 2023). This may be due to our study including a substantial number of fluxes measured at temperatures of −6 to −0.5 °C, which were insufficiently represented in previous datasets.

Figure 6 (a) CO₂ flux as a function of snowpack basal temperature, with exponential model fit for this study and similar studies in winter Arctic environments. Shaded area represents the 95 % regression parameter confidence interval around this study's exponential model (black line). (b) CH₄ flux as a function of soil-snow interface temperature. Soil temperature and flux uncertainty propagated from measurements and flux calculations are shown by the grey error bars.

Methane fluxes have less distinct relationship with soil surface temperature than CO₂ fluxes. T1-Shrub, T2-Mixed and T3-Tussock microsites were all below −6 °C and showed very little CH₄ emissions, with two high emitting microsites (up to 0.58 mgC m⁻² d⁻¹ or 0.56 nmol m⁻² s⁻¹) in the coldest transect T3. All fluxes measured in T4-Forest and had basal temperatures were above −5 °C and showed weak net CH₄ sinks between 0 and −0.1 mgC m⁻² d⁻¹ (−0.10 nmol m⁻² s⁻¹). Although within the same temperature range as T1-Shrub and T2-Mixed (−8 °C < T < −4.5 °C), T5-Inlet had the highest CH₄ emissions with three profiles displaying uptake of the same magnitude as T4-Tussock. T5-Inlet had the most diverse surface vegetation and snow depths across its sampling transect, with the first 100 m of the transect situated on a well-drained and wind scoured eastward facing slope covered with dwarf birch, and the second 100 m along the lake inlet with taller alder shrubs, tussocks and several micro drainage channels which are generally saturated with water across the shoulder and growing seasons.

3.4 Sampling method analysis

We separately calculated $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ from only basal and near surface, and from snow basal and average ambient air samples to assess advantages and need for snow gas profiles. Both subsampling strategies tightly predicted fluxes similar to those fluxes calculated using multiple evenly spaced snowpack samples (Fig. 7). Both subsampling methods showed near 1:1 ratio for $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ with R²> 0.99 and RMSE two orders of magnitude smaller than $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and one order of magnitude smaller than $F_{{\mathrm{CH}}_{\mathrm{4}}}$. As two-sample and one-sample methods show similar linear relationships, concentration gradient and consequently vertical fluxes have been shown to be almost entirely driven by the snowpack basal concentration.

Figure 7 Linear model fitted to $F_{{\mathrm{CO}}_{\mathrm{2}}}$ (a, c) and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ (b, d) calculated from snowpack basal and near surface profile (top) samples, and from basal profile samples and average ambient air (bottom). Dotted lines show the 1:1 ratio line and shaded areas show the linear model 95 % confidence intervals.

Transect subsampling comparisons show that $F_{{\mathrm{CO}}_{\mathrm{2}}}$ Md_n variability remains within 90 % of the fully sampled dataset Md_full with n>9 in transects T1-Shrub, T4-Forest, and T5-Inlet, and n > 16 for T3-Tussock (Fig. 8). T3-Tussock subsampled $F_{{\mathrm{CO}}_{\mathrm{2}}}$ Md_n did not stabilise within subsamples suggesting surface cover variability within T3 may require a larger sampling size than we conducted (n > 25). $F_{{\mathrm{CH}}_{\mathrm{4}}}$ transect Md_n required more subsamples to capture surface cover scale variability in Md_full, with n> 12 for T2-Mixed and T4-Forest, while other transects required n > 20 when comparing means but did not reach a stable plateau within N samples when comparing medians. TVC site scale analysis (black lines Fig. 8a and b) shows that $F_{{\mathrm{CO}}_{\mathrm{2}}}$ were within 90 % of Md_full with transect n> 12 for medians and n> 13 for means, while CH₄ failed to be adequately reproduced with subsamples below N.

Figure 8 Surface cover subsampling analysis. Horizontal axes show the number of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ subsamples (Md_n) randomly selected from total number of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ measurements from each surface cover transect picked out of the total number of samples. Vertical axes show the proportion of times the ratio of medians of subsamples and all measurements in a transect (Md_n $/$ Md_full) were within 19.4 % Md_full. The thicker black line shows the same exercise combining all transects.

4 Discussion

4.1 Transect scale CO₂ and CH₄ fluxes

We consistently observed positive (soil-to-atmosphere) CO₂ emissions ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) at all microsites with significantly different rates between surface cover transects. Instantaneous $F_{{\mathrm{CO}}_{\mathrm{2}}}$ were similar to those derived using syringe extraction methods by Mavrovic et al. (2023) at TVC during the late 2022 and 2023 winters, and comparable to pan-Arctic measurements using comparable in-situ concentration measurements or a wide range of other techniques including eddy covariance and flux chambers (Natali et al., 2019; Pirk et al., 2016; Webb et al., 2016). However, our study resolved a paucity of Arctic $F_{{\mathrm{CO}}_{\mathrm{2}}}$ at soil temperatures between −6 and 0 °C. Filling this measurement gap reduced error in exponential model fits by 0.01 and 0.004 gC m⁻² d⁻¹ for Natali et al. (2019) and Mavrovic et al. (2023), respectively (Fig. 6a). A wide basal temperature range and lower $F_{{\mathrm{CO}}_{\mathrm{2}}}$ uncertainty at temperatures below −6 °C, constrains exponential $F_{{\mathrm{CO}}_{\mathrm{2}}}$ response to soil temperature, indicating $F_{{\mathrm{CO}}_{\mathrm{2}}}$ tends to zero at lower temperatures observed at TVC ($<-\mathrm{9}$ °C).

Measured winter CO₂ flux magnitudes at TVC are consistent with early observations from Arctic Alaska showing sustained cold-season emissions under snow cover (Fahnestock et al., 1998, 1999; Jones et al., 1999b; Welker et al., 2000). Subsequent work has demonstrated that snow depth, thermal insulation, and active-layer dynamics exert strong controls on winter respiration and the age of respired carbon (Lupascu et al., 2018; Pedron et al., 2023; Schimel et al., 2006; Sullivan et al., 2008). Measurements at TVC extend these findings by resolving fine-scale spatial variability in winter fluxes and by evaluating how simplified sampling strategies capture this variability under dry snowpack conditions.

Significant differences in the magnitudes of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ between transects coincide with spatial patterns of basal temperatures, indicating CO₂ emissions from each surface cover are dominated by microbial respiration in frozen soils (Pedron et al., 2022). Consistently high emitting microsites in T4-Forest are typical of fluxes observed in permafrost free northern boreal forests with similar snow depths (Mavrovic et al., 2023). Near zero basal temperatures at T4-Forest suggest a portion of the vertical soil column might remain at or near zero-curtain conditions, a narrow temperature range (0.05 °C) in which changes in heat energy alters the phase of water without changing soil temperature (Zona et al., 2016; Lupascu et al., 2018; Devoie et al., 2022). $F_{{\mathrm{CO}}_{\mathrm{2}}}$ outliers (Fig. 4) at low surface temperatures (basal temperature $<-\mathrm{15}$ °C), may have resulted from wintertime microbial activity in localised warmer soil pockets, or slow release of trapped CO₂ produced during autumn freeze-up similarly to those found in other tundra sites at similar latitudes (Pedron et al., 2022). Talik formations have been reported in dry tundra environments by Kim et al. (2021) at similar latitudes (Cambridge Bay, NU Canada), but have yet to be verified in TVC.

Consistent net uptake measured in the T4-Forest transect between −6 and 0 °C reveals that high-latitude forested soils can function as significant net sinks even under frozen conditions, representing an integrated balance of diffusive transport and microbial methanogenic and methanotrophic turnover within soils (Mast et al., 1998). Few measurements of wintertime CH₄ uptake specific to land surface types have been reported in high-latitude tundra environments, mainly due to the historical difficulty of resolving low-magnitude fluxes through Arctic snowpacks. Recent studies have begun to address this gap, with winter CH₄ uptake rates in northern closed-crown boreal forest uplands recently reported at −0.43 ± 0.34 to −0.47 ± 0.26 mg C m⁻² d⁻¹ (Mavrovic et al., 2025). For comparison, our observed winter forest uptake (−0.04 to 0.08 mg C m⁻² d⁻¹) represents a persistent, albeit diminished sink relative to peak growing season consumption. At TVC, automated chambers have recorded summer CH₄ uptake rates for lichen and shrub communities between −0.37 and −0.44 mg C m⁻² d⁻¹ (converted from −0.49 to −0.59 mg CH₄ m⁻² d⁻¹), suggesting winter forest uptake magnitudes are roughly one order of magnitude lower than peak summer rates in the same region (Voigt et al., 2023).

Landscape-scale tundra winter CH₄ sinks have also been observed in Svalbard, Greenland, and Norway, where soil organic layers are less spatially pervasive than in TVC (D'Imperio et al., 2023; Donateo et al., 2025; Treat et al., 2018). In Svalbard's dry upland tundra, median CH₄ fluxes remain negative throughout the winter, ranging from −0.18 to −0.37 mg C m⁻² d⁻¹ (reported as −0.17 to −0.36 nmol m⁻² s⁻¹), while notable net uptake events during the freezing shoulder season can reach −0.72 mg C m⁻² d⁻¹ (reported as −0.69 nmol m⁻² s⁻¹) (Donateo et al., 2025). Growing season and start-of-cold-season CH₄ uptake measured at TVC across drier landforms, including upland tundra and polygonal rims, further suggest that methanotrophic activity persists long after plant senescence (Ivanova et al., 2025; Voigt et al., 2023).The magnitude of this winter sink is particularly important when contrasted with comparatively wetter tundra sites in Northern Alaska, where summer emissions from even relatively dry polygon rims can range from 4.5 to 9.6 mg C m⁻² d⁻¹ (McEwing et al., 2015). Comparisons between different pan-Arctic land covers highlight that while mid-winter CH₄ uptake may be a small component of the annual budget, persistence of aerobic oxidation in well-drained forested and upland soils provides a non-negligible offset to winter CH₄ emissions from saturated hotspots.

Near-zero basal temperatures at T4-Forest suggest that much of the near-surface soil profile likely remained at zero-curtain temperature conditions, When soil temperatures are poised near 0 °C due to release of latent heat during the phase transition of water to ice, allowing liquid water content and ice to coexist in soil pore space (Devoie et al., 2022; Tucker, 2014; Outcalt et al., 1990). Remaining liquid soil water over winter sustains microbial metabolic activity in thin films surrounding soil particles even at temperatures as low as −18 °C (Dutch et al., 2024; Elberling and Brandt, 2003; Sullivan et al., 2008).

Northern boreal forest CH₄ uptake rates have been reported to be within the same order of magnitude as those found in T4-Forest, where soils also remained in zero-curtain conditions throughout winter (Lee et al., 2023; Mavrovic et al., 2025). Evidence for microbial CH₄ consumption within the T4-Forest is provided by the significant inverse relationship between $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ (Fig. 5). This inverse relationship aligns with global patterns where winter CO₂ emissions explain nearly 96 % of the variability in CH₄ oxidation, effectively serving as a proxy for the biological activation and oxygen diffusion required for substrate consumption (Voigt et al., 2023; Kim et al., 2024). High-affinity methanotrophs likely maintain this activity in the aerobic upper soil layers during mid-winter; however, the preferential coupling of basal temperatures with CO₂ rather than CH₄ flux magnitudes suggests that CH₄ consumption may also take place in deeper, relatively warmer soil layers (Donateo et al., 2025; Mikan et al., 2002). Further variability in CH₄ uptake magnitudes could also be explained by oxygen, and the liquid portion of soil moisture available in upper soil horizons have been identified as important regulators of uptake (D'Imperio et al., 2023; Lee et al., 2023; Rößger et al., 2022).

However, soil temperature measurements taken at the soil-snow interface may be insufficient to infer liquid water availability across the upper soil layers. As the active layer typically freezes faster from the surface downward due to large soil-atmosphere thermal gradients, interface temperatures may drop below ideal temperatures for microbial activity, while deeper carbon-rich layers remain in a zero-curtain state for weeks or months longer (Pedron et al., 2022; Yi et al., 2020). Relying solely on surface temperatures contributes to the cold-soil temperature biases in carbon fluxes observed in terrestrial biosphere models such as CLM5.0, which often fail to resolve the extended duration of these subsurface thermal refugia and resulting non-growing season carbon losses (Damseaux et al., 2025; Dutch et al., 2024; Tao et al., 2021).

We also measured CH₄ uptake in locations with notably cooler soils in T5-Inlet (temperature $<-\mathrm{6}$ °C), which also contained the highest sitewide CH₄ emitting microsites. The T5-Inlet transect crosses a lake inlet area that is typically saturated or submerged during the autumn shoulder season, forming a patchwork of dry grassy hummocks, tall alder shrubs, and 10–40 cm deep drainage channels often permeating the active layer. While soil surface temperatures in T5-Inlet did not differ significantly from those in T1-Shrub, greater heterogeneity in snow depths and elevated spatial variability in soil water saturation in T5-Inlet may account for enhanced positive $F_{{\mathrm{CH}}_{\mathrm{4}}}$ and suppressed $F_{{\mathrm{CO}}_{\mathrm{2}}}$ typical of anoxic environments. Conditions for preferentially anoxic microbial activity are likely driven by prolonged thermal retention in saturated microsites during winter due to increased latent heat release during freeze-up, similarly documented in hydrologically comparable “wet” tundra systems in northern Alaska (Blanc-Betes et al., 2025).

Alternatively observed net CH₄ uptake in saturated landcovers may be limited by physical encapsulation of soil layers by ice, where basal ice or thick frozen soil layers creates an impermeable barrier occluding efflux of CH₄ produced in deeper, unfrozen soil horizons (Elberling and Brandt, 2003; van Bochove et al., 2001). Such diffusion-limited transport, exacerbated by the occlusion of soil pores by ice, can decouple subsurface CH₄ production from soil-atmosphere exchange, potentially resulting in apparent uptake signals if atmospheric CH₄ is consumed in the uppermost aerated soil or if production is physically trapped at depth (Mavrovic et al., 2025; Treat et al., 2018; Schimel et al., 2006). Slower CH₄ diffusion through ice-filled pores in frozen soils can further mask the temperature dependence of winter CH₄ fluxes, suggesting that surface net uptake may reflect a lack of transport from deeper source layers rather than high methanotrophic potential (Mavrovic et al., 2025). Consequently, the fine-scale variability at T5-Inlet likely reflects a gradient where some microsites maintain active methanogenesis while others are governed by the passive trapping of greenhouse gases within the soil profile.

Further efforts to partition $F_{{\mathrm{CH}}_{\mathrm{4}}}$ into its autotrophic and heterotrophic respiration components through isotopic analysis (e.g., Pirk et al., 2016) or via metagenomic sequencing may help identify the relative importance of soil surface temperature, liquid moisture availability and nutrient availability on CH₄ flux, and allow estimation of CH₄ specific temperature response curves under methanogenic or methanotrophic prevailing conditions. Identifying sampling locations prior to snowfall and optimizing repeated flux measurements to capture more variability in vegetation, soils, and hydrology, throughout the growing and autumn shoulder seasons would also allow process level understanding of abiotic and biotic drivers of flux variability, providing a more comprehensive dataset to upscale fluxes to landscape and regional spatial extents.

4.2 Profile linearity and flux uncertainty

Overall GHG concentrations throughout snowpacks were linear in nature (Fig. 3). Surprisingly, denser upper snowpack wind crust layers did not appear to limit CO₂ and CH₄ diffusion rates. Consistently significant linear fits across transects independent of snow depth (except T3-Tussock, where fluxes were too small) demonstrate that the snowpack acted as a uniform porous material (i.e., independent of snow structure). Concentrations at distinct vertical sampling points concentration profiles were reproducible across replicate profiles and consistent between full-profile and two-point flux estimates (Fig. 7), indicating that extraction did not draw down or homogenise air from different layers or ambient air. As GHG concentrations follow a linear trend from basal to air above snowpacks, microsite level concentration gradients are controlled by their respective basal concentrations. Future GHG emissions could therefore be estimated using basal concentration variability within landscapes.

Robust agreement between full-profile fluxes and estimates derived from basal and near-surface concentrations was evident under the predominantly dry, cold snowpack conditions encountered at TVC. Where snowpacks lack extensive ice layers and experience limited mid-winter melt-freeze cycles, molecular diffusion will remain the dominant transport mechanism. Melt-refreeze events promoting the formation of spatially coherent ice lenses within snowpacks may cause diffusive transport to deviate from linear assumptions and decouple subsurface microbial production from surface exchange (Elberling and Brandt, 2003; Mavrovic et al., 2023). Soil stratigraphic barriers or basal snowpack ice lenses can trap gases produced in deep soil layers, potentially leading to zero-gradient fluxes within the overlying snow despite ongoing soil metabolic activity

Non-linear departures in snowpack gas concentration profiles, particularly homogenization of gas concentrations observed in basal depth hoar layers, introduce a known potential for gradient underestimation (Dutch et al., 2022). Such concentration gradient (dC $/$ dz) underestimations, particularly in low carbon flux environments such as the T3-Tussock transect, may reflect a reduced concentration change over the sampled depth, resulting in $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ magnitudes that are lower than the actual soil metabolic flux rates or mechanical release. Wind-induced advection (wind pumping) can further perturb gradients in upper snow layers, although these effects are often transient and confined to periods of high wind speeds (Jones et al., 1999a; Mavrovic et al., 2025), although those effects were negated in this study by sampling while wind speeds remained below 3.5 m s⁻¹.

Previous studies have also identified non-linear profiles in snowpack GHG distributions when impermeable soil layers were present (van Bochove et al., 2001). In such cases, using snowpack basal concentrations would likely overestimate instantaneous flux magnitudes, and would affect estimates in the timing of emissions rather than seasonal budget estimates. While profiles at TVC were predominantly linear, significant gas trapping below frozen soil barriers can artificially increase concentration gradients once released through thermal or mechanical cracking, or conversely, mask subsurface production entirely (Jones et al., 1999a). Such diffusion-limited transport decoupling suggests that measurements taken strictly within the snowpack may periodically reflect physical encapsulation rather than active biological oxidation (Elberling and Brandt, 2003; van Bochove et al., 2001).

The application of Fickian diffusion models to wintertime greenhouse gas exchange assumes linear concentration gradients through the snowpack. Uncertainty bounds (Figs. 5 and 6) incorporate the standard error associated with linear regression of concentration versus depth, accounting for variance introduced by non-linear yet significant vertical profile concentration gradients. Robustly quantification of uncertainty for $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ from both instrument measurement error and linear model uncertainty was then made possible. High resolution snow density profile measurements (SMP) and high precision in situ gas measurements using a portable IRGA allowed a tighter relationship between $F_{{\mathrm{CO}}_{\mathrm{2}}}$ magnitudes and basal temperatures than comparable studies and gave statistical confidence on reported $F_{{\mathrm{CH}}_{\mathrm{4}}}$ magnitudes in a low flux environment. Some profiles only became linear at 100 to 300 mm above the ground surface with concentrations homogenised below these depths. Profiles where CO₂ and CH₄ concentrations were more homogenous in the depth hoar layer likely caused gradient underestimation thus leading to lower $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$magnitudes. Conversely, bulk snowpack density and temperature values derived from the entire snow column could lead to $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$overestimation due to an underestimation of effective snow porosity and tortuosity, making the snowpack less resistant to vertical molecular diffusion. In either case, we recalculated $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ using per-microsite range of densities and air temperatures and found resulting instantaneous flux variability to be within reported propagated flux uncertainty.

We found CO₂ and CH₄ concentration measurement error from the portable IRGA to have negligible effect on $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ uncertainties, which were mostly controlled by snowpack gas concentration gradient confidence intervals. Flux uncertainties are often reported but rarely well constrained, whether from microsite scale measurements due to low sample size, or from integrated landscape scale flux measurements such as eddy covariance, where systematic errors are difficult to separate from random measurement uncertainty (Levy et al., 2022). The importance of accurate uncertainty quantification in greenhouse gas flux measurements is underscored by (Fisher et al., 2014), who highlighted the significant uncertainties in carbon cycle models for the Alaskan Arctic. Similarly, (Björkman et al., 2010) demonstrated that different methods for estimating winter CO₂ effluxes can yield results that vary by up to two orders of magnitude, emphasizing the need for robust error analysis and uncertainty calculation. Measurement error propagation to $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ uncertainties presented here are nevertheless within similar magnitudes as snowpack gradient method derived fluxes reported by (Mavrovic et al., 2023, 2024).

4.3 Spatial variability and scalability

Understanding the distribution of spatial variance in CO₂ and CH₄ soil-atmosphere fluxes across different land surface types is vital to enable up-scaling to the landscape scale, where flux uncertainties must account for surface cover aggregation error. Past studies reporting winter fluxes have either compromised spatial coverage for temporally continuous data (Zhu et al., 2014; Seok et al., 2009), or report spatially aggregated landscape scale fluxes that pose challenges to partitioning into their component surface covers (Levy et al., 2022; Pirk et al., 2017). Fast, instantaneous CO₂ and CH₄ concentration measurements using a roving, high-precision portable IRGA and SMP allowed robust quantification of spatial variability within surface cover types. This is particularly relevant to upscaling small-scale fluxes from tens of metres to scales typically used in Terrestrial Biosphere Models (TBMs), where finer-scaled data is aggregated by averaging discrete data points to much larger grid cells (Levy et al., 2022).

Measurements from TVC can directly inform TBMs by providing high-resolution empirical bounds for temperature sensitivity (Q₁₀) of cold-season biogeochemical processes as a primary source of uncertainty in Arctic carbon simulations (Dutch et al., 2024). The relationship between late-winter soil temperature, snow depth, and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ suggests that deterministic models using grid-cell mean snow depths likely underestimate winter activity by failing to capture the disproportionate contribution of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ from warmer frozen soils under deeper snow. Our findings suggest that bulk snowpack diffusivity remains a robust approximation for models, provided that refined snow thermal conductivity parameterizations are implemented to correct the cold soil temperature biases that currently inhibit accurate flux simulations. A need for revised temperature-response curves in Arctic tundra extends to CH₄ dynamics; where persistent uptake observed in the forested transects indicates that current TBM methane production and oxidation schemes require better low temperature calibration to account for microbial oxidation occurring at soil temperatures below 0 °C. Following the process-based evaluation by Dutch et al. (2024), which highlighted how default model configurations often simulate near-zero winter fluxes at TVC due to overly restrictive soil moisture thresholds (Ψ_min) and respiration temperature responses, a transition is recommended toward sub-pixel distribution-based parameterizations (e.g., probability density functions) that can represent sub-grid heterogeneity in snow insulation and soil thermal regimes.

Errors in the estimation of spatial variability using methods such as Bayesian Kriging can be constrained by informing the number of measurements needed to quantify spatial variability within homogenous surface covers (Levy et al., 2022). A subsampling experiment (Figure 8), demonstrated that at TVC the homogenous surface cover scale $F_{{\mathrm{CO}}_{\mathrm{2}}}$ variance can generally be captured with n >=12 sampling sites, while larger sample sizes were needed for $F_{{\mathrm{CH}}_{\mathrm{4}}}$ (n > 25). Both median and average $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ at T3-Tussock failed to be reproduced with fewer subsamples because fluxes were either one order of magnitude smaller than those found in other transects or were not significantly different to zero. T5-Inlet covered both a well-drained low snow slope and a flat high snow and historically wet drainage area, exhibiting both high CH₄ emitting and low uptake microsites. Although closely grouped spatially (within 150 m), T5-Inlet's more heterogeneous snow depths and soil moisture (as well as associated nutrient availability), in contrast to other transects most likely explained why $F_{{\mathrm{CH}}_{\mathrm{4}}}$ average and median were not well reproduced with subsamples. This suggests future scaling attempts should consider a data informed method, such as multivariate spatial interpolation (and prior information) to divide up land surface for spatial flux upscaling.

Vertical snow structure had little effect on concentration profile linearity (and thus $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$), and surface cover scale flux could confidently be derived from concentrations near the soil-snow and air-snow interfaces (Fig. 7). We show no significant difference between calculating soil-atmosphere fluxes from basal snowpack CO₂ and CH₄ concentrations and ambient air above the snowpack, compared to calculations from a larger number of equally spaced samples throughout the snowpack. Together these findings reduce the measurement burden needed for rapid in situ instantaneous flux measurement, opening up potential for landscape-scale $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$ winter budgets at larger spatial scales using remotely sensed snow properties while expanding measurement spatial coverage improving quantification of variability over a wider set of surface covers (Dutch et al., 2024; Ivanova et al., 2026; Virkkala et al., 2022).

5 Conclusion

This study quantified wintertime CO₂ and CH₄ fluxes across contrasting tundra surface covers in the western Canadian Arctic, linking snowpack structure and soil surface temperatures to 10 m spatial variability in winter carbon fluxes. Our results provide valuable guidance for refining sampling protocols and provide key observational constraints for improving terrestrial biosphere model representations of Arctic winter carbon dynamics under ongoing climatic warming. Key findings presented in this work are:

• Carbon dioxide emissions were evident throughout TVC, but at significantly different rates within and between surface cover transects, primarily driven by basal soil temperature and local snow depth.

• Methane uptake was consistent throughout the forested transect (T4-Forest), bacterial methanotrophic activity in upper soil layers.

• Vertical snow structure has little effect on concentration profile linearity (and thus $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and $F_{{\mathrm{CH}}_{\mathrm{4}}}$). Reliable flux estimates using the snow gradient diffusion method were obtained using only near-surface and basal gas concentration measurements, reducing the need for full-depth gas profiles in future studies.

• Small-scale heterogeneity (10 m spacing) in snow depth, soil thermal regime, and vegetation exerts strong controls on wintertime CO₂ and CH₄ flux dynamics.

Appendix A

Figure A1 Snow free photos of vegetative cover at measurement transects. (a) T1-Shrub: dwarf birch shrub dominated with mosses and lichen; (b) T2-Mixed: mix of dwarf shrubs; lichens, and tussocks, (c) T3-Tussock: tussock field; (d) T4-Forest: black and white spruce sparsely forested area; (e) T5-Inlet: tussocks, alder shrubs in a lake inlet area, dwarf birch on neighbouring slope. Photo credits: Gabriel Hould Gosselin (a, d, e) and Kathryn Bennett (b, c).