Environmental controls of winter soil carbon dioxide fluxes in boreal and tundra environments

The carbon cycle in Arctic–boreal regions (ABRs) is an important component of the planetary carbon balance, with growing concerns about the consequences of ABR warming for the global climate system. The greatest uncertainty in annual carbon dioxide (CO₂) budgets exists during winter, primarily due to challenges with data availability and limited spatial coverage in measurements. The goal of this study was to determine the main environmental controls of winter CO₂ fluxes in ABRs over a latitudinal gradient (45^∘ to 69^∘ N) featuring four different ecosystem types: closed-crown coniferous boreal forest, open-crown coniferous boreal forest, erect-shrub tundra, and prostrate-shrub tundra. CO₂ fluxes calculated using a snowpack diffusion gradient method (n=560) ranged from 0 to 1.05 g C m² d⁻¹. To assess the dominant environmental controls governing CO₂ fluxes, a random forest machine learning approach was used. We identified soil temperature as the main control of winter CO₂ fluxes with 68 % of relative model importance, except when soil liquid water occurred during 0 ^∘C curtain conditions (i.e., T_soil≈0 ^∘C and liquid water coexist with ice in soil pores). Under zero-curtain conditions, liquid water content became the main control of CO₂ fluxes with 87 % of relative model importance. We observed exponential regressions between CO₂ fluxes and soil temperature in fully frozen soils (RMSE=0.024 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; 70.3 % of mean $F_{{\mathrm{CO}}_{\mathrm{2}}}$) and soils around the freezing point (RMSE=0.286 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; 112.4 % of mean $F_{{\mathrm{CO}}_{\mathrm{2}}}$). $F_{{\mathrm{CO}}_{\mathrm{2}}}$ increases more rapidly with T_soil around the freezing point than at T_soil<5 ^∘C. In zero-curtain conditions, the strongest regression was found with soil liquid water content (RMSE=0.137 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; 49.1 % of mean $F_{{\mathrm{CO}}_{\mathrm{2}}}$). This study shows the role of several variables in the spatio-temporal variability in CO₂ fluxes in ABRs during winter and highlights that the complex vegetation–snow–soil interactions in northern environments must be considered when studying what drives the spatial variability in soil carbon emissions during winter.

1 Introduction

Carbon stocks and fluxes in the Arctic and boreal biomes (hereafter called Arctic–boreal regions; ABRs) constitute large components of the planetary carbon balance (Tarnocai et al., 2009; van Huissteden and Dolman, 2012; Carreiras et al., 2017). ABRs store substantial quantities of carbon due to inherently slow decomposition rates, largely attributable to cold temperatures (Ravn et al., 2020). ABRs are warming up to 4 times faster than the rest of the planet with potential feedbacks to the global climate system (Derksen et al., 2019; Rantanen et al., 2022). Although ongoing warming of ABRs has the potential to lengthen growing seasons, enhance plant growth, and increase aboveground carbon storage (Sturm et al., 2005; McMahon et al., 2010), the growing-season vegetation response is variable and complex (Myers-Smith et al., 2020). Warmer air and soil temperatures enhance production and release of carbon dioxide (CO₂) from ecosystem respiration, comprising heterotrophic respiration by microbes decomposing soil organic matter, and autotrophic respiration by above- and belowground plant components (Bond-Lamberty and Thomson, 2010). The release of previously frozen carbon stocks is particularly important in regions undergoing permafrost thaw (ground completely frozen for at least 2 consecutive years) (Schuur et al., 2015; Natali et al., 2021; Miner et al., 2022). If increases in ecosystem respiration exceed those of photosynthetic CO₂ uptake from enhanced plant growth, ABRs may shift from a weak net CO₂ sink to a net CO₂ source, thereby generating a potentially non-negligible positive feedback to the global climate system (Hayes et al., 2011; Gauthier et al., 2015; Natali et al., 2019; Bruhwiler et al., 2021; Virkkala et al., 2021; Braghiere et al., 2023).

During winter months in ABRs, landscapes are generally snow-covered and photosynthesis is considered negligible, and therefore winter CO₂ fluxes derive primarily from soil respiration (Christiansen et al., 2012; Webb et al., 2016). It is expected that complex soil–vegetation–snow interactions will lead to regional and local variability in soil respiration rates across ABRs because of relationships between vegetation types, snow cover, soil properties, soil moisture, and soil temperature (Gouttevin et al., 2012; Busseau et al., 2017; Loranty et al., 2018; Grünberg et al., 2020; Royer et al., 2021). Higher soil temperatures promote microbial activity and increase CO₂ production from soil organic matter decomposition during winter (Natali et al., 2019). A snowpack acts as an important thermal insulative layer for the soil during winter, keeping soils warmer than the ambient air (Domine et al., 2016b). Vegetation affects snow properties by increasing snow depth where wind trapping occurs (Callaghan et al., 2011a, b; Busseau et al., 2017), decreasing snow density and thermal conductivity around shrubs (Gouttevin et al., 2012; Domine et al., 2015, 2016a), decreasing albedo due to protruding branches (Ménard et al., 2012), and causing earlier spring snowmelt due to vegetation thermal conductivity (Wilcox et al., 2019; Kropp et al., 2022). However, Dominé et al. (2022) showed that shrub branches within the snowpack can contribute to mid-winter soil cooling by conducting temperature through the snowpack. Hence, the complex vegetation–snow–soil interactions in northern environments must be considered when studying what drives the spatial variability in soil carbon emissions during winter. Soil microbial activity can also be limited by lack or saturation of available water, meaning that higher amounts of available soil liquid water content (LWC) should allow higher heterotrophic respiration rates by increasing soil microbial activity as long as the soil environment is not anaerobic (Linn and Doran, 1984; Knowles et al., 2015). Anaerobic soil conditions are usually found in fully water saturated soils.

High uncertainties in winter ABR CO₂ exchange between the ground surface and atmosphere are in part due to limited data availability because of difficulties in accessing these vast, remote regions and the harsh winter conditions creating technical challenges for CO₂ fluxes measurements (Natali et al., 2019; Virkkala et al., 2022). Methods currently available to measure wintertime CO₂ fluxes include (1) the eddy covariance technique (Baldocchi, 2003), (2) chamber measurements under or above the snowpack (McDowell et al., 2000), and (3) snowpack gradient diffusion methods (Sommerfeld et al., 1993). Each of these has their advantages and limitations. The eddy covariance (EC) technique exploits the atmosphere's turbulent nature to estimate net CO₂ fluxes at high temporal resolution without environmental disturbance (Baldocchi et al., 2001; Pastorello et al., 2020). Data gaps are common during the ABR winter since the EC equipment is energy-intensive and prone to failure in low temperatures. In addition, solar power supply systems are limited by low sunlight (Jentzsch et al., 2021; Pallandt et al., 2022). Furthermore, the EC equipment is stationary and covers a large footprint (250–3000 m). In contrast, plot-scale chamber techniques for measuring CO₂ fluxes are portable methods with a small footprint (<1 m) (Subke et al., 2021; Maier et al., 2022). Chambers can be used either above the snowpack or directly on the ground. Placing a chamber on the snowpack does not provide a direct measurement of soil CO₂ fluxes due to CO₂ retention and lateral diffusion within snowpacks, generally creating a negative bias and uncertainties linking the snow and/or atmosphere fluxes to soil fluxes (McDowell et al., 2000; Björkman et al., 2010; Webb et al., 2016). Chambers can also be placed directly on the ground by excavating the snow cover (Elberling, 2007), providing a direct measurement of soil CO₂ fluxes that is, however, prone to a positive bias generated by a tunnel effect due to the snow excavation (McDowell et al., 2000; Björkman et al., 2010). Unavoidable snow cover disturbance also reduces the possibility of revisiting locations for temporal surveys because the soil thermal regime is altered by the snow disturbance. Alternatively, permanent chambers can be installed before the first snowfall, but this disturbs the state of the ground and snow cover around the chamber (Webb et al., 2016). The snowpack diffusion gradient method uses snow porosity and tortuosity to estimate CO₂ fluxes from the gas concentration gradient along a vertical snow profile including ambient air above the snowpack (Sommerfeld et al., 1993; Pirk et al., 2016; Kim et al., 2019). In this study, the snowpack diffusion gradient method will be used to evaluate the spatial variability in CO₂ fluxes in ABRs because of its portability and minimal environmental disturbance.

The goal of this study was to determine the main environmental controls of winter CO₂ fluxes in ABRs. A total of 560 snowpack diffusion gradient measurements were made over a latitudinal gradient of four different ecosystem types common in ABRs in Canada: closed-crown coniferous boreal forest, open-crown coniferous boreal forest, erect-shrub tundra, and prostrate-shrub tundra. Spatio-temporal measurements of snowpack CO₂ diffusion gradients were performed at several locations in the four sites during the 2020–2021 and 2021–2022 winters (December to May). Firstly, a random forest (RF) machine learning analysis was used to evaluate the relative importance of the following environmental variables known to exert control over winter CO₂ fluxes: soil temperature, soil LWC, vegetation type, snow water equivalent, snow depth, and several snow-density-related measurements. Secondly, the response of and uncertainty in winter CO₂ fluxes to the most impactful environment variables determined by the RF model were quantified through regression analysis.

2 Method

2.1 Study sites

To cover different vegetation types and a wide range of soil temperature (T_soil) regimes and snow conditions found in ABRs, four study sites were selected across Canada (Fig. 1 and Table 1). Each site represents a specific ecosystem type (Royer et al., 2021), and vegetation types within each of those ecosystems were determined using vegetation maps specific to each site. Cambridge Bay (CB), situated on Victoria Island in the Canadian Arctic Archipelago was the northernmost site located in the Arctic tundra, dominated by lichen and prostrate-shrub tundra. Ponomarenko et al. (2019) generated a detailed ecotype map of the Arctic tundra biome present in the CB study area. Here, these ecotypes were grouped by water availability into three tundra vegetation types from which the sampling locations (S) were selected: dry (S=94), sub-hydric (S=24), and hydric (S=110). Trail Valley Creek (TVC), Northwest Territories, situated just north of the treeline in the transitional zone between the boreal and Arctic biomes close to the Mackenzie Delta, is dominated by erect-shrub tundra with remaining tree patches (Martin et al., 2022). Grünberg et al. (2020) produced a vegetation map of the TVC study area using airborne orthophotos, vegetation height, and field observations from which seven vegetation types and landforms were identified: lichen (S=68), tussock (S=21), dwarf shrub (S=19), tall shrub (S=26), polygon (S=21), riparian shrub (S=17), and black spruce tree patch (S=18). Havikpak Creek (HPC) is located just south of the treeline, at about 50 km south of TVC in an open-crown black-spruce-dominated forest constituting the only type of vegetation present (S=30) (Krogh et al., 2017). Montmorency Forest (MM) is the southernmost site located in a closed-crown balsam-fir-dominated boreal forest constituting the only type of vegetation present (S=110) (Barry et al., 1988). The CB, TVC, and HPC sites are underlain by continuous permafrost, while the MM site is permafrost-free.

Figure 1 Study site locations in Canada. The Arctic biome is delimited following the Conservation of Arctic Flora and Fauna (CAFF) working group of the Arctic Council, and the boreal biome is delimited following Potapov et al. (2008). Permafrost extent (Brown et al., 2002) is estimated in percent area: continuous (>90 %–100 %), discontinuous (>50 %–90 %), sporadic (10 %–50 %), and isolated patches (<10 %).

Table 1 Study sites with the number of sampling locations in Canada and the number of CO₂ flux measurements (N) for each site.

2.2 Snowpack diffusion gradient method

2.2.1 Theoretical framework for CO₂ flux calculation

During winter in ABRs, soil respiration produces CO₂ below the snowpack. Consequently, a vertical CO₂ diffusion gradient is maintained through the snowpack ($\mathrm{d}\left[{\mathrm{CO}}_{\mathrm{2}}\right]/\mathrm{d}z$; g C m⁻⁴), with CO₂ concentration ([CO₂]; g C m⁻³) decreasing with snow height from the soil surface (z; m) (Jones et al., 1999). Hereafter, [CO₂] is expressed in grams of carbon per cubic meter but units of concentration could also be expressed in relative units (i.e., parts per million – ppm) using the ideal gas law. The snowpack diffusion gradient method uses $\mathrm{d}\left[{\mathrm{CO}}_{\mathrm{2}}\right]/\mathrm{d}z$ within the snowpack and Fick's first law for gas diffusion through porous media to estimate CO₂ fluxes ($F_{{\mathrm{CO}}_{\mathrm{2}}}$; $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) (Sommerfeld et al., 1993; Zhu et al., 2014):

$$\begin{array}{}\text{(1)}& F_{{\mathrm{CO}}_{\mathrm{2}}}=-\mathit{\phi }\mathit{\tau }{D}_{\mathrm{a}}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{CO}}_{\mathrm{2}}\right]}{\mathrm{d}z}},\end{array}$$

where φ represents the porosity of the snow medium, τ its tortuosity, and D_a the air diffusion coefficient of the diffused gas in square meters per day. The porosity of snow can be assessed from its density (Kinar and Pomeroy, 2015):

$$\begin{array}{}\text{(2)}& \mathit{\phi }=\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\text{snow}}}{{\mathit{\rho }}_{\text{ice}}}}+\mathit{\theta }\cdot \left({\displaystyle \frac{{\mathit{\rho }}_{\text{water}}}{{\mathit{\rho }}_{\text{ice}}}}-\mathrm{1}\right),\end{array}$$

where ρ represents the density of snow, water, and pure ice (${\mathit{\rho }}_{\text{ice}}=-\mathrm{0.0001}\cdot T_{\text{ice}}+\mathrm{0.9168}$ with T_ice as ice temperature in ^∘C and ρ_ice in g cm⁻³; Harvey et al., 2017) and θ is the snow liquid water content. The tortuosity is strongly correlated with porosity. Du Plessis and Masliyah (1991) established the following relationship:

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

Tortuosity can also be approximated as $\mathit{\tau }\approx {\mathit{\phi }}^{\mathrm{1}/\mathrm{3}}$ (Millington, 1959; Mast et al., 1998). $\mathrm{d}\left[{\mathrm{CO}}_{\mathrm{2}}\right]/\mathrm{d}z$ is obtained by measuring the [CO₂] vertical profile at various snow depths. Standard air diffusion coefficients of CO₂ (unit: m² d⁻¹) are available in the literature but must be corrected for temperature and pressure (Marrero and Mason, 1972; Massman, 1998):

$$\begin{array}{}\text{(4)}& D_{\mathrm{a}}=\mathrm{0.2020}\cdot {\left({\displaystyle \frac{T}{T_{\mathrm{o}}}}\right)}^{\mathrm{1.590}}\cdot e^{-\frac{\mathrm{0.3738}}{T/T_{\mathrm{o}}}},\end{array}$$

where T is the air temperature and T_o is the freezing point in kelvins. The diffusion gradient method assumes that gas fluxes are the result of simple, linear, gradient-induced diffusion in uniform porosity through snow cover (McDowell et al., 2000). A snowpack with strongly heterogeneous density (i.e., vertical stratification) can induce a bias when gas flow is altered by dense layers or ice crusts, typically leading to $F_{{\mathrm{CO}}_{\mathrm{2}}}$ overestimation (Seok et al., 2009). Such layers were rarely found in our study sites. The diffusion gradient assumption also does not hold when strong wind events occur, decreasing snowpack CO₂ concentration through wind-pumping and inducing a negative bias in CO₂ fluxes (Seok et al., 2009). Consequently, $\mathrm{d}\left[{\mathrm{CO}}_{\mathrm{2}}\right]/\mathrm{d}z$ was not measured in days following a strong wind event. Monitoring of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ at a few sampling locations did not show any relationship between $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and wind speed or atmospheric pressure (e.g., Fig. A1).

2.2.2 Data collection

All data were collected during the 2020–2021 and 2021–2022 winters between December and May (Table 1). The CO₂ concentration gradient was measured by collecting gas samples at various depths in the snowpack. Each gradient profile consisted of five gas samples collected at (1) 5 cm above the snowpack (ambient air), (2) 5 cm depth below the snowpack surface, (3) $\mathrm{1}/\mathrm{3}$ of total snow depth, (4) $\mathrm{2}/\mathrm{3}$ of total snow depth, and (5) the soil–snow interface. Gas present in snow pores was collected with a thin, hollow, stainless-steel rod (50–120 cm long, 4 mm outer diameter, and 2 mm inner diameter), starting with gas samples in the upper snowpack and then pushing the sampling rod downward to collect gas samples deeper in the snowpack to minimize snow disturbance (Fig. 2a). Gas was collected in a 60 mL syringe (Air-Tite Luer Lock, Virginia Beach, Virginia) connected to the rod via a three-way valve. Gases were transferred into 12 mL hermetic glass vials (Labco Exetainer^®, Labco Ltd., Lampeter, UK), which were sent to the Université du Québec à Trois-Rivières laboratory to be measured with a gas analyzer to obtain CO₂ concentrations. At each site, several sampling locations were selected to cover the maximum range of vegetation types and snowpack characteristics, covering areas of 0.05–22.5 km². At each sampling location, two to four replicate profiles were measured at 50 cm spacing to test the repeatability of the sampling. A minimal spacing of 57.5 cm was required between sampling positions, since it corresponds to the radius of the 60 mL sampling volume of each gas sample, based on a snow density range of 100–650 kg m⁻³.

Figure 2 (a) Gas sampling equipment for the CO₂ concentration gradient measurement. (b) Typical snow depth profile of an Arctic snowpack (picture from Trail Valley Creek close to a tree patch).

For typical Arctic snowpacks, samples at $\mathrm{1}/\mathrm{3}$ depth are usually in wind slabs, the dense and cohesive surface snow layer formed by strong Arctic winds. Samples at $\mathrm{2}/\mathrm{3}$ depth are usually in depth hoar, the lower snow layer with low density and cohesion formed by a strong temperature gradient driving vertical vapor flux through the snowpack (Fig. 2b). Typically, boreal snowpacks are deeper than in Arctic tundra and display a more continuous vertical stratification with increasing snow density at the bottom of the snowpack. In HPC, snowpack depths were 40–80 cm in March, while snowpack depths at MM were 100–200 cm (Fig. A2). For comparison, by March, snowpacks at CB were 10–75 cm deep and 15–150 cm at TVC.

Once the gas samples were collected, a vertical profile of snow and soil properties was measured to calculate the CO₂ air diffusion coefficient from the snow temperature, as well as snow porosity and tortuosity from snow density. Snow properties were measured every 5 cm, including snow temperature (Snowmetrics digital thermometer; Fort Collins, Colorado; $\mathrm{1}/\mathrm{10}{}^{\circ }$ resolution), snow density (Snowmetrics digital scale, 100 and 250 cm³ snow cutters used to weigh snow samples; σ(ρ_snow)≈9 %; Proksch et al., 2016), snow liquid water content (hand test from Fierz et al., 2009), and snow stratigraphy. Examples of snow vertical stratification along with CO₂ concentration measurements can be found in Appendix A (Fig. A3). T_soil was measured at 1 cm depth under the soil–snow interface as it was not possible to go deeper in frozen soil and no permanent sensors were installed (Snowmetrics digital thermometer; Fort Collins, Colorado; $\mathrm{1}/\mathrm{10}{}^{\circ }$ resolution); three measurements of T_soil were averaged. Snow depth measurements were done with a ruler graduated every 1 cm (σ(d_snow)≈0.5 cm).

The CO₂ concentration of 86 % of gas samples was measured using a LI-COR LI-7810 ${\mathrm{CH}}_{\mathrm{4}}/{\mathrm{CO}}_{\mathrm{2}}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ Trace Gas Analyzer (LI-COR Biosciences, Lincoln, Nebraska; σ<1 %; N=483). The gas samples were passed through an open loop along a continuous flow of a 200 ppm CO₂ calibration gas (Linde Canada, Ottawa, Ontario, Canada). Based on a calibration curve using 0, 400, and 1000 ppm CO₂ calibration gases (Linde Canada), the CO₂ concentration of gas samples was calculated (detailed protocol: https://www.licor.com/documents/xst0ld9jozfby78bmpdqi9i7rmjjjjmg, last access: 9 December 2023).

Randomly distributed gas samples collected during the 2020–2021 winter were analyzed with a Picarro G2201-i CRDS gas analyzer (Picarro, Santa Clara, California; σ<0.1 %; N=26) to validate the method used with the LI-7810 to determine CO₂ concentration. CO₂ concentrations estimated from the LI-7810 and Picarro G2201-i gas analyzers were not significantly different in their concentration range and distribution (Fig. A4; R²=0.92). At TVC in March 2022, a portable LI-850 ${\mathrm{CO}}_{\mathrm{2}}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ Gas Analyzer was used (σ<1.5 %; N=38), allowing for CO₂ concentrations to be measured on the same day as sample collection (avoiding the need for bottling and transportation). CO₂ concentrations estimated from the LI-7810 and LI-850 gas analyzers were not significantly different in their concentration range and distribution (Fig. A4b; R²=0.82).

2.2.3 Evaluation of CO₂ flux uncertainties

An uncertainty assessment was conducted to evaluate CO₂ flux precision based on the snowpack diffusion gradient method. The uncertainty assessment focuses on random errors, as systematic errors are discussed at the end of Sect. 2.2.1. From sampling to flux estimation, several steps could add uncertainty to the results. Uncertainties can be subdivided into four sources: gas concentration estimates; gas transfer, transport, and storage; evaluation of the snowpack $\mathrm{d}\left[{\mathrm{CO}}_{\mathrm{2}}\right]/\mathrm{d}z$; and snow pit measurements (i.e., snow density and temperature). Gas concentration uncertainties were evaluated from the gas analyzer precision as assessed by the manufacturer and tested using calibration gases. Six CO₂ reference gases of 400 ppm were bottled during two different field campaigns and were processed among the gas samples from the snowpack to ensure the transfer, transport, and storage protocol did not lead to sample contamination. The $\mathrm{d}\left[{\mathrm{CO}}_{\mathrm{2}}\right]/\mathrm{d}z$ uncertainties were evaluated with the standard deviation from the coefficient of determination ($\mathit{\sigma }=\sqrt{(\mathrm{1}-R^{\mathrm{2}})/(N-\mathrm{1})}$; Bowley, 1928). $F_{{\mathrm{CO}}_{\mathrm{2}}}$ uncertainty was estimated by propagation of the uncertainties in $\mathrm{d}\left[{\mathrm{CO}}_{\mathrm{2}}\right]/\mathrm{d}z$ and snow density using Eq. (1) (Taylor, 1997). The uncertainty in ρ_snow was fixed at 9 % (Proksch et al., 2016), while the uncertainty in $\mathrm{d}\left[{\mathrm{CO}}_{\mathrm{2}}\right]/\mathrm{d}z$ was estimated based on the root mean squared error of the linear regression for each snowpack concentration gradient measurement.

2.3 Soil volumetric liquid water content at the Montmorency Forest site

Conditions of a 0 ^∘C curtain exist when the soil temperature is around freezing point (0 ^∘C) and a mix of ice and liquid water coexist in the soil pore space because the phase transition between water and ice is slowed due to latent heat (Outcalt et al., 1990). Hence, liquid water content (LWC; m³ m⁻³) and ice fractions can be used as a freezing–thawing indicator during the zero-curtain period. The MM study sites were equipped with TEROS 12 soil moisture sensors (METER Group) at 5 cm depth. LWC was only monitored at the MM site since it was the only site where T_soil in upper layers remained around 0 ^∘C for the whole winter, allowing the presence of liquid water in the soil throughout winter. The Zhang et al. (2010) empirical soil liquid water and ice mixing model was used to calculate soil liquid water content (m_uw) (Eqs. 5 to 8). LWC was estimated to be negligible at the CB, TVC, and HPC sites since T_soil was between −5 and −25 ^∘C. The model from Zhang et al. (2010) supports LWC, at T_soil colder than −5 ^∘C, being negligible.

$$\begin{array}{}\text{(5)}& {\displaystyle }\text{LWC}=a\cdot {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{b}}}{{\mathit{\rho }}_{\mathrm{w}}}}\cdot {\left|T_{\text{soil}}\right|}^{-b},\text{(6)}& {\displaystyle }\begin{array}{r}\ln a=\mathrm{0.5519}\cdot \ln \text{SSA}+\mathrm{0.2618},\\ \ln b=-\mathrm{0.264}\cdot \ln \text{SSA}+\mathrm{0.3711},\end{array}\end{array}$$

where ρ_w and ρ_b (g cm⁻³) represent liquid water and soil bulk density respectively, T_soil (^∘C) represents soil temperature, and SSA (m⁻¹) represents soil particles' specific surface area described by Sepaskhah et al. (2010).

$$\begin{array}{}\text{(7)}& {\displaystyle }\text{SSA}=\mathrm{3.89}\cdot d_{\mathrm{g}}^{-\mathrm{0.905}},\text{(8)}& {\displaystyle }\ln d_{\mathrm{g}}=f_{\mathrm{c}}\cdot \ln M_{\mathrm{c}}+f_{\text{si}}\cdot \ln M_{\text{si}}+f_{\text{sa}}\cdot \ln M_{\text{sa}},\end{array}$$

where d_g represents the soil geometric mean particle-size diameter (mm) and f and M represent soil component fraction and mean particle-size diameter (mm). Soil components are clay (M_c=0.001 mm), silt (M_si=0.026 mm), and sand (M_sa=1.025 mm). Soil bulk density and gravimetry was evaluated using a soil sampling protocol similar to the National Forest Inventory protocol (CFI, 2008). Undisturbed soil samples were collected in three homogenous horizons of a soil profile using 400 cm² cores. Volumetric soil samples were dried (103 ^∘C) and weighed to determine bulk density. Gravimetric samples were used to determine sand (%, 50–2000 µm), silt (%, 2–50 µm), clay (%, <2 µm), and organic content (g kg⁻¹). The soil texture was determined by the hydrometer method (Bouyoucos, 1962), whereas the organic content was determined with a LECO organic analysis instrument (LECO Corporation, St. Joseph, Michigan).

2.4 Random forest algorithm

Random forest (RF) is an ensemble machine learning method based on a multitude of decision trees (Breiman, 2001). Each decision tree of our RF model (scikit-learn 1.2.1 library from Python 3.10.3) is trained on a random subset of environmental variables drawn from the dataset input: T_soil, LWC, vegetation type, snow water equivalent, snow depth, snow mean density, snow maximum density, snow porosity, snow tortuosity, wind slab thickness (if present), and wind slab fraction relative to total snow depth (if present). Each decision tree generates a $F_{{\mathrm{CO}}_{\mathrm{2}}}$ prediction, and the overall RF prediction is the average of all prediction trees. A strength of the RF algorithm is that it performs well even when input variables are correlated with each other (Liaw and Wiener, 2002; Strobl et al., 2008; Kibtia et al., 2020). Our RF model was composed of 500 fully decomposed decision trees. Our dataset was randomly divided into a training subset (75 %) and a testing subset (25 %), preserving the relative distribution between vegetation types. Our RF model performance was assessed using the coefficient of determination (R²), explained variance, and mean absolute error. We used our RF model to identify the relative importance of winter CO₂ flux predictors. Relative importance of each environmental variable was computed with the permutation method, i.e., alternatively removing variables from the RF model and evaluating the performance decrease, which was measured via the coefficient of determination.

3 Results

3.1 CO₂ flux uncertainties

Evaluation of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ precision showed that the two main sources of uncertainty are associated with snow density measurements, in agreement with Sommerfeld et al. (1996), and with $\mathrm{d}\left[{\mathrm{CO}}_{\mathrm{2}}\right]/\mathrm{d}z$ linear regression (mean R²=0.790 (σ=0.236) for $F_{{\mathrm{CO}}_{\mathrm{2}}}\ge \mathrm{0.01}$ $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; N=398) (Table A1). Snow density uncertainty (σ(ρ_snow)≈9 %) impacted snow porosity and tortuosity in Eq. (1). From the linear fit of Fig. 3, the average $F_{{\mathrm{CO}}_{\mathrm{2}}}$ uncertainty can be estimated at 19.4 %, which provides sufficient accuracy to observe the impact of environmental variables on winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$.

Figure 3 CO₂ flux ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) uncertainty relationship to $F_{{\mathrm{CO}}_{\mathrm{2}}}$ for four study sites and two winters: 2020–2021 and 2021–2022. Specifications of the linear fit can be found in the upper left. The data dot color indicates the study site, and its symbol (i.e., circle or x-shaped) indicates the winter during which it was collected.

The overall [CO₂] precision of around 1 % shows that the measurement technique is not a main source of uncertainty in $F_{{\mathrm{CO}}_{\mathrm{2}}}$ estimates. Gas concentration estimations from the LI-7810 have a precision of 0.88 % at 400 ppm according to the manufacturer. The stability of the [CO₂] measurement was evaluated over 169 measurements, displaying a standard deviation of 0.09 %. The LI-7810 was further tested using a 400 ppm calibration gas with a 1 % [CO₂] precision (Linde Canada). A linear calibration fit equation was used to estimate [CO₂] of small gas samples, using three calibration gases (200, 400, 1000 ppm) plus the theoretical zero intercept. Average uncertainty in the linear regression was 0.76 % over six calibration runs with a standard deviation of 0.15 %. The average accuracy of the reference [CO₂] bottled among the gas samples from the snowpack was 1.11 %.

3.2 Spatio-temporal variability in winter CO₂ fluxes associated with abiotic controls

The RF model determined T_soil and LWC to be the two main predictors of winter CO₂ fluxes. We found two temperature and LWC regimes of winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ (Fig. 4). The first regime was when the soil was frozen with T_soil<0 ^∘C and LWC<0.2 m³ m⁻³, leading to $F_{{\mathrm{CO}}_{\mathrm{2}}}$ being mainly controlled by T_soil. The second regime was when LWC>0.2 and <0.42 m³ m⁻³ but with a fraction of its water in the form of ice (zero-curtain condition), causing LWC to be the main control of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ instead of T_soil. While the first regime mostly corresponds to Arctic study sites, the second regime only includes one study site (MM) located in the southern boreal forest. Therefore, conclusions from the second regime should be less generalized than those from the first regime. Subsequent evaluation focused on the response of winter CO₂ fluxes to T_soil and LWC using exponential regressions in order to better understand the role of these two variables in winter CO₂ fluxes.

Figure 4 Random forest (RF) performance and variable relative importance. Variables used are soil temperature (T_soil), soil liquid water content (LWC), snow density (ρ_snow), snow water equivalent (SWE), and vegetation type. (a) The first iteration integrated the complete dataset, and (b) the second iteration only integrated the Montmorency Forest dataset with LWC>0 m³ m⁻³. The values displayed by the bar plot are the mean variable relative importance over 100 permutations, while the error bars are the standard deviation.

3.2.1 Variable importance determined by the random forest model

T_soil was the $F_{{\mathrm{CO}}_{\mathrm{2}}}$ predictor with the highest relative importance (68 %) when using the complete dataset (Fig. 4a), followed by LWC (17 %). Snowpack characteristics, ρ_snow (11 %) and snow water equivalent (SWE) (2 %), had a lower relative importance in the RF model. Contrary to what might be expected, the vegetation type had near-negligible relative importance (1 %) in $F_{{\mathrm{CO}}_{\mathrm{2}}}$ prediction. The RF model was developed starting with all environmental variables available: T_soil, LWC, vegetation type, SWE, snow depth, mean ρ_snow, max ρ_snow, φ, τ, wind slab fraction, and wind slab thickness. Although the correlation of several snow parameters did not decrease the RF model performance, snow parameters impacted the assessment of variable relative importance by splitting the relative importance between the correlated variables. Consequently, variables with lower importance and with no significant impacts on the RF performance were progressively removed. The two selected snow parameters that had significant impact were SWE and ρ_snow. MM was the only site where soil LWC was present, enabling the assessment of the relative importance of this variable. When using only data from MM in the RF model (Fig. 4b), the relative importance of T_soil (12 %) for $F_{{\mathrm{CO}}_{\mathrm{2}}}$ was lower than with all combined datasets since T_soil was near 0 ^∘C for all measurements. At MM, LWC becomes the main predictor (87 %) of $F_{{\mathrm{CO}}_{\mathrm{2}}}$, while ρ_snow importance drops (2 %) and SWE importance remains similar (<1 %).

Figure 5 CO₂ flux ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) as a function of surface soil temperature (T_soil) for T_soil<0 ^∘C. An exponential regression was fitted with the data (black line) and compared to the exponential regression by Natali et al. (2019) from an external dataset (blue line).

3.2.2 Soil temperature

Figures 5 and 6 show the relationship between winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and T_soil. Figure 5 focuses on T_soil<0 ^∘C from CB, TVC, HPC, and MM. An exponential regression was used to evaluate the relationship between T_soil and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ estimates (RMSE=0.024 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). $F_{{\mathrm{CO}}_{\mathrm{2}}}$ values at MM when T_soil<0 ^∘C and LWC <0.2 m³ m⁻³ were included in this graph because they are more strongly correlated to T_soil than LWC (see Sect. 3.2.3). Note that the low number of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ measurements with T_soil between −6 and −0.5 ^∘C restrict the capacity to evaluate the regression within this range. Using the exponential regression of Natali et al. (2019), we obtained an RMSE of 0.030 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. The regression of Natali et al. (2019) generally shows an overestimation of fluxes for $T_{\text{soil}}<-\mathrm{5}$ ^∘C but an underestimation for T_soil>5 ^∘C when compared to our exponential regression. The systematic bias between our dataset and the regression of Natali et al. (2019) is minimal (mean bias = −0.0025 $\mathrm{g}\mathrm{C}{\mathrm{d}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$). We also observed the isolated occurrence of comparably large winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ of up to 0.36 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ at temperatures below −10 ^∘C (Fig. 5). These measurements of high $F_{{\mathrm{CO}}_{\mathrm{2}}}$ at low temperature seem to be genuine since the repeatability was verified over the three sampling profiles performed at each site. Nevertheless, we were not able to explain these strong $F_{{\mathrm{CO}}_{\mathrm{2}}}$ fluxes and no environmental variables measured in our study could be linked to those occurrences. It has been suggested that gas bursts during autumn freeze-up in permafrost environments might be due to gas compression by ice formation and ground cracking (Pirk et al., 2015). This hypothesis can be considered to explain the high $F_{{\mathrm{CO}}_{\mathrm{2}}}$ observed in this study, although the high $F_{{\mathrm{CO}}_{\mathrm{2}}}$ observed occurred at a near-surface T_soil between −25 and −10 ^∘C, so the freeze-up would have to occur at lower depths in the soil. Figure 6 displays the higher winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ from MM where T_air values are higher and the important snowpack insulation keeps the soil at temperatures of around 0 ^∘C through the entire winter. $F_{{\mathrm{CO}}_{\mathrm{2}}}$ increases more rapidly with T_soil around freezing point than at T_soil<5 ^∘C, which is shown by the higher temperature-dependency parameter (B=2.82 ${}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$) of the MM site exponential regression (RMSE=0.286 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) compared to the exponential regression of Fig. 5 (B=0.18 ${}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$). This discrepancy in temperature dependency creates a discontinuity between the measurements at T_soil<5 ^∘C and T_soil≈0 ^∘C that did not allow for a continuous temperature-dependency regression across all the study sites. The lower RMSE of the exponential regression of Fig. 5 (RMSE=0.024 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; 70.3 % of mean $F_{{\mathrm{CO}}_{\mathrm{2}}}$) compared to the exponential regression of the MM site (RMSE=0.286 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; 112.4 % of mean $F_{{\mathrm{CO}}_{\mathrm{2}}}$) might be due to the impact of soil LWC at the MM site (see Sect. 3.2.3).

Figure 6 CO₂ flux ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) as a function of soil temperature (T_soil) at the Montmorency Forest study sites where soil liquid water content (LWC) was greater than 0 m³ m⁻³ throughout winter. An exponential regression was fitted to the data (black line).

3.2.3 Soil liquid water content

The relationship between LWC and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ during winter at MM (RMSE=0.137 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; 49.1 % of mean $F_{{\mathrm{CO}}_{\mathrm{2}}}$) was stronger than between T_soil and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ (RMSE=0.286 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; 112.4 % of mean $F_{{\mathrm{CO}}_{\mathrm{2}}}$), when excluding the sampling location that contained a thick organic soil layer with very high soil moisture due to its location near the bottom of a microtopographic depression (Fig. 7). Other MM sampling locations with a thin organic layer shared a similar soil composition dominated by mineral soils. The strong correlation between LWC and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ was mostly observed at LWC>0.2 and <0.42 m³ m⁻³. The plateau observed in Fig. 7 indicates that T_soil might be a better predictor than LWC at LWC<0.2 m³ m⁻³.

Figure 7 CO₂ flux ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) as a function of soil volumetric liquid water content (LWC) at the Montmorency Forest study site. An exponential regression was fitted to the data (black line), excluding the thick-organic-layer site (red markers).

3.2.4 Vegetation types

Figure 8 shows winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ across the four study sites for different vegetation types. Since CB vegetation is mostly prostrate-shrub tundra, CB ecosystems were regrouped by water availability. On average, higher winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ values at CB were observed in environments experiencing wetter conditions during the growing season. At TVC, several vegetation and land cover types are present. $F_{{\mathrm{CO}}_{\mathrm{2}}}$ values from MM were higher than for the other sites. Higher $F_{{\mathrm{CO}}_{\mathrm{2}}}$ can be explained by warmer mean annual average temperature, a deeper snowpack, and winter T_soil around 0 ^∘C (see Sect. 3.4).

Figure 8 Boxplot of CO₂ flux ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) across 12 vegetation types and 4 sites. $F_{{\mathrm{CO}}_{\mathrm{2}}}$ values from Montmorency Forest (MM) are on a separate scale because they are much higher than the colder environments (b). Cambridge Bay (CB) sites are ordered by increasing water availability, and Trail Valley Creek (TVC) sites are ordered by increasing mean soil surface temperature in March 2021 and 2022 (a). Havikpak Creek (HPC) and MM were composed of a single vegetation type (a). Outliers were defined as $F_{{\mathrm{CO}}_{\mathrm{2}}}>\text{Q3}+\mathrm{1.5}$ IQR, where Q3 is the third quartile and IQR the interquartile range. Outliers are out of the y-axis range for the dry tundra (four), sub-hydric tundra (four), hydric tundra (four) and lichen (three). The outliers can be found in Fig. 5.

Vegetation type was not identified as a strong predictor of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ by the RF model. Nonetheless, we observed differences in the mean and range of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ for the various vegetation types probed in this study. This might be due to the strong correlation between vegetation type and soil temperature (Fig. A5), as well as relationships between vegetation and soil type, including soil organic matter content and soil pore size. The RF algorithm showed vegetation type relative importance increased to 42 % when T_soil was removed from the environmental variables, although the removal of T_soil decreased RF performance substantially (R²=0.40). Therefore, vegetation could be used as a proxy variable for T_soil if the latter is not available to predict $F_{{\mathrm{CO}}_{\mathrm{2}}}$, but with poorer results.

4 Discussion

4.1 Controls of winter CO₂ fluxes

The RF model predictors' relative importance showed that during winter, T_soil emerged as the dominant predictor of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ when the soil was frozen. Nevertheless, in the closed-crown boreal forest site (i.e., MM) where zero-curtain conditions persisted throughout winter, soil LWC took precedence as the dominant predictor as there was minimal variation in T_soil under these conditions. Our results confirm the strong winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ dependency on T_soil shown by Natali et al. (2019), although we observed fluxes lower than reported by Natali et al. (2019) at $T_{\text{soil}}<-\mathrm{5}$ ^∘C and mostly higher fluxes at $T_{\text{soil}}>-\mathrm{5}$ ^∘C. Considering the two regressions of the relationship between T_soil and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ have large uncertainties attached to them, the difference between them falls inside the uncertainty margin (Fig. 5). It should be noted that the Natali et al. (2019) regression was obtained using $F_{{\mathrm{CO}}_{\mathrm{2}}}$ estimates from several methods including eddy covariance, chamber methods, and snowpack diffusion measurements, whereas our study exclusively uses the latter. Several studies have shown bias between the different measurement methods; eddy covariance and soil chamber methods displayed positive biases when compared to snowpack diffusion measurements (McDowell et al., 2000; Björkman et al., 2010; Webb et al., 2016), while the snow chamber displayed negative biases when compared to the snowpack diffusion measurements (McDowell et al., 2000). It should be remembered that the T_soil used in the latter study refers to near-surface temperature; T_soil at greater depths may vary and affect the correlation with $F_{{\mathrm{CO}}_{\mathrm{2}}}$.

4.2 Zero-curtain conditions

Soil LWC was observed only at the MM site, where T_soil was around 0 ^∘C throughout winter. In zero-curtain conditions, LWC was shown to become the dominant control of winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$, while T_soil importance diminished. It should be noted that it would be ill-advised to generalize the relationship between soil LWC and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ as it is only based on data from one study site, and it cannot be ruled out that this relationship is site-specific depending on soil and vegetation composition. Nevertheless, our study highlighted the important impact of LWC on $F_{{\mathrm{CO}}_{\mathrm{2}}}$ around the soil freezing point when there is a mixed state of ice and free water in soils. When the soil is under zero-curtain conditions, the latent heat governs the ice and liquid water ratio in the soil (Devoie et al., 2022). Hence, LWC and ice fractions can be used as a freezing–thawing indicator during the zero-curtain period and help better quantify the $F_{{\mathrm{CO}}_{\mathrm{2}}}$ fluxes in boreal forest environments where zero-curtain conditions prevail (Prince et al., 2019). This result is particularly important in ABRs since the duration and frequency of zero-curtain periods are expected to increase in a warming climate (Yi et al., 2015, 2019; Tao et al., 2021). Further research on winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ in zero-curtain conditions should investigate different sites to assess if the relationship between $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and soil LWC is site-specific or dependent on soil properties. It should be noted that one of the measurement locations at MM displayed low $F_{{\mathrm{CO}}_{\mathrm{2}}}$ despite its LWC being the highest of all sites. The soil composition of this site consisted of a thick (>30 cm) soil organic top layer, whereas all other measurements were done at sites with thinner (3–10 cm) organic layers on top of mineral soil. It is well known that anaerobic conditions created by high soil moisture (at least >50 %) constrain soil CO₂ respiration rates during the growing season because many microorganisms require oxygen for organic matter decomposition, which they lack if soil pores are filled with water (Linn and Doran, 1984; Davidson and Janssens, 2006).

4.3 Snowpack importance

Our study shows that abiotic variables related to T_soil, LWC, and physical snowpack properties explain the majority of variance in winter CO₂ fluxes. It should be noted that we did not incorporate variables related to temporal dynamics such as the previous days' soil temperature and LWC, which have been shown by Harel et al. (2023) to be of importance during the growing season. However, winter soil variables are not expected to be as dynamic as during the growing season because of the snowpack insulating properties. The RF model showed that SWE and mean snow density were the snow characteristics that provided the greatest improvement of the RF model, although to a lesser degree than T_soil and LWC. The importance of snow characteristics in $F_{{\mathrm{CO}}_{\mathrm{2}}}$ is linked with the strong correlation to T_soil (Domine et al., 2016b; Pedron et al., 2023), although a snapshot of snow conditions provides limited abilities to infer T_soil, as shown in Slater et al. (2017). Snow properties' temporal information is required to predict the impact of snow insulation on T_soil, with the most important period being in the autumn freeze-up when air temperature decreases below the freezing point. Snow characteristics are closely linked to topography (Meloche et al., 2021) and thus soil wetness and soil carbon content (Gouttevin et al., 2012). Regarding the snowpack diffusion gradient method, the snowpack is used to estimate winter CO₂ fluxes. An average snow density was used to estimate snow porosity and tortuosity used in CO₂ flux calculations (Eq. 1), which does not consider the vertical stratification of the snowpack. However, the diffusion gradient remained linear despite vertical stratification in snow density (e.g., Fig. A3 where the average ratio between the standard deviation and mean of $D_{\text{air}}\cdot \mathit{\phi }\cdot \mathit{\tau }$ is around 10 %), which points toward a minimal impact of this assumption on our results.

4.4 Soil biogeochemistry

The unexplained variance (16 %) suggests that winter CO₂ fluxes might have been controlled by other environmental variables such as soil physical–chemical properties regulating soil biogeochemistry and soil redox conditions, which were neither addressed nor measured in this study. CO₂ production is governed by the availability and quality of labile C compounds regulating the decomposition of soil organic matter (Michaelson and Ping, 2005; Wang et al., 2011) and by the activity and composition of the soil microbial community (Monson et al., 2006). Soil type and structure, for example the thickness of the organic layer, soil pore size distribution, and soil pH, may be further strong controls on CO₂ production (Steponavičienė et al., 2022; Yli-Halla et al., 2022). All these variables vary widely across the heterogeneous tundra terrain (Virtanen and Ek, 2014), where small-scale moisture, vegetation, and soil conditions occur among hummock and inter-hummock depressions (Wilcox et al., 2019). Further analysis is required to understand the role of physical–chemical soil properties in $F_{{\mathrm{CO}}_{\mathrm{2}}}$ during winter.

4.5 Relevance for terrestrial biosphere models

Large uncertainties remain in terrestrial biosphere models used to estimate CO₂ fluxes in ABRs (Fisher et al., 2014; Tei and Sugimoto, 2020; Birch et al., 2021; Virkkala et al., 2021), especially regarding the respiratory release of CO₂ via soil respiration (the sum of heterotrophic respiration and belowground autotrophic respiration) during winter (Natali et al., 2019). The limited number of observational data available has restricted model improvements, testing, and evaluation (Virkkala et al., 2022). Modeling the ABR carbon cycle is critical for climate projections since a warmer climate should lead to higher T_soil, thus increasing ABR winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ (Mellander et al., 2007; Throop et al., 2012; Wieder et al., 2019). Several terrestrial biosphere models are currently in use (Fisher et al., 2022), such as CLM (Community Land Model; Lawrence et al., 2019) and CLASSIC (Canadian Land Surface Scheme Including Biogeochemical Cycles; Melton et al., 2020; Seiler et al., 2021). The $F_{{\mathrm{CO}}_{\mathrm{2}}}$ relationships to T_soil and LWC observed in this study could be used to inform terrestrial biosphere models through the parametrization of winter soil respiration sensitivity to soil temperature (e.g., Q10 temperature coefficient) and LWC in zero-curtain conditions. Our study shows that permanent installation of the snow gradient method (Seok et al., 2009; Zhu et al., 2014; Graham and Risk, 2018) would be suitable for gathering the temporal non-growing-season CO₂ fluxes in ABRs required to fully test terrestrial biosphere models.

5 Conclusion

Our study showed that T_soil is the main control of winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ at T_soil<0 ^∘C in ABRs. The relative importance analysis of our RF model showed that T_soil was the main predictor of $F_{{\mathrm{CO}}_{\mathrm{2}}}$, followed by LWC. However, we found that at our site maintaining zero-curtain conditions throughout winter, LWC becomes the main control of winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$. We observed non-negligible winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ that may partially offset growing-season CO₂ uptake in ABRs. Consequently, winter $F_{{\mathrm{CO}}_{\mathrm{2}}}$ must be properly estimated in terrestrial biosphere models and climate models. Additionally, future research should focus on linking the effects of abiotic variables on $F_{{\mathrm{CO}}_{\mathrm{2}}}$ during winter, as we determined here, with soil biogeochemistry, microbial functioning, and vegetation.

Appendix A

Table A1 Uncertainty sources in $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and their uncertainty. [CO₂] precision was evaluated at a concentration of 400 ppm.

Figure A1 CO₂ fluxes ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) at a sampling location in the Trail Valley Creek erect-shrub tundra (lichen) between 19 March and 27 March 2022. Atmospheric pressure and wind speed were obtained from Environment and Climate Change Canada's Meteorological Service of Canada meteorological station at Trail Valley Creek (https://climate.weather.gc.ca/historical_data/search_historic_data_e.html, last access: 21 November 2023). The date format is YYYY-MM-DD.

Figure A2 Violin plot of the snow depth range of sites where $F_{{\mathrm{CO}}_{\mathrm{2}}}$ was estimated. The black stripes inside the violins represent data points. The study sites are Cambridge Bay (CB), Havikpak Creek (HPC), Trail Valley Creek (TVC), and Montmorency Forest (MM). CB sites are ordered by increasing hydricity, and TVC sites are ordered by increasing mean soil surface temperature in March 2021 and 2022.

Figure A3 Examples of snow diffusion coefficient ($D_{\text{snow}}=D_{\text{air}}\cdot \mathit{\phi }\cdot \mathit{\tau }$) vertical stratification and CO₂ concentration ([CO₂]) gradient measurements as a function of snow height (h_snow) from the ground level. The coefficient of determination (R²), [CO₂] gradient (m), and y-axis intercept (b) for the linear regressions on the [CO₂] gradient measurements are provided. The ratio between D_snow standard deviation (σ(D_s)) and average ($\stackrel{\mathrm{‾}}{{\mathit{D}}_{\mathbf{s}}}$) is provided in percent. The data come from (a) Montmorency Forest balsam fir closed-crown coniferous boreal forest on 26 February 2021, (b) Cambridge Bay prostrate-shrub tundra (hydric tundra: hydric sedge fen) on 15 April 2022, (c) Trail Valley Creek erect-shrub tundra (lichen) on 26 March 2022, and (d) Havikpak Creek black spruce open-crown coniferous boreal forest on 16 March 2022.

Figure A4 Comparison of winter CO₂ flux calculated from CO₂ concentration estimated with different gas analyzers. The LI-7810 gas analyzer was used as the reference and is compared to a Picarro G2201-i and LI-850. In the Arctic biome (a), the correlation coefficient is 0.924 for the Picarro instrument and 0.821 for the LI-850. In the boreal biome (b), the correlation coefficient is 0.929 for the Picarro instrument.

Figure A5 Violin plot of the soil temperature (T_soil) range of sites where $F_{{\mathrm{CO}}_{\mathrm{2}}}$ was estimated. The black stripes inside the violins represent data points. CB sites are ordered by increasing hydricity, and TVC sites are ordered by increasing soil surface temperature in March 2021 and 2022.