**5. Conclusions**

Based on in situ soil temperature data measured at three TP sites (BJ, QOMS, and NADORS), we determined the hourly, daily, and monthly soil apparent thermal diffusivity values of the 0.0 m to 0.20 m layer for 2014–2016 by using a conduction–convection method combined with DHR. The hourly, daily, and monthly k values of the 0.0 m to 0.20 m layer were obtained. The hourly and daily k values ranged from 0.3 × 10−<sup>6</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> to 1.9 × 10−<sup>6</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> at BJ, and from 1.0 × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> to 4.0 × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> at QOMS and NADORS. The monthly k ranged from 0.4( ±0.0) × 10−<sup>6</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> to 1.1( ±0.2) × 10−<sup>6</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> at BJ, from 1.7( ±0.0) × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> to 3.3( ±0.2) × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> at QOMS, and from 2.1( ±0.3) × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> to 3.1( ±0.1) × 10−<sup>7</sup> m<sup>2</sup> s<sup>−</sup><sup>1</sup> at NADORS. The results suggested that k was not constant over a day, and k showed seasonal variations. The variations of k with θ appeared to be roughly similar for unfrozen soil at these sites and years, namely, k increased sharply before it reached a peak value as θ increased, and then it tended to be stable or varied slightly with further increases in θ. The correlation coefficients (r) between k and θ ranged from 0.37 to 0.80, and 0.80 to 0.92 on hourly and monthly timescales, respectively. However, the relationship between k and θ changed when soil temperature was below 0 ◦C. Our results also suggested that the k and NDVI values were significantly related on monthly and annual timescales, with r ranging from 0.73 to 0.93. These results broaden our understanding of the relationship between in situ k and θ. The presented values of k at various timescales can be used as soil parameters when modeling land–atmosphere interactions at these TP regions.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/rs14174238/s1, Figure S1. The variations of soil temperature at (a) BJ, (b) QOMS, (c) NADORS, respectively. Figure S2. The amplitude (A, ◦C) and phase ( Φ, rad) of soil temperature 2015–2016 at BJ. Figure S3. The amplitude (A, ◦C) and phase ( Φ, rad) of soil temperature 2015–2016 at QOMS. Figure S4. The amplitude (A, ◦C) and phase ( Φ, rad) of soil temperature 2015–2016 at NADORS. Figure S5. The variation of soil apparent thermal diffusivity (k, m<sup>2</sup> s<sup>−</sup>1) with soil moisture (θ, m<sup>3</sup> m<sup>−</sup>3) on a daily timescale in 2014 (in the 1st column), 2015 (in the 2nd column) and 2016 (in the 3rd column) at (a–c) BJ, (d–f) QOMS, and (g–i) NADORS, respectively.

**Author Contributions:** B.T.: formal analysis, writing—original draft, writing—review and editing, funding acquisition. H.X.: formal analysis, writing—review and editing. R.H.: conceptualization, writing—review and editing. L.B.: conceptualization, supervision. J.G.: supervision, funding acquisition, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Second Tibetan Plateau Scientific Expedition and Research program (STEP) (2019QZKK0102), the Natural Science Foundation of China (U2142209), and the China Postdoctoral Science Foundation (2021M703558).

**Data Availability Statement:** The data that support the findings of this study are available from the Science Data Bank (https://doi.org/10.11922/sciencedb.00103; Ma et al. 2020 [52]) and additionally at the National Tibetan Plateau Data Center (https://doi.org/10.11888/Meteoro.tpdc.270910).

**Acknowledgments:** We appreciate the access to the NASA datasets. The authors appreciate all of the hard work done by researchers attending the Tibetan Plateau (TP) Experiment.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

#### **Appendix A. Determination of Soil Temperature Amplitude and Phase with DHR**

Gordon et al. [57] described how to extract soil amplitude and phase information from a soil temperature time series with the Dynamic Harmonic Regression (DHR). DHR, a simplification of the unobserved component model, has the form as follows:

$$y\_t = T\_t + \mathbb{C}\_t + \mathfrak{e}\_t \tag{A1}$$

where *yt* is the observed soil temperature time series, *Tt* is a trend or zero-frequency component, *Ct* is a cyclical component, and *et* is an irregular, white-noise component [50]. The *Ct* is modeled as a sum of the fundamental signal and its associated harmonics, as follows:

$$\mathbf{C}\_{l} = \sum\_{i=1}^{N} \left[ a\_{i,t} \cos(\omega\_{i}t) + b\_{i,t} \sin(\omega\_{i}t) \right] \tag{A2}$$

where *ai*,*<sup>t</sup>* and *bi*,*<sup>t</sup>* are stochastic time-varying parameters and *ωi* (I = 1:N) are the fundamental frequency and its harmonics (*<sup>ω</sup>*1 ) up to the Nyquist frequency (*<sup>ω</sup>N* ). This DHR model is a non-stationary extension of the discrete Fourier transform, where the amplitude (A) and phase (Φ) of the soil temperature for each time series vary with time. Identification of the time-varying parameters is achieved in a stochastic state formulation using two-step Kalman filtering and fixed-interval smoothing [50].

After obtaining the time-varying parameters, the A and Φ of any harmonic component at discrete time can be calculated by the following equations:

$$A\_{i,t} = \sqrt{a\_{i,t}^2 + b\_{i,t}^2} \tag{A3}$$

$$\Phi\_{i,t} = \tan^{-1}(a\_{i,t}/b\_{i,t})\tag{A4}$$

where *Ai*,*<sup>t</sup>* and *<sup>Φ</sup>i*,*<sup>t</sup>* are the amplitude and phase for the component with frequency *ωi* at time *t*, respectively.
