**4. Discussion**

### *4.1. The Variations of Soil Apparent Thermal Diffusivity*

With k methods based on the solution of the soil heat transfer equation, normally only daily or longer timescale k values were obtained in previous studies. In this study, we used a conduction–convection method combined with DHR to obtain hourly k. The daily and monthly k were also provided for 2014–2016 at the three sites (Figures 4 and 5 and Table 3). Figure 4 visually shows that in a day, k is not necessarily constant over a day, and it can change drastically when the soil become wetter suddenly. The knowledge of k with higher temporal resolution may have grea<sup>t</sup> implications in improving soil heat and water models. The issue of underestimating soil temperature at night with the conduction-convection method [24] may be largely resolved by using hourly k values instead of daily k values as inputs. Besides, higher temporal resolution k may help to improve the modeling of permafrost distributions [11].

Figure 4 also indicates that k had obvious seasonal variations at BJ and QOMS, while it did not vary greatly during wetting at NADORS. The soil types for both QOMS and NADORS were sandy and gravel, but θ at NADORS during wetting (> 0.10 m<sup>3</sup> m<sup>−</sup>3) was greater than that at QOMS (Figure 2b,c). Previous studies suggested that k variations during wetting depended on the magnitude of θ, and k was insensitive to changes in θ when θ reached certain thresholds (e.g., 0.1~0.2 m<sup>3</sup> m<sup>−</sup><sup>3</sup> for sand soils, [16]), which could explain why k did not vary much during wetting at NADORS.

The k values at QOMS and NADORS were much less than those at BJ (Figures 4 and 5). In addition to the soil moisture (Figure 2), the soil bulk density at QOMS and NADORS was expected to differ from that at BJ, since their soil types were sand and gravel, whereas BJ's soil type was sandy silt loam. Previous studies indicated that k varied with bulk density as well as soil moisture [17,18,59–61]. Therefore, the combined effects of the soil texture, moisture and bulk density resulted in the relative magnitude of k at the three sites. At BJ, the k values in 2014 were much larger than those at 2015 and 2016 (Figure 6a–c). The main reason may be due to differences in soil bulk density and soil moisture content.

With the conduction–convection method, Gao [26] determined the daily k ranging from 0.1 × 10−<sup>6</sup> to 2.0 × 10−<sup>6</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> at BJ during DOYs 195 to 258, 1998, which was within the range of k at BJ for 2014–2016 in this study. Additionally, with the conduction–convection method, Zhou et al. [28] determined daily k based on soil temperatures measured at 0.8 m and 3.2 m depths in 39 weather sites in the TP during 1980 to 2001. They reported that the magnitude of k in most areas of the TP was 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup>1, and relatively high k values were obtained at the central and eastern parts of the plateau with an order of magnitude of 10−<sup>6</sup> m<sup>2</sup> s<sup>−</sup>1. Our results that k at BJ are larger than those at QOMS and NADORS are broadly consistent with the findings of Zhou et al. [28]. The soil thermal properties were determined at a cold semi-desert site on the western TP for about two years by Wang et al. [48]. They calculated k as the ratio of λ/Cv, which were determined by soil temperature and heat fluxes measured at two layers according to the method proposed by Zhang and Huang [49]. Their daily k values ranged from 3.0 × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> to 9.0 × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup>1, which was larger than our k values at QOMS and NADORS, and less than those at BJ. The differences in k were attributed to the differences in soil texture, moisture, and the method of determining k. Even with a well-calibrated soil heat flux sensor, it is difficult to measure the soil heat flux accurately [62–65], because the soil heat and moisture fluxes are disturbed [66]. Heat flux plates measure only sensible heat as it moves past the plate by means of the temperature gradient, which exists across the plate. Latent heat, which is hidden in the evaporative process, is not detected [67]. Therefore, cautions should be exercised when determining soil thermal properties using soil heat flux plate data.

Using monthly soil temperature data measured in the 0.2 m to 3.2 m soil layer at the northern permafrost regions, Zhu et al. [11] obtained k with a conduction method, ranging from 0.2 × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> to 2.0 × 10−<sup>6</sup> m<sup>2</sup> s<sup>−</sup>1. The minimum k values occurred at sites with large soil organic carbon (SOC, 60~70 kg C m<sup>−</sup>3), and they were less than our minimum results. A possible reason is that there is more SOC at the northern permafrost regions than at our sites, and the k of SOC is an order of magnitude smaller than that of typical soil minerals [68].

Few studies provided long-term k using a conduction–convection method on different timescales over the TP. The long-term k values obtained in this study on different timescales can be used as input for land surface models over this region.

#### *4.2. The Relationship between Soil Apparent Thermal Diffusivity and Soil Moisture*

Normally, k increases rapidly with increasing θ to a maximum and then decreases with further increases in θ. This is explained by the fact that the heat capacity increases linearly with θ, whereas the thermal conductivity experiences its most rapid rise at low θ, leading to the ratio of thermal conductivity to heat capacity to have an internal maximum as a function of θ [29,69]. The increase rate of k with θ differs for different soil textures. A laboratory experiment showed that for sandy soils, k increased rapidly with θ when θ was less than 0.1~0.2 m<sup>3</sup> m<sup>−</sup>3, and then it remained stable or slightly decreased; while the variation of k with θ was smaller for silty and clay-textured soils [16].

Figure 6 shows that although the probability distributions of k and θ differed, especially at different sites, the variation trends of unfrozen soil k versus θ were roughly similar, which was consistent with previous studies reported based on laboratory measured data under controlled conditions (e.g., [16,17,29,59]) and in situ data (e.g., [22,25,70]). We should be aware that this trend of k versus θ reported with laboratory measurements generally only applies when the soil temperature is room temperature, since measurements in the laboratory are usually conducted at room temperature. However, if some ice is present in a soil layer, the k versus θ relationship is expected to deviate from this trend, since ice has a higher k value than water (1.1 vs. 0.14 × 10−<sup>6</sup> m<sup>2</sup> s<sup>−</sup>1, [68]). That is why there are some "outliers" in the trend of k versus θ as mentioned above, i.e., k varied greatly within a narrow range of θ, as the blue point when θ < 0.1 m<sup>3</sup> m<sup>−</sup><sup>3</sup> shown in Figure 6a,c,g–i.

About 2/3 blue points in Figure 6c (when θ ranged from 0.06 m<sup>3</sup> m<sup>−</sup><sup>3</sup> to 0.075 m<sup>3</sup> m<sup>−</sup>3) appeared on DOYs 332–364. To elucidate this phenomenon, the variations of soil temperature, θ and k over DOYs 300–365, 2016, at BJ are shown in Figure 9. One can see that on DOYs 332–364 soil liquid water content, θ in the 0.05–0.20 m layer showed a downward trend, and as the soil temperature dropped to the freezing point, more and more ice was expected to form. At the same time, θ at the depth of 0.10 m varied slightly, while k tended to increase greatly. Therefore, one can see that k had a wide range from 0.6 × 10−<sup>6</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> to 0.95 × 10−<sup>6</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> with a narrow range of θ in Figure 6c. Note that almost no water loss is assumed since DOY 300 here, as soil water evaporation is generally small when soil temperature is below freezing. Therefore, the decrease in liquid water content in the soil layer is expected to be a result of increasing ice content.

For Figure 6i, the "outliers" (when θ < 0.054 m<sup>3</sup> m<sup>−</sup><sup>3</sup> and k ranged from 1.4 × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> to 3.4 × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup>1) appeared on DOYs 1–63, 2016 at NADORS. The variations of soil temperature, θ and k over DOYs 1–70, 2016, at NADORS are shown in Figure 10. In contrast to Figure 9, the soil temperature tended to increase from DOYs 1 to 63, resulting in an increasing amount of time that soil temperature was above freezing (Figure 10a). We expect that ice content in the 0–0.20 m soil layer tended to decrease over time, although θ appeared to vary slightly. Figure 10 also shows that as the ice content decreased with increasing soil temperature, k tended to decrease, eventually varying from approximately 1.4 × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> to 3.4 × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup>1. Therefore, one can see that k varied greatly with a narrow range of θ in Figure 6i. Figures 9 and 10 explains why large k changes occurred over a small range of θ in Figure 6a,c,g–i.

**Figure 9.** The (**a**) soil temperature (Ts, ◦C) at the depths of 0.05 m, 0.10 m and 0.20 m, and soil apparent thermal diffusivity (k, m<sup>2</sup> s<sup>−</sup>1), and (**b**) soil moisture (θ, m<sup>3</sup> m<sup>−</sup>3) at the depths of 0.05 m, 0.10 m and 0.20 m over DOYs 300–365, 2016, at BJ. The labels of k are marked in cyan when the soil amplitude at the 0.2 m depth is less than 0.5 ◦C.

Compared to BJ and NADORS, there were no apparent "outliers" at QOMS (Figure 6d–f). The soil at QOMS was driest among the three sites, and θ decreased to the minimum value (<0.03 m<sup>3</sup> m<sup>−</sup>3) in September from the summer (Figure 2b), and small θ lasted until winter. Some laboratory experiments indicated that the unfrozen water content of freezing soil was largely controlled by the initial volumetric water content [71], which may prove why there was almost no ice after September at QOMS. Ochsner and Baker [20] presented some in situ measurements of soil thermal properties across a full range of temperatures encountered in freezing and thawing soil, and the measurement and model both showed that for temperatures between −5 ◦C and 0 ◦C, soil thermal properties were strongly temperature dependent. They explained that temperature dependence was primarily the result of latent heat transfer processes when water underwent a phase change. Therefore, although soil temperature changed from positive to negative at QOMS in autumn, the "frozen" soil had a lesser effect on k since little ice was produced. Thus, few "outliers" existed during the cold months at QOMS (Figure 6d–f).

To indirectly investigate the effect of soil moisture on k from the vegetation aspect, the relationship between monthly k and NDVI was also examined (Figure 8), which is similar to the relationship between k to θ (Figure 7). The reason may be that NDVI is closely related to θ at these sites. The r coefficients between NDVI and θ ranged from 0.73 to 0.92 on a monthly timescale, and from 0.69 to 0.82 on an annual timescale (Table 4). Compared to the other two sites, more vegetation was present in BJ, so the NDVI of BJ represented more vegetation information than bare soil. Therefore, the correlation between k and NDVI was weaker than that between k and θ at BJ, while closer to that between k and θ at the other two sites, as shown in Figures 7 and 8.

With in situ soil temperature, we determined k by the conduction–convection method for 2014–2016 at three TP sites and examined the relationship between k and θ. Our results indicated that k had a clear relationship with θ for unfrozen soil, but the relationship changed when the soil temperature was less than 0 ◦C and the initial θ was not too small. These findings broaden our understanding of the relationship between in situ k and θ.

**Figure 10.** The (**a**) soil temperature (Ts, ◦C) at the depths of 0.0 m and 0.20 m, and soil apparent thermal diffusivity (k, m<sup>2</sup> s<sup>−</sup>1), and (**b**) soil moisture (θ, m<sup>3</sup> m<sup>−</sup>3) at the depths of 0.0 m, 0.10 m and 0.20 m over DOYs 0–70, 2016, at NADORS. The labels of k are marked in cyan when the soil amplitude at the 0.2 m depth is less than 0.1 ◦C.

**Table 4.** Correlation coefficient (r) of the monthly soil moisture and NDVI for 2014-2016 at the three sites.


1\*\* *p*-value < 0.01, 2 \* *p*-value < 0.05.

### *4.3. Limitations*

In this study, long-term k was determined for frozen and unfrozen conditions. Due to the uncertainty of the k method, several k values were removed during the transition periods with soil thawing/freezing when soil temperature variations were low, as mentioned in Section 2.4. Another method or sensor is needed to determine k during soil thawing/freezing periods.

Our results indicated a potential effect of ice content on the relationship between k and θ, while the amount of ice content was inferred by combining the variations of soil temperature and soil moisture, rather than direct measurements. The ice content measured in situ is vital to quantify the relationship between k and soil moisture (including liquid and ice phases). The thermo-TDR sensor is a candidate for in situ measurement of both liquid and ice contents, since its performance was satisfactory in laboratory experiments reported in previous studies (e.g., [72–74]). In addition, thermo-TDR was also used in the field to measure soil thermal properties during thawing and freezing (e.g., [20]). By using actively heated fiber Bragg grating (AH-FBG) sensors, Wu et al. [75] measured the ice content of frozen soil in laboratory. The AH-FBG sensor integrates the functions of active heating and temperature measurement, which can accurately detect the thermal response of frozen soil [75]. We recommend using soil thermocouples and thermo-TDR sensors or only AH-FBG sensors for soil temperature, water content (liquid and ice phase) and thermal property measurements over multiple thawing and freezing cycles to more deeply explore the time variations of k and its relationship with water content.

In this study, we did not examine the effects of other soil factors (e.g., soil texture, SOC) on k due to a lack of data. Zhu et al. [11] suggested that SOC is the dominate factor (among soil texture, bulk density, moisture, and SOC) controlling the variability of diffusivity at 200 sites in high latitude regions, and k is a strong predictor for simulated permafrost extent. In future investigations, additional soil factors should be included in the study of long-term variations of soil thermal properties.
