On the Arctic near-surface permafrost and climate sensitivities to soil and snow model formulations in climate models

Abstract

This study investigates the sensitivity of the Canadian Regional Climate Model (CRCM5) simulated near surface permafrost and its climate interactions to soil and snow formulations. In particular, sensitivities to the depth of the soil column, inclusion of organic soils and modified snow conductivity formulation are investigated. The impact of these modifications are first assessed in offline simulations performed with the Canadian Land Surface Scheme (CLASS), which is the land surface scheme used in CRCM5, when driven by ERA-40/ERA-Interim for the 1957–2008 period. Analysis of CLASS simulations shows major improvements in the simulated permafrost extent, particularly with a deeper soil column. Inclusion of organic soil decreased the summer ground heat flux and therefore the summer soil temperatures, leading to improvements in the simulated active layer thickness (ALT). The impact of the new snow thermal conductivity formulation is moderate compared to the effect of organic soils, but reduces the cold biases in winter soil temperatures. CRCM5 experiments revealed similar sensitivities to soil depth, organic soil and snow conductivity changes as with the offline simulations. Significant changes are noted in the land–atmosphere interactions, through modified energy and moisture partitioning at the surface resulting from the inclusion of the organic soils. The inter-annual variability of the ALT shows larger sensitivities to summer temperatures for mineral soil while experiments including organic soils show increased sensitivities to annual temperatures. The ALT trends in the CRCM5 are similar to the observed values, despite the overestimation of ALT associated with a warm bias in the CRCM5 climate.

1 Introduction

Permafrost covers approximately one quarter of all the exposed land of the continental northern hemisphere (Zhang et al. 1999). Many recent studies (Lemke et al. 2007; Osterkamp 2007; Frauenfeld and Zhang 2011) suggest increases in soil temperatures in these permafrost regions, particularly during the last few decades. These increases in soil temperatures are reflected in the increasing active layer thickness (ALT)—defined as the maximum annual thaw depth—and permafrost degradation. The permafrost degradation can alter significantly the soil structure and hydrology of the region, lead to the formation of thermokarst lakes and/or drainage of existing lakes, and change wetlands and vegetation cover (Hinzman et al. 2005). Permafrost degradation is also of concern from the viewpoint of the large amounts of carbon (~1,672 Pg of carbon; Tarnocai et al. 2009) that are currently sequestered in permafrost. Since the majority of the microbial decomposition occurs in the seasonally thawed active layer, increases in the ALT and permafrost thaw will increase decomposition of the soil organic matter, resulting in the release of carbon dioxide and/or methane to the atmosphere, depending on the type of decomposition, i.e. aerobic or anaerobic. This release of greenhouse gases, associated with permafrost thaw, will act as a positive feedback to climate warming (Zimov et al. 2006) estimated between 2.8 and 7.8 °C by the end of the century (Solomon et al. 2007) without the inclusion of this potential positive feedback. MacDougall et al. (2012) estimated that the climate warming from the permafrost carbon feedback could result in an additional warming in the 0.23–0.27 °C range by the end of the twenty-first century. Most climate models still do not consider all these above factors adequately, which contribute to large uncertainties in the projected climate changes (Schaefer et al. 2012).

The observed changes in soil temperatures, ALT and near-surface permafrost are not only related to higher temperatures, but also to changes in snow cover extent and duration (Lemke et al. 2007; Frauenfeld and Zhang 2011). Most of the permafrost studies available to date are based on land surface model (LSM) simulations driven by observed meteorological data (Oelke et al. 2004; Dankers et al. 2011) or by outputs from Global Climate Models (Sushama et al. 2007; Lawrence et al. 2008; Schaefer et al. 2011; Koven et al. 2011; Burke et al. 2013). These offline LSM simulations have provided important insights related to the evolution of permafrost, but do not represent the two-way feedbacks between the land and the atmosphere. Efforts to simulate permafrost interactively in climate models, both global and regional, is currently an active area of research. Studies by Smerdon and Stieglitz (2006) and Alexeev et al. (2007) documented limitations of climate models in simulating near-surface permafrost. They report that realistic simulation of soil temperatures in cold permafrost regions using LSMs with a zero-flux bottom boundary condition requires a total soil column depth of at least 30 m. Furthermore, Nicolsky et al. (2007), Lawrence et al. (2008) and Dankers et al. (2011) demonstrated the need to include soil organic carbon in LSMs for realistic simulations of soil thermal and moisture regimes and therefore ALT and permafrost extent, which in turn is important for realistic surface energy and water partitioning at the surface (Lawrence and Slater 2008; Rinke et al. 2008). Indeed, organic material acts as an insulator, with its low thermal conductivity and relatively high heat content. The implementation of soil organic carbon in LSMs therefore decreases the summer soil temperatures, decreasing the annual thaw depth compared to mineral soil formulations, increasing the area covered by near-surface permafrost in models for the same surface air temperatures (Lawrence and Slater 2008; Rinke et al. 2008; Dankers et al. 2011). Organic material also influences the soil moisture. In nature, most areas where permafrost is present are characterized by nearly saturated sub-surface conditions with a drier surface layer (Hinzman et al. 1991). Strong variability is observed in the soil moisture of the surface layer due to the efficient transport of water deeper in the soil column. This enhanced downward transport is due to increased porosity, high hydraulic conductivity and weak suction of organic material (Quinton and Gray 2003).

The simulated near-surface atmospheric fields in climate models will be sensitive to the parameterization method adopted for organic soils. Lawrence and Slater (2008) and Rinke et al. (2008) obtained different sensitivities of near-surface atmosphere to the implementation of a soil organic material in their respective models. Lawrence and Slater (2008) showed large increase in the sensible heat flux over most of the high-latitude regions of the northern hemisphere. These changes led to increased 2 m-air temperature, deeper and drier atmospheric boundary layer and reduced occurrence of low-level clouds. On the contrary, Rinke et al. (2008) described large increases in the latent heat flux, thereby cooling the surface temperature and enhancing low-level clouds. Given the large uncertainties in the atmospheric response, there is a need for more focused studies to identify the main reasons behind these uncertainties related to the representation of organic material in climate models.

The objectives of the present study are twofold. The first objective is to assess the sensitivity of the Pan-Arctic soil temperature and moisture regimes to the depth of the soil model, to the implementation of a soil organic carbon parameterization and to a modification of the snow thermal conductivity in the Canadian Land Surface Scheme (CLASS) through a series of offline experiments. This allows an evaluation of the Pan-Arctic representation of permafrost conditions excluding the complex land–atmosphere interactions and feedbacks. The second objective is to study the impact of the LSM improvements, particularly the organic soil parameterization, in coupled land–atmosphere experiments using the fifth generation of the Canadian Regional Climate Model (CRCM5) on the simulated permafrost, surface climate and land–atmosphere interactions.

The outline of this paper is as follows. Section 2 provides a brief description of CLASS and CRCM5. Model configurations and the soil organic carbon parameterization are presented in Sect. 3. Section 4 details the observational datasets and reanalysis used in model evaluation. Analysis of the offline simulations performed with CLASS followed by that of CRCM5 simulations are presented in Sect. 5. Summary and conclusions are presented in Sect. 6.

2 Model description

The offline LSM simulations presented in this study are performed with CLASS (Verseghy 1991, 2008; Verseghy et al. 1993), and the regional climate model used is a developmental version of the CRCM5 (Zadra et al. 2008; Martynov et al. 2013). One might note that the CRCM5 uses CLASS as its land surface scheme. A brief description of CLASS and CRCM5 follows.

2.1 The Canadian Land Surface Scheme

The basic prognostic variables in CLASS version 3.5 (Verseghy 2008) consist of the temperatures and the liquid and frozen moisture contents of the soil layers; the mass, temperature, density and albedo of the snow pack; the temperature and intercepted rain and snow on the vegetation canopy; the temperature and depth of ponded water on the soil surface; and an empirical vegetation growth index. At each time step, CLASS calculates the characteristics of the vegetation canopy on the basis of the vegetation types present over the modeled area. In a pre-processing step, 23 surface and vegetation types with assigned background values of parameters such as albedo, roughness length, annual maximum and minimum leaf area index, rooting depth, etc. are aggregated over four main vegetation categories in CLASS: needleleaf trees, broadleaf trees, crops, and grass (Verseghy 2008). Relevant to this study, tundra is represented as grass while the boreal forest and taïga shows a gradient from dense forested areas towards open forest and shrubs. Adjustments of the vegetation parameters allow realistic representation of variations in vegetation density.

CLASS is particularly suited for permafrost studies due to its flexible soil vertical formulation. Each grid cell corresponds to a single soil profile of sand, clay and bedrock at the grid resolution computed independently, with no lateral heat and soil moisture transfers between adjacent soil columns. CLASS uses sand and clay concentrations derived from Wilson and Henderson-Sellers (1985). A spatially varying depth to bedrock database is derived from Webb et al. (2000). If the depth to bedrock occurs within a soil layer rather than at the interface between two layers, CLASS assigns the specified soil characteristics to the fraction of the layer above bedrock, and values corresponding to solid rock to the portion below. Thermal conductivity of rock is the same as for sand, i.e. 2.5 W m⁻¹ K⁻¹, with zero porosity. Figure 1b presents the depth to the bedrock over the study domain.

Fig. 1

a CRCM topography (m); b depth to bedrock (m). c Soil organic concentration (kg m⁻²) from IGBP-DIS for the first 100 cm of soil. d Number of soil organic layers considered in Off_OM47, Off_OMSC47, C_OM47 and C_OMSC47. Red regions represent deep peatlands where Letts et al. (2000) parameterization is used while green (blue) regions represent grid points where 10 cm (30 cm) of organic soil are used

Snow in CLASS is modelled as a single variable depth layer. The effective thermal conductivity of snow, λ_s, is determined from the snow density, ρ _s , assumed constant with depth, as described in Mellor (1977):

$$\uplambda_{\text{s}} = 2.576 \times 10^{ - 6} \,\uprho_{\text{s}}^{2} + 0.0074$$

(1)

where λ_s, is in W m⁻¹ K⁻¹ and ρ _s is in kg m⁻³. To represent snow ageing, the snow density increases exponentially with time from a fresh snow value of 100–300 kg m⁻³, according to an expression derived from the field measurements of Longley (1960) and Gold (1958). CLASS treats refreezing of percolating melt-water or rain, which can lead to increases in snow density. After snowfall, ρ _s is recalculated as the weighted average of the previous density and that of the new snow.

2.2 The Canadian Regional Climate Model

The CRCM5 (Zadra et al. 2008; Martynov et al. 2013) is based on a limited-area version of the Global Environment Multiscale (GEM) model used for Numerical Weather Prediction at Environment Canada (Côté et al. 1998). GEM employs semi-Lagrangian transport and (quasi) fully implicit marching scheme. The following parameterizations are used in CRCM5: deep convection following Kain and Fritsch (1990), shallow convection based on a transient version of Kuo (1965) scheme (Bélair et al. 2005), large-scale condensation (Sundqvist et al. 1989), correlated-K solar and terrestrial radiations (Li and Barker 2005), subgrid-scale orographic gravity-wave drag (McFarlane 1987), low-level orographic blocking (Zadra et al. 2003), and turbulent kinetic energy closure in the planetary boundary layer and vertical diffusion (Benoit et al. 1989; Delage and Girard 1992; Delage 1997). CRCM5 simulations are performed with 56 atmospheric levels, with the model top at 0.1 hPa.

3 Model parameters and experiments

3.1 Soil carbon data and parameterization

In CLASS, high concentration of soil organic carbon can be represented with the peatland parameterization of Letts et al. (2000). This parameterization assumes that the entire soil column consists of organic material or peat. Three organic soil/peat classes: fibric, hemic and sapric, are considered for layers one, two and three and below to account for the effect of compaction of the organic material, with variations in the hydraulic properties of the uppermost 0.5 m of organic soil (Letts et al. 2000). Bedrock can be located at any depth in the soil column.

To account for organic material present outside of deep peatlands areas, a simple parameterization of soil organic carbon (SOC) has been introduced. This parameterization was implemented to improve the representation of moderate concentration of organic material and thereby its influence on the surface fluxes and the thermal and hydraulic properties of the soil. This parameterization redistributes the estimated SOC from the Global Soil Data Task Group of the International Geosphere-Biosphere Programme Data and Information System (IGBP-DIS) (Fig. 1c) assuming that high concentration of organic material is located at the surface, and the concentration decreases rapidly with depth. This assumption is reasonable since the SOC estimation from IGBP-DIS is representative of the first 100 cm of the soil column. SOC from IGBP-DIS is interpolated from its original 1° × 1° grid to the model grid. Vertical layers of CLASS are “filled” with organic carbon from the surface down until the observed soil carbon content is depleted. At the moment, no fractional percentage of soil organic carbon is allowed within a soil layer, resulting in either 100 % organic carbon content or 100 % mineral composition. For the vertical levels filled with organic material, the thermal and hydraulic properties follow Letts et al. (2000). The net effect is a 10 cm organic layer for most of the high latitudes and 30 cm for certain grid points located in Scandinavia and in the Ob River Valley (Fig. 1d).

The above implementation method is similar to that used by Rinke et al. (2008) where they assume that soil layers are entirely composed of SOC. Rinke et al. (2008) used three different organic classification: lichen, peat and moss, with parameters from Beringer et al. (2001) spatially distributed based on vegetation cover while in this study we distribute a single SOC category, with vertically varying parameters, based on observed values from the IGBP-DIS data. The introduction of 10 cm of pure organic soils at the surface is similar to other studies of Lawrence and Slater (2008) and Rinke et al. (2008) although CLASS only possess one vertical layer over that depth (see Sect. 3.2). Implementation of SOC in the Community Land Surface Scheme (CLM; Lawrence and Slater 2008) and in the Joint UK Land Environment Simulator (JULES; Dankers et al. 2011) uses a redistribution of the observed SOC assuming vertically varying profiles and allows fractional SOC to be present at any soil layer. All methods have important limitations caused by the spatial heterogeneity of the soil carbon distribution. In reality, SOC is mostly accumulating in valleys and wetlands while ridges usually show limited SOC concentration.

3.2 Model experiments

All experiments are performed over a Pan-Arctic domain, at 0.5° horizontal resolution on a rotated latitude-longitude grid (Fig. 1). To evaluate the impact of various formulations on the simulated near-surface permafrost, four pairs of experiments are performed. Each pair consists of an offline CLASS experiment and a matching CRCM5 experiment using identical CLASS formulations (see Table 1).

Table 1 List of experiments performed in this study

To evaluate the sensitivity of the soil thermal and moisture regimes to the depth of the soil column, two soil column configurations—shallow and deep—are considered. The first pair of simulations, Off_Mine6 and C_Mine6, uses a shallow configuration with six soil layers that are 0.1, 0.2, 0.3, 0.5, 0.9 and 1.5 m thick for a total depth of 3.5 m. The second pair of simulations, Off_Mine47 and C_Mine47, uses the deeper soil column where 41 additional layers of 1.5 m thickness were added to the shallow configuration for a total of 47 levels for a total depth of 65 m. These two pairs of simulation using only mineral soils allow an evaluation of the model sensitivity solely to the vertical extension of the soil.

To assess the impact of soil organic carbon, the third pair of simulations (Off_OM47 and C_OM47) includes both the peatland parametrization (Letts et al. 2000) and the additional SOC parameterization (Sect. 3.1) using a deep soil column. Finally, the fourth pair of simulations (Off_OMSC47 and C_OMSC47), allows to evaluate the model sensitivity to changes in the snow density—snow thermal conductivity relation by using the quadratic relation derived by Sturm et al. (1997) (Eq. 2):

$$\begin{aligned} \uplambda_{\text{s}} & = 0.023 + 0.234\,\uprho_{\text{s}} \quad {\text{for}}\, \left\{ {\uprho_{\text{s}} < 0.156} \right\} \\ \uplambda_{\text{s}} & = 0. 1 3 8- 1.0 1\,\uprho_{\text{s}} + 3. 2 3 3\,\uprho_{\text{s}}^{2} \quad {\text{for}} \,\left\{ {0. 1 5 6\le \uprho_{\text{s}} \le \, 0. 6} \right\} \\ \end{aligned}$$

(2)

where λ_s, is in W m⁻¹ K⁻¹ and ρ_s is in g cm⁻³. This relation leads to overall reduced snow conductivity compared to the Mellor (1977) formulation.

The offline CLASS simulations are driven by atmospheric forcings, i.e. precipitation, downward solar and longwave radiation, near-surface air-temperature and specific humidity, surface wind speed and air pressure, from ERA-40 reanalysis (Uppala et al. 2005) for the 1957–1994 period and from ERA-Interim data (Dee et al. 2011) for the 1995–2008 period. ERA-40 and ERA-Interim reanalysis were spatially interpolated from their original 2.5° and 1.5° grids, respectively, to the model grid at 0.5°. Both reanalysis were available at 6-hourly time interval, which were linearly interpolated to the model 30-min time step. No temperature or precipitation corrections were applied to the reanalysis data. For each configuration, soil temperature and moisture content were initialized using the respective fields from 200-year spinup runs using atmospheric conditions from a climatology derived from ERA-40 over the 1970–1999 period.

The CRCM5 simulations are driven at their lateral boundaries by ERA-40 and ERA-Interim reanalysis. All coupled experiments were spun-up for an additional 20 years using repeated 1,957 atmospheric lateral boundary conditions from ERA-40. This additional spinup was executed to assure that soil temperature and moisture in the first 20 m of the soil column were in equilibrium with the CRCM5 climate before executing hindcast experiments (no figure).

4 Data for model evaluation

4.1 Permafrost extent

The observed permafrost extent is derived from the International Permafrost Association (IPA) map (Brown et al. 1998). IPA classifies permafrost into four different categories based on the areal extent of permafrost: 90–100 % coverage is considered continuous permafrost, 50–90 % discontinuous, 10–50 % sporadic and less than 10 % isolated. The data from IPA represents an estimate of the permafrost extent valid for the second half of the twentieth century (Burke et al. 2013). The original IPA map was reported to the model grid at 0.5° resolution using the dominant near-surface permafrost category over the model grid cell. Given the model resolution, sporadic and isolated permafrost regions cannot be adequately captured by the simulations. For this reason, simulated permafrost extent is compared to continuous and discontinuous permafrost regions.

4.2 Active layer thickness

The observed ALT dataset available from the Circumpolar Active Layer Monitoring program (CALM) (Brown et al. 2003; Brown 1998) is used to evaluate the simulated ALTs. The Arctic part of the CALM network consists of 79 stations mainly located in the arctic and sub-arctic lowlands over period from 1990 to present. Three methods are mainly used to determine the ALT: by mechanical probing, thaw tubes; or by inferring the thaw depth from ground temperature measurements recorded by thermistors. More details about the data and methodology are available on the CALM website (http://www.udel.edu/Geography/calm/index.html). Comparison are performed for the CALM sites located north of 60°N, with an altitude difference of less than 150 m with the altitude of the model representative grid cell. This way, stations where altitude differences could lead to significant mismatch in the surface climate between the CALM site and the representative model grid are eliminated. As observed values are not available for all sites for the entire period, average observed and modeled values for the 1990–2008 period are compared.

4.3 Russian historical soil temperature dataset

The in situ soil temperature observations over Russia (Zhang et al. 2001) is used for evaluating simulated soil temperatures. This dataset is a collection of monthly soil temperatures measured at meteorological stations for 13 different depths from 0.02 to 3.2 m using bent stem thermometers, extraction thermometers, and electrical resistance thermistors. The stations are located in different climatic regions of Russia, between 35°E and 180°E, providing useful large-scale soil information for the model evaluation. Although the original dataset covers the period from 1882 to 1990, data coverage in space and time suffers from numerous gaps. Data concentration is relatively high in the 1980–1990 period and will constitute the reference period for comparison with model results. Measurements were generally performed over bare soil without a surface organic layer, therefore potentially overestimating the seasonal cycle of soil temperature at depth compared to areas covered by organic material (Burke et al. 2013; Gilichinsky et al. 1998). For all observed soil levels, stations located North of 55°N with at least one value available per month were selected. Monthly averages were computed over the 10-year period. All observed vertical levels falling within a model layer are averaged to maximize the number of stations for the comparison. Then, similar to ALT comparison, soil temperature comparisons are performed only for those observation stations with differences in altitude less than 150 m with the model representative grid cell.

4.4 Surface air temperature and precipitation

The model evaluation of 2 m-air temperature were performed against ERA-Interim (Dee et al. 2011) at 1.5° horizontal resolution and two station-based gridded observational datasets: the global gridded data from the University of Delaware (Udel; Willmott and Matsuura 1995) and from the Climatic Research Unit version 3.1 (CRU3.1; Mitchell and Jones 2005). These datasets combine weather station records and use spatial interpolation methods to produce monthly means at 0.5° horizontal resolution over land.

As a complement to the precipitation data from ERA-Interim and UDel, data from the Global Precipitation Climatology Center (GPCC; Schneider et al. 2011) is also used in this study. The version 6 of GPCC includes data from global stations to provide gridded monthly means of precipitation at 0.5° horizontal resolution for the 1951 to present period. No corrections were applied to any of the station-based observational datasests for the undercatch of solid precipitation in the Arctic.

4.5 Snow water equivalent

Two datasets of snow water equivalent (SWE) are used for model evaluation: the Global Snow Monitoring for Climate Research version 1.2 (GlobSnow; Luojus et al. 2010) and the Canadian Meteorological Center snow analysis (CMC SWE; Brown and Brasnett 2010). The GlobSnow product is derived from a combination of ground based data and satellite microwave radiometer-based measurements. Due to the nature of the radiometer observations, the SWE product is reliable over areas with seasonal dry snow cover. Areas with sporadic wet snow or a thin snow layer are not reliably detected and typically not present in the SWE product. The CMC SWE dataset consists of Northern Hemisphere snow depth analysis. Snow depth data is obtained from surface synoptic observations, meteorological aviation reports, and special aviation reports. Monthly averages and climatologies of snow depth and estimated SWE are provided, where SWE was estimated using a density look-up table.

4.6 Former Soviet Union hydrological snow surveys

This dataset contains observations made using measuring rods and snow balance at sites throughout the Former Soviet Union between 1966 and 1990. These observations include snow depth and snow water equivalent measured over a nearby snow course transect. The transect snow data used in this study are the spatial average of 100–200 individual measuring points taken on the 10th, 20th, and final day of the month across the typical terrain of the region (Krenke 1998).

For all snow transect data, stations located north of 55°N with at least 75 % of valid data over the 1980–1990 period were selected. Snow data comparison is performed only for those observation stations with differences in altitude less than 150 m with the model representative grid cell. To differentiate between open land and forested areas, a selection is performed on the stations data based on the dominant vegetation type present over the representative model grid cell. By applying these criteria, a total of 351 stations are used, 176 located over open land and 175 over forested areas, mainly located in the Western Part of Siberia. The monthly averages are computed over the 10-year period.

5 Results

5.1 Offline results

5.1.1 Permafrost extent

Figure 2 shows the observed (Brown et al. 1998) and CLASS simulated near-surface permafrost extent North of 45°N for each offline experiment. To be defined as near-surface permafrost, the temperature of at least one soil layer in the top 5.0 m must remain below 0 °C for 24 consecutive months.

Fig. 2

a Observed permafrost extent (continuous, discontinuous, sporadic and isolated) from the International Permafrost Association (IPA) (Brown et al. 1998); b–e Modeled permafrost extent and ALTs for the offline simulations for the 1990–2008 period

The shallow 6-layer configuration, Off_Mine6, greatly underestimates the permafrost extent with 1.28 × 10⁶ km² (Table 2), with ALT below 2.0 m limited to the Canadian Arctic Archipelago and the northernmost regions of Siberia. The model fails to capture the observed extent of permafrost, primarily due to the zero-flux boundary condition at 3.5 m, which leads to the overestimation of the simulated annual cycle of soil temperatures, as shown by Smerdon and Stieglitz (2006), letting all six soil layers thaw regularly in summer.

Table 2 Permafrost extent (×10⁶ km²) located North of 45°N and averaged active layer thickness (m) over 1990–2008 period

Results from the Off_Mine47 experiment show that increasing the soil depth to 65 m substantially increases the near-surface permafrost extent. Comparison between the observed (13.25 × 10⁶ km²) and simulated (12.27 × 10⁶ km²) permafrost extent for continuous and discontinuous permafrost regions shows relatively good agreement, with an average simulated ALT of 2.2 m.

The incorporation of organic soils (Off_OM47) slightly increases the areal extent of continuous and discontinuous permafrost by 0.33 × 10⁶ km² to a total of 12.60 × 10⁶ km². The major impact is a large decrease in the average ALT (1.7 m) over the continuous and discontinuous permafrost regions. Since the ALT mainly represents the summer temperatures, the introduction of soil carbon in CLASS lowers summertime temperatures by insulating more effectively the deeper soil layers from the atmosphere (see Sect. 5.1.3).

Combining SOC with snow conductivity from Sturm et al. (1997) (Off_OMSC47) decreases near-surface permafrost extent compared to the SOC experiment (Off_OM47) (Fig. 2). Permafrost is lost from the southern edge but also in the continuous and discontinuous region, reaching 11.48 million km² compared to 12.6 × 10⁶ km² for Off_OM47 (Table 2). Average ALT over the continuous and discontinuous region also increases to 2.16 m, similar to the 2.21 m obtained for Off_Mine47. As will be demonstrated in the following sections, the decreased snow conductivity tends to warm winter soil temperatures through increased insulation relative to the other experiments.

5.1.2 Active layer thickness

The simulated ALTs are overestimated compared to the observed values at CALM sites, with averaged biases ranging from 39 to 123 cm (Fig. 3). Due to the underestimated representation of near-surface permafrost in Off_Mine6 (Fig. 2), the ALT comparisons with the CALM sites are limited to few sites located in the very cold climate of the Canadian Arctic Archipelago, explaining the smaller bias value of this experiment. The deep mineral, Off_Mine47, and deep organic with modified snow conductivity Off_OMSC47 show large average biases in simulated ALTs, 113 and 123 cm respectively, while the deep organic run Off_OM47 shows relatively smaller average bias of 71 cm. The larger biases of the Off_OMSC47 experiment is explained, as discussed in greater details in following sections, by warmer winter soil temperatures resulting in increased vulnerability to thaw. Careful evaluation of the ALT shows larger overestimation of the ALT for CALM sites located inland, mostly along the Mackenzie River and Alaska where overestimation increases from coastal sites to those inland (not shown).

Fig. 3

Averaged 1990–2008 observed versus modeled annual maximum ALTs for CLASS offline simulations: Off_Mine6 (red); Off_Mine47 (green); Off_OM47 (blue); Off_OMSC47 (black). For all experiments, the number in brackets indicates the number of stations used for comparison

It should be noted that CLASS is run at relatively coarse horizontal resolution (0.5° × 0.5°), while observations are mainly point-scale or representative of a much smaller area (1 km²), and therefore not very representative of large area means, in particular for complex terrain (Dankers et al. 2011; Oelke et al. 2003). Moreover, grid-averaged soil properties and coarse-resolution atmospheric forcings might also introduce biases in the comparison between observed and simulated ALTs. Therefore, comparison between observed and simulated ALTs should be viewed as an indication of the model performance rather than a true in situ validation. Nevertheless, results clearly show a general tendency to overestimate maximum thaw depth for all stations.

5.1.3 Soil temperatures

Figure 4 shows a comparison of the mean annual soil temperatures from the Russian historical temperature dataset at 0.2, 0.85 and 2.75 m corresponding to the midpoints of the 2nd, 4th and 6th model soil layers. All simulations show relatively good agreement for regions with milder climate, while they all underestimate the mean annual soil temperatures at all levels for stations located in colder regions. Nevertheless, the simulation using the deep configuration with SOC and modified snow conductivity (Off_OMSC47) shows smaller biases compared to other experiments, less than −4.5 °C, at all soil levels for stations between 35 and 100°E, while larger biases are noted for grid cells located between 100 and 190°E (Table 3).

Fig. 4

(Left) Simulated and observed mean annual soil temperatures for stations within the region bounded by 20–190°E and 55–90°N. (Right) Mean annual cycle of observed (grey) and simulated soil temperature: Off_Mine6 (red); Off_Mine47 (green); Off_OM47 (blue) and Off_OMSC47 (black). Dashed line on top right plot represents 2 m-air temperature for UDel (grey) and that from ERA-Interim (black) used to drive CLASS offline

Table 3 Annual mean soil temperature biases relative to longitudinal position of the stations

Figure 4 also presents the annual cycle of soil temperatures, which clearly show a strong wintertime cold bias. For all the stations (35–190°E), CLASS underestimates winter temperatures down to −8 °C in February at 0.2 m depth for all runs except Off_OMSC47 which tends to have smaller biases. Surface air temperature used to drive the model shows no significant difference with the gridded dataset from UDel suggesting that differences are most likely due to the model formulation or a dry bias in the input precipitations. Indeed, all experiments show a negative bias in SWE over Central Siberia (not shown) as a result of the underestimation of precipitation in the reanalysis products. The reduced cold bias of the Off_OMSC47 experiment is a result of the reduced snow thermal conductivity from Sturm et al. (1997) allowing a more effective insulation of the soil column by the snow cover, reducing the soil cooling during winter, as discussed in the next section. Summer soil temperature is underestimated compared with observations and the introduction of the SOC further increases the near-surface cold biases. It must be noted that the top organic material is removed from the observational sites, most likely increasing summer soil temperature and amplifying the annual cycle in observations (Gilichinsky et al. 1998).

The shallow configuration Off_Mine6 tends to overestimate the annual cycle amplitude at 2.75 m as a result of the zero-flux bottom boundary conditions imposed at 3.5 m. This confirms the necessity to use deeper soil column configuration especially for high latitudes.

Generally, CLASS captures the annual cycle of the soil temperature reasonably well despite the systematic cold bias in winter. The introduction of SOC does impact the summertime temperatures, while modified snow thermal conductivity plays a major role in the improvement of wintertime soil temperatures.

5.1.4 Impact of soil organic carbon over Northwest Siberia

In this section, we evaluate the changes in the modelled soil temperatures over Northwest Siberia, specifically the region bounded by 55–90°E and 55–75°N, as this region is characterized by deeper permeable soil column (Fig. 1b) and high concentration of organic carbon represented by deep peatlands or by one or two layers of SOC (Fig. 1d). Furthermore, soil temperatures are relatively well simulated for this region compared to observations (Table 3). The baseline experiment for comparison is the Off_Mine47 experiments, allowing to better present the effects of SOC and snow conductivity.

Figure 5 presents the monthly averaged soil temperatures from the surface down to 3.5 m. The shallow mineral experiment Off_Mine6 does not show large differences for the first 0.6 m (three layers) in the annual cycle compared to Off_Mine47, but the amplitude of the annual cycle is overestimated at deeper levels due to the zero-flux bottom boundary condition imposed at 3.5 m, as shown previously in Fig. 4.

Fig. 5

Averaged mean monthly soil temperature for the 1990–2008 period for Northwest Siberia (55–90°N, 60–90°E; blue shaded region in the left central panel), for various CLASS offline simulations and differences relative to the deep mineral configuration (Off_Mine47). No vertical interpolation of the temperature is done therefore results are presented on model levels

The deep SOC experiment, Off_OM47, shows decreased summer soil temperature with reduced penetration of the 0 °C isotherm, thus explaining the decrease in the simulated ALT over the region (Fig. 2). Compared to the mineral formulation, this experiment shows lower summer maximum temperatures in July with differences of −5.4 and −5.6 °C at 0.2 and 0.45 m respectively. This is due to the reduced summertime heat exchanges between the atmosphere and the soil, as shown by reduced ground heat flux (Fig. 6), due to the lower thermal conductivity of the SOC and increased soil water content (addressed in the next section).

Fig. 6

Monthly averages of ground heat flux at the soil/snow (atmosphere) interface for Northwestern Siberia (60–90°E, 55–75°N) for the 1990–2008 period for: Off_Mine6 (red), Off_Mine47 (green), Off_OM47 (blue) and Off_OMSC47 (black)

Similar to the Off_OM47, summer temperatures in Off_OMSC47 are lower compared to the mineral formulation. The temperature differences are smaller with values of −4.23 and −4.25 °C in July at 0.2 and 0.45 m respectively. Wintertime temperatures are significantly higher in this experiment compared to the other deep configuration experiments due to the decreased snow thermal conductivity. In February, the soil temperatures are warmer by +5.03 and +4.25 °C at 0.2 and 0.85 m respectively. The increased winter temperature is caused by reduced heat flux from the soil to the atmosphere through the snow pack (Fig. 6), a direct consequence of the reduced snow thermal conductivity of Sturm et al. (1997) compared to Mellor (1977). The warmer winter soil temperatures in Off_OMSC47 have an influence on the summer temperatures as well. The warmer soil temperatures in spring combined with very similar ground heat flux in spring and summer (Fig. 6) compared to Off_OM47 leads to a deeper penetration of the 0 °C isotherm (Fig. 5), explaining the larger ALT in this particular experiment (Fig. 3). These results demonstrate an important sensitivity of the simulated annual cycle of the soil temperatures to the formulation of the snow thermal conductivity-snow density relation in CLASS.

5.1.5 Changes in hydrology and surface soil variables over Northwest Siberia

The implementation of SOC has a direct effect on the hydrology associated with the increased porosity and soil hydraulic conductivity. This section evaluates the changes in the simulated hydrology over Northwest Siberia. Deep peatlands are present mostly over the southern part of the region while over the northern part, mainly one SOC layer is added (Fig. 1d). The number of permeable soil layers also varies geographically from 3.5 m (six layers) in the southern part to 0.6 m (three layers) in the coastal area (Fig. 1b).

Figure 7 presents the average annual cycle of the main surface hydrological variables for the 1990–2008 period over Northwest Siberia. Soil liquid and frozen water contents increase in the first soil layer in the SOC experiments as a direct consequence of the larger soil porosity of the organic soils (~0.9) compared to mineral soils (<0.49). The first layer water content, liquid and solid combined, is smaller in summer due to increased evapotranspiration.

Fig. 7

Annual cycle of a liquid water content and b frozen water content for the first soil layer, c surface runoff, d sub-surface runoff, e SWE and f evaporation, over Northwestern Siberia (60–90°E, 55–75°N) computed over the 1990–2008 period for CLASS offline simulations: Off_Mine6 (red), Off_Mine47 (green), Off_OM47 (blue) and Off_OMSC47 (black). Numbers on top of the subpanels represent the annual average values for Off_Mine6, Off_Mine47, Off_OM47 and Off_OMSC47 in order. Grey lines snow observational data of SWE for GlobSnow (full line) and CMC SWE (dashed)

The peak in surface runoff occurs in May and June for all experiments and is related to the springtime snowmelt. The simulated differences in the surface runoff and the drainage are subject to many factors: the soil temperatures, the presence of permafrost, the presence of peatlands and changes in evapotranspiration. In CLASS, drainage occurs solely at the bottom of the permeable soil depth, varying between 0.6 and 3.5 m over the region (Fig. 1b). In CLASS, drainage is inhibited by the presence of a frozen soil layer above the depth to bedrock, thus increasing (decreasing) the potential for surface runoff (drainage). In this situation, the soil column behaves similarly to a “bucket” with soil moisture content close to saturation, except for the first two layers where evapotranspiration is active to remove ground water (not shown). For all cases with near-surface permafrost, the degree of saturation generally high, particularly when ALT is above the bedrock, causing runoff when 100 % saturation is reached. Soil moisture loss in this case occurs only through evapotranspiration (not shown). This situation is more present for Off_OM47 where most of the region does not allow drainage due to shallower ALTs (Fig. 2). The presence of shallow frozen ground therefore explains the minimum (maximum) drainage (surface runoff) this particular experiment. The Off_OMSC47 experiment, with milder soil temperatures, earlier thaw and deeper ALT allows a higher fraction of the melt water to percolate into the soil and increases drainage, especially for grid points where only one or two SOC layers are added (not shown).

5.2 Results from CRCM5 experiments

While the offline simulations discussed in the previous section helped assess the sensitivity of soil thermal and moisture regimes to the different modifications in the LSM, it is essential to use coupled land–atmosphere simulations to study the interactions and sensitivity of surface climate to changes in soil and snow formulations.

5.2.1 General model evaluation

The simulation of permafrost and soil conditions in a coupled land–atmosphere model is dependent on the quality of the simulated surface climate. Hence, a general evaluation of the CRCM5 performance against reanalysis and observational databases is presented for important variables influencing the soil: the 2 m-air temperature, precipitation and snow properties. The model evaluation is presented for the deep mineral configuration (C_Mine47).

Figure 8 compares the simulated 2 m-air temperature for winter (DJF) and summer (JJA) with ERA-Interim, UDel and CRU3.1. Summer temperatures are relatively well simulated and biases are less than 4 °C when compared to ERA-Interim over most of the domain. Compared to UDel, the model has a general cold bias over most of North America and over Central and Eastern Siberia. Comparison with CRU3.1 dataset shows different patterns in the temperature differences, with larger areas with positive differences (e.g. Kolyma Mountain Range). The differences in the two station-based climatological datasets, UDel and CRU3.1, can be explained by the low density of the observational network in high-latitudes and the differences in the interpolation methods used. One could see these differences between observed datasets as an approximate measure of the observational uncertainties.

Fig. 8

Comparison of averaged 2 m-air temperature (°C) for DJF (left) and JJA (right) for C_Mine47 (1st row) and differences with ERA-Interim (2nd row), UDel dataset (3rd row) and CRU TS3.1 dataset (4th row), for the 1990–2008 period

Winter temperature shows systematic cold bias over Western Russia, interior Alaska and over the Rocky Mountains. Outside these regions, warm biases are widespread, with biases of 4–6 °C around Hudson Bay and over the Central Siberian Upland. Maximum warm biases are located in the Lena River Valley and in the surrounding mountainous regions (6–12 °C). Although biases are smaller in comparison with ERA-Interim, the existence of similar patterns in the comparison with all datasets, despite their uncertainties, lead to the conclusion that the CRCM5 shows a systematic warm bias over these regions.

Figure 9 shows the simulated winter and summer precipitation fields and its comparison with ERA-Interim, UDel and GPCC datasets. Winter precipitation shows minimum precipitation located over the Arctic Ocean extending to the Canadian Arctic Archipelago and to Siberia, in good agreement with observations. Comparisons with ERA-Interim, UDel and GPCC show general overestimation of winter precipitation over mountainous regions and largely over North America. A dry bias over West Siberian Plains is also present and is consistent in comparison with all datasets. One should note that uncertainties in wintertime precipitation are large mainly because of the undercatch of solid precipitation by ground-based observation stations, a systematic bias for which neither UDel nor GPCC have been corrected. One hypothesis for the underestimation of the precipitation over West Siberia is the underestimation of the Icelandic Low system in CRCM5 (not shown) that would limit the moisture transport from the Nordic Sea region towards Western Siberia.

Fig. 9

Comparison of averaged precipitation (mm month⁻¹) for DJF (left) and JJA (right) for C_Mine47 (1st row) and differences with ERA-Interim (2nd row), UDel dataset (3rd row) and GPCC dataset (4th row), for the 1990–2008 period

Simulated summer precipitation shows the important signature of complex topography over the whole domain with maximum values located over the Pacific Coast of North America and over Siberian Mountains of the Altai, Central Siberian Upland, the Kolyma and Stanovoy Range. This signature of complex topography is hardly visible in station-based observational datasets, most likely due to the sparseness of the stations network, leading to the systematic wet biases in those regions as seen in the difference fields. While the comparison with both station-based datasets show coherent biases (mostly wet) over the high latitudes of the continental North America, comparison with ERA-Interim shows a very different pattern with large dry bias over Alaska decreasing in magnitude towards Hudson Bay. These differences again reflect the large uncertainties in the measurement of precipitation and the representativeness of gridded datasets derived from stations compared to reanalysis products. In that sense, one might consider solely the regions where both station-based datasets and reanalysis biases are coherent as signs of significant biases in the CRCM5. Two regions, the West Siberian Plains and the vicinity of Hudson Bay distinctly appear as areas where the simulated CRCM5 precipitation might suffer from a systematic dry precipitation bias.

Figure 10 compares the average SWE from November to March over 1990–2008 period against the gridded datasets from GlobSnow and from the CMC. SWE data over complex topography is not represented in GlobSnow due to its low reliability (Luojus et al. 2010) and are masked out in Fig. 10. Simulated SWE is underestimated over Eastern Russia compared to GlobSnow data, in agreement with the winter precipitation bias presented in Fig. 9. Simulated SWE over Alaska, Yukon and the Northwest Territories in Canada also show some underestimation although the comparison with precipitation datasets showed a wet bias in the CRCM5. Compared to the CMC SWE, model results show an overestimation of SWE over most of the domain except for an underestimation over the Central Siberian Upland. The discrepancies in the snow data are another illustration of the observational uncertainties leading to difficult model evaluation over the Arctic region.

Fig. 10

Comparison of average SWE (cm) for NDJFM for: a C_Mine47, b GlobSnow, c the difference C_Mine47 − GlobSnow, d CMC SWE analysis, e the difference C_Mine47 − CMC. Mountainous areas are masked in GlobSnow dataset due to insufficient data and high uncertainties; therefore no comparison is made with modeled snow mass over these areas

As a complementary snow analysis, Fig. 11 compares the annual cycle of the snow density averaged for the selected stations form the Former Soviet Union Snow Survey with results form the C_MINE47 experiment over the 1980–1990 period. The CRCM5 shows a shorter snow season compared to observations. A total of 351 stations (175 over forested areas and 176 over open land) were used for model evaluation. On average, the snow density (Fig. 11) appears to be well simulated by the CRCM5 during the October–March period, but significant differences can be noted during the beginning and ending of the snow cover period. Generally, the simulated snow density for all experiments is less variable compared to that observed; this is expected as the observations are in situ, while the simulated density correspond to averages for grid cells that are 50 km × 50 km in size. Separate analysis for forested and open land (not shown) revealed underestimation of the snow density through most of the snow cover period over open land, which might be caused by unrepresented snow processes such as wind compaction and deep hoar development—processes important for the tundra (Schaefer et al. 2009)—in the model.

Fig. 11

Annual cycle of the average snow depth (cm; top panel), snow mass (kg m⁻²; middle panel) and snow density (kg m⁻³; bottom panel) for open land and forested observation stations (grey) and for the C_Mine47 experiment (red). Solid lines represent the average and dashed lines delimits the plus and minus one standard deviation

In general, the CRCM5 shows reasonable skill in reproducing the climatic means over the 1990–2008 period. Albeit some strong biases in winter temperatures over Eastern Siberia, temperature and precipitation can generally be considered in good agreement with the observations and within the range of observational errors. The SWE shows large differences between the observational datasets with larger differences over regions of complex topography. Average snow density is well reproduced by the CRCM5 but shows a systematic underestimation of the snow variability compared to station observations.

5.2.2 Permafrost extent

Figure 12 shows CRCM5 simulated permafrost extent and ALT North of 45°N averaged over the 1990–2008 period. Similar to the offline simulations, the shallow configuration C_Mine6 underestimates the permafrost extent with only 350,000 km², and ALTs below 2.0 m limited to the northernmost latitudes of Siberia and over the Canadian Arctic Archipelago.

Fig. 12

a Observed permafrost extent (continuous, discontinuous, sporadic and isolated) from the International Permafrost Association (IPA) (Brown et al. 1998); b–e Modelled permafrost extent and ALTs for the CRCM5 simulations for the 1990–2008 period

The deep mineral configuration (C_Mine47) shows substantially more near-surface permafrost compared to the shallow configuration with a total coverage of 10.5 × 10⁶ km². The matching offline experiment Off_Mine47 shows larger permafrost extent (12.3 × 10⁶ km²) with the southern permafrost limit extending further south, mostly on the Siberian side. This reduced extent is combined with generally deeper averaged ALT over the continuous and discontinuous permafrost regions in C_Mine47 compared to Off_Mine47, 3.08 versus 2.21 m respectively (Table 2). This increased ALT in CRCM5 is likely related to the winter warm biases compared to ERA-Interim (Fig. 8).

Compared to C_Mine47, the C_OM47 and C_OMSC47 experiments show increased permafrost extent and reduced ALTs over the continuous and discontinuous regions. Summertime changes in soil temperature resulting from the implementation of SOC (figure not shown) shows very similar response to offline simulations (Fig. 5) with maximum cooling near the surface. Comparison of soil temperatures and ALTs for coupled experiments with observed data are similar to those for the offline simulations and are therefore not presented here.

5.2.3 Surface energy balance

Offline experiments showed that the implementation of SOC changes the summer surface energy partitioning by decreasing the ground heat flux (Figs. 5, 6). In the CRCM5, the decreased ground heat flux is compensated by increases in the surface turbulent fluxes. The surface energy partitioning between the latent and sensible heat fluxes can potentially have a large impact on the surface climate in coupled models. Two previous studies addressed the impact of SOC on near-surface climate and showed conflicting results. Lawrence and Slater (2008) showed a large increase in the sensible heat flux over the latent heat flux, causing increased 2 m-air temperatures, and a deeper and dryer atmospheric boundary layer, which decreased the low-level clouds. Rinke et al. (2008) on the other hand showed a large increase in the latent heat flux, a decrease in the 2 m-air temperature and an increase in the low-level cloud cover. Large uncertainties in the atmospheric response in coupled models exists and shows high sensitivity to the implementation technique and parameterization.

Figure 13 presents the differences in the June–July–August latent heat flux (LHF), sensible heat flux (SHF) and 2-m air temperature (T2M), averaged over the 1990–2008 period, for the C_OM47 and C_Mine47 experiments. The noted increase in the energy redistribution towards the atmosphere is a direct consequence of the decreased ground heat flux (no figure), similar to results from offline experiments (Fig. 6). The maximum increases in the LHF are located in coastal regions while the increase in the SHF are mainly located inland, mostly over Central and Eastern Siberia. Most regions where differences exist in both fluxes are statistically significant at a 95 % confidence level (Fig. 13). T2M changes show a general decrease (increase) over regions where the latent (sensible) heat flux increases dominates, in good agreement with the results from both Rinke et al. (2008) and Lawrence et al. (2008) but with higher spatial heterogeneity in the signals. Despite the significant changes in the surface turbulent fluxes, the statistically significant temperature changes are mostly present outside of the continuous permafrost region (Fig. 12).

Fig. 13

(Top) Comparison of summer (June–July–August) climatology of (left) sensible heat flux, (middle) latent heat flux and (right) 2 m-air temperature between C_OM47 and C_Mine47 over the 1990–2008 period. (Bottom) Statistical significance of the differences using a t test

Figure 14 presents the mean 1990–2008 annual cycle of the LHF and SHF along with the T2M over two main vegetation categories present over the Arctic region; forested areas (needleleaf trees) and tundra (grasslands). Besides the differences in vegetation properties, another important difference between the two categories is their geographic distribution and therefore climate, with forested areas with milder climate mostly present in the southern part of the domain, and tundra with colder climate covering the higher latitudes and coastal regions of the Arctic Ocean (Fig. 14c).

Fig. 14

Average annual cycle of a sensible heat flux, b latent heat flux (W m⁻²) and c 2 m-air temperature (°C) for needle leaf trees (full lines) and grasslands (dashed lines) between 55–75ºN and 45–270ºE for the 1990–2008 period for simulations Off_Mine6 (red), Off_Mine47 (green), Off_OM47 (blue) and Off_OMSC47 (black). d Annual cycle of soil saturation for the first two soil layers: 0–10 cm (full lines) 10–30 cm (dashed lines) for needle leaf trees. Identical color codes are used to designate experiments

For both vegetation classes, the surface energy fluxes (Fig. 14) are separated in two distinct groups, with generally larger fluxes noted for the two experiments using the SOC parameterization compared to the two experiments using mineral soil. This separation is simply the result of the decreased ground heat flux in the summer for the experiments using SOC, increasing the surface energy redistribution towards the atmosphere.

For the forested areas, the SHF and LHF show differences in their annual cycle. The LHF (Fig. 14a) shows larger increases during early summer for the two SOC experiments. Differences in LHF are maximum in June, rapidly decreasing to reach similar values to the mineral experiments in August. The SHF, on the other hand, is similar between all four experiments until May. In June and July, the soil organic experiments show large increase in the SHF (Fig. 14b). Over tundra, the LHF increase is generally larger compared to the changes in the SHF in good agreement with spatial distribution seen in Fig. 13. The differences between the SOC and the mineral experiments are maximum in June by up to 9.4 W m⁻² while the LHF reaches maximum values in July. The SHF increases are generally smaller, as can be noted on Fig. 13 limited to a maximal 6.8 W m⁻² in June.

To understand the surface energy partitioning, one must study the availability of soil water for evapotranspiration. Firstly, CLASS allows transpiration to occur for air temperatures above 0 °C, provided soil liquid water is available to the roots in the soil. This explains the delay in the increase of the LHF component over grasslands compared to the forested areas, directly related to the colder surface climate. Secondly, results for offline experiments (and similarly for CRCM5) showed increased soil water content in the first soil layer in experiments using the SOC parameterization (Fig. 7), largely caused by the increased soil porosity. Figure 14d presents the soil saturation for the upper two layers for the needleleaf vegetation category. For the organic soil experiments, the soil saturation reaches maximum value in May for the first layer, rapidly decreasing to reach values near the retention capacity of the organic soils (0.27) in July and August. Therefore, water available for evapotranspiration is higher in May and June compared to July explaining the larger increase in the LHF over that period while SHF shows the largest differences in July, when the surface soil layer is drier.

The net effect of the introduction of soil organic carbon is a drying of the first soil layer combined with an increase in the soil saturation in deeper layers, where the depth to bedrock allows more than one permeable layer. The downward displacement of soil water is caused by the larger hydraulic conductivity and low suction of the organic soils, a result also obtained by Lawrence and Slater (2008). Furthermore, in the organic soils parameterization, if the liquid water content of the first soil layer is above the retention capacity, CLASS creates a vertical gradient within this layer imposing the retention capacity at the surface increasing to reach saturation at the bottom of the soil layer (Verseghy 2008). This formulation therefore presents a relatively dry surface to the atmosphere, limiting direct evaporation from the top soil layer, explaining the combination of increased LHF and SHF.

The CRCM5 shows moderate atmospheric response to the increase in the summer surface turbulent heat fluxes (Fig. 14c). The summer 2 m-air temperature changes are limited to maximal values within ±2 °C over both vegetation categories. Spatially, changes in 2 m-air temperatures shows summertime cooling over North America and Western Russia reaching changes of −0.5 to −1.5 °C, while Eastern Russia warms between 0.5 and 1 °C (Fig. 13).

In summary, CRCM5 results show a more moderate response of the surface climate to the implementation of SOC compared to the studies of Lawrence and Slater (2008) and Rinke et al. (2008). Increases in LHF and SHF are determined by the availability of water in the near-surface layers explaining the seasonality of the changes in the surface turbulent fluxes, especially over the forested areas. The limited influence of the SOC implementation on the surface and atmospheric boundary layer climate is likely caused by the similar increases in both fluxes, having little impact on the boundary layer stability, illustrating the atmospheric models sensitivity to the implementation of SOC.

5.2.4 Sensitivity of simulated ALT to atmospheric parameters

In this section, we evaluate the simulated ALT trends for the 1960–2008 period and its sensitivity to atmospheric variables. Figure 15 presents the ALT trends for the three simulations using the deep soil configuration for those points that retain near-surface permafrost until the end of 2008. The experiment using mineral soils, C_Mine47, has the largest trend values over most of the domain, compared to experiments using SOC. One might note signs of permafrost degradation over the simulation period for areas near the southern limit of the permafrost region (Fig. 15). The extent of this area is larger in C_Mine47 compared to the experiments using SOC. For grid points where the ALT trend is statistically significant at 90 % confidence level, the average ALT trend in C_Mine47 is 12.3 cm/decade, a larger value compared to the ~8 cm/decade obtained by Oelke et al. (2004). Maximum values are mainly located over Northwest and Eastern Siberia, regions where the average ALT is overestimated compared to other modelling experiments using offline LSMs (Oelke et al. 2004; Burke et al. 2013). The experiment C_OM47 shows smaller trends in ALT amongst all experiments, 6.7 cm/decade, a consequence of the colder winter soil temperatures and the effective insulation of the SOC in summer. Nevertheless, the trend is statistically significant over most of the continuous and discontinuous permafrost regions except for Central Siberia, where all experiments show non-significant trends, in good agreement with Oelke et al. (2004) and Burke et al. (2013). The C_OMSC47 experiment shows significant trends over most of the domain, with maximum trends located in Southern and Eastern Siberia and over North America from Alaska to the Hudson Bay with an average trend of 9.9 cm/decade over grid points where the trend is significant at 90 % confidence level.

Fig. 15

(Top) Trends in ALT for the 1960–2008 period for C_Mine47 (left), C_OM47 (middle) and C_OMSC47 (right). Grey regions represent grid cells where permafrost is not present for the entire simulation period. (Bottom) Statistical significance of the trends defined using the Mann–Kendall test

To understand the sensitivity of the ALT to atmospheric variables, an analysis of the relation between ALT and atmospheric variables directly influencing the soil thermal regime was performed over areas of significant ATL trends. To identify the most important atmospheric variables having an impact on the ALT, a trend analysis and correlation between atmospheric parameters and ALT were investigated. Statistically significant correlations at 95 % confidence level are found between the ALT and the following atmospheric parameters: 2 m-air temperatures (T2M), the degree-day thawing (DDT) and freezing (DDF) indexes, and the length of both thawing and freezing season. The DDT (freezing) index is defined as the sum of the above-zero (sub-zero) daily average 2 m-air temperatures from October to September. Other parameters, such as SWE, net shortwave radiation at the surface, annual and seasonal precipitations show no significant correlation with the ALT.

Figure 16 shows the relation between the ALT anomaly and selected atmospheric parameters: T2M, DDT and DDF. Over the 1960–2008 period, the average annual 2 m-air temperature warms by 0.34, 0.29 and 0.27 °C/decade for experiments C_Mine47, C_OM47 and C_OMSC47 respectively. The differences in the magnitude of the trends between the experiments are due to the different spatial distribution of the grid points where the ALT trends are significant (Fig. 14b). This increase in temperature is significantly correlated with increases in the ALTs (Fig. 16) with values of 0.68, 0.85 and 0.80 for C_Mine47, C_OM47 and C_OMSC47 respectively.

Fig. 16

Relation between ALT departure and average annual 2 m-air temperature (top), degree-day thawing index (middle) and degree-day freezing index (bottom) for experiments: C_Mine47 (green), C_OM47 (blue) and C_OMSC47 (black) over the respective regions where ALT trends are significant at a 90 % confidence level

Since the ALT—defined as the annual maximal thaw depth—occurs at the end of the summer, the thawing index should be a better indicator than the annual average 2 m-air temperatures. The DDT shows an increase over the simulated period with values of 28.8, 30.1 and 24.2 °C day/decade for C_Mine47, C_OM47 and C_OMSC47 respectively, coherent with a warming of the surface climate over the Arctic. The ALT and the DDT in C_Mine47 experiment shows the highest correlation amongst the three experiments (0.91). The other simulations using SOC are also highly correlated with the DDT but with smaller values: 0.81 and 0.79 for C_OM47 and C_OMSC47, respectively.

The DDF shows larger trends compared to DDT with values of −97.2, −76.1 and −68.1 °C day/decade. The DDF and the ALT are negatively correlated, significant at 95 % confidence level, with values of −0.44, −0.74 and −0.64 for C_Mine47, C_OM47 and C_OMSC47 respectively.

Although all experiments agree on the relations between a deepening of the ALT with the warming of the surface climate over the 1960–2008 period, some differences related to the soil configuration are noted. Firstly, the mineral soil experiment, C_Mine47 shows maximal correlation of the ALT with the DDT index amongst the experiments, while minimal correlations are found with T2M and DDF for that particular experiment. This shows that the mineral configuration is most sensitive to the summer temperatures. According to the larger ground heat flux in the mineral experiment compared to the SOC experiments (Fig. 6), results presented in this section suggest a short-term sensitivity (or “memory”) of the soil to the air temperature anomalies of the same summer with limited sensitivity to the previous winter or summer temperatures. On the other hand, the SOC experiment C_OM47 shows higher correlation with the annual temperatures, showing larger sensitivity of this particular experiment to both the summer and winter temperatures. This is partly due to the isolative effect of the SOC reducing the impact of summer temperatures while the larger snow conductivity, compared to C_OMSC47, shows relatively important impact of the winter air temperature changes on the soil column temperature. One might note that the high correlations between the temperature-based atmospheric parameters and the ALT on a yearly basis shows similarities to the observational study of Frauenfeld and Zhang (2011) over Russian stations, where antecedent conditions from the previous year do not appear to play a major role in affecting the subsequent near-surface soil conditions.

Results from the CRCM5 experiments showed overestimated ALT and trends over most of the Pan-Arctic domain compared to observations (CALM) and other numerical experiments using LSM (Oelke et al. 2004; Burke et al. 2013). Two hypotheses are formulated to explain the overestimated ALT and trends: (1) the warm bias in the CRCM5 surface climate (Fig. 8) and (2) an accelerated warming in the CRCM5 over the simulated period compared to observations. Figure 16 compares the time series of the DDF (Fig. 17a) and DDT (Fig. 17b) indexes from UDel, CRU3.1 to simulated results for areas where the three deep soil configuration experiments shows significant trends in the ALT, mostly over Eastern Siberia and Alaska (Fig. 15). Due to the limited observations of daily temperatures in the Arctic, estimates of the DDF and DDT indexes were computed from the monthly temperatures from UDel and CRU3.1. This method is known to introduce some uncertainties in the estimates of the indexes between 10 and 15 % over the continuous and discontinuous permafrost regions (Frauenfeld et al. 2007). Compared to CRU3.1 data, the simulated DDF is underestimated while in relatively good agreement with the UDel data. Inversely, simulated DDT is overestimated compared to CRU3.1 while underestimated compared to UDel. In general, the CRCM5 experiments are between CRU3.1 and UDel estimates, showing that CRCM5 is within the observational uncertainties but tends to be biased warm, in good agreement with the simulated biases noted over Eastern Siberia and North America (Fig. 8). Statistically significant differences at 95 % confidence level exist between the simulated and observational estimates, which highlight the CRCM5 warm bias over the region, especially when compared to CRU3.1, noted to be colder than UDel, especially over Eastern Siberia (Fig. 8).

Fig. 17

(Left) Time series of degree-day freeze (a) and degree-day thaw (b) over grid points where all experiments shows statistically significant ALT trends at 90 % confidence level. Experiments C_Mine47 (green), C_OM47 (blue), C_OMSC47 (black) and observations from UDel (cyan) and CRU (magenta) are presented. Normalized time series and statistically significant linear trends are presented for DDF (c) and DDT (d)

Standardizing the time series of DDF (Fig. 17c) and DDT (Fig. 17d) with respect to their long-term mean and standard deviation, removing any biases, allows estimating and comparing the trends. Standardized trends and interannual variability for DDF and DDT are almost identical for CRU3.1 and UDel. The CRCM5 experiments tend to overestimate both trends and the inter-annual variability of the indices. A t test was performed over the standardized data that showed that the observed and simulated distributions are not significantly different at 95 % confidence level. Therefore we can conclude that the warming trends are well reproduced by the CRCM5 although the warming is slightly accelerated in the model.

In summary, the simulated ALTs in CRCM5 experiments are overestimated compared to observations and other studies using offline LSMs and significantly correlated to surface air temperatures and freezing and thawing indexes. The snow and precipitation show little correlation with ALTs and likely play a secondary role in the evolution of the soil temperatures. Despite the warm bias noted in the CRCM5 surface air temperatures, the warming trends over the 1960–2008 period are not significantly different from the observations, suggesting that the overestimated ALTs are likely related to that warm bias.

6 Summary and conclusions

The warming observed in the Arctic in the recent decades and the warming suggested by coupled general circulation models for the twenty-first century could have an important impact on the thermal state of the Arctic soils. Although direct observations are invaluable sources for monitoring the Arctic Climate, the limitations in the available data both spatially and in time strongly supports the modelling approach to get continuous data for large-scale studies of the Arctic. Despite an increasing number of studies addressing the changes in the thermal state of the Arctic permafrost, the LSMs generally require improvements to their numerical formulation to adequately represent the soil organic carbon and its thermal and hydrological properties. Moreover, a limited number of studies actually focus on the land–atmosphere interactions resulting from the implementation of SOC and its potential feedbacks on the surface climate. These two different aspects are addressed in this study.

The first objective of this study was to understand the sensitivity of simulated Arctic soil temperature and moisture regimes, particularly for near-surface permafrost and ALT, to soil layer configuration and soil organic carbon, using offline simulations with CLASS. In agreement with the work of Smerdon and Stieglitz (2006), the shallow soil configuration showed overestimated annual cycle of soil temperatures, leading to a large underestimation of the permafrost extent in the Arctic. Other modeling groups showed better permafrost extent using shallow (~3.5 m) configuration, e.g. Lawrence et al. (2008) and Dankers et al. (2011), using the Community Land Model (CLM) and the Joint UK Land Environment Simulator (JULES), respectively. Though the above studies had similar total soil column depths to that of the shallow configuration considered here, there are important differences with respect to the depth to bedrock. In this study, the depth to bedrock is based on the dataset from Webb et al. (2000), while Dankers et al. (2011) and Lawrence et al. (2008) assumed permeable soil composed of sand, silt, clay and organic material for the entire simulated soil column. Since a large part of the study domain has less than 1.0 m of permeable soil according to Webb et al. (2000), the thermal conductivity is larger in our experiments due to the presence of shallow bedrock and reduced thermal inertia due to the limited soil water content and phase changes in the soil column. The combination of these effects most likely leads to larger thaw depths compared to a fully permeable soil column.

Deepening the soil column depth to 65 m improved the permafrost extent although the ALT was overestimated compared to observations and to other modeling studies (Oelke et al. 2003, 2004; Lawrence et al. 2008; Dankers et al. 2011; Burke et al. 2013). The implementation of SOC greatly reduced the ALTs and summer soil temperatures over the Pan-Arctic domain compared to the experiment using only mineral soils. By combining the SOC and the decreased snow thermal conductivity of Sturm et al. (1997), soil temperatures show colder summertime temperatures due to the SOC, and warmer wintertime temperatures as a consequence of the reduced heat flux from the soil towards the atmosphere through the snow pack. The LSM CLASS therefore showed an important sensitivity of the soil temperatures to the snow density-snow thermal conductivity relation that should be investigated further.

The hydrological response to the SOC implementation revealed complex interactions between the surface runoff and drainage to the permeable soil depth, presence of deep peatlands and near-surface permafrost. In general, simulation with colder temperatures and shallower ALTs showed increased (decreased) surface runoff (drainage), particularly over areas of high soil carbon concentrations and peatlands.

CRCM5 experiments showed similar sensitivity to soil column depth and SOC implementation as the offline experiments. Simulated near-surface permafrost extent showed similar response although the CRCM5 deep mineral configuration showed reduced extent and larger ALT values compared to its offline counterpart. These warmer soil temperatures are caused by the warm bias in the CRCM5 climatology. As for offline experiments, significant decreases in simulated ALTs result from the implementation of SOC.

The summer surface turbulent fluxes in CRCM5 show important increases caused by the SOC implementation, as a direct response to the decreased ground heat flux and energy redistribution towards the atmosphere. While Lawrence and Slater (2008) and Rinke et al. (2008) obtained opposite and significant signals in the surface turbulent heat fluxes, results from this study show moderate response with very limited impact on the surface climate and atmospheric boundary layer. This is caused by the similar increases in both sensible and latent heat fluxes in the SOC experiments, controlled by the availability of near-surface soil moisture and saturation level. Only limited regions of the Pan-Arctic domain show significant changes to the 2 m-air temperature. These results emphasized the large sensitivity of atmospheric models to SOC formulation.

The warming of the near-surface climate in the CRCM5 experiments over the 1960–2008 period is responsible for the positive and significant trends in the simulated ALT. The mineral soil experiment showed the largest trends and permafrost degradation amongst the deep soil experiments. The large sensitivity of that particular experiment to the increasing DDT index and its high correlation with the ALT suggests that changing summer temperatures dominate the ALT sensitivity. Although the experiments using SOC also showed significant correlations between the ALT trends and DDT index, increased correlations were also found with the annually averaged 2 m-air temperature and the DDT index. These higher correlations show that the SOC experiments are less sensitive to the summer temperatures and that increasing temperatures during winter also plays an important role.

Analysis of the trends in DDT and DDF showed that the CRCM5 surface climate trends are not significantly different from that of CRU3.1 and UDel. The overestimated ALTs and trends in the model are then most likely linked to the near-surface warm biases in the model.

Despite the presence of shortcomings in CLASS formulation and biases in the CRCM5 surface climatology, substantial improvements in the representation of soil thermal and moisture regimes, hence near-surface permafrost, were achieved in this study by using a deeper soil column, implementing soil organic carbon and by the modification of the snow thermal conductivity—snow density relation. Based on the present study, it appears important to improve the representation of snow in CLASS. For example, Burke et al. (2013) showed that a more realistic multi-layer snow scheme significantly reduced the winter cold bias in soil temperatures previously observed in JULES (Dankers et al. 2011). Although CLASS snow model reproduces the average snow conditions over the Arctic, it also underestimates the spatial variability of the snow variables. The snow density, especially over the tundra, could be improved by adding some snow processes actually not represented in CLASS, such as the effect of wind compaction and deep hoar development. The study by Schaefer et al. (2009) showed clear improvements in their snow representation, with largest impacts on the soil temperatures over the tundra, by adding snow classes and associated processes to their multi-layer snow model.