Snow-Covered Soil Temperature Retrieval in Canadian Arctic Permafrost Areas, Using a Land Surface Scheme Informed with Satellite Remote Sensing Data

Abstract

High-latitude areas are very sensitive to global warming, which has significant impacts on soil temperatures and associated processes governing permafrost evolution. This study aims to improve first-layer soil temperature retrievals during winter. This key surface state variable is strongly affected by snow’s geophysical properties and their associated uncertainties (e.g., thermal conductivity) in land surface climate models. We used infrared MODIS land-surface temperatures (LST) and Advanced Microwave Scanning Radiometer for EOS (AMSR-E) brightness temperatures (Tb) at 10.7 and 18.7 GHz to constrain the Canadian Land Surface Scheme (CLASS), driven by meteorological reanalysis data and coupled with a simple radiative transfer model. The Tb polarization ratio (horizontal/vertical) at 10.7 GHz was selected to improve snowpack density, which is linked to the thermal conductivity representation in the model. Referencing meteorological station soil temperature measurements, we validated the approach at four different sites in the North American tundra over a period of up to 8 years. Results show that the proposed method improves simulations of the soil temperature under snow (Tg) by 64% when using remote sensing (RS) data to constrain the model, compared to model outputs without satellite data information. The root mean square error (RMSE) between measured and simulated Tg under the snow ranges from 1.8 to 3.5 K when using RS data. Improved temporal monitoring of the soil thermal state, along with changes in snow properties, will improve our understanding of the various processes governing soil biological, hydrological, and permafrost evolution.

1. Introduction

The thermal insulation capacity of snow cover, which is well-known, limits winter cooling of underlying soil in Northern regions, thereby influencing the permafrost thermal regime [1]. Permafrost now covers 25% of the landmass of the Northern Hemisphere, and recent studies have shown that about 50% of near-surface permafrost (the top few meters) could thaw by 2050, and 90% by 2100. Sustained by constant warming, this thaw could produce a very significant positive feedback effect (see [2]). Since 70% of northern land areas are covered by snow for more than six months of the year [3], the state of snowpack evolution is one of the key variables governing the evolution of the active layer of permafrost [4,5,6]. Current climate models usually have simple representations of snow (i.e., simplified snow microstructure and density profiles) that lead to erroneous snowpack thermal conductivity [5,7,8,9,10,11]. Moreover, the significant observed shrub expansion in the arctic [12] also modifies snowpack state, which in turn could alter heat transfer mechanisms [6,8,13]. Snowpack thermal conductivity changes under climate warming will affect soil temperatures and permafrost active layer evolution processes, as well as energy and biogeochemical fluxes between land and the atmosphere [14].

Several approaches have been put forward to retrieve soil surface temperatures directly from microwave radiometry, such as those proposed by Holmes et al. [15], Royer and Poirier [16] and Prigent and Rossow [17]; however, none can retrieve soil temperatures under snow. Model-based estimates that include a snow model to better differentiate land surface temperatures (LST) from soil temperatures (Tg) are thus needed to improve retrieval accuracy. Wang et al. [18] reviewed different land surface models (LSM) to retrieve Tg, but they are highly dependent on meteorological driver data, as outlined by Brun et al. [19]. Holmes et al. [20] used microwave observations to retrieve Tg in snow-free conditions only, while Langer et al. [21] proposed satellite-based modeling, combining thermal observations from MODIS and snow data (Globsnow, [22]) to retrieve soil temperature. However, lingering uncertainties remain in the Globsnow dataset, where significant biases were observed in the North (see the discussion in [23]). These important biases in snow mass prediction induce inaccuracies in Tg retrievals. Kohn and Royer [24] previously advanced an approach that uses a complex multilayer physical snow model constrained by microwave brightness temperatures, but which is computationally heavy to implement over large areas.

In this paper, we explore the potential of using satellite-based remote sensing (RS) data to improve the seasonal evolution of the thermal conductivity of the snowpack (snow insulation) in a simple 1-D offline land-surface scheme (the Canadian Land Surface Scheme, or CLASS) [25]. The aim is to improve the soil temperature (Tg) retrievals under snow. We hypothesize that the main uncertainty in the Tg estimates arises from erroneous snowpack microstructure representation in the model, rather than from the soil conductivity parameterization.

2. Materials and Methods

2.1. The Canadian Land Surface Scheme-Specific Surface Area Model

The Canadian Land Surface Scheme (CLASS) 1-D Version 3.5 [25] in an off-line column version was driven by reanalysis of meteorological data [26]. CLASS is used operationally in the Canadian global circulation model [27] and the Canadian Regional Climate Model (CRCM, [28]). CLASS computes the surface energy budget as a function of soil and land-cover types. The standard operational configuration for CLASS consists of three soil layers of 0.10 m, 0.25 m, and 3.75 m in thickness. However, in this study, to allow variations of the deeper layer (3.75 m) and to better account for the effect of phase changes in the soil, the soil was stratified over 47 layers (47 calculation nodes) from 0.1 m to 1.5 m thick, over a total depth of 65 m [29]. In CLASS, the top soil layer is fixed at 0.1 m, which is thicker than the effective soil layer for microwave emission depth [30]. Nevertheless, we assume a correspondence between both layers (see Section 5). CLASS includes prognostic equations for energy and water conservation for the 47 soil layers, while a thermally and hydrologically distinct snowpack is treated as an additional variable-depth layer (see below). In order to represent subgrid-scale variability in a simple manner, CLASS adopts a “pseudo-mosaic” approach and divides the land fraction of a given grid cell into four sub-areas: bare soil, vegetation, snow over bare soil, and snow with vegetation. In this study, for all considered arctic sites, we only used homogeneous grid cells with vegetation corresponding to the grass plant functional type, defined by its structural attributes, including albedo, leaf area index (LAI), roughness length, canopy mass, and rooting depth.

CLASS includes a one-layer snow model described in detail by Bartlett et al. [31] and Brown et al. [32]. The thermal conductivity of snow was calculated from snow density using the empirical relationship described in [33]. However, the CLASS snow model does not simulate the snow grain metamorphism that is needed for the assimilation of passive microwave brightness temperature (Tb). This motivated the implementation of snow grain metamorphism within CLASS by adding an offline snow specific surface area (SSA) model [34]. This offline multilayer model (CLASS-SSA) simulates the decrease of SSA based on snow age, snow temperature, and temperature gradient under dry snow conditions [35], while considering the liquid water content of the snowpack for wet snow metamorphism [36].

The meteorological forcing data were extracted from the National Centers for Environmental Prediction (NCEP) North American Regional Reanalysis (NARR) data [37], and consist of total precipitation mass occurring between two time steps (3 h; kg m⁻²), 2 m air temperature (°C or K) and humidity (%), 10 m wind speed (m s⁻¹), and incoming surface shortwave (SWd) (W m⁻²) and longwave radiation (LWd) (W m⁻²). Meteorological data from the NARR nearest neighbor pixel of each site was used to drive the model. In this study, CLASS-SSA was run at a time step of 30 min, and since NARR provides data for a three-hour time step, all variables were interpolated to a 30-min time step, except precipitation, which was kept on three-hour intervals. A one-year spin-up period was repeated three times, providing three years of spin-up in total, to initialize the model.

2.2. The Helsinki University of Technology Microwave Radiative Transfer Model

Tb was simulated using the one-layer Helsinki University of Technology (HUT) radiative transfer model [38]. The attenuation coefficient of snow was calculated using an empirical equation dependent on the frequency and the effective grain size of the snowpack. The radiative transfer of this simplified one-layer model is based on a two-flux approximation. CLASS-SSA outputs (snow temperature and density and snow grain effective radius derived from simulated SSA and soil temperature) were used as inputs to HUT. The optical radius, calculated from SSA values simulated by CLASS-SSA, was computed to effective grain size (required snow grain input into HUT) by applying a correction factor of 3.7, put forward by Roy et al. [34]. In this study, we only used Tb at 11 and 19 GHz frequencies, given their lower sensitivity to snow grain size scattering (see [39]). Thus, the simplified approach based on CLASS-SSA–HUT models has a small impact on the calculated Tb [34]. When the simulated CLASS soil temperature was below zero, we used soil permittivity values of frozen soil [40], otherwise the soil Dobson permittivities were used [38]. Hereafter, the combined CLASS-SSA models are referred to as CLASS for simplicity.

2.3. Reference Study Sites

Four arctic sites where continuous in-situ Tg measurements are available across Canada and the United States were chosen as reference: North Slope, Alaska; Inuvik, Northwest Territories (NWT); Daring Lake, NWT; and Salluit, Nunavik, Québec (Figure 1 and

Table 1. Details of study sites (see Figure 1): mean temperatures (T) (summer: June–August; winter: December–February), total precipitation, and snow water equivalent (SWE).

Sites	Lat., Long., Alt.	Mean Annual T (°C)	Mean Summer T (°C)	Mean Winter T (°C)	Annual Precip. (mm)	Max. SWE (mm)	Period of Analysis
North Slope	70.27N, 148.88W, 15 m	−11	14	−21	330	84.3	October 2007
Inuvik	68.43N, 133.33W, 81 m	−7.2	10	−27	280	115.2	November 2003
Daring Lake	64.52N, 111.54W, 459 m	−6	8	−29	300	136.1	August 2007
Salluit	61.87N, 75.23W, 312 m	−7.8	6	−22	300	236.2	October 2006

). The meteorological characteristics of these sites (provided by Environment and Climate Change Canada) are listed in Table 1. Mean maximum snow water equivalent (SWE) values were derived from the Liston and Hiemstra [3] reconstruction for 1979–2008, showing a significant low-to-high longitudinal gradient from west to east. The percentages of water bodies within the microwave 25 × 25 km² satellite pixel chosen for these sites are (1) North Slope: 2%; (2) Inuvik: 8%: (3) Daring Lake: 10%; and (4) Salluit: 4%. Given its relatively low coverage, the water fraction is assumed to have limited impact on soil temperature retrieval.

2.4. Satellite Data

We used two types of satellite data to constrain the retrieval in this study: microwave Tb from AMSR-E and the thermal infrared LST from AQUA and TERRA MODIS data. The Tb data at 10.7 GHz (“11 GHz”) and 18.7 GHz (“19 GHz”) at the horizontal (H) and vertical (V) polarizations were extracted from the National Snow and Ice Data Center (NSIDC) EASE-Grid 25 × 25 km² product [42]. Both ascendant (13:30) and descendent (01:30) orbits were used. The atmospheric contributions to the satellite Tb were removed by using the millimeter-wave propagation model [43] implemented in the HUT model [38]. The atmospheric model was driven by air temperature and the moisture of the atmospheric layers above the surface from the 29 North American regional reanalysis [37] atmospheric layers.

The MODIS land-surface temperature (LST) data were taken from the NASA 1-km resolution—the daily L3-V6 product (MODIS-Terra MOD11A1 and MODIS-Aqua MYD11A1 Land Surface Temperature products, [44]). This LST product was retrieved from the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC) at the United States Geological Survey (USGS)/Earth Resources Observation and Science (EROS) Center (Sioux Falls, South Dakota, United States). The data is corrected for atmospheric effects and for surface emissivity variation. We resampled the LST data by averaging the LST values at the EASE-Grid 25 km resolution.

3. Inversion Approach

The basic approach used to retrieve Tg was to constrain CLASS simulations with coincident satellite observations. The different temperature variables considered in this study are defined in Figure 2: land surface (skin) temperature (Ts), brightness temperature at 11 and 19 GHz (Tb), and soil temperature (Tg(i) at each layer i). The successive steps for the retrieval procedure are described below (Section 3.1). For the summer period (without snow), we used satellite data to scale CLASS air temperature at 2 m (T_air) and soil moisture simulations by adjusting respectively the summer NARR T_air and precipitation (drivers). For the winter period, we used satellite data to scale CLASS T_air and snow depth by adjusting the winter NARR T_air and precipitation (drivers), respectively. The snow cover density simulations were then adjusted against the Tb polarization ratio at 11 GHz (see Section 3.2).

3.1. Adjustment of Meteorological Driving Data

North American Regional Reanalysis (NARR) meteorological data suffers spatial bias that must be corrected for in local analyses (see [45]). Here, we only adjusted NARR 2-m T_air and the amount of precipitation. Both were corrected by minimizing the root mean square error (RMSE) between the CLASS simulations and the measured satellite-based Ts and Tb. For the summer period, Ts is the surface soil temperature, while for winter, Ts is the snow surface temperature. Also, for summer, Tb is used for fitting the soil moisture, while for winter, it is used for fitting the snow depth—both parameters are governed by precipitation, with T_air > 0 or < 0 °C, respectively.

For the T_air adjustment, we minimized the simulated CLASS Ts against a coincident-measured MODIS LST at four passes (TERRA ascendant (descendant) orbit at 22:30 (10:30) local time; and Aqua ascendant (descendant) orbit at 13:30 (01:30) local time). We adjusted the T_air(h) at each hour (h) using a correcting factor (β_air) such as

$${\mathrm{Tair}}_{\mathrm{corr}}\left(\mathrm{h}\right)={\mathrm{Tair}}_{\max }-{\mathsf{\beta }}_{\mathrm{air}}·\left(1-\cos \frac{\mathsf{\pi }·\left(\mathrm{h}-{\mathrm{t}}_{\max }\right)}{\mathrm{T}}\right)·\left({\mathrm{Tair}}_{\max }-{\mathrm{Tair}}_{\min }\right)$$

(1)

where Tair_corr(h), Tair_max and Tair_min are respectively the corrected, maximum and minimum NARR air temperatures when T = 24 h, and t_max is the time of solar maximum. The factor β_air was optimized by minimizing the RMSE between measured MODIS LST and simulated Ts. Note that the assumption in Equation (1) of sinusoidal behavior about local solar noon could slightly be biased by a possible shift in the daily Tair_max from solar noon by advection, which sometimes occurs at high latitudes and in winter when solar forcing is weak [46]. However, because Tair_corr(h) is then corrected against MODIS observations, such a bias does not significantly impact the results.

Moreover, as the down-welling longwave radiation (LWd) is linked to T_air, we corrected the NARR LWd with the adjusted NARR T_air using a simple linear regression, since LWd_NARR and T_airNARR are well correlated, with R² from 0.5 in summer and 0.82 in winter (not shown). The equation of correction was defined by:

$${\mathrm{LWd}}_{\mathrm{corr}}={\mathrm{LWd}}_{\mathrm{NARR}}-\mathrm{p}\Delta \mathrm{Tair}$$

(2)

where LWd and LWd_NARR are respectively the down-welling longwave radiation corrected and from NARR data; p is the correction factor (i.e., the slope of the regression between LWd_NARR and T_airNARR); and ΔT_air is the difference between T_aircorr and T_airNARR.

For the precipitation amount adjustment, Tb values were simulated at 11 and 19 GHz, using the CLASS outputs as inputs to HUT. For the summer period (T_air > 0 °C), we hypothesize that measured Tb at 11 and 19 GHz are sensitive to the surface soil moisture (SM) [47,48,49], which we assumed to be mainly linked to precipitation. The sensitivity of Tb to SM increases as frequency decreases. The most sensitive AMSR-E bands are thus the X (11 GHz) and Ku (19 GHz) bands. We did not use the C (6 GHz) band, which is contaminated by radio frequency interference, nor the Ka (37 GHz) band, which is known to be less sensitive to SM (higher frequency) and also noisier, because of the vegetation above the soil (see the review [47]). Note that the dual-polarization measurements (Tb H/V), sometimes used for SM products [47], would have probably improved the sensitivity to SM, since the polarization ratio attenuates the physical temperature effects. Figure 3 illustrates the sensitivity of Tb with V of 10.7 and 18.7 GHz to the soil moisture, expressed as %vol. In this example, the non-linear variations are of the order of −1.2 and −0.8 K/% SM, respectively, for 10.7 and 18.7 GHz.

For the winter period, solid precipitation data (deemed to occur when T_air < 0 °C) were adjusted using Tb to correct snow depth. Because of the strong sensitivity of high frequencies (typically 37 GHz or higher) to snow grain size and density (see the review papers from [39,51,52,53,54]), we used both 10.7 and 18.7 GHz channels to optimize the snow depth by adjusting the solid precipitation from the initial NARR values. Figure 4 illustrates typical brightness temperature variations with the snow depth, showing a weak but significant sensitivity of snow depth, with variations of the order of 3 to 4 K/m at 18.7 GHz, and of 3 to 2 K/m at 10.7 GHz (ΔTb/ΔHsnow).

Despite the lower sensitivity of Tb-V 11 and 19 GHz to snow depths, which leads to possible biases in snow depth correction [56], the minimization procedure was carried out using a nonlinear, least-square fitting algorithm to minimize the RMSE between modelled and satellite-based Tb. We thus adjusted the precipitation by a factor Pr_c:

$${\mathrm{Precipitation}}_{\mathrm{corr}}={\mathrm{Pr}}_{\mathrm{c}}{\mathrm{Precipitation}}_{\mathrm{NARR}}$$

(3)

where Pr_c is simply a correcting factor in term of percentage of the initial NARR precipitation amount.

In general, the results of the minimization, with a residual RMSE on Tb 19 V GHz of 2–3 K, lead to a correction factor β_air (Equation (1)) on the order of 0.6–0.75, as well as a reduction of precipitation by 40% to 50% (Pr_c, Equation (3)) [56].

The snowpack insulation effect was corrected in the second processing step (see Section 3.2).

3.2. Snow Density Module Adjustment and Analyzed Variables

The last processing step for retrieving Tg was to constrain the snow density simulated by CLASS using microwave brightness temperature at 11 GHz. The snow density can be linked to the polarization ratio at 11 GHz (PR11 = Tb_11H/Tb_11V) using the same polynomial equation over a large range of snow depths (Figure 5). Assuming that density remains the main uncertainty in the snowpack state (since snow depth and temperature were indirectly adjusted through precipitation and air temperature), this relationship allows the adjustment of the mean snowpack density seasonal evolution using AMSR-E PR11 observations.

Figure 6 shows an example of the daily variations of the measured AMSR-E PR11 compared to simulated PR11 using the outputs of CLASS simulations (snow depth, temperature, and density) for four consecutive winters at the Salluit site (see Figure 1). The observed differences in winter can be attributed to underestimated snow density simulations (red line in Figure 6).

Snow density values from CLASS were adjusted using the observed PR11 in an iterative minimization procedure. We used the Levenberg–Marquardt least squared error minimization algorithm [57], which shows rapid convergence (a maximum of 10 iterations on snow density values to reach a minimum threshold when ΔPR (simulated–measured) is equal to 0.005, corresponding to a variation of snow density of ±20 Kg m⁻³; see example in Figure 7). Using the Sturm et al. [33] relationship between snow effective thermal conductivity (Keff, W m⁻¹ K⁻¹) and density, Keff values were corrected following the derived adjusted snow density values (Figure 7). The new CLASS run, i.e., the snow energy budget and heat transfer within the snowpack with this new value of snow density and the newly derived thermal conductivity, results in a new estimate of the temperature of the first soil layer below the snow.

In summary, the analysis was focused on the simulated variables summarized in

Table 2. Definition of simulated variables used for the analysis.

Variable (See Figure 2)	Acronym	Definition
Surface Temperature	BO Ts	Ts simulated by CLASS, before optimization
	AO Ts	Ts simulated by CLASS, after optimization using RS data (MODIS LST and AMSR-E Tb) (meteorological adjustment)
First layer soil temperature	BO Tg	Tg simulated by CLASS, before optimization
	AO Tg	Tg simulated by CLASS, after optimization using RS data (MODIS LST and AMSR-E Tb) (meteorological adjustment)
	AC Tg	Tg simulated by CLASS, after optimization using RS data (MODIS LST and AMSR-E Tb) and after correction for snow density (using PR 11 GHz)

, and compared to in-situ measurements.

4. Results

4.1. Summer Soil Temperature Estimates

The summer analysis is used here only to test and calibrate the approach and the efficiency of the simplified T_air and soil moisture adjustment procedure. Results for the Inuvik site are depicted in Figure 8, where a good convergence of the simulated Tb 19 GHz toward the measured Tb (RMSE of 2.7 K and bias of −1.7 K) is observed. Simulations were improved significantly at each site.

Table 3. RMSE and Bias (in K) at the four tundra sites for the available summer time series for the surface temperature (Ts) and soil-temperature (Tg). BO and AO refer respectively to the values before and after optimization procedures (meteorological adjustment) (see definitions in Table 2).

Summer periods (July–August)
Sites	# of Years		BO Ts	AO Ts	BO Tg	AO Tg
Salluit	5	RMSE	5.8	3.4	7.2	3.6
		Bias	2.4	−2.2	3.2	1.9
North Slope	4	RMSE	4.1	2.1	4.3	1.8
		Bias	0.9	−1.4	0.1	−0.2
Inuvik	8	RMSE	5.2	1.9	3.7	2.9
		Bias	1.1	−0.5	1.7	1.4
Daring Lake	2	RMSE	4.1	2.2	9	3
		Bias	−1	−0.7	3.3	1.2
Total	19	RMSE	4.8	2.4	6.1	2.8
		Bias	0.9	−1.2	2.1	1.1

summarizes the results for each studied site over several years. For the 19 summer periods analyzed (four sites and all the years used between 2003 and 2011), results show a reduction by a factor of 2 of the RSME and the bias between measured and simulated Tg, after local adjustment of driving meteorological variables.

4.2. Snow-Covered Soil Temperature Retrieval in Winter

The winter analysis results for the Inuvik site are highlighted in Figure 9. This example shows the improvement between the raw simulated Tg without RS data—i.e., before the optimization (CLASS BO Tg, the black dotted line in Figure 9c)—compared to the simulated Tg after correcting simulated snow density (CLASS AC Tg, the red line in Figure 9c). Figure 9c clearly shows that the initial simulated soil temperature was too warm, due to an underestimated snow density that led to an overestimation of snow insulation. The adjustment of snow density, using the optimization of simulated PR11 (HUT) with satellite-based measured PR11 (Figure 9d), was able to increase the density and in turn, increase the thermal conductivity (Figure 9e, red line). This higher thermal conductivity reduces the soil temperature under the snow and improves the agreement with soil temperature measurements at the meteorological station.

Table 4. RMSE and bias (in K) at the four tundra sites for the studied winter time series, for the surface temperature (Ts) and soil temperature (Tg). BO refers to the values before the optimization procedure (meteorological adjustment) and AC to the values after correction of the snow density (see definitions in Table 2).

Winter (December–January–February)
Sites	# of Years		BO Ts	AC Ts	BO Tg	AC Tg
Salluit	5	RMSE	7.3	3.5	7.2	3.3
		Bias	1.1	0.2	2.7	0.9
North Slope	4	RMSE	6.8	2.3	7.9	3.1
		Bias	0.7	−1.7	2.3	1.4
Inuvik	8	RMSE	6.8	2.9	7.0	2.0
		Bias	0.8	−1.4	2.2	0.5
Daring Lake	2	RMSE	6.3	1.9	7.7	2.5
		Bias	2.0	−0.4	3.4	0.8
Total	19	RMSE	6.8	2.7	7.5	2.7
		Bias	1.2	−0.8	2.7	0.9

summarizes the results for each studied site, for both the snow surface temperature (Ts), and the snow-covered soil temperature (Tg). The overall results (19 winter periods) show a reduction of 64% in the RSME between the measured and simulated Tg when adjusting with RS data. The mean bias is also significantly reduced from 2.7 K to 0.9 K.

4.3. Time Series Analysis

Figure 10 shows the results of corrected CLASS Tg variations (cyan line) in Inuvik over eight winter periods, compared to in-situ measurements (blue line) and Tg variations derived from NARR data (black line). We show that CLASS-corrected Tg leads to a good representation of soil temperature under the snowpack along the whole period, showing different conditions of temperatures and snow precipitation. In this specific case, Tg values from NARR are too cold, close to Tair, and are probably associated with an underestimated snow insulation effect. It is also suggested that the mean winter Tg is not linked to snow depth in this case (see the discussion from [5]), highlighting the need to better consider the snow insulating effect on Tg. The cumulation of Tg underestimations/overestimations can have significant impacts on the permafrost evolution (see Section 5).

5. Discussion

5.1. Method Limits

The objective of this study was to demonstrate the usefulness of the polarization ratio at 11 GHz to improve simulated snow density through the winter, using a simple one-layer snow model in a climate-model land surface scheme (CLASS, in this case), in order to improve soil temperature retrievals in winter. It is clear that the adjustment methodology used is simple and could be improved with state-of-the-art assimilation schemes (see [58,59,60,61,62,63]). These studies demonstrate possible improvements in snow property retrievals, generally SWE and/or snow depth, using ensemble Kalman or particle filter assimilation schemes. However, despite the simple and computationally light method developed in the present study, we clearly show that the PR11 GHz used with a snow radiative transfer model (HUT) allows constraining snow density and soil temperature under snow in a land surface scheme.

Moreover, the first soil layer in CLASS is fixed at 10 cm, which is substantially thicker than the 11 GHz emission depth (from approximately 5 cm for a dry soil to 0.5 cm for wet soil, [64,65]). The soil temperature values derived from CLASS informed with RS data (AC Tg) should be seen as an effective retrieved temperature. The magnitude of the difference before and after correction (a factor of two; see Table 4) shows that this issue does not invalidate the observed improvements. Future improvements could be made using a heat transfer soil model with a higher vertical resolution near the surface.

One also should note that the proposed approach must be adapted if the surface is covered by vegetation such as shrubs. Vegetation in tundra and taiga zones generates a trapping effect on snow, leading to greater snow depth [66] as well as a modified snowpack state [6,67] that could alter the Tg variations. Because shrub expansion is expected to increase under observed warming [12], the derived snow–vegetation interaction processes should be taken into account in the future in LSM. In addition, for shrubs with a height above the snow cover, vegetation scattering and emission contribution on the microwave signal have to be accounted for in the HUT model using a vegetation radiative transfer model [68].

5.2. Degree-Day Index Comparison

The observed differences between soil temperatures derived from optimized CLASS outputs informed with RS and NARR data, for example, could have a significant impact on permafrost evolution assessment. Figure 11 shows the cumulative effect of winter AC Tg (the frozen degree-day index (FDDI), the sum of negative surface soil temperatures over the year). The values were averaged over the studied period at each site. The NARR-cumulated soil temperatures (FDDI-Ntg, the cross in Figure 11) are significantly underestimated compared to in-situ measurements (FDDI-Stg) and to those retrieved from CLASS (FDDI-Ctg; the proposed approach). The latter is well correlated to in-situ measurements. A significant offset of 1500 to 2000 degree-days appears between these datasets. It is hypothesized that the observed negative bias in temperatures from NARR FDDI is likely a result of underestimated insulating capacity of the snowpack (too-low snow depth or density, for example). These results show that providing RS data to a land surface scheme for improving snow, and in turn soil temperature, could lead to an improved quantification of the evolution of the active layer in permafrost regions.

6. Conclusions

This paper presents a method that uses the polarization ratio from the satellite microwave brightness temperatures at 11 GHz (PR11 H/V) to improve Arctic snow density winter simulations from a simplified snow model in a land surface scheme (LSM) (the CLASS model, in this case). Improved snow density leads to a better estimation of soil temperature under the snow (Tg). Arctic snow has the peculiarity of being very dense near the surface, due to frequent snow-blowing events and sustained cold temperatures, while the bottom of the snowpack is typically less dense, with thick depth hoar layers formed through temperature gradient metamorphism [7]. This leads to a combination of high thermal conductivity for the upper layers and low conductivity for the bottom layers of the snowpack. The result of this combination can significantly alter snowpack-insulating properties, as these layers develop through the winter [67]. No LSM is currently able to accurately simulate such a specific density stratigraphy. The snow models generally produce the reverse density stratigraphy, with low density at the surface and high snow density at the bottom, because they are mostly based on compaction and ignore the water vapor fluxes that lead to depth hoar formation (see [67]). Simulations with CLASS LSM clearly show typical overestimates of soil temperature under the snow. We show that the model can be well corrected by adjusting the mean snow density over the winter (improvement of Tg by a factor of 2) using the PR11 H/V as reference, without the use of any in-situ measurements. In comparison, for the four tundra sites across Canada studied here, the NARR reanalysis dataset shows underestimated (too cold) soil temperatures under the snow, with values close to air temperatures, very likely due to underestimated snowpack insulation.

The examples shown in this paper outline the potential of microwave remote sensing to improve Tg estimates under the snow over the Arctic, which can strongly impact the soil temperature evolution, with significant consequences for permafrost, hydrology (runoff), and soil biochemistry (microbial activities).