Explaining CO₂ fluctuations observed in snowpacks

Winter soil carbon dioxide (CO₂) respiration is a significant and understudied component of the global carbon (C) cycle. Winter soil CO₂ fluxes can be surprisingly variable, owing to physical factors such as snowpack properties and wind. This study aimed to quantify the effects of advective transport of CO₂ in soil–snow systems on the subdiurnal to diurnal (hours to days) timescale, use an enhanced diffusion model to replicate the effects of CO₂ concentration depletions from persistent winds, and use a model–measure pairing to effectively explore what is happening in the field. We took continuous measurements of CO₂ concentration gradients and meteorological data at a site in the Cape Breton Highlands of Nova Scotia, Canada, to determine the relationship between wind speeds and CO₂ levels in snowpacks. We adapted a soil CO₂ diffusion model for the soil–snow system and simulated stepwise changes in transport rate over a broad range of plausible synthetic cases. The goal was to mimic the changes we observed in CO₂ snowpack concentration to help elucidate the mechanisms (diffusion, advection) responsible for observed variations. On subdiurnal to diurnal timescales with varying winds and constant snow levels, a strong negative relationship between wind speed and CO₂ concentration within the snowpack was often identified. Modelling clearly demonstrated that diffusion alone was unable to replicate the high-frequency CO₂ fluctuations, but simulations using above-atmospheric snowpack diffusivities (simulating advective transport within the snowpack) reproduced snow CO₂ changes of the observed magnitude and speed. This confirmed that wind-induced ventilation contributed to episodic pulsed emissions from the snow surface and to suppressed snowpack concentrations. This study improves our understanding of winter CO₂ dynamics to aid in continued quantification of the annual global C cycle and demonstrates a preference for continuous wintertime CO₂ flux measurement systems.

1 Introduction

The global soil carbon (C) pool stores three times the amount of C as the atmosphere. Organic C reserves in high-latitude soil are disproportionately affected by anthropogenic climate change (IPCC, 2013). Careful assessment of the soil C pool and corresponding fluxes in these often snow-covered, high-latitude regions is critical for understanding the future global C cycle, as increasing global temperatures are likely to stimulate soil CO₂ emissions (Raich et al., 2002).

Cold and wet conditions, like snow cover, pose challenges for measuring wintertime carbon dioxide (CO₂) fluxes (Liptzin et al., 2009), leading many studies to focus on ecosystem respiration during the growing season. For instance, seasonal variation in soil CO₂ fluxes is not always discussed in meta-analyses of global soil C studies, whether or not wintertime measurements were included in individual studies (Scharlemann et al., 2014). Despite this skewed focus, soil CO₂ is still produced throughout the winter, even at −7 ^∘C (Brooks et al., 1996; Coxson and Parkinson, 1987; Flanagan and Bunnell, 1980). In some cases, the insulating snowpack can prevent soils from freezing completely, stimulating soil CO₂ emissions (Grogan and Jonasson, 2006; Larsen et al., 2007; Monson et al., 2006). Further, snow is a porous medium where soil CO₂ emissions easily pools, complicating measurement techniques. There has been an observed decrease in northern hemispheric snow cover and an earlier onset of spring melt since the 1950s as a result of climate change (Dyer and Mote, 2006; IPCC, 2013). Dyer and Mote (2006) indicated that these changes in snow cover are associated with increasing air temperatures and variations in precipitation. Additionally, increasing air temperature results in increased water vapour in the atmosphere (an increase of approximately 7 % in water-holding capacity of air per 1 ^∘C warming), generating more intense precipitation events, including snow storms (Trenberth, 2011). Despite decreases, snow covers 44–53 % of northern hemispheric land area during winter months (Barry, 1992). With the complex interplay between changes in precipitation, temperature, snow cover, and CO₂ emissions in recent and future decades, winter soil CO₂ measurements are important for accurate estimates of annual CO₂ soil respiration (Fahnestock et al., 1999).

There are several methods of measuring CO₂ fluxes through snowpacks including the snowpack gradient technique, chamber methods, and eddy covariance (EC). The snowpack gradient technique is a commonly used technique and, based on Fick's first law of diffusion, uses CO₂ concentration measurements through a vertical profile from the soil to the snowpack surface to calculate flux (McDowell et al., 2000; Seok et al., 2009). This technique minimizes disturbance to the snowpack when compared with the chamber method and does not require homogenous terrain, as for EC. However, the snowpack gradient technique requires many assumptions and cannot easily account for advective transport of CO₂ through snowpacks (McDowell et al., 2000; Seok et al., 2009). Measurement frequencies of wintertime CO₂ fluxes in past gradient studies have ranged widely, from only twice per winter to half-hourly (Liptzin et al., 2009). Measurements of wintertime CO₂ fluxes recorded at a higher frequency (half-hourly) have shown that wintertime CO₂ fluxes can be surprisingly variable, depending more on transport of CO₂ than on microbial variation (Bowling et al., 2009; Seok et al., 2009). This variability presents a problem because it obfuscates any biological sensitivity to environmental drivers. Under what conditions does the soil microbial community thrive over winter? This is difficult to determine when observed variations are caused by abiotic factors. For example, Seok et al. (2009) observed patterns of high temporal variability in wintertime subniveal CO₂ flux, ranging from 0 to 1.2 µmol m⁻² s⁻¹ during a period of relatively steady soil conditions (temperature, moisture) below 0 ^∘C. Steady soil conditions therefore rule out a microbial driving force when variable fluxes were observed. As for advective transport, it does not increase production of CO₂ in soils, but changes the rate of exchange (Bowling and Massman, 2011).

Although we accept the assumption that CO₂ production occurs in snow-covered soils, there are methodological limitations for quantifying this CO₂ production. Transport of this CO₂ out of soils into the overlying media, whether snow or open air, is driven by two main mechanisms: diffusion and advection (also known as bulk flow or mass flow) (Janssens et al., 2001; Roland et al., 2015). The mode of this gas transport through snowpacks affects the timing and magnitude of CO₂ release to the atmosphere and will potentially create significant lags between the times of CO₂ production and emission. Under calm conditions, it is generally accepted that trace gases are transported out of soils and through snowpacks into the overlying atmosphere via diffusion. Explained by Fick's first law, the background theory of diffusion assumes that trace gas transport out of soils or through a snowpack occurs vertically, with the magnitude of fluxes determined by the concentration gradient (Seok et al., 2009). Advective transport from wind, however, can also affect the transport of trace gases such as CO₂ through porous media like soil and snow (Janssens et al., 2001; Kelley et al., 1968).

Studies are increasingly showing that this non-diffusive (advective) mass transport (e.g. wind) through snow is significant, and must be taken into consideration (Bowling and Massman, 2011; Rains et al., 2016), while considering the appropriate timescale. Advective transport of trace gases through naturally permeable media occurs due to variations in atmospheric pressure at the surface and have been studied on both high-frequency timescales (seconds to minutes; Massman et al., 1995) and low-frequency timescales (barometric, Bowling and Massman, 2011). These natural advective flows are ubiquitous and should also be considered on the mid-range timescale of hours to days (Rains et al., 2016). Bowling and Massman (2011) make it clear that wind pumping in the snowpack enhances outward rates of transport. They measured slower bulk air velocities in snow, which fell within the range of 10⁻³ to 10⁻² m s⁻¹, implying that the contribution of advection to trace gas transport through snowpacks was smaller than that of diffusion. Modelling results from Massman et al. (1997) indicate that advective transport can either enhance or diminish fluxes by a wide range of 1.5 to 25 %, and so further studies with field experiment components are required. A more recent study by Bowling and Massman (2011) found enhanced transport of CO₂ beyond diffusive transport by up to 40 % in the short term and 8 to 11 % when considering the snow-covered season as a whole. The net combined effect of advective and diffusive transport in snow environments on CO₂ and other trace gas transport is considered to be an enhancement to diffusive transport. These studies that investigated advective influence on CO₂ transport in snow systems encouraged further study in this area, and so we intended to help fill this gap with our study. To do so, we investigated the mid-range timescale of the re-establishment of consistent CO₂ concentration gradients in the snowpack after a wind-induced disturbance using both field and modelling methods.

Our overarching objective was to help overcome the methodological limitations of quantifying wintertime CO₂ production. Specifically, we aimed to quantify the effects of advective transport of CO₂ in soil–snow systems on the subdiurnal to diurnal (hours to days) timescale and to mechanistically describe these behaviours using a one-dimensional advective–diffusive model adapted for the soil–snow–atmosphere system.

2 Methods

2.1 Continuous automated field monitoring

The primary motivation for establishing our field stations was to determine the relationship between wind speed, snowpack ventilation, and snowpack CO₂ concentration. The site selected is on a plateau in a recovering boreal system at North Mountain, Nova Scotia, Canada, in the Cape Breton Highlands National Park. Wintertime snow patterns at North Mountain allow for snowpacks of up to 3 m, with the last of the snow melting in May or June, depending on the timing and amount of snow in a given year. Average annual air temperature at North Mountain is 5.1 ^∘C (1999–2013). Average winter air temperature is −6.1 ^∘C (January–March, 1999–2013). An insulating snowpack is often established before soils have a chance to freeze completely. Therefore, soils often remain above 0 ^∘C throughout the winter, and over-winter CO₂ production from these soils is very likely (Grogan and Jonasson, 2006; Larsen et al., 2007; Monson et al., 2006), as soils produce CO₂ down to −7 ^∘C (Brooks et al., 1996; Coxson and Parkinson, 1987; Flanagan and Bunnell, 1980). Average annual wind speed is 17.3 km h⁻¹, with highest wind speeds in the winter (20.7 km h⁻¹, January–March, 1999–2013). High winds and variable meteorological conditions (intense snow squalls, freeze–thaw cycles) create varying snow depths within close proximity (tens to hundreds of metres).

Two measurement stations were installed 60 m apart at North Mountain in late 2013, with data collection from 12 November 2013 to 26 March 2014 and 15 April to 29 April 2015. The sites are referred to as NM1 (North Mountain 1: 46^∘49^′7.41${}^{\prime \prime }$ N, 60^∘40^′20.16${}^{\prime \prime }$ W) and NM2 (North Mountain 2: 46^∘49^′9.15${}^{\prime \prime }$ N, 60^∘40^′18.67${}^{\prime \prime }$ W). The key environmental difference between the two sites was the predictably differing snow depth. At each of the two stations, CO₂ concentration through the snow profile was measured at three depths (0, 50, and 125 cm from the soil surface) using Vaisala CARBOCAP^® Carbon Dioxide Probe GMP343 sensors. A Campbell Scientific CR3000 data logger was used at NM1 and a Campbell Scientific CR1000 data logger was used at NM2 to control the instrumentation, recording values every 30 min. To save power and to minimize potential heating impacts, the GMP343 sensors were turned on for 10 min preceding measurement, a measurement was taken averaged over 1 min, and then the sensors were turned off for the remainder of the 30 min interval. Optics heaters of the GMP343 sensors were kept off entirely, as there was a very limited risk of condensation formation in the relatively constant temperature environment of a snowpack. This further reduced potential sensor heat from < 3.5 W (optics heaters on) to < 1 W (optics heaters off). Together, turning the GMP343 sensors off regularly and keeping the optics heaters off at all times minimized any small potential heating impacts of the sensors. Data were collected from the data loggers at the end of the winter. One BP Solar 50 W solar panel and one Discover D12550 12 V battery were used to power each of the two stations. Snow depth was measured at both stations using SR50A Sonic Ranging Campbell Scientific sensors. A Young Wind Monitor (model 05103) anemometer measured wind speed at NM1. Figure 1 gives the general structure of these stations.

Figure 1 Schematic of initial (2014) CO₂ monitoring stations (NM1, NM2) at North Mountain, Cape Breton. Snowpack CO₂ sensors were at 0, 50, and 125 cm within the snowpack (diagram not to scale).

To enhance the field campaign, adjustments were made to the NM2 station for winter 2015 by adding additional CO₂ measurements throughout the vertical profile. Specific measurements recorded at NM2 include CO₂ concentration at 5 cm depth in the soil and soil surface and at 25, 50, 75, and 100 cm above the soil surface (in the snowpack). We continued to record ambient air CO₂ concentration, wind speed, and snow depth. Measurement recording frequency for all measurements was adjusted to hourly for 2015. The profiler system for the enhanced concentration profile experiment contained two Eosense eosGP (dual-channel nondispersive infrared) sensors to measure CO₂ concentrations for select time periods over the 2015 winter. A pump within the station enclosure extracted air samples from the various sampling locations via flexible nylon tubing, carrying the air to the sensor. In-snow and in-soil terminal ends of nylon tubing sampled from 550 mL PVC tubes that had openings covered with high-density polyethylene membranes to exclude liquid water. Data extracted from winter 2015 for analysis ranged from 15 April to 29 April 2015.

2.2 Field data analysis

To examine the degree of concentration decrease after wind ventilation started, we focused on periods in which the likelihood of steady-state gas transport was maximized (initial winter 2014 experiment). This is an assumption of the snowpack gradient technique, and we assumed that disturbance to the snowpack, including snowfall, results in deviations from steady state (McDowell et al., 2000). We extracted data for time periods during which snow depth had not changed more than several centimetres in the previous 3 days, meaning that there had been no melt or appreciable new snow. To do this, we took the rolling 4 h mean of the snow depth values and found the difference between each set of consecutive snow depth values. We retained the values for which the difference of the rolling mean was < 0.001 m. We conducted regression analyses of CO₂ concentration at the three depths and the corresponding wind speeds during these steady-state periods. The ideal situation (the best set of environmental conditions for which a strong negative correlation could be found) was satisfied when winds increased slowly, then abated several hours later (and vice versa). In order to select data where characteristic response patterns of concentration depletion with increasing wind were present, data were further filtered to satisfy the following conditions: (1) the relationship produced a slope < 0, i.e. there was a negative relationship between the two variables, and (2) R² ≥ 0.1. Any relationships that had a strength of < 0.1 were discarded to eliminate weak relationships that may have occurred due to highly turbulent winds, overly short-term winds, overly persistent winds, or other mechanisms that would have resulted in significant complexity. Mean R² values were calculated, divided by site (NM1 and NM2) and height within snowpack (0, 50, and 125 cm). Our data filtering technique was biased towards selecting periods of steady state and negative correlations between wind speed and CO₂ concentration. While the criteria seem demanding, in practice they were less restrictive than one might expect, and nearly one-fifth of all the measured data passed these filters and were included in the final analysis. With the use of our filtering process, our analysis does not represent an estimate of CO₂ flux during the snow-covered period.

We inspected the enhanced concentration profile experiment data (winter 2015) as a time series to analyze the effect of changing wind speed on CO₂ concentration at various levels within the snowpack. To quantify the effect of wind on CO₂ snowpack concentration, we identified the time periods when an abrupt increase in wind speed resulted in a rapid decrease in CO₂ concentration. These time periods were then used to determine the rate at which CO₂ decreased with an increase in wind speed. This was done in order to directly compare the field data with the modelled CO₂ data (see Sect. 2.4).

2.3 Model development and sensitivity testing

We developed a model to explore the control of three parameters on the CO₂ dynamics of a soil–snow system: soil diffusivity, snow diffusivity after initialization (advective wind intensity), and snow depth. The goal of this model was to use a diffusive model to mimic advective wind events through a snowpack. A previously existing multilayer 1-D soil diffusion model (Nickerson and Risk, 2009) was adapted for the soil–snow system. The exchange of CO₂ between layers was determined by Fick's first law, which assumes that gas transport through a diffusive medium is controlled by the concentration gradient, and occurs vertically. Fick's first law is given as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{F}_{{\mathrm{CO}}_{\mathrm{2}}}=-D_{{\mathrm{CO}}_{\mathrm{2}}}\left({\displaystyle \frac{\partial C_{{\mathrm{CO}}_{\mathrm{2}}}}{\partial z}}\right),\end{array}$$

where $F_{{\mathrm{CO}}_{\mathrm{2}}}$ is CO₂ flux (µmol m⁻² s⁻¹), $D_{{\mathrm{CO}}_{\mathrm{2}}}$ is CO₂ diffusivity within the snowpack (m⁻² s⁻¹), and $\frac{\partial C_{{\mathrm{CO}}_{\mathrm{2}}}}{\partial z}$ is the CO₂ concentration gradient of the snowpack (µmol m⁻³). The diffusivity of CO₂ within the snowpack can be calculated empirically using snowpack porosity (based on density), tortuosity, the diffusion coefficient of the specific gas under standard temperature and pressure, ambient pressure, and snowpack temperature (Seok et al., 2009). We tested a range of diffusivities (soil and snow), along with snow depth, but for simplicity we did not test ranges for individual parameters that are used to calculate diffusivity (e.g. snowpack porosity, tortuosity).

Table 1 A 1-D soil CO₂ diffusion model was adapted for the soil–snow system. The model simulated step changes in transport rate over a broad range of plausible synthetic cases. Soil diffusivity ranged logarithmically, whereas snow diffusivity and snow depth ranged linearly. We ran the model with all possible permutations of these parameters.

The purpose of the induced change of an increased snowpack CO₂ diffusivity was to mimic observed changes in CO₂ flux and snowpack concentration. Specifically, the induced increase in snowpack CO₂ diffusivity was used to simulate an advective wind event within a diffusive model. With Atlantic Computational Excellence Network (ACEnet) high-performance computers, we used model runs to explore the control of each of the three parameters on the CO₂ dynamics of the soil–snow system. The three parameters investigated were soil diffusivity (m² s⁻¹), snow diffusivity at step change (m² s⁻¹), and snow depth (cm). The tested range for each of the parameters is given in Table 1.

We initialized the model using a linear CO₂ concentration profile through the layers, determined by soil CO₂ diffusivity, layer height, and atmospheric CO₂ concentration (set at 380 ppm). Each model simulation began with the system in equilibrium state, which means storage flux was set to 1 µmol m² s⁻¹. We define storage flux here as the change in CO₂ storage in the snowpack, analogous to the exchange of CO₂ between the snowpack and the atmosphere. Varying numbers of snow layers were added on top of the 100 cm of modelled soil layers with the following distinctions: (1) we assumed that snow has a higher porosity than the underlying soil, so the snow layer diffusivities were always set to a value higher than the soil layers, and (2) we assumed that snow does not produce CO₂, so we removed CO₂ production from the snow layers.

Initial condition snow diffusivity was held constant at $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ for all simulations. Since snow diffusivity encompasses porosity, density, and tortuosity, these parameters also remained constant for initial conditions for all simulations: we assumed a homogeneous snowpack and did not test a range of snow diffusivities for initial conditions. To mimic a range of wind events, after initialization, we tested a range of snow diffusivities. Our 10 test values for this snow diffusivity, mimicking advective “wind events”, ranged linearly from $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{6}}$ (equal to the snow diffusivity at initial conditions) to $\mathrm{9.08}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ (approximately the diffusivity of CO₂ in air) (Table 1). We tested a plausible range of soil diffusivity and snow depth values (parameters used for initializing), though these remained unchanged through the “wind event” in each simulation. We tested a range of soil diffusivities to mimic a range of CO₂ emission rates out of the soil into the overlying snowpack. A range of snow depths was tested to mimic the natural environment that we tested in the field.

Figure 2 shows an example of the apparent storage flux and corresponding change in snowpack CO₂ concentration at every 10 cm, with an induced change in CO₂ snowpack diffusivity, which was the mechanism used to mimic an advective “wind event”. In summary, to simulate how the modelled diffusive system responds to an advective wind event, the model simulated induced changes in transport rate (snow diffusivity) within the snowpack over a range of plausible synthetic base cases (soil diffusivity and snow depth). We ran the model with all possible permutations of the three parameters given in Table 1.

Figure 2 An instantaneous change in snowpack diffusivity after model initialization mimics advection. Panel (a) shows the modelled storage flux with an induced change in snowpack CO₂ diffusivity. Panel (b) shows the corresponding change in snowpack CO₂ concentration at every 10 cm. Soil diffusivity is $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹, stepped snow diffusivity is $\mathrm{9.08}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹, and snow depth is 60 cm. The soil–atmosphere arrow indicates depths within the 60 cm snowpack: highest modelled CO₂ concentrations occur at the soil–snow interface, whereas lowest modelled CO₂ concentrations occur at the snow–atmosphere interface, before and after the advective “wind event”.

It is very likely that lateral CO₂ flux occurs within the snowpacks at our field sites, especially with the presence of wind slabs, sun crusts, and ice lenses at the sites. These features are unaccounted for in our modelling, as modelling lateral CO₂ transport through a snowpack in addition to vertical transport would require a 3-D model. Our overall objective with this model was to observe and understand the differences in diffusive and advective transport through snowpacks. As such, we refrained from overcomplicating the 1-D model (e.g. Fick's second law of diffusion).

For sensitivity analysis, we calculated fractional change. Each post-wind event CO₂ value was compared to a CO₂ value under the same conditions as if a wind event had not occurred:

$$\begin{array}{}\text{(2)}& {\displaystyle }\text{fractional change}=\left|{\displaystyle \frac{w-n}{n}}\right|,\end{array}$$

where w is a post-wind event and n is an event under no elevated wind conditions.

Figure 3 (a) Time series of wind speed and CO₂ concentration at 50 cm above the ground within the snowpack from 06:30 on 2 January 2014 to 14:00 on 3 January 2014 at NM1. Average snow depth at NM1 over this time period was 124 cm. (b) The corresponding linear regression of CO₂ concentration versus average wind speed (R²=0.70, P<0.001).

2.4 Field–model comparisons

In order to properly compare the field and modelled data, we determined the rate at which modelled CO₂ responded to the induced “wind events”. This refers to the change in CO₂ concentration over time (ppm s⁻¹) as a result of the change in snowpack CO₂ diffusivity after initialization. Of the modelled data, we considered only scenarios with a soil diffusivity of $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹. Additionally, only “low wind” and “high wind” events were considered, which had induced snow diffusivities of $\mathrm{1.72}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{9.08}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹, respectively. Output included CO₂ concentration at every 10 cm within the modelled environment (both soil and snow). For field–model comparison purposes, we only considered the CO₂ concentration of the topmost layer of snow.

We processed the enhanced concentration profile experiment data (winter 2015: 15 April–29 April) by calculating the rate of change of CO₂ concentration (ppm) per unit time (seconds) after a noticeable wind event.

3 Results

3.1 Snowpack CO₂ concentration profile experiment

Initial field campaigns (2014) (Graham, 2018) showed a relationship between wind speed and CO₂ concentration within the snowpack at NM1 and NM2. Wind speed sometimes had a very strong effect on CO₂ concentration within the snowpack (Figs. 3 and 4).

Trace amounts of snow at NM1 and NM2 began accumulating at the beginning of data collection (11 November 2013), with appreciable (> 25 cm) snowfall at both stations occurring on 15 December 2013 and remaining through the winter. Maximum snow depth at NM1 was 188 cm (26 March 2014), whereas maximum snow depth at NM2 was 137 cm (4 January 2014).

Figure 4 (a) Time series of wind speed and CO₂ concentration at 125 cm above the ground from 06:30 on 2 January 2014 to 14:00 on 3 January 2014 at NM1. Average snow depth at NM1 over this time period was 124 cm. Therefore, these CO₂ values were a good representation of the snow–air interface. (b) The corresponding linear regression of CO₂ concentration versus average wind speed (R²=0.53, P<0.001).

There was a negative correlation between average wind speed and CO₂ concentration 50 cm above the ground, an example of which can be seen in Fig. 3a. During this period of 31.5 h, snowpack CO₂ concentration at this height above soil ranged from 587 to 965 ppm. Wind speeds over this same time period ranged from 3.2 to 31.1 km h⁻¹. The corresponding linear regression (Fig. 3b) shows the effect that average wind speed exerted on CO₂ concentration (R²=0.70, P<0.001). As wind speed increased, CO₂ concentration decreased at a rate of 14.4 ppm km⁻¹ h.

Figure 4 shows measurements at NM1 over the same time period from 125 cm above ground. These CO₂ values were very close to predicted atmospheric concentrations, as the average snow depth over this time period at NM1 was 124 cm, very near the measurement height. The closeness of the measurement height to the snow surface indicates these values were likely a good representation of the CO₂ concentration at the snow–air interface. Despite increased atmospheric mixing, average wind speed exerted good control over CO₂ concentration (Fig. 4a). This result is reinforced with the corresponding linear regression (Fig. 4b; R²=0.53, P<0.001). As wind speed increased, CO₂ concentration decreased at a rate of 1.57 ppm km⁻¹ h.

Table 2 Summary of regression analysis between CO₂ concentration within the snowpack and wind speed. Data were filtered to satisfy the following conditions: (1) snow cover was considered to be at equilibrium, (2) the relationship produced a slope < 0, and (3) R² ≥ 0.1. N is the number of time periods that satisfy all three conditions. Each time period covered a minimum of 6 h. The y intercept is the mean CO₂ concentration when wind speed is  0 km h⁻¹. Slope is the mean change in CO₂ concentration with a 1 km h⁻¹ increase in wind speed. R² is the mean strength of the relationship between CO₂ concentration in the snowpack and mean wind speed. n is the mean number of half-hourly observations within each N. Duration is the mean duration of N. Instrumentation error for the NM1 0 cm CO₂ probe prevented data collection at that height. Instrumentation error for the NM1 0 cm CO2 probe prevented data collection at that height, and so data at that height is not available (NA).

We conducted a regression analysis of CO₂ concentration versus average wind speed for filtered data for winter 2014 (11 November 2013 to 26 March 2014), as per the three conditions specified in the Methods section. From this summary table (Table 2), there were some identifiable trends with the increasing height of CO₂ concentration measurement. With the increase from 50 to 125 cm at NM1 and 0 to 125 cm at NM2, there was a decrease in the y intercept, which was the mean predicted value of CO₂ concentration when average wind speed was 0 km h⁻¹. Additionally, the average slope of individual regressions became flatter with an increase in measurement height. Finally, the strength of the relationship (R²) decreased with an increase in measurement height (towards the open air). Instrumentation error for the NM1 0 cm CO₂ probe prevented data collection at that height.

The measurements that satisfied all conditions accounted for an average of 15.1 % of the data collected at a given station (NM1, NM2) and height in the snowpack (0, 50, 125 cm). This value does not represent an estimate of the CO₂ flux during the snow-covered period, since we used a biased filtering process to identify wind events during periods of steady snow cover.

Figure 5 Time series of enhanced profiler experiment (winter 2015) CO₂ concentrations throughout the snowpack and wind speed at NM2 over 2 weeks during late winter 2015 (15 April–29 April). Measurements were recorded hourly.

3.2 Enhanced concentration profile experiment

We collected CO₂ concentration profile data at the enhanced NM2 station from 16:00 on 15 April 2015 to 11:00 on 29 April 2015, a total of 331 uninterrupted hours (Fig. 5) (Graham, 2018). Average snow depth over this time period was 157 cm, ranging from 149 to 167 cm. Average air temperature was −1.4 ^∘C, ranging from −8.6 to 7.6 ^∘C.

Figure 5 shows a time series of CO₂ concentration throughout the snowpack (0, 25, 50, 75, and 100 cm from the ground), and mean wind speed. There was considerable variability in snowpack CO₂ concentration and wind speed over the 2-week period, with snowpack CO₂ values ranging from 357 to 4161 ppm and wind speeds ranging from 0.0 to 34.0 km h⁻¹. Average wind speed over the 2-week period was 13.5 km h⁻¹.

Average CO₂ concentration decreased with increasing proximity to the atmosphere: 1244, 1076, 1007, 886, and 867 ppm at 0, 25, 50, 75, and 100 cm, respectively. Average atmospheric CO₂ concentration over this sampling period was relatively constant at 512 ppm. For some time periods between 15 April and 29 April 2015, there was a slight negative correlation between wind speed and snowpack CO₂ concentration (Fig. 5), but this was not tested using the methodology of testing the winter 2014 data.

Figure 6 (a) Modelled results at top of snowpack shown as the fraction of CO₂ concentration depleted from the snowpack. (b) Modelled results at bottom of snowpack shown as the fraction of CO₂ concentration depleted from the snowpack. (c) Modelled storage flux, shown as factor increase in short-term CO₂ flux. Scenarios at equilibrium (8 days) were incalculable (not shown), as there was no change in CO₂ concentration once equilibrium was reached. (d) Modelled CO₂ at the topmost soil layer, shown as the fraction of CO₂ concentration depleted from the snowpack. There was very minimal effect on the fraction of CO₂ depletion for immediate scenarios (10 min), and so there is significant overlap of the two 10 min scenarios (light wind and severe wind).

3.3 Modelling

Figure 6 shows results from sensitivity testing of an enhanced diffusion model used to simulate advection and the effect of several parameters as deviations from a base case (Table 1). Model activity was investigated at the following layers: the topmost snow layer (CO₂ concentration in Fig. 6a and storage flux out of the top of the layer in Fig. 6c), the bottommost snow layer (CO₂ concentration in Fig. 6b), and the topmost soil layer (CO₂ concentration in Fig. 6d).

Results are shown as fractional depletion of CO₂ concentration in the snowpack (Fig. 6a, b, d), and factor increase in short-term CO₂ storage flux (Fig. 6c). Of the three parameters (soil diffusivity, snow diffusivity mimicking advection, and snow depth), soil diffusivity had negligible control on layers involving snow (Fig. 6a, b, and c) and is therefore not represented in those panels. Soil diffusivity showed some control on the modelled soil layer (Fig. 6d).

We also considered time when analyzing the modelled data to investigate how CO₂ concentration is affected during the “wind event” recovery period as the system works its way towards equilibrium (immediate change) and once the modelled system had recovered to an equilibrium state. Equilibrium refers to no change in the modelled storage flux or when storage flux had returned to the initialized condition of 1 µmol m⁻² s⁻¹. The two time “scenarios” considered were (1) at 10 min and (2) at 8 days following the simulated “wind event”. The 10 min scenario represented “immediately following a wind event” and the 8-day scenario represented “once equilibrium had been reached”.

In the modelled topmost layer of snow (Fig. 6a), the maximum fraction to which CO₂ concentration was depleted was 0.39 once equilibrium was reached after a severe wind event. Snow depth had no effect on CO₂ depletion for either equilibrium scenarios at the top of the snowpack. For scenarios immediately following a wind event, severe winds had a greater effect on the fraction of CO₂ depleted, but this effect decreased with increasing snow depth (approaching no CO₂ depletion).

CO₂ concentration at the bottommost layer of snow (Fig. 6b) behaved similarly to the CO₂ concentration in the topmost layer. Depletions at the bottom of the snowpack were up to 2 times that of the depletions at the top of the snowpack (maximum fraction of 0.81 once equilibrium was reached after a severe wind event, with 100 cm of snow). Scenarios that immediately followed a wind event showed that severe winds had a greater effect on CO₂ depletion, although this decreased with increasing snow depth, reaching a minimum fraction of 0.06 at 100 cm.

Storage flux from the top of the snowpack into the modelled atmosphere is shown as a factor increase in short-term CO₂ flux (Fig. 6c). Scenarios at equilibrium (at 8 days post-event) are not shown, as there was no change in CO₂ concentration once equilibrium was reached. Of the scenarios that immediately follow a wind event, light and severe winds had similar effects on factor increase with 20 cm of snow: a factor of 0.53 (light wind) and a factor of 0.25 (severe wind). With increasing snow depth, severe winds showed a much greater fractional increase (9.92) in storage flux than light winds (1.15).

At the topmost soil layer (Fig. 6d), CO₂ concentration was affected by soil diffusivity and unaffected by snow depth. With increasing soil diffusivity at equilibrium, a greater fraction of CO₂ was depleted from the soil layer. Severe winds depleted a greater fraction than light winds. There was essentially no effect on the fraction of CO₂ depletion immediately following wind events (at 10 min post-event) of any severity, and therefore there is significant overlap of the two 10 min lines in Fig. 6d.

4 Discussion

4.1 Wind causes short-lived advective anomalies

Findings of the initial snowpack CO₂ concentration profile experiment showed that there was a negative correlation between wind (advective) events and the CO₂ concentration in a snowpack, on a timescale of hours to days. This was clear from specific examples (Figs. 3 and 4) as well as from the overall summary of linear regressions performed between CO₂ snowpack concentration and wind speed (Table 2). However, this was not continuous over the entire winter and was only true under particular conditions where filtering criteria were satisfied. The balance of the datasets that did not meet criteria was simply noisy with visible but weak trends. These time periods that did not meet the criteria may have resulted from the presence of vertical density variations (wind slabs, ice lenses) within the snowpacks at our field sites, plausibly causing lateral CO₂ flux. In addition to finding a negative correlation between wind events and CO₂ concentration within the snowpacks, analysis of data from the first experiment showed that there was a CO₂ concentration gradient throughout the snowpack, with highest concentrations closest to the soil and lowest concentrations closest to the atmosphere. This was consistent with previous literature, which indicates that within the porous medium locations closer to the source of production of the trace gas (e.g. CO₂) correspond with higher concentrations of that trace gas (Seok et al., 2009).

This work reinforced earlier observations of depleted CO₂ concentrations in field datasets (Seok et al., 2009), although we did not measure or calculate CO₂ storage flux directly in the field at the snow surface. However, we inferred that sporadic changes in snow–atmospheric flux would have been present from the large decreases in concentration. Positive storage fluxes were balanced by negative storage fluxes following wind events. It is important to consider concentration gradients to help with our understanding of the underlying physical processes of CO₂ transport through snowpacks.

As the measurements taken at each of the snowpack heights at each of the stations satisfied all specific conditions for an average of 15.1 % of the time analyzed, we can conclude that advection showed some control over snow CO₂ transport for this location for the equivalent of 20.4 days during the 135-day period in 2014 (12 November 2013 to 26 March 2014). This value did not represent the percentage of annual flux during the snow-covered season (Liptzin et al., 2009), though it did confirm that advective transport needed to be taken into account when studying snowpack CO₂ transport. It also gives an indication of how much data were eliminated for analysis, biasing our results.

The enhanced concentration profile experimental data reinforced the results of the initial findings and added CO₂ concentration measurements throughout the snowpack, increasing the total in-snow measurements from three to five. This gave us a clearer indication of how the CO₂ concentration gradient behaved, even without taking snow properties into account. These data covered the late winter period, so ice layers within the snowpack were likely present. Despite this, the wind seemed to have an effect on CO₂ snowpack concentrations, even at 0 cm with a snowpack of 157 cm.

Table 3 Summary table of change in modelled CO₂ concentration per second at 1, 2, 4, 6, and 24 h since the wind event (step change in modelled snowpack diffusivity) at the topmost layer in the model. Snow depths of 20, 60, and 100 cm are shown, along with lowest and highest simulated wind speeds.

Some authors have used turbulent atmospheric pressure pumping to explain anomalous CO₂ storage fluxes but have often focused this work on shorter, high-frequency timescales of seconds to minutes (Massman et al., 1995). On the longer, low-frequency range of the timescale, Bowling and Massman (2011) and Massman et al. (1995) mentioned the importance of synoptic-scale changes in atmospheric pressure. Additionally, Rains et al. (2016) showed that changes in wind speed at multiple hour frequencies (greater than 10 h) could be more effective than atmospheric pressure pumping when explaining changes in snowpack CO₂ concentration. These processes of different timescales and different mechanisms (atmospheric pressure changes, wind) affect CO₂ concentration gradients and fluxes measured with Fick's law by ventilating diffusive media, like snowpacks. The ventilation, no matter the timescale, affects the CO₂ concentration gradient by mixing atmospheric air into diffusive media where CO₂ typically pools, thereby affecting the CO₂ flux from the top of the snowpack. Our work showed how persistent wind and an enhanced diffusive profile controlled CO₂ concentration and fluxes across timescales of hours to days, in the mid-range between very high-frequency pressure pumping and low-frequency barometric pressure effects. The low-frequency synoptic processes occur on a longer timescale than the wind depletion events discussed in this study, though would be present here as well and would likely contribute to some of the variability (Robinson and Sextro, 1997; Tsang and Narasimhan, 1992). With more longer continuous wintertime CO₂ records, similar to this one, it may be possible to extricate these synoptic process periodicities in addition to the mid-range frequencies we investigated.

Table 4 Summary table of change in actual CO₂ concentration per second for four events in April 2015 when a decrease in CO₂ concentration corresponded to an increase in wind speed. CO₂ concentration was measured in the snowpack at 100 cm from the ground. Rate of change of CO₂ concentration, snow depth, start time, end time, range of CO₂, and range of wind speed are given in the table.

Figure 7 Conceptual diagram of diffusive versus storage flux. (a) Diffusive flux through a snowpack, with CO₂ originating from soils and consistently passing through a diffusive medium into the atmosphere as a result of a concentration gradient. Small arrows indicate low levels of diffusive flux that are prevalent and constant through time. (b) Storage flux through a snowpack, with CO₂ originating from soils, pooling in a diffusive medium, and then released to the atmosphere at a higher rate (than diffusive flux) following a high wind event, which has ventilated the top of the diffusive medium and steepened the concentration gradient. The larger arrow indicates the higher rate and lower frequency of storage flux out of snowpacks when compared with diffusive flux.

4.2 A diffusive model can help explain advective questions

The 1-D diffusional transport model and enhanced diffusion approach was able to replicate the CO₂ depletions seen in the field in this experiment, as well as those in previous observations (Seok et al., 2009) and in other plausible situations. Advective events were created with induced increases in snowpack diffusivity after model initialization, which worked well to mimic wind events.

In general, when snowpack diffusivity was instantaneously increased in this diffusive transport model, we observed rapid changes in the snowpack CO₂ concentration, CO₂ storage flux, and soil CO₂ concentration, similar findings to Bowling and Massman (2011). This effectively simulated advective events observed in the field. According to this model, severe wind events always produced more dramatic results than light wind events in terms of both rate of change (flux) and overall concentration change.

This modelling work showed that we can simplify the impacts of sustained advection on CO₂ in a soil–snow system to an effective diffusion problem. This approach was less complicated than other models that use the diffusive–advective coupled solution approach.

4.3 Field–model comparisons

To determine the applicability of the model to real-world scenarios, we compared our field and model results. To do so, we calculated the rate of change of CO₂ concentration (ppm) per unit time (seconds) after a wind event for both the modelled wind events and the field wind events (using the 2015 enhanced experiment). Figure 5, which displays a time series of CO₂ concentrations and wind speed over 2 weeks in late April 2015, shows that despite similarly variable wind conditions, snowpack CO₂ concentrations throughout the first week vary less than the CO₂ concentrations observed in the second week. The lack of variation in the first week could be due to a variety of reasons, including the composition of the snowpack or other meteorological conditions like temperature or humidity. Despite the variation through the 2-week period, it was still possible to discern change in CO₂ concentration after a wind event (Table 4).

Table 3 summarizes the calculated rates of change of modelled CO₂ concentration at varying snow depths, at low and high simulated wind speeds (induced change in snow diffusivity), and at various times since the modelled wind event. All of these modelled measurements were taken from the topmost snow layer. Table 4 shows a similar summary for four wind events in the field in April 2015. All of these CO₂ field measurements were taken at 100 cm from the ground within the snowpack, which was the in-snow measurement farthest from the ground and closest to the atmosphere at the time.

Change in modelled CO₂ concentration per second (Table 3) did not align perfectly with the change in field CO₂ concentration per second (Table 4) after a wind event. However, the rates of change in the field events ($-\mathrm{0.07},-\mathrm{0.04},-\mathrm{0.20},-\mathrm{0.04}$ ppm s⁻¹) were of approximately the same order of magnitude as the rates of change in the modelled events (ranging from 0.00 to −2.08 ppm s⁻¹). This indicated that the model was able to mimic advective events with some accuracy. Though it may be possible, as in other studies (Latimer and Risk, 2016), to apply an iterative procedure to our model with the conditions we observed in the field (e.g. initial CO₂ concentration), we deemed that to be unnecessary. This is because our primary goal was to properly illustrate the underlying physics of CO₂ transport through snowpacks. As such, matching the model conditions exactly to the field conditions was not required.

This study showed the importance of continuous monitoring of CO₂ concentrations and fluxes from soils through snowpacks. Similarly, Webb et al. (2016) and Rains et al. (2016) highlighted the non-growing season contributions to annual CO₂ flux. Webb et al. (2016) showed that different wintertime measurement methods at one Alaskan site resulted in a 4-fold range in CO₂ loss. The EC method showed the highest fluxes, as more CO₂ was released under windy conditions and the EC method was able to measure fluxes in turbulent conditions (Webb et al., 2016). Rains et al. (2016) noted that there are benefits and disadvantages to the EC, flux gradient, and chamber methodologies for non-growing season soil CO₂ flux measurements and that accurate parameterization of advective transport through snowpacks is important, regardless of the methodology. This is particularly true because of the likelihood that CO₂ flux through snowpacks is often underestimated (Rains et al., 2016). Accompanying these findings, we agreed that infrequent measurement can lead to significant error in the annual C budget of various ecosystems once inaccurate values are scaled up (Fig. 7). The effects of advection on these soil–snow systems can lead to variability in storage flux, as effective diffusion is closely related to wind. Snowpack depth, density, and layering will also affect the timing and amounts of CO₂ storage flux from these systems. We recommend that future studies utilize continuous CO₂ monitoring methods and consider the advective effects of wind in order to capture the uncertainties of soil CO₂ emissions in snow-covered ecosystems.

5 Conclusions

Although this study was conducted at one site over two winters, the findings have implications for measuring wintertime CO₂ fluxes in snow-covered environments. This is important for continued careful assessment of the soil C pool and fluxes of these snow-covered regions, which are experiencing increasing temperatures and variations in precipitation patterns.

As seen from the fieldwork in winters 2014 and 2015, advective transport by wind is important for CO₂ concentration (and therefore flux) through a soil–snow profile. Additionally, this process can be simulated with some accuracy by a model of enhanced diffusion. In both field and model cases, we observed how sustained winds could deplete CO₂ concentration in the snowpack and create storage flux outward to the atmosphere. During the re-equilibration phase, fluxes across the snow–air interface would have been depressed, as most of the production contributed initially to pore space storage. This process of buildup and release occurs with regularity in snow profiles and is likely more severe in snowpacks than in soil, which has lower permeability and is therefore less vulnerable to wind invasion.

Transport lags are the main effect of diffusion and advection. Measurements such as EC, which can be made above the snow profile with speed, are at an advantage for detecting storage flux events. While useful for total accounting purposes, EC records may not be effective in determining specific overwinter biological soil CO₂ production. For this, sensors within or at the base of the snowpack would also be needed, allowing the results to quantify soil–snow fluxes or concentration gradients within the first few centimetres of snow. Additionally, in situ sensors are typically cheaper and can be more easily and frequently deployed than EC methods. Alternatively, the model used here, which accurately simulated gas transport physics, could be applied through an inversion scheme to determine microbial changes in CO₂ production by removing the effects of snow gas transport.

This study shows snow profile CO₂ depletions that exist on timescales of hours to days. Putting this knowledge into practice would help to improve our understanding of global winter soil CO₂ release because it improves our efforts to quantify winter fluxes. As a start, we recommend that researchers approach winter data like they do summer data, which means that they should consider using continuous automated approaches for wintertime CO₂ flux observations, as done in this study. We also recommend close collaboration between the modelling community and soil field scientists. This will ensure that available physical models are being effectively used for stripping flux data of transport-related artefacts, thereby isolating soil biological behaviour.