Snow thermal conductivity controls future winter carbon emissions in shrub–tundra

The Arctic winter is disproportionately vulnerable to climate warming and approximately 1700 Gt of carbon stored in high-latitude permafrost ecosystems is at risk of degradation in the future due to enhanced microbial activity. Few studies have been directed at high-latitude cold season land–atmosphere processes and it is suggested that the contribution of winter season greenhouse gas (GHG) fluxes to the annual carbon budget may have been underestimated. Snow, acting as a thermal blanket, influences the Arctic soil temperatures during winter and parameters such as snow effective thermal conductivity (K_eff) are not well constrained in land surface models, which impacts our ability to accurately simulate wintertime soil carbon emissions. To address this, we investigated the impacts of implementing a K_eff parameterisation more suitable to Arctic snowpacks into the Community Land Model (CLM5.0). A point-model version of CLM5.0 forced by an ensemble of NA-CORDEX (North American Coordinated Regional Downscaling Experiment) future climate realisations (RCP 4.5 and 8.5) indicates that median winter CO₂ emissions will have more than tripled by the end of the century (2066–2096) under RCP 8.5. Implementing the refined K_eff parameterisation increases simulated winter CO₂ in the latter half of the century (2066–2096) by 130 % and CH₄ flux by 50 % under RCP 8.5 compared to the widely used default K_eff parameterisation. The influence of snow K_eff parameterisation within CLM5.0 on future simulated CO₂ and CH₄ is at least as significant, if not more so, than climate variability from a range of NA-CORDEX projections to 2100. The average difference in refined K_eff compared with the default K_eff raises minimum winter soil temperatures by 4–7 °C by the end of the century under RCP 4.5 and 8.5. Furthermore, CLM5.0 simulations using the refined K_eff show an extension of the early winter (September–October) zero-curtain period, by nearly a month under RCP 8.5. Consequently, recent increases in both zero-curtain duration and winter CO₂ emissions appear set to continue to 2100. Modelled winter soil temperatures and carbon emissions further highlight the importance of climate mitigation in preventing a significant increase in winter carbon emissions from the Arctic in the future.

1 Introduction

It is estimated that ∼ 1700 Gt of carbon is stored in the permafrost ecosystems of northern latitudes (Miner et al., 2022), which accounts for half of the global soil carbon storage (Hugelius et al., 2014) and is at risk of degradation due to climate change (Koven et al., 2013). The near-surface permafrost temperature has been increasing 1.1 °C decade⁻¹ since 1987 (AMAP, 2017), which acts to reduce the near-surface permafrost extent and increase active layer thickness (Intergovernmental Panel on Climate, 2023); such loss of permafrost drives the further release of carbon in the form of CO₂ during winter. Arctic winters are estimated to warm by 4.8 °C by 2100 compared to 2.2 °C in the summer over the same period (Webb et al., 2016; Christensen et al., 2013). Climate warming has enhanced the warming of the soil at Arctic sites in the recent past (AMAP, 2017; Ednie and Smith, 2015) and future warming will enhance the microbial decomposition of soil organic matter, driving the subsequent release of carbon dioxide (CO₂) and methane (CH₄) (Natali et al., 2019; Schuur et al., 2015). Furthermore, the winter carbon flux is crucial to the annual carbon budget, as demonstrated through flux measurements from a tundra site, where approximately 60 % of CO₂ uptake in the growing season was then lost in the non-growing period, and CH₄ emissions during this period accounted for 30 % of the annual budget (Kittler et al., 2017). Arctic tundra ecosystems, once considered to be carbon sinks, are increasingly acting as carbon sources due to elevated winter carbon emissions driven by rising temperatures and deeper snow cover (Pongracz et al., 2021; Fahnestock et al., 1999; Welker et al., 2000; Schimel et al., 2004). Future winter emissions will likely offset growing season uptake under Representative Concentration Pathways (RCP) 4.5 (medium global emission scenario) and 8.5 (high global emission scenario) (Natali et al., 2019). Additionally, the Arctic winter is longer than the growing season (the latter being typically 2–5 months in length), and therefore may represent a substantial input of carbon to the atmosphere that is not accurately represented in current climate models nor quantified accurately in global carbon budgets (Treat et al., 2024). Assessing simulations of heat and gas fluxes to and from subnivean soils is therefore critical in understanding how well models, such as the Community Land Model (CLM5.0; Lawrence et al., 2019), can be expected to simulate future carbon cycling in the Arctic.

Arctic snow is a key determinant of ground temperature and plays a major role in the wider hydrological and ecological Arctic system (Callaghan et al., 2011). The hydrological cycle in the Arctic is projected to intensify throughout the 21st century and new model estimates from the Coupled Model Intercomparison Project 6 (CMIP6) support a future transition from a snow- to rain-dominated Arctic in the summer and autumn by 2100 (McCrystall et al., 2021). Despite future temperatures still being cold enough to enhance snowfall, winter rainfall also intensifies, with projected increases even matching that of snow (Bintanja and Andry, 2017). Snow is a key determinant of ground temperature and the freeze–thaw state through its insulating properties (Bigalke and Walsh, 2022). The prospect of a rain-dominated Arctic has implications for soil temperature, through changes in the snow cover duration (Mudryk et al., 2021; Mudryk et al., 2020), insulative and structural snow properties (Cohen et al., 2015; Serreze et al., 2021), and increased soil moisture (Trenberth, 1998; Cohen et al., 2015). In the recent past (1979–2018) snowpack properties have featured significant increases in the spring snow bulk density (May and June), a downward trend in the snow cover duration, and an upward trend in wet snow, all of which have resulted in a soil surface temperature increase of +0.41 K decade⁻¹ in winter (Royer et al., 2021). Such temperature increases, induced by snowpack evolution may contribute to increased soil heterotrophic respiration and carbon release (CO₂ and CH₄) in future winters (Natali et al., 2019).

The “shoulder” seasons, defined here as the early (after the summer plant senescence has ended) and late (before the snowmelt is complete) winter period (Natali et al., 2019), are a critical period in the context of the annual Arctic carbon budget, especially in terms of CH₄, where zero-curtain (when soil temperatures are poised near 0 °C, such as in early winter) emissions alone contribute around ∼ 20 % of the annual budget (Zona et al., 2016). In terms of CO₂, historic early winter season emission rates have increased 73 % ± 11 % since 1975 (October–December) (Commane et al., 2017) in response to a warming climate, demonstrating the vulnerability of this period to future warming. To date Earth system models (ESMs) have poorly simulated respiration in the early winter shoulder season (Commane et al., 2017), which is a key component of the annual Arctic carbon budget (Kittler et al., 2017). This study addresses this knowledge gap by examining how projected changes in the shoulder season air temperature, precipitation, snow thermal conductivity, and soil temperature influence the cold season carbon dynamics (Figs. 2–5).

Land surface models (LSMs) such as CLM5.0 struggle to accurately reproduce cold season CH₄ and CO₂ emissions due to uncertain numerical representation of mechanisms controlling heterotrophic respiration (Zona et al., 2016; Commane et al., 2017; Tao et al., 2021; Natali et al., 2019). One such mechanism is the mediation of winter soil temperatures by the insulating effect of seasonal snow. The default parameterisation of snow thermal conductivity (K_eff, $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$; Jordan, 1991) in CLM5.0 produces winter soil temperatures that are colder than observations (Dutch et al., 2024), which may contribute to an underestimation of the soil respiration (SR) in the early winter (Commane et al., 2017). Sturm et al. (1997) presented an alternative parameterisation for K_eff derived from extensive Arctic snow measurements (n = 488) that was implemented into CLM5.0 by Dutch et al. (2024). The Sturm et al. (1997) parameterisation is more relevant for Arctic snow types compared to Jordan (1991), as it accounts for snow layers such as basal depth hoar and wind slab within the snowpack and produced more accurate simulations of subnivean soil temperature (Dutch et al., 2024). Further, the CLM5.0 default soil moisture threshold for decomposition (Ψ_min = −2 MPa) is too high to permit sub-zero SR and this has been identified as a limitation in winter simulations (Tao et al., 2021; Dutch et al., 2024). Similarly, the CLM5.0 default settings of Q₁₀ (1.5) and Q₁₀ch4 (1.3), which dictate respiratory responses to changes in temperature, are too low for Arctic tundra environments (Müller et al., 2015; Dutch et al., 2024). These parameters Ψ_min, Q₁₀, and Q₁₀ch4, alongside K_eff, require adjustment to realistically simulate soil temperature (TSOI), SR, and methane flux (FCH4) under cold season conditions. Implementing the Sturm et al. (1997) parameterisation reduces cold bias in simulated soil temperatures by two-thirds. Further, while the default CLM5.0 produces negligible winter net ecosystem exchange (NEE), combining the Sturm et al. (1997) parameterisation with a mid-range value for Ψ_min (−20 MPa) produces winter NEE values consistent with field observations (Dutch et al., 2024; Fig. 5). Our study applies Q₁₀ and Ψ_min values suitable for tundra soil as defined by Dutch et al. (2024) to investigate the influence of different K_eff parameterisations on future projections of soil temperature and carbon emissions.

This study assesses the sensitivity of SR (i.e. soil CO₂ flux to the atmosphere during winter) and CH₄ flux to K_eff, Ψ_min, Q₁₀, and Q₁₀ch4 within CLM5.0, and presents future projections of carbon emissions for Trail Valley Creek (TVC) to 2100. Firstly, future changes in Arctic meteorology at TVC are analysed using an ensemble (n = 33) of future NA-CORDEX (North American Coordinated Regional Downscaling Experiment) projections (McGinnis and Mearns, 2021), representing varying global climate model–regional climate model (GCM–RCM) combinations; RCP scenarios (4.5; n = 6 and 8.5; n = 27), and RCM resolutions (0.22°; n = 11, and 0.44°; n = 22). The NA-CORDEX ensemble is then applied to a point model of the default version CLM5.0 (“CORDEX-Jordan”) to simulate present (2016–2046) and future (2066–2096) snow water equivalent (SWE), 10 cm ground temperature (GT10), and soil moisture (SOILLIQ) and resulting winter SR and CH₄ fluxes. The ensemble of experiments is then re-run, with the representation of K_eff changed from Jordan (1991) to Sturm et al. (1997) to assess the impact on the soil temperature; this ensemble is named “CORDEX-Sturm”. Finally, a plausible parameter space of Ψ_min, Q₁₀, and Q₁₀ch4 for Arctic environments is applied to each experiment grouping (CORDEX-Jordan and CORDEX-Sturm) to assess the effect of assumptions governing the moisture- and temperature-related controls on SR, resulting in an overall set of 396 experiments to quantify the effects of parametric and future climate uncertainty on carbon emissions at TVC.

2 Methodology

2.1 Study location

Point scale CLM5.0 simulations were produced for TVC (68°45^′ N, 133°30^′ W), a research basin on the boundary of the boreal-tundra ecotone ∼ 50 km north of Inuvik in the Inuvialuit Settlement Region of western Inuit Nunangat in the Western Canadian Arctic. The basin lies within the Inuvialuit Settlement Region near the Inuvik-Tuktoyaktuk Highway. Vegetation consists of shrub–tundra and is dominated by mineral earth hummocks. The mean annual air temperature from 1999–2018 averaged −7.9 °C (Grünberg et al., 2019), with average seasonal variations of +10 °C in summer and −25 °C in winter (1995) (Marsh et al., 2002). Air temperatures increased 1.1 °C decade⁻¹ during the 1990–2018 period, with the strongest warming observed in May (2.8 °C decade⁻¹) (Grünberg et al., 2019). The snow depth at TVC is typically < 50 cm, with deeper packs associated with variations in topography and vegetation (King et al., 2018). The snow is often redistributed by wind, resulting in scour and snow drift features (Thompson et al., 2016). Snow cover typically lasts for 8 months (October–May) and the underlying permafrost is between 350–500 m in thickness, with an active layer of 0.5–1.0 m (Wilcox et al., 2019; Dutch et al., 2024). Tundra environments occupy ∼ 31 % of the Canadian Arctic (Bliss and Matveyeva, 1992; Quinton et al., 2000) and well-drained, shrub- and lichen-covered uplands, such as those at TVC, comprise 80 % of the Arctic-boreal region (Voigt et al., 2023). TVC has long lasting concurrent hydrometeorological and eddy covariance (EC) carbon flux data since 2013 (Tutton et al., 2024), which can be used to bias-correct future meteorological simulations (see Sect. 2.2).

2.2 Meteorological forcing data and simulations using CLM5.0

The NA-CORDEX meteorological forcing data (incoming shortwave and longwave radiation, precipitation, humidity, wind speed, air pressure, and min and max temperature) for the grid cell closest to TVC provide CLM5.0 with daily meteorological conditions to 2100. This forcing dataset was generated from a suite of 33 bias-corrected GCM–RCM combinations for North America (McGinnis and Mearns, 2021) run under full transient conditions with a historical period spanning 1950–2005 and with scenarios RCPs 4.5 and 8.5 for 2006–2100 at 0.22 and 0.44° resolution. To minimise the influence of climate model bias on model results, the meteorological forcing data from the climate models were bias-adjusted to match the statistical characteristics of the meteorological observations (Tutton et al., 2024) in the overlapping 1992–2022 period following Cannon (2018) and Cannon et al. (2022). The CLM5.0 was run in “point mode” (a 0.1° × 0.1° grid cell) and centred at the location of TVC as per Dutch et al. (2021). A model spinup is required to stabilise the carbon pools within CLM5.0 prior to model initiation; model simulations were spun up for 512 years and full details of these spinup conditions are outlined in Dutch et al. (2024).

2.3 The CLM5.0 parameter assessment

K_eff describes the rate of heat transfer from the atmosphere through the snow to underlying soils and as such K_eff influences the heterotrophic respiration rates and subsequent production of CO₂ and CH₄. The CLM5.0 default K_eff parameterisation using Jordan (1991) overestimates K_eff by a factor of 3 for tundra snow (Dutch et al., 2021). In the CLM5.0, K_eff is parameterised as a function of snow density (ρ, kg m⁻³) (Eq. 1), which is calculated from relative proportions of ice mass (m_i) and liquid water mass (m_lw), weighted by the snow cover fraction (F_sno) multiplied by the height of the snow layer (h_sl) for each grid cell. K_eff can then be calculated using coefficients in Eq. (2) (Jordan, 1991), which utilises the thermal conductivity of air and ice (K_ice−K_air), or Eq. (3) (Sturm et al., 1997):

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{\rho }={\displaystyle \frac{m_{\mathrm{i}}\times m_{\text{lw}}}{F_{\text{sno}}\times h_{\text{sl}}}}\text{(2)}& {\displaystyle }\begin{array}{rl}{K}_{\text{eff}}=& K_{\text{air}}+\left(\right((\mathrm{7.75}\times {\mathrm{10}}^{-\mathrm{5}}\times \mathit{\rho })\\ & +(\mathrm{1.105}\times {\mathrm{10}}^{-\mathrm{6}}\times {\mathit{\rho }}^{\mathrm{2}})\left)\right(K_{\text{ice}}-K_{\text{air}}\left)\right)\end{array}\text{(3)}& {\displaystyle }\left\{\begin{array}{l}{K}_{\text{eff}}=\mathrm{0.138}-\mathrm{1.01}\mathit{\rho }+\mathrm{3.233}{\mathit{\rho }}^{\mathrm{2}},\mathrm{0.156}\le \mathit{\rho }\le \mathrm{0.6}\\ K_{\text{eff}}=\mathrm{0.023}+\mathrm{0.234}\mathit{\rho },\mathit{\rho }<\mathrm{0.156}\end{array}\right.\end{array}$$

The use of the snow thermal conductivity parameterisation following Sturm et al. (1997) improved soil temperature simulations in CLM5.0 at TVC (Dutch et al., 2021) and at other sites with additional LSMs (Royer et al., 2021). Here, the Sturm et al. (1997) K_eff parameterisation was applied to future CLM5.0 simulations and compared with the default parameterisation which follows Jordan (1991). CLM5.0 uses Ψ_min to describe the soil moisture threshold for respiration (MPa). Without adjustments to Ψ_min, CLM5.0 simulates near-zero SR for the majority of the snow-covered, non-growing season, which is not consistent with field observations (Dutch et al., 2024). Reducing Ψ_min facilitates respiration in moisture-limited soils (Tao et al., 2021). Ψ_min is applied in CLM5.0, as described by Andren and Paustian (1987):

$$\begin{array}{}\text{(4)}& r_{\mathrm{W}}={\sum }_{j=\mathrm{1}}^{\mathrm{5}}\left\{\begin{array}{cc}\mathrm{0}& \text{for }{\mathrm{\Psi }}_j<{\mathrm{\Psi }}_{\text{min}}\\ \frac{\log ({\mathrm{\Psi }}_{\text{min}}/{\mathrm{\Psi }}_j)}{\log ({\mathrm{\Psi }}_{\text{min}}/{\mathrm{\Psi }}_{\text{max}})}{w}_{\text{soil},j}& \text{for }{\mathrm{\Psi }}_{\text{min}}<{\mathrm{\Psi }}_j<{\mathrm{\Psi }}_{\text{max}}\\ \mathrm{1}& \text{for }{\mathrm{\Psi }}_j>{\mathrm{\Psi }}_{\text{max}}\end{array}\right\}\end{array}$$

where the rate scalar for soil water potential (r_W) is calculated from Ψ_j that is the soil water potential in soil layer j, and ${\mathrm{\Psi }}_{\text{min}}/{\mathrm{\Psi }}_{\text{max}}$ as upper and lower limits. The default values of ${\mathrm{\Psi }}_{\text{min}}/{\mathrm{\Psi }}_{\text{max}}$ for CLM5.0 are −2 and −0.002 MPa, respectively.

As well as K_eff and Ψ_min, we also focus on the temperature soil decomposition modifier r_T, which is parameterised using a Q₁₀ function which describes the dependence of biological metabolic processes, such as respiration, on temperature. r_T is calculated as follows:

$$\begin{array}{}\text{(5)}& r_{\mathrm{T}}=Q_{\mathrm{10}}^{\left(\frac{T_j-T_{\text{ref}}}{\mathrm{10}}\right)}\end{array}$$

where T_j is the temperature of soil layer j and T_ref is a reference temperature with a default value of 25 °C. CLM5.0 uses a fixed Q₁₀ of 1.5; however, this is unrealistic for Arctic ecosystems, which typically have an average Q₁₀ of 5.4 (Chen et al., 2020). A range of Ψ_min (“minpsi_ hr”, −2 to −20 MPa) and Q₁₀ (“q10_ hr”, “froz_ q10”, 1.5 to 7.5) was used to assess the sensitivity of CLM5.0 for future simulations to 2100. A parameter “q10ch4” within the methane module of CLM5.0 which controls the temperature dependency of CH₄ production is available and was tested from its default value 1.33 to a plausible extreme of 4 (Müller et al., 2015; Riley et al., 2011) with a reference temperature of 22 °C. We implemented a broad range of Q₁₀ and q10ch4 (1.5–7.5) to capture variability in the temperature sensitivity of SR associated with differences in carbon pool lability (Fierer et al., 2005; Yan et al., 2017). The full range of tested parameters is presented in Fig. A1.

A major focus of this study is the winter season when snow is on the ground (i.e. the snow-covered non-growing season). We define this period as the time when all model ensemble members (RCP 4.5 n = 6, RCP 8.5 n = 27) agree that the SWE is > 5 mm. Simulations of SR and FCH4 are filtered by these constraints so the analysis in Figs. 5 and 6 is focused only on carbon fluxes across a common snow-covered, non-growing season in all scenarios and forcing datasets. The resulting data are therefore weighted towards a shorter snow cover period as ensemble members capture a range of snow cover durations. Periods of the winter season that are critical to driving the timing and rates of respiration are described as follows. The “zero-curtain period” is defined as days where the soil temperature falls between −0.75 and 0.75 °C (Tao et al., 2021; Zona et al., 2016) while “shoulder seasons” describes annual transitional periods of snowmelt (April–May) and plant senescence to freeze up (September–October) (Shogren et al., 2020). Testing for statistically significant differences between distributions of future meteorological conditions and median winter carbon fluxes was evaluated using the non-parametric Kolmogorov–Smirnov (K-S) test.

Figure 1 Simulated median monthly precipitation and temperature over two 30-year time periods: 2016–2046 (black) and 2066–2096 (red) under RCP 4.5 and RCP 8.5, for TVC, using meteorological data from an ensemble of 33 NA-CORDEX GCM–RCM combinations (RCP 4.5 n = 6, RCP 8.5 n = 27). The shaded areas represent the 25th and 75th percentiles (the variability between individual ensemble members for precipitation is presented in Fig. A1).

3 Results

3.1 NA-CORDEX meteorological forcing data

Bias-corrected NA-CORDEX model ensemble projections show an increase in future winter air temperatures at TVC, which intensifies under RCP 8.5, with a 9 °C increase from the present to the future in median January air temperatures (Fig. 1a and b). The timing and rate of air temperature increase in the spring is comparable between the present and the future, whereas the rate of cooling after the summer maximum decreases in the future under RCP 8.5 (a decrease of 6.8 °C per month in the present and 5.7 °C per month in the future for July–December). Under RCP 8.5, significant increases in precipitation are projected from July–December, the largest increase occurring in August from 41 to 50 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$ (two sample K-S test, test statistic: D = 0.16, p < 0.05, Fig. 1d). In shoulder season months (April, May, September, and October), projected air temperatures increase on average by 1.6 °C under RCP 4.5 and 3.7 °C under RCP 8.5. Warmer future air temperatures influence the phase of future precipitation with more precipitation falling as rain than snow, especially in transitional shoulder seasons between summer and winter (Fig. 2).

Figure 2 The half-violin plots show the frequency distribution of precipitation events as a function of air temperature for TVC in April, May, September, and October under RCP 4.5 (a, c, e, g) and RCP 8.5 (b, d, f, h), from an ensemble of 33 NA-CORDEX GCM–GCM combinations (RCP 4.5 n = 6, RCP 8.5 n = 27). The white violins with black outlines represent 2016–2046 and red violins 2066–2096. The height of each violin represents the frequency of precipitation events at a given temperature. The solid and dashed black lines at 0 and 2 °C show transitional temperatures between snow and rain, where CLM5.0 treats the transition in the precipitation phase as a linear ramp. The inset boxplots show the monthly total precipitation for the two 30-year periods.

Changes in air temperature in May and October cause a phase shift (snow to rain) in the total monthly precipitation under RCP 8.5 from the present to the future (Fig. 2d and h). In May, in the present (2016–2046) 50 % of precipitation events in the current period fall within 0–20 °C (as rain or as a mix of snow and rain) compared to 87 % in the future (2066–2096). In October, only 2 % of precipitation events fall within 0–20 °C compared to 40 % in the future. Furthermore, in May, future (red) peak precipitation frequency occurs at 2.7 °C compared to −0.6 °C at present (black), and in October at −0.5 °C compared to −4.6 °C at present (Fig. 2) under RCP 8.5, which indicates a shoulder season precipitation shift from majority snow to a mix of rain/snow. By contrast, differences between the current and future total precipitation are small (Fig. 2a–h, insets), and increases from the present to the future are significant only in September (two sample K-S test: 31 to 40 mm, D = 0.17, p < 0.05) and October (two sample K-S test: 23 to 28 mm, D = 0.12, p < 0.05) under RCP 8.5. As such, changes in the magnitude and phase of precipitation in the shoulder seasons may influence future snow-season length, snow thermal conductivity (due to an increase in melt conditions), and the subsequent temperature and moisture of soils.

3.2 CLM5.0 simulations

3.2.1 Environmental conditions

An increase in future rainfall in September impacts both the timing of winter onset of snow accumulation and the magnitude of peak SWE later in the winter (Fig. 3). The duration of winter snow cover is shorter in the future where future spring melt-out occurs between 18 to 25 d earlier in both RCP scenarios, while the timing of the initial snow accumulation remains similar (± 4–6 d). The projected future peak annual SWE shows an increase of 14 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$ under RCP 4.5 and a decrease of 6 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$ under RCP 8.5 and early spring snowmelt is coincident with earlier increases in soil temperature in mid-late spring (Fig. 3a–d). Compared to CORDEX-Jordan, snow thermal conductivity using CORDEX-Sturm is lower on average by 0.07 and 0.14 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ under RCP 4.5 and RCP 8.5, respectively.

Figure 3 The CLM5.0 simulated median daily SWE (a, b), soil liquid water (12 cm) content (c, d), snow thermal conductivity (e, f), and 10 cm soil temperature (g, h) over two 30-year time periods: 2016–2046 (black) and 2066–2096 (red) under RCP 4.5 and RCP 8.5, for TVC, using input meteorological data from an ensemble of 33 RCM–GCM combinations (RCP 4.5 n = 6, RCP 8.5 n = 27). The solid and dashed lines show the ensemble median values for the CORDEX-Jordan and CORDEX-Sturm experiments, respectively.

Using CORDEX-Sturm, the median minimum winter soil temperature is 4–7 °C warmer by the end of the century (2066–2096) compared to CORDEX-Jordan under both RCP scenarios (Fig. 3g and h). Warmer air temperatures predicted for the end of this century (Fig. 1a and b) will warm the soil, and snow cover will then insulate the soil from the colder atmosphere, maintaining these elevated soil temperatures throughout the winter (Fig. 3g and h).

Due to earlier snowmelt, soils are predicted to thaw earlier in the year under both RCP scenarios (RCP 4.5 = 19 d, RCP 8.5 = 26 d; Fig. 3e and f). The timing of soil moisture increase in May and June remains similar in the present (2016–2046) between RCP 4.5 and 8.5 but occurs 20–30 d earlier in the future (2046–2096). By comparison, the timing of soil moisture decrease in September and October remains similar at Julian day 257–267 between the two RCP scenarios and 30-year time periods. The modelled soil moisture plateaus at 10 kg m⁻² from October–May and reaches saturation 39.167 kg m⁻² from late May–early October.

Figure 4 The CLM5.0 simulated median daily soil temperature with depth over two 30-year time periods: 2016–2046 and 2066–2096 under RCP 4.5 (a, b) and 8.5 (c, d). The first row plots show CLM5.0 under the default K_eff parameters (Jordan, 1991) and the second row plots show an alternative parameterisation for snow thermal conductivity (Sturm et al., 1997). The solid black lines indicate 0.75 °C and the dashed −0.75 °C, indicating the zero-curtain period between them.

Evidence of a late winter (February–March) zero-curtain period is minimal across all but one scenario (Fig. 4h); however, soil temperatures at depths up to 10 m suggest that the zero-curtain duration from August–December increases by 10 and 36 d in the future under CORDEX-Jordan and CORDEX-Sturm, respectively (Fig. B1). In the future under RCP 8.5, autumn shoulder season (September–October) near-zero temperatures penetrate up to 6 m deeper into the soil column compared with the present day in the CORDEX-Sturm simulations (Fig. 4h). An earlier spring snowmelt, increasing winter soil temperatures, and longer zero-curtain duration combine to create the conditions for increased cumulative seasonal carbon emissions in the future at TVC.

3.2.2 CO₂ and CH₄ emissions

Cumulative snow-covered CO₂ and CH₄ fluxes are projected to increase under both RCP scenarios (Fig. 5), with the highest CO₂ increase produced for CORDEX-Sturm (34.2 g CO₂ m⁻²; from 21.1 to 55.3 g CO₂ m⁻²) under RCP 8.5. An earlier onset of snow in 2066–2096 facilitates earlier emissions of CO₂ and CH₄ during the autumn under RCP 4.5 (12–13 d post 1 September; Fig. 5a and c) compared with RCP 8.5 (41–47 d post 1 September; Fig. 5b and d). A delayed snow onset reduces overall cumulative carbon emissions during the snow cover period under RCP 8.5 (Fig. 5b and d) compared with RCP 4.5 (Fig. 5a and c). By the start of November, cumulative CO₂ emissions reach 50 %–90 % of their winter totals, highlighting the importance of autumn shoulder season emissions to the winter CO₂ budget. Under RCP 4.5 there is minimal increase in cumulative CO₂ after mid-late November (after day 80) with increases of 0.02 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and 0.00015 $\mathrm{g}{\mathrm{CH}}_{\mathrm{4}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Comparatively, under RCP 8.5, higher levels of emissions continue deeper into the winter, with average increases of 0.070 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and 0.003 $\mathrm{g}{\mathrm{CH}}_{\mathrm{4}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (day 80 to day 240) (Fig. 5). Applying the Sturm et al. (1997) snow thermal conductivity scheme within CLM5.0 increases the total accumulation of winter CO₂ at TVC from 2016–2100 by 50 %–150 % and CH₄ by 50 %–74 % across both RCP scenarios compared to CORDEX-Jordan. Fig. D1 presents a comparison between cumulative emissions limited to when snow is on the ground as per Fig. 5 and those unconstrained by snow cover. The cumulative emissions not limited to the snow-covered period are 3–4 times higher than those constrained by snow and highlight the impact of the early winter period on accumulated winter soil carbon emissions.

Figure 5 The CLM5.0 simulated cumulative CO₂ and CH₄ flux to the atmosphere during snow cover periods over two 30-year time periods: 2016–2046 (black) and 2066–2096 (red) under RCP 4.5 and RCP 8.5, for TVC, using the meteorological forcing from an ensemble of 33 RCM–GCM combinations (RCP 4.5 n = 6, RCP 8.5 n = 27). The solid lines show the median values from the CORDEX-Jordan ensemble and the dashed lines show values for the CORDEX-Sturm ensemble. The variable plot line lengths are indicative of changes in the snow-season length.

Figure 6 The CLM5.0 simulated median SR for the winter season, comparing the CORDEX-Jordan ensemble with CORDEX-Sturm for two RCP scenarios (RCP 4.5 n = 6, RCP 8.5 n = 27). The shaded areas represent the 75th and 25th percentiles and represent the ensemble distribution. The ensemble median and uncertainty calculations incorporate a range of Ψ_min, Q₁₀, and Q₁₀ch4 (Fig. E1).

From 2016–2046 to 2066–2096, the median winter CO₂ flux increases by 18 % (0.11–0.13 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) under RCP 4.5 and more than triples (0.08 to 0.3 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) under RCP 8.5, using CORDEX-Sturm. CORDEX-Sturm significantly increases the winter CO₂ emissions compared with CORDEX-Jordan, by 57 % (2016–2046) to 62 % (2066–2096) in RCP 4.5 and 100 % (2016–2046) to 130 % (2066–2096) in RCP 8.5 (two sample K-S test: D = 0.26, p < 0.01, Fig. 6). Likewise, under RCP 8.5 CORDEX-Sturm significantly increases the median CH₄ emissions by 50 % during 2066–2096 compared to CORDEX-Jordan (2.5$\times {\mathrm{10}}^{-\mathrm{4}}$ to 3.8$\times {\mathrm{10}}^{-\mathrm{4}}$ $\mathrm{g}{\mathrm{CH}}_{\mathrm{4}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) (two sample K-S test: D = 0.13, p < 0.01, Fig. 6), further showing the impact of snow representation on the simulation of CH₄ from Arctic soils. Between RCP 4.5 and 8.5, the median CO₂ emissions more than double (0.13 to 0.30 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) for the last 30 years of the century using CORDEX-Sturm. Simulations with reduced Ψ_min (−20 MPa) consistently stimulate higher carbon fluxes under RCP 4.5 and 8.5 compared to simulations using the default Ψ_min value (−2 MPa) (Fig. E1). This is particularly evident under CORDEX-Sturm, which shows average increases of 0.13–0.16 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and 0.0014–0.0018 $\mathrm{g}{\mathrm{CH}}_{\mathrm{4}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ from 2016–2100. CLM5.0 simulations indicate that warmer soil conditions induced by climate warming under RCP 8.5 intensify the winter carbon emissions from TVC soils in the future.

4 Discussion

4.1 Future temperature, precipitation, and snow cover

NA-CORDEX forecasts suggest rising future winter air temperatures (Fig. 1) which may increase winter soil heterotrophic respiration and contribute significantly to winter carbon emissions, offsetting growing season carbon uptake in the future (Natali et al., 2019; Commane et al., 2017). The mobilisation of “legacy” carbon (Pedron et al., 2023) within the Arctic permafrost is of concern in a warming climate (Schuur et al., 2015) and warmer winters will enhance future climate change through carbon output from the decomposition of soil carbon. Understanding changes in future Arctic air temperature and meteorology is key in determining the influence of a warming world on Arctic winter carbon emissions and vice versa. The NA-CORDEX ensemble simulations demonstrate a shift from a snow- to rain-dominated Arctic in the future, reported elsewhere in the literature (Mccrystall et al., 2021), which is intensified under RCP 8.5 (Fig. 2). The shift in precipitation phase in September and October will lead to a change in snow properties, particularly the effective conductivity of the snowpack. Additionally, melt conditions and rain-on-snow events in the initial stage of the snow cover season provide energy for melt through the condensation of water vapour onto the snowpack (Harr, 1981; Mazurkiewicz et al., 2008). An increase in the precipitation fraction falling as rain is therefore a likely driver of earlier snowmelt in CLM5.0 simulations (Fig. 3). Future changes in precipitation have key, yet uncertain, implications for the Arctic surface and soil energy balance with regard to snow cover depth, timing, and duration, which are major factors controlling the temperature variability in the upper 3 m soil (Callaghan et al., 2011).

4.2 Future soil carbon emissions

The CLM5.0 parameterisations for K_eff, Q₁₀, and Ψ_min explored by Dutch et al. (2024) were found to be highly suitable for representing winter soil temperatures and carbon fluxes under present-day conditions. This alignment between observations and simulations increases confidence in the model's ability to simulate future Arctic soil processes through to 2100. Our simulations show that by the end of the century (2066–2096), CLM5.0 simulated median winter soil CO₂ emissions under RCP 4.5 are predicted to be less than half of those under RCP 8.5, showing that projected winter carbon emissions are tied to the magnitude of global temperature increase. Soil temperatures are linked to changes in the atmosphere by the thermal conductivity of the snow layer. Under RCP 8.5, cumulative winter carbon emissions from soils at TVC are projected to increase in the future despite a reduction in the snow cover duration. A reduced K_eff (0.07–0.14 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), as introduced by CORDEX-Sturm (Fig. 3), increases cumulative winter CO₂ and CH₄ emissions by 50 %–150 % for 2016–2100 compared with CORDEX-Jordan under both RCP 4.5 and 8.5. By comparison, Natali et al. (2019) projected a 41 % increase in CO₂ emissions by 2100 compared to the present. Additionally, under RCP 8.5, the magnitude of the influence of CORDEX-Sturm on winter carbon fluxes is comparable to the uncertainty associated with future climate projections, and snow-related processes – snow cover extent and duration (Natali et al., 2019) and snow thermal conductivity (this study) – emerge as common critical drivers of this uncertainty.

The response of carbon fluxes from permafrost zones is highly sensitive to hydrological change and an increase in CO₂ emissions post-2050 (Fig. 6) is indicative of soil drying where, as temperatures increase, soil CH₄ production is strongly suppressed (Lawrence et al., 2015). The zero-curtain period is an important period for facilitating cold season emissions from tundra ecosystems (Tao et al., 2021) and an increasing proportion of the soil at or around 0 °C presents a risk of elevated carbon emissions from Arctic soils in the future. The seasonality of soil temperature is critical in controlling winter CO₂ emissions from soils, particularly the length of time at which soil remains at or near 0 °C at the beginning of winter (Schimel et al., 2006). Future CO₂ and CH₄ emissions show greater seasonal variability under CORDEX-Sturm compared with CORDEX-Jordan, particularly under RCP 8.5, which suggests that soil moisture and thermal dynamics are more sensitive to snow cover in the CORDEX-Sturm configuration (Fig. E1). Persistent carbon emissions throughout the winter period are partly attributable to a permanently unfrozen active soil layer (Zona et al., 2016) and its influence in the future will be intensified by a longer zero-curtain period. Simulations performed in this study indicate an increase in both the duration and depth of the early winter zero-curtain period under future climate conditions (Fig. 4). Under CORDEX-Sturm, the projected early winter zero-curtain period extends up to 26 d longer and reaches depths up to 6 m deeper than under CORDEX-Jordan, further highlighting the impact of snow representation on simulated soil temperatures. Increases in the duration of the zero-curtain period is concurrent with measured borehole data from 2006–2015 at an Arctic tundra site, which show an increase of up to 20 d, with large emissions of CO₂ between September and December in years with a longer zero-curtain period (Euskirchen et al., 2017; Larson et al., 2021). A deeper active layer (Aalto et al., 2018) and an increase in unfrozen soil (Schaefer and Jafarov, 2015; Natali et al., 2019; Elberling and Brandt, 2003) as a result of a longer and deeper zero-curtain period will increase SR in the future, accelerate the degradation of permafrost and likely mobilise the legacy carbon from Arctic soils.

Future changes in the shoulder season air temperature, soil moisture, and snow cover control both CO₂ and CH₄ due to moisture limitation, oxygen limitation, and soil temperature. A major fraction of cold season CH₄ emissions occur in the early winter shoulder season, particularly during the zero-curtain stage (Zona et al., 2016) and interannual variability in both CO₂ and CH₄ depends largely on this period (Kittler et al., 2017). Soil moisture fluctuations (Fig. 3) critically impact the rates of microbial decomposition of organic matter via methanogenesis, which favours wet, anoxic environments. CH₄ emissions are closely tied to soil moisture in the upper 30 cm of soil and are closely correlated with soil moisture fluctuations during the soil freeze-in period (Sturtevant et al., 2012). In the early winter (September–October) the soil profile is saturated, which drives anoxic conditions favoured by methanogens, which increases with depth (Arndt et al., 2020). Soil moisture increases in late September between RCP 4.5 and 8.5 (Fig. 3) resulting in wetter soil, alongside increased zero-curtain period at depth (Fig. 4) under RCP 8.5, present a risk of increased CH₄ emissions during this early winter period in the future at TVC. Such soil moisture increases, by slowing the soil freezing process have implications for the extension of elevated CH₄ emissions longer into the winter season (Sturtevant et al., 2012). Understanding the potential impacts of climate warming on both the early and late winter shoulder seasons is therefore key for assessing the risk of elevated winter carbon emissions from tundra soils. As well as persistent winter carbon emissions, “pulses” of CO₂ and CH₄ have been observed during both the early and late winter shoulder seasons (Raz-Yaseef et al., 2017; Mastepanov et al., 2013) where trapped gases are released as soils fluctuate between freeze and thaw. Soil emission pulses, which are not simulated by CLM5.0, are enhanced by an increase in rain-on-snow events. Such events accelerate soil warming in the spring, making soils more susceptible to cracking and rapid gas release (Raz-Yaseef et al., 2017). A shift towards a rain-dominated Arctic (Fig. 2) may therefore increase the abundance of spring gas pulses from soils in both frequency and distribution, increasing cumulative carbon emissions annually. Further work is required to better understand the impact of these pulses on annual carbon budgets and their representation within LSMs.

5 Conclusions

Warming air and soil temperatures are facilitating higher rates of heterotrophic respiration in CLM5.0 simulations for TVC soils, resulting in an expected increase in winter carbon output in the future under RCP scenarios 4.5 and 8.5. Such winter emissions outpace those modelled by Natali et al. (2019) and therefore reinforce the vulnerability of Arctic soils and the carbon stored within to climate warming. Soil temperature simulations show a lengthening of the zero-curtain period in the future as well as deeper penetration of near-zero temperatures into the soil column which will mobilise legacy carbon in the future. A longer zero-curtain period, leading to a persistently unfrozen active layer of soil, where higher rates of soil respiration persist for longer into the snow cover period, poses a risk of significantly increased cumulative carbon output in the future from tundra soils.

The modelled winter carbon emissions demonstrate the importance of climate mitigation in preventing a significant increase in the Arctic winter carbon budget. The use of parameterisation for K_eff from Sturm et al. (1997) reduces soil temperature biases and presents significantly higher CO₂ production from 2016–2100 compared to the default parameters from Jordan (1991). An improved representation of K_eff exacerbates already increasing CO₂ emissions from the Arctic tundra caused by a warming climate and the influence of the selected K_eff parameterisation is shown to be as important, if not more, than the variability in the future climate on simulated carbon emissions from Arctic tundra. Projected CO₂ and CH₄ emissions are highly sensitive to parameters K_eff, Ψ_min, Q₁₀, and Q₁₀ch4, which govern SR. We find that lower Ψ_min consistently increases cold season carbon fluxes and higher Q₁₀ suppresses them, which aligns with the findings of Dutch et al. (2024). Implementing the Sturm et al. (1997) K_eff parameterisation increases the sensitivity of the modelled carbon emissions to Ψ_min and Q₁₀ compared with the default K_eff parameterisation. Further work should aim to improve representations of snow thermal conductivity in LSMs while increasing the spatial coverage of future simulations to allow a more holistic outlook for Arctic soil carbon emissions and better inform climate mitigation strategies and carbon budgets.

Appendix A: NA-CORDEX future meteorology

Figure A1 The NA-CORDEX ensemble median monthly precipitation (January–December) split by individual ensemble member (RCP 4.5 n = 6, RCP 8.5 n = 27).

Appendix B: Zero-curtain duration

The CLM5.0 simulations indicate minimal changes in the late winter zero-curtain duration but a notable increase in the early winter zero-curtain duration. The early winter has been identified as a key period in facilitating carbon emissions from Arctic tundra, particularly CH₄ (Tao et al., 2021). The variation in future zero-curtain duration from August to December is presented below.

Figure B1 The August–December zero-curtain duration with soil depth for the CORDEX-Jordan (solid black and red lines) and CORDEX-Sturm (dashed black and red lines) for 2016–2046 and 2066–2096.

To further constrain the CLM5.0 parameter uncertainties, the values of Ψ_min and Q₁₀ were adjusted alongside two parameterisations for K_eff (Jordan and Sturm). A range of realistic Ψ_min and Q₁₀ values were chosen as per Dutch et al. (2024) alongside a range of Q₁₀ch4.

Appendix C: CLM5.0 parameter assessment space

The full combination of applying changes in K_eff, Ψ_min, Q₁₀, and Q₁₀ch4 as well as the variability provided by the CORDEX ensemble is presented in the main document. Figure E1 shows the upper and lower bounds of the simulated CO₂ and CH₄ when parameter adjustments are applied, with Q₁₀ = 7.5, Ψ_min = −2, Q₁₀ch4 = 4.0 being the lower extreme and Q₁₀ = 1.5, Ψ_min = −20 MPa, Q₁₀ch4 = 1.3 being the upper extreme.

Figure C1 Schematic showing the range of Q₁₀ and Ψ_min values applied to CLM5.0 using the CORDEX-Jordan (blue) and CORDEX-Sturm (red) for the modelling experiment.

Appendix D: Cumulative soil carbon emissions as constrained by “snow-down” and irrespective of snow

Figure D1 The CLM5.0 simulated cumulative CO₂ and CH₄ flux to the atmosphere constrained by the snow-down period (red, black) and irrespective of snow (green, blue) over two 30-year time periods: 2016–2046 (black, green) and 2066–2096 (red, blue) under RCP 4.5 and RCP 8.5, for TVC, using the meteorological forcing from an ensemble of 33 RCM–GCM combinations (RCP 4.5 n = 6, RCP 8.5 n = 27). The solid lines show the median values from the CORDEX-Jordan ensemble and the dashed lines show the values for the CORDEX-Sturm ensemble. The variable plot line lengths are indicative of changes in the snow-season length.

Appendix E: Cumulative soil carbon emissions as constrained by “snow-down” and irrespective of snow

Figure E1 The CLM5.0 simulated median SR for the winter season, comparing the CORDEX-Jordan (blue) ensemble with CORDEX-Sturm (red) for two RCP scenarios (RCP 4.5 n = 6, RCP 8.5 n = 27). The shaded areas represent the 75th and 25th percentiles and represent the CORDEX ensemble distribution. Each plot shows the extremes of Q₁₀, Ψ_min, and Q₁₀ch4 from the chosen parameter values seen in Fig. C1, where the lower end Ψ_min = −2, Q₁₀ = 7.5, Q₁₀ch4 = 4 suppresses carbon output and the upper end Ψ_min = −20 MPa, Q₁₀ = 1.5, Q₁₀ch4 = 1.3 stimulates the carbon output. The average difference between the displayed parameters (solid versus dashed lines, Δ) is included on each subplot.

Appendix F: Frequency distribution of precipitation events in March, June, August, and September 2016–2100

Figure F1 The half-violin plots show the frequency distribution of precipitation events as a function of air temperature for TVC in March, June, August, and November under RCP 4.5 (a, c, e, g) and RCP 8.5 (b, d, f, h) over two 30-year time periods: from an ensemble of 33 NA-CORDEX GCM–GCM combinations (RCP 4.5 n = 6, RCP 8.5 n = 27). The white violins with black outlines represent 2016–2046 (black) and the red violins 2066–2096 (red) under RCP 4.5 and RCP 8.5, for TVC, using the meteorological data from an ensemble of 33 NA-CORDEX GCM–GCM combinations (RCP 4.5 n = 6, RCP 8.5 n = 27). The height of each violin represents the frequency of precipitation events at a given temperature. The solid and dashed black lines at 0 and 2 °C show the transitional temperatures between snow and rain, where CLM5.0 treats the transition in precipitation phase as a linear ramp. The inset boxplots show the monthly total precipitation for the two 30-year periods.