The importance of interactions between snow, permafrost and vegetation dynamics in affecting terrestrial carbon balance in circumpolar regions

The LPJ-STM model (Jiang et al 2016) is developed by incorporating a soil thermal model (STM) (Jiang et al 2012) into the standard LPJ model (Sitch et al 2003, Gerten et al 2004). It has been applied to simulate large-scale vegetation distribution, soil thermal regime and carbon cycling (Jiang et al 2012, 2016). In LPJ-STM, surface air temperature is treated as the surface soil temperature to simulate the soil thermal regime. In this study, we replace this simple scheme by considering the influences of snow and changing plant canopy due to vegetation dynamics on surface soil temperature dynamics. The model is then used to evaluate the impacts of vegetation-permafrost interactions on carbon dynamics in the region.

2.1.1. Snow scheme

Following Lyu and Zhuang (2018), we treat snow thickness and snow thermal conductivity explicitly in a snow-soil continuum. Snow thickness is estimated from precipitation and air temperature (Pongracz et al 2021). We assume that the snow bottom temperature is the same as soil surface temperature, which is calculated as the following when snow exits:

$$\begin{equation}{T_{{\text{soil}}\;{\text{surf}}}} = {\text{ }}\frac{{{T_{{\text{air}}}} + r * {T_{{\text{soilbase}}}}}}{{1 + r}}\end{equation} \tag{ 1 }$$

where T_air is surface air temperature, T_soilbase is soil bottom temperature. r is calculated as:

$$\begin{equation}r = {\text{ }}\frac{{{K_{{\text{soil}}}}*{Z_{{\text{snow}}}}}}{{{K_{{\text{snow}}}}*{\text{ }}{Z_{{\text{soil}}}}}}\end{equation} \tag{ 2 }$$

where K_soil is soil thermal conductivity (Wm⁻¹K⁻¹), K_snow is snow thermal conductivity (Wm⁻¹K⁻¹), Z_snow is snow depth (m), and Z_soil is soil column depth (m). More details can be found in Lyu and Zhuang (2018).

2.1.2. Canopy energy balance scheme

Following Gu (1998) and Gu et al (1999), we treat the canopy and the soil as two different layers (figure S1). The energy balance of both canopy and soil surface is modeled as:

$$\begin{equation}{\text{Q}} = {\text{SH}} + {\text{LH}} + {\text{SW}} + {\text{LW}}\end{equation} \tag{ 3 }$$

where Q is the energy budget, SH is the sensible heat flux between the canopy (soil) and the air, LH is the latent heat flux, SW is the short-wave radiation absorbed by the canopy (soil) and LW is the long-wave radiation of the canopy (soil). The sensible heat flux is calculated based on the temperature and the near-surface wind velocity. The latent heat flux is estimated with the temperature and the soil water content. The short-wave radiation and long wave radiation are estimated with the leaf area index (LAI). More details can be found in the supporting information.

In this study, monthly mean surface air temperature, cloudiness, precipitation, number of wet days and atmospheric CO₂ concentrations are used to drive the original model. For revised model simulations, in addition to those data used for original model simulation, near-surface wind velocity and relative humidity are also used to drive the revised models. Historical climate data for the period of 1959–2021 are from ERA5 (Hersbach et al 2020). Three climate scenarios from IPCC representative concentration pathways (RCPs): RCP2.6, RCP4.5 and RCP8.5 (Jones et al 2011) are used for the simulations from 2022 to 2100. To eliminate the errors between the ERA5 dataset and RCP datasets, we modify the RCP data based on the monthly climatology from 2006 to 2021. We first calculate the monthly climatology for each grid, then the ratio between the two datasets for each variable. For the temperature, we calculate the difference. Finally, we modify the RCP datasets by multiplying this ratio, but for temperature, plus the difference.

The models are spun-up for 1000 years before 1959 with cyclic climate data from 1959 to 1988 to achieve an equilibrium state. For 1959–2100, the historical data (ERA5, 1959–2021), followed by each scenario (2022–2100), are used to drive the simulation at a spatial resolution of 0.5° × 0.5° (total 25 063 grid cells) for the circumpolar region (45–90° N). Since both snow and canopy can influence the soil surface temperature, we organize three sets of simulations: the original LPJ-STM, the modified model that only contains the scheme of the snow (LPJ-STM-S) and the modified model that contains both the schemes of snow and canopy (LPJ-STM-SC). By comparing these results between LPJ-STM-S and the original LPJ-STM, we can estimate the snow influences on vegetation, permafrost, and carbon dynamics. The differences between LPJ-STM-S and LPJ-STM-SC simulations reveal the impacts of the canopy.

The model is calibrated using the observed soil temperature from AmeriFlux for 20 sites. Table S1 in supporting information gives the parameters for model improvement. The calibration is conducted with Model-Independent Parameter Estimation (PEST; Doherty 2004, Doherty et al 2022) which is a software using Bayesian method to calibrate the model. It requires iterating the model simulations to find the optimal values for the parameters. Figure 1 shows the results from three selected sites. The LPJ-STM-SC simulated soil surface temperatures agree well with the observations for all the three sites (Amiro 2016, Desai 2022, Ueyama et al 2022) while the air temperatures show relatively large differences to the observations. The RMSEs between the simulated and observed soil surface temperature (RMSE_s) are all much smaller than the RMSEs between the air temperature and the observation (RMSE_a) for all the sites. The seasonal fluctuations get dampened with higher simulated temperature in winter (mainly due to the snow) and lower temperature in summer (mainly due to the canopy). LPJ-STM-SC estimates the soil surface temperature better than LPJ-STM that uses surface air temperature as upper boundary conditions for soil thermal model.

Zoom In Zoom Out Reset image size

Figure 2 shows the differences of the simulated monthly mean soil surface temperature due to snow effects (in January) and canopy effects (in July) for the period of 2001–2010 (historical) and 2091–2100 (RCP2.6, RCP4.5 and RCP8.5). For historical results, we can see that in January (winter), the soil surface temperature is generally much warmer than the air temperature due to the insulation of the snow, especially for the East Europe and Northeast Canada where the warming can be up to 10 °C. During July (summer), the simulated soil surface temperature is consistently lower than the air temperature due to the influences of the canopy, with a magnitude of 0.5 °C–4 °C. For the future, the temperature differences show similar patterns under each warming scenario. The snow's warming effect is still quite obvious except for the West Europe area under RCP8.5 since this region is relatively warmer than other areas (Minobe et al 2008) and has less snow cover.

Zoom In Zoom Out Reset image size

Deep soil (20 cm) temperature differences exhibit similar patterns to the soil surface (figure S2) with warming effect of the snow in winter and cooling effect of the canopy in summer for the 2000s. The magnitude of the changes is lower since the deep soil is less impacted by the air and land surface processes.

LPJ-STM simulates vegetation dynamics based on the growth, competition and disturbances. Hence, we can assess the influences of the soil temperature change on vegetation distribution. In this study, we examined six plant function types (PFTs) for the circumpolar region (table S1). Here, we compare the differences of three dominant PFTs. LPJ-STM-S estimates less boreal needle-leaf evergreen tree (BNE) and more boreal broad-leaf summer-green tree (BBS) and C3 grass (figures 3(a)–(c)) coverage than the original LPJ-STM. The loss of BNE coverage (5.0 × 10⁵ km² in 2020) is first taken by the grass, then gradually by the BBS (2.3 × 10⁵ km² in 2020), since the warmer environment due to the snow favors the growth of the deciduous broad-leaf trees. Similarly, the canopy's cooling effect can lead to more BNE (up to 2.5 × 10⁵ km² during 1990–2000) and less BBS (1.1 × 10⁵ km² in 2020) coverage (figures 3(d)–(f)). These results demonstrate that the soil thermal regime can directly impact the vegetation distribution, in turn, affect the carbon dynamics.

Zoom In Zoom Out Reset image size

Compared with LPJ-STM, LPJ-STM-S estimates lower present-day NPP (figure 3(a)). Since the snow can significantly increase the winter soil surface temperature, the root respiration will increase when snow exits, resulting in the lower NPP (about −0.33 Pg C/year) (figure 4(a)). Also, the increased soil surface temperature leads to higher R_H (figure 4(b)) (about 1.21 Pg C/year). Thus, the NEP estimated by LPJ-STM-S is less than LPJ-STM (about −1.54 Pg C/year), with more carbon released to the atmosphere due to the influences of the snow. LPJ-STM-S also estimates lower vegetation carbon (figure S3(a)) and lower soil carbon (figure S3(b)). The differences get larger with the time goes on and the largest absolute values are about 15 and 45 Pg C (2010–2020), respectively.

Zoom In Zoom Out Reset image size

Different from the snow effects, LPJ-STM-SC estimates higher NPP (figure 4(d)) than LPJ-STM-S (about 0.24 Pg C/year) since the canopy will decrease the soil surface temperature, especially in summer (figure S3). The lower soil surface temperature will lead to lower soil respiration and higher NPP. Similarly, lower estimated R_H (about −0.76 Pg C/year) by LPJ-STM-SC (figure 4(e)) can be found for the period of 1959–2021. In this way, the NEP estimated by LPJ-STM-SC is higher than LPJ-STM-S (about 1.00 Pg C/year) with more carbon sequestered in the ecosystem. The estimated vegetation carbon and soil carbon pools are higher, with peak values of about 8 and 30 Pg C (2010–2020), respectively.

Under warming scenarios, LPJ-STM-S still estimates lower NPP than LPJ-STM and the differences between the two models are generally larger than that in the historical period (about −0.5 Pg C/year). Changes in R_H show similar patterns to the historical period but the differences gradually decreased (from 1.0 to 0.5 Pg C/year for RCP2.6 and RCP4.5). Especially, for RCP8.5, the R_H differences are around 0 at the end of the century. As the difference of NPP and R_H, the NEP estimated by LPJ-STM-S is lower under all the scenarios, indicating the impact of snow will release more carbon into the atmosphere in the future. In addition, the differences between the two models show a trend of decreasing, especially under RCP8.5 where the differences of NEP are close to 0 at the end of the century. Similarly, the differences of vegetation carbon remain negative, but the difference decreases, especially under the RCP8.5 (figure S3(a)). The soil carbon shows a consistent change under scenarios with more and more carbon released from the soil (figure S3(b)).

For the future, LPJ-STM-SC also estimates higher NPP and lower R_H than LPJ-STM-S, especially under the RCP8.5. Changes in R_H generally becomes weak over time and oscillates around 0 at the end of the century. As a result, the NEP differences are larger than 0 for all the scenarios, suggesting more carbon is sequestered in ecosystems due to the cooling effects of the canopy. LPJ-STM-SC also estimates higher vegetation carbon (figure S3(c)) and higher soil carbon (figure S3(d)) under all the scenarios while the differences of vegetation carbon get smaller from 2020 to 2080. The differences of the soil carbon keep increasing under all the scenarios, indicating there is a constant cooling effect of the canopy.