Soil temperature estimation at different depths over the central Tibetan Plateau integrating multiple Digital Earth observations and geo-computing

ABSTRACT

Soil temperature (ST) plays a critical role in ecosystems. Monitoring high-resolution ST profiles remains challenging due to the inherent heterogeneity of ST in space, time, and depth. To address this challenge, in this study, we integrate remote sensing techniques and deep learning methods to retrieve spatiotemporal continuous ST from 2011 to 2020 at four depths (5, 10, 20, and 40 cm) over the central Tibetan Plateau (TP). Landsat and MODIS observations were fused to obtain land surface properties at different resolutions. The fused variables were integrated with in-situ ST and soil moisture (SM) measurements to estimate spatiotemporally continuous (0.0005°, 0.0025°, and 0.0125°) ST profile through deep learning-based training models. The deep belief network (DBN) was applied to estimate ST: (a) from layer to layer (LW-DBN), and (b) from land surface properties directly (Direct-DBN). The ten-fold cross-validation indicates that both approaches achieve promising results (R² > 0.836, MAE < 2.152 °C), and the Direct-DBN outperformed the LW-DBN at all spatial scales and depths. ST retrieval at deeper depths and coarser resolutions tend to have better monitoring accuracy. Further analysis indicates an increment of ST at all four depths over the study period, which provides valuable insights into global warming.

KEYWORDS:

• Soil temperature profile
• deep learning
• Tibetan Plateau

1. Introduction

The global temperature tends to increase, as reported by Intergovernmental Panel on Climate Change (IPCC). The rising temperature will accelerate sea-level rise, increase floods, worsen droughts, and exacerbate heat waves (Kohn and Royer 2010; S. Li and Xiao 2022; Shukla et al. 2022; Taylor, Myrold, and Bottomley 2019). As one of the essential components of the Earth’s surface energy budget and the assessment of global climate change, soil temperature (ST) plays a vital role in climatology, agriculture, hydrology, and the environment. It is also one of the critical factors affecting the emission of greenhouse gases (GHG) and numerous ecosystems (Choler 2018; Kohn and Royer 2010; S. Li et al. 2022; Li and Xiao 2022; Shukla et al. 2022; Taylor, Myrold, and Bottomley 2019; Yan et al. 2018; Zheng, Hunt, and Running 1993). ST at different depths are crucial agrometeorological indicators for ecological modeling and precision agricultural activities (Huang et al. 2020).

Temperature is crucial for the rate of chemical reactions and the activity of soil organisms, and ST profile estimation is significant to understand the ecosystem's responses to climate change. An increase in ST over permafrost regions may affect active layer depth, deplete soil organic carbon and increase carbon dioxide emissions. Numerous techniques have been proposed to measure the ST profile (Li et al. 2022). In-situ measurements are still the most convincing way to provide accurate ST profile monitoring at different depths, however, field experiment is costly, time-consuming, and spatially discontinuous. Furthermore, the uneven distribution of observing stations makes it hard to truly represent the spatiotemporal dynamics of ST, which is crucial to land surface process, ecosystem response analysis, sustainable management, and agricultural applications (Alizamir et al. 2020; Dharssi et al. 2013; Huang et al. 2020; Zhan et al. 2014). Accurately measuring the spatiotemporal distribution of ST helps understand the ecological environment change and contributes to precision agriculture management (Alizamir et al. 2020; Huang et al. 2020; Qin et al. 2020).

Consequently, numerous models have been proposed for ST profile estimation (Dharssi et al. 2013). Models for ST estimation at different depths can be mainly classified into physical models and statistical models (Huang, Li, and Lu 2008; Zhan et al. 2014). 1) Physical models, also called process-based models, are based on the interactions between the atmosphere, soil, and energy conservation, however, the use application of physical models is limited by the large input requirements from soil physical constants to climatic variables (Araghi et al. 2019; Zoras, Dimoudi, and Kosmopoulos 2012). 2) Statistical models directly estimate the relationship between ST and related variables with a high explanatory degree, while poor regional generality and collinearity are the major issues in the large-scale application of statistical models (Lobell and Burke 2010; Moriondo, Maselli, and Bindi 2007; Roberts et al. 2017). Therefore, to improve the accuracy of ST retrieval models and promote the application on a larger scale, more effective approaches are urgently required to process and analyze the existing multi-source data and their potential for ST monitoring at different depths. Machine learning and deep learning models have been introduced and are widely explored to assist ST prediction with higher accuracy and regional applicability, such as random forest (RF), extreme learning machine (ELM), and artificial neural networks (ANN) (Feng et al. 2019; Samadianfard et al. 2018; Sanikhani et al. 2018).

ST profile mapping at high resolutions is significant due to the spatiotemporal heterogeneity of temperature, and the demand has increased over recent decades. However, the methods discussed above are all point-based ST monitoring techniques. Traditional interpolation methods, such as the kriging and the inverse density weighting (IDW) methods, generate spatial continuous observations but are limited by the distribution and density of in-situ observing stations and are affected by the complexity of different environmental conditions (Shen et al. 2020). Remote sensing has been a promising technique to conquer this challenge based on the strong relationship between ST and thermal remote sensing measurements, which has been verified by numerous studies (Huang, Li, and Lu 2008; Xu et al. 2020; Zhan et al. 2014). Integrating thermal remote sensing products and models provides practical alternatives to obtain spatiotemporal continuous ST estimation with reasonable accuracy. Lugassi, Ben-Dor, and Eshel (2010) used information in soil spectroscopy to accurately reconstruct spatial distributions of surface ST. Kohn and Royer (2010) used passive microwave remote sensing to estimate ST under snowpack with lower error than ST derived from the same model but without the application of remote sensing. Zhan et al. (2014) utilized four thermal measurements to reconstruct the diurnal cycle of ST with an error of approximately 1.5 K, showing that thermal remote sensing can be applied to estimate ST efficiently. Especially the yielded credible thermal remote sensing land surface temperature (LST) products have been widely utilized to estimate ST (Huang, Li, and Lu 2008; Xu et al. 2020; Zhu et al. 2018).

However, the methods discussed above are all about surface ST retrieval through remote sensing techniques. Due to the limited penetration of relevant bands, ST beneath surface soil cannot be retrieved through remote sensing directly. It is desirable to combine model simulations with remote sensing observations and apply other land surface properties as input variables for ST monitoring at different depths. The comprehensive analysis of daily ST trends at predetermined depths can improve the quantitative understanding of the sustainable management of ecosystems (Araghi, Mousavi-Baygi, and Adamowski 2017).

The Tibetan Plateau (TP), known as the ‘Third Pole of the Earth’, plays a critical role in the Asian monsoon system, land-atmosphere interactions, and global eco-safety (Hu et al. 2019; Tan, Jia, and Wang 2021; Zhao et al. 2021). TP is one of the most sensitive regions to climate change, and has the most significant areas of permafrost terrain in the mid-latitude and low-latitude regions (Li et al. 2015; Liu and Chen 2000; Feng, Tang and Wang 1998; Tan, Jia, and Wang 2021). ST variations over permafrost regions reflect changes in the atmosphere and the Earth's surface, and it also indicates frozen ground status over TP (Li et al. 2015; Zhang, Li, and Pang 2008; Zou et al. 2017). Moreover, increasing ST can enhance freeze–thaw cycles and decrease the duration of seasonal freeze–thaw periods that may influence the stability of the permafrost (Ming et al. 2022). Now that permafrost on TP is degrading, it is essential to understand ST variations over TP and develop robust models to accurately estimate spatiotemporal continuous ST profile. However, few studies targeting on monitoring ST over TP at different depths with high resolution. In recent years, the permafrost in the TP has undergone degradation as temperatures have risen significantly. The change in active layer depth has altered the heat flow and energy balance within the soil. This has led to more complex energy exchanges between the land surface and atmosphere, making accurate estimation of surface temperatures more challenging (Cao et al. 2020; Qin et al. 2020; Yang et al. 2020). Research targeting on TP mainly utilized in-situ ST observations, model assimilations and reanalysis products. In-situ observations are point based, model assimilations and reanalysis products are with coarse spatial resolutions (Cao et al. 2020; Li et al. 2022; Yang et al. 2020). ST’s variation over TP is sensitive and the dynamic of ST varies in different regions, ST monitoring with high resolutions is essential to climate change and earth system science (Qin et al. 2020).

This study proposed deep learning-based approaches to provide spatiotemporal estimations of daily ST with high spatial resolutions at four depths (5, 10, 20, and 40 cm) over the central TP by integrating in-situ measurements and various remote sensing observations, which has not been fulfilled by previous studies and is significant to climate change and permafrost science: 1) The Land Remote Sensing Satellite (Landsat) and the Moderate Resolution Imaging Spectroradiometer (MODIS) bands with corresponding wavelengths were fused to obtain daily reflectance and radiance at different scales (0.0005°, 0.0025°, 0.0125°). 2) The fused reflectance/ radiance were applied to retrieve LST and vegetation index (VI) at three spatial resolutions, which were further integrated with factors acquired from the Global Land Surface Satellite (GLASS) and interpolated SM profile to the DBN model to estimate ST at various depths. 3) The ST beneath the surface was retrieved through the land surface properties directly and from the ST of the previous layer and SM at the current layer separately. 4) The estimated ST profile from two approaches at three spatial scales was validated by in-situ measurements from the multiscale Soil Moisture and Temperature Monitoring Network on the central TP (CTP-SMTMN) network (Chen et al. 2017).

2. Study area and data

2.1. Study area

The study area is over the central of Qinghai-Tibetan Plateau, as shown in Figure 1 (Chen et al. 2017). The Qinghai-Tibetan Plateau, located on the eastern Eurasian continent, is the highest and the most extensive plateau in the world, with an average altitude of more than 4000 m (Liu et al. 2017; Y. Wang et al. 2022). The study was conducted with a latitude ranging from 31° N to 32° N, a longitude ranging from 91°30′ E to 92°30′ E, and an average elevation of 4500 m above, as shown in Figure 1(b). This area is located in the Naqu basin and is characterized by low biomass, high soil moisture dynamic range, and a typical freeze–thaw cycle (Zhao et al. 2013). The dominant land cover is the alpine meadow; silt and sand are the primary components of the soil texture (Qin et al. 2013; Zhao et al. 2013). The air pressure is low, approximately 70 kPa, due to the high altitude (Yu et al. 2020). The study area is with a typical semiarid monsoon climate. The precipitation is usually concentrated from May to September and is limited during winter; the annual precipitation is approximately 400 mm (Su et al. 2011; Zhao et al. 2021). There were 38 stations in total applied in this study, as shown in Figure 1(c).

Figure 1. Study area.

2.2. Datasets

In-situ measurements, Landsat 7, Landsat 8, MODIS, and the GLASS product were integrated to generate spatiotemporal continuous high-resolution ST at different depths.

2.2.1. In-situ observations

Daily in-situ ST and SM measurements at four depths (5, 10, 20, and 40 cm) for each station of the CTP-SMTMN network were collected over the study area from January 1st, 2011, to December 31st, 2020. The daily in-situ ST at different depths was applied to train, calibrate and validate the proposed model through Matlab software, 3/4 of the measurements were applied to train the model, and 1/4 were used for validation. Figure 2 illustrates the ST distribution at each depth. The kriging interpolation method was utilized to generate the spatial continuous SM at corresponding depths through the in-situ SM measurements.

Figure 2. Temperature distribution at different depths (05, 10, 20, 40 cm).

2.2.2. Landsat

The Landsat program has been continually monitoring the Earth since 1972, and up to now, with nine satellites. Landsat 7 (2011.1 ∼ 2020.12) and Landsat 8 (2013.4 ∼ 2020.12), which cover the study period, were selected in this study. The reflectance of the blue band, red band, and near-infrared (NIR) band and radiance of thermal infrared (TIR) band provided by the Enhanced Thematic Mapper Plus (EMT+) onboard Landsat 7; reflectance of blue band, red band, and NIR band provided by the Operational Land Imager (OLI), and radiance of TIR band provided by the Thermal Infrared Sensor (TIRS) onboard Landsat 8 were chosen, respectively. Both Landsat 7 and Landsat 8 are with 16-day temporal resolution. From January 2011 to April 2013, the 16-day Landsat 7 observations were fused with MODIS corresponding bands. And from April 2013 to December 2020, with the availability of Landsat 8 products, the combination of Landsat 7 and Landsat 8 provided land surface observations with 8-day temporal resolution over the study area, which were fused with MODIS products for ST retrieval.

2.2.3. Modis

The MODIS land products applied in this study were the daily Surface Reflectance Product (MOD09), the daily Land Surface Temperature and Emissivity Product (MOD11), and the 8-day Global Evapotranspiration (ET) Product (MOD16) with a spatial resolution of 1000, 1000, and 500 m, respectively. The reflectance of the blue band, red band, and NIR band were extracted from MOD09 and then fused with Landsat observations with the corresponding channels. The radiance of the thermal band was collected from MOD11 and fused with thermal radiance obtained by Landsat.

Detailed information about the bandwidths and the spatial resolutions of the corresponding channels for Landsat 7, Landsat 8, and MODIS that were analyzed are shown in Table 1. Before fusing, MODIS and Landsat observations were resampled to the same three scales (0.0005°, 0.0025°, 0.0125°).

Table 1. Bandwidth and spatial resolution of MODIS and Landsat corresponding bands.

	Bandwidth (μm)			Resolution (m)
	Landsat 7	Landsat 8	MODIS	Landsat 7	Landsat 8	MODIS
Blue	0.45 ∼ 0.52	0.45 ∼ 0.51	0.46 ∼ 0.48	30	30	500
Red	0.63 ∼ 0.69	0.64 ∼ 0.67	0.62 ∼ 0.67	30	30	250
NIR	0.77 ∼ 0.90	0.85 ∼ 0.88	0.84 ∼ 0.88	30	30	500
TIR	10.40 ∼ 12.50	10.60 ∼ 11.19	10.78 ∼ 11.28	60	100	1000

2.2.4. Glass

The GLASS is developed based on multiple spatiotemporal continuous land surface satellite remote sensing products. Compared with other land surface observations, the GLASS products have no gaps or missing values and are highly accurate (Zhao et al. 2013). Considering there is no significant variation of the canopy within a short period, the 8-day GLASS leaf area index (LAI) product, the fraction of absorbed photosynthetically active radiation (FAPAR) product, and the net primary production (NPP) product were applied in this study with spatial resolutions of 250, 250, and 500 m, respectively. The GLASS products were samples to three scales (0.0005°, 0.0025°, 0.0125°) for ST retrieval.

3. Methods

In this study, daily ST profiles at four depths (5, 10, 20, 40 cm) were estimated by integrating in-situ measurements, Landsat, and MODIS land surface products through deep learning techniques for three spatial scales (0.0005°, 0.0025°, 0.0125°), as illustrated in Figure 3. The red line indicates the direct deep learning approach that estimated ST at all depths from land surface properties directly, and the blue line indicates the layer-wise approach that estimates ST at the deeper depth from the ST of the upper depth layer by layer.

Figure 3. Flow chart of this study.

The research is mainly with three parts: 1) the corresponding bands (blue, red, NIR, and TIR) of MODIS and Landsat with similar wavelengths are fused through an improved method to obtain daily land surface observations at three spatial scales; 2) the downscaled blue, red and NIR bands were applied to calculated VI, TIR band was used to retrieve LST through single-channel method for further ST estimation; 3) the downscaled LST and VI were integrated with MODIS ET product and GLASS NPP, FAPAR products to monitor daily ST at four depths.

3.1. MODIS and Landsat data fusion

Satellite remote sensing provides land surface observations at high spatial or temporal resolutions; however, no sensing system currently offers global coverage observations with high spatial and temporal resolutions. Due to the high spatiotemporal heterogeneity of ST, land surface properties for ST monitoring with high spatial and temporal resolutions are significant, especially for regional studies. Since the channels for ST monitoring applied in this study of Landsat and MODIS are highly consistent, the study area can be treated homogeneously with grassland coverage, and Landsat and MODIS can roughly be considered acquired at the same time. Landsat and MODIS corresponding channels (blue, red, NIR, and TIR) were fused through an improved fusion method, which has been applied to VI, LST, air temperature, and ST monitoring with promising results (Shen et al. 2020; Xu et al. 2018; Zhang et al. 2022). And for pixels with missing values, the value of the same pixel on the next nearest observing time was applied to replace the invalid one to obtain more valid datasets for further analysis.

The theory of the fusion methods is based on the relationship between the observations of reflectance/ radiance from MODIS and Landsat, which can be expressed as: (1) $\begin{array}{c}L=m\cdot M+n\end{array}$(1) where $L$ and $M$ represent the reflectance/ radiance of Landsat and MODIS, respectively; $m$ and $n$ are the regression parameters. And for a specific pixel given location and time, the relationship can be described as: (2) $\begin{array}{c}L(x,y,t)=m\cdot M(x,y,t)+n\end{array}$(2) where ($x,y$) represents a given location and $t$ is the acquisition date. The relationship between Landsat observations of reflectance/ radiance and those of MODIS of a specific pixel at observed time $t_0$ and predicted time $t_p$ can be described as: (3) $\begin{array}{c}L(x,y,t_0)=m\cdot M(x,y,t_0)+n\end{array}$(3) (4) $\begin{array}{c}L(x,y,t_p)=m\cdot M(x,y,t_p)+n\end{array}$(4)

Through Equations (3) and (4), the reflectance/ radiance with good spatial resolution at the predicted time $t_p$ can be calculated as: (5) $\begin{array}{c}L(x,y,t_p)=L(x,y,t_p)+m\cdot [M(x,y,t_p)-M(x,y,t_0)]\end{array}$(5) Before the mission of Landsat 8, the reflectance/ radiance of Landsat 7 16-day blue, red NIR, and TIR bands were integrated with MODIS corresponding bands from January 2011 to April 2013 to obtain land surface measurements with good resolutions. With the availability of Landsat 8 products from April 2013, combined Landsat 8-day observations from April 2013 to December 2020 were integrated with daily MODIS products to generate reflectance/ radiance of blue, red, NIR, and TIR bands at the spatial resolutions of 0.0005°, 0.0025° and 0.0125°.

3.2 LST monitoring

The accuracy of LST calculation is crucial for ST estimation. The single channel method is applied to calculate LST through the downscaled land surface observations, which can be expressed as: (6) $LST={\displaystyle \frac{T_b}{1+\left(\lambda \cdot {\displaystyle \frac{T_b}{\rho }}\right)\cdot \ln \epsilon }}$(6) (7) $\begin{array}{c}\rho ={\displaystyle \frac{hc}{\sigma }}\end{array}$(7) where $T_b$ represents the brightness temperatures of the TIR band, which can be calculated as Equation (8); λ is the wavelength of the emitted radiance; $\epsilon$ is the land surface emissivity, which can be calculated as Equation (9); $h$ is Plank’s constant (6.626 × 10⁻³⁴ J Sec), $c$ is the velocity of light (2.998 × 10⁸ m/sec), and $\sigma$ is Stefan Boltzmann’s constant (5.670 × 10⁻⁸ Wm⁻²K⁻⁴). (8) $T_b={\displaystyle \frac{K_2}{\ln \left({\displaystyle \frac{K_1}{L_{\lambda }}+1}\right)}}$(8) where $L_{\lambda }$ is the at-sensor radiance of the TIR band, $K_1$ and $K_2$ are calibration constant, which is 666.09 and 1282.71 in this study, respectively. (9) $\begin{array}{c}\epsilon =\{\begin{array}{c}{a}_{\epsilon }\cdot R_{red}+b_{\epsilon },NDVI<0.2\\ {\epsilon }_v{P}_v+{\epsilon }_s\cdot (1-P_v)+C,0.2\le NDVI\le 0.5\\ {\epsilon }_v+C,NDVI>0.5\end{array}\end{array}$(9) where $a_{\epsilon }$ and $b_{\epsilon }$ are empirical parameters, ${\epsilon }_v$ (0.986) and ${\epsilon }_s$ (0.967) are emissivity values of vegetation and bare soil, $NDVI$ is the normalized difference vegetation index, $P_v$ is the vegetation fraction, and C is the value of surface roughness effects. $NDVI$, $P_v$ and C in Equation (9) can be calculated as: (10) $NDVI={\displaystyle \frac{R_{red}-R_{nir}}{R_{red}+R_{nir}}}$(10) (11) $\begin{array}{c}{P}_v={\left({\displaystyle \frac{NDVI-NDV{I}_{min}}{NDV{I}_{max}-NDV{I}_{min}}}\right)}^2\end{array}$(11) (12) $C=(1-{\epsilon }_s)\cdot {\epsilon }_v\cdot 0.55\cdot (1-P_v)$(12)

3.3. ST retrieval at different depths

As the second generation of neural networks and a branch of machine learning, the deep learning method was applied in this study to simulate the non-linear relationship between ST and multi-source data and to estimate the spatiotemporal ST profile. The deep belief network (DBN) proposed by Hinton (2009) is a generative model that can identify data features, classify data, and generate data according to the maximum probability by training the weight between neurons.

To improve model performances for ST monitoring at various depths, we introduced the downscaled blue band, red band, and NIR band to generate the Enhanced Vegetation Index (EVI), which has been applied to ST retrieval with convincing results at the spatial resolutions of 0.0005°, 0.0025°, and 0.0125° as inputs of this model (Xu et al. 2020). EVI can be calculated as follows: (13) $EVI=2.5\times {\displaystyle \frac{R_{nir}-R_{red}}{R_{nir}+6\cdot R_{red}-7.5\cdot R_{blue}+1}}$(13) where $R_{nir}$, $R_{red}$ and $R_{blue}$ represent reflectance of the NIR band, red band, and blue band, respectively.

In this study, DBN is utilized to estimate ST at multiple depths combining the downscaled LST and EVI from the previous procedures, MODIS ET products, and GLASS LAI/ NPP/ FAPAR products. The MODIS and GLASS products were resampled to the same spatial resolutions as LST for ST retrieval. The in-situ ST measurements from the CTP-SMTMN were used to verify and validate the model. Specifically, we applied two different approaches to estimate ST at different depths: 1) directly estimate ST (the Direct-DBN) at four depths through the remote sensing land surface measurements and the interpolated SM measurements at corresponding depth as illustrated by Equation (13); and 2) a layer-wise approach (the LW-DBN), which estimates surface ST at a depth of 5 cm first and then use the upper layer ST estimations with interpolated SM measurements of the next depth to estimate ST at the next depth, one depth at a time, as illustrated by Equation (15). DBN is one of the most commonly used deep learning models and has been widely utilized in remote sensing for monitoring land surface properties. In the applied DBN architecture, the output of one training layer becomes the input to the next training layer, and the training data is generated by determining the maximum probability weights between neurons. Matlab is utilized to develop the DBN algorithm through functions provided by the Matlab deep learning tools Deep Neural Network (Tanaka, 2023). ST and input factors were first normalized to improve the model’s learning ability. Then, they were separated into training group and testing group for model training. After that, unfold DBN to NN (Neural Network), add the output layer, and tune the network. Finally, ST was monitored through denormalization. Through this iterative process, the full neural network is constructed. The grid search method is applied in this study to tune the hyper-parameters, for both the Direct-DBN and the LW-DBN, we utilized a 4-layer DBN model with two hidden layers (the neurons of each hidden layer are 50 and 20), and with ‘number of epochs’ of 130, ‘batch _ size’ of 120, ‘momentum’ of 0.2, ‘learning rate’ of 0.2. (14) $\begin{array}{c}S{T}_n=f_D^n(LST,EVI,ET,LAI,NPP,FAPAR,S{M}_n)\end{array}$(14) where $f_D^n$ represents the Direct-DBN training model for ST monitoring at depth $n$ through land surface properties,$LST$ is land surface temperature,$EVI$ is enhanced vegetation index, $ET$ is evapotranspiration, $LAI$ is leaf vegetation index,$NPP$ is net primary production, $FAPAR$ is fraction of absorbed photosynthetically active radiation, $SM$ is soil moisture; the $n$ value is 1, 2, 3, and 4, which indicates 5, 10, 20, and 40 cm, respectively. (15) $\begin{array}{c}S{T}_n=\{\begin{array}{c}{f}_D^n(LST,EVI,ET,LAI,NPP,FAPAR,S{M}_n),n=1\\ f_{LW}^n(S{T}_{n-1},S{M}_n),n>1\end{array}\end{array}$(15) where $f_{LW}^n$ represents the layer-wise approach for deep ST monitoring integrating the ST estimation of the upper layer and SM of the current layer, $LST,EVI,ET,LAI,NPP,FAPAR,S{M}_n$ are with the same meaning as Equation (14).

The performances of the methodologies proposed in this study were evaluated using the coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE) through in-situ measurement from the CTP-SMTMN. (16) $R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(\widehat{y_i}-y_i)}^2}{{\sum }_{i=1}^n{(\widehat{y_i}-\overline{y})}^2}}$(16) (17) $RMSE=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^n{(\widehat{y_i}-y_i)}^2}{n}}}$(17) (18) $MAE={\displaystyle \frac{{\sum }_{i=1}^n|\widehat{y_i}-y_i|}{n}}$(18) where $\widehat{y_i}$ is the retrieved ST for sample i, $y_i$ is the in-situ ST for the same sample, $\overline{y}$ is the average ST observations, and $n$ is the number of observations.

4. Results

4.1. ST monitoring at different depths and scales

ST at four depths (5, 10, 20, and 40 cm) over the central of TP from 2011 to 2020 for three spatial scales (0.0005°, 0.0025°, and 0.0125°) was estimated through deep learning methodology integrating in-situ measurements and various remote sensing products. In-situ ST measurements were applied to calibrate and validate the approaches proposed in this study. The LW-DBN and the Direct-DBN models’ overall performances were evaluated using R², RMSE, and MAE, as shown in Tables 2–4. For the surface ST retrieval (5 cm), both the LW-DBN and the Direct-DBN integrated land surface products, including the downscaled LST, MODIS products, GLASS products, and the interpolated surface SM measurements. For the other three layers (10 cm, 20 cm, and 40 cm), the LW-DBN applied the ST estimation of the upper layer and SM measurements of the current layer. At the same time, the Direct-DBN integrated the land surface products and the interpolated SM of the current layer to monitor ST at different layers directly.

Table 2. R² of two approaches at different depths and scales.

	0.0005°		0.0025°		0.0125°
Approach	LW	Direct	LW	Direct	LW	Direct
5 cm	0.851		0.855		0.865
10 cm	0.839	0.847	0.842	0.849	0.855	0.865
20 cm	0.859	0.866	0.861	0.870	0.866	0.876
40 cm	0.836	0.846	0.837	0.846	0.843	0.859

Table 3. RMSE (${}^{\circ }\mathrm{C}$) of two approaches at different depths and scales.

	0.0005°		0.0025°		0.0125°
Approach	LW	Direct	LW	Direct	LW	Direct
5 cm	2.813		2.772		2.670
10 cm	2.933	2.738	2.752	2.591	2.649	2.437
20 cm	2.678	2.453	2.589	2.322	2.306	2.270
40 cm	2.660	2.449	2.457	2.219	2.273	2.139

Table 4. MAE (${}^{\circ }\mathrm{C}$) of two approaches at different depths and scales.

	0.0005°		0.0025°		0.0125°
Approach	LW	Direct	LW	Direct	LW	Direct
5 cm	2.152		2.101		2.024
10 cm	2.130	2.124	2.114	2.093	2.019	1.962
20 cm	1.803	1.745	1.785	1.714	1.878	1.843
40 cm	1.714	1.686	1.694	1.650	1.707	1.634

From the statistical results, both the LW-DBN and the Direct-DBN have a performance of R² greater than 0.836 and MAE smaller than 2.152${}^{\circ }\mathrm{C}$. In previous studies, the layer-wise approach is the most widely applied method for ST retrieval at different depths. However, the Direct-DBN, which used the land surface properties to monitor ST at specific depths directly, shows better performance than the LW-DBN according to the statistical results (higher R², lower RMSE, and lower MAE) at all three scales and four depths. For both of the two approaches, the retrieved ST profile with the coarser spatial resolution is slightly better than estimations with higher resolutions (0.0125° > 0.0025° > 0.0005°); for ST at different depths, RMSE and MAE decrease with the increase of layer depth (5 cm > 10 cm > 20 cm > 40 cm), while there is no significant correlation between R² and the depth of the layer. The following analysis of the results is produced mainly based on the Direct-DBN at 0.0025° spatial resolution.

4.2. Overall performance of the Direct-DBN

Figure 4 illustrates comparisons between the in-situ ST measurements from the CTP-SMTMN network and ST estimations from the Direct-DBN approach. Figure 4(a) shows the scatter plot between the observed and predicted ST profile at four depths (5, 10, 20, and 40 cm) with three spatial resolutions (0.0005°, 0.0025°, and 0.0125°) and dots located along the red line indicates the predicted value is equal to the in-situ ST measurements. Figure 4(b) is the histogram for the bias of predicted ST at four depths and three scales. Both the scatter plots and histogram indicate that the models developed in this study were valid and trustworthy, with the distribution of dots along the 1:1 red line for the scatter plot and the bias accumulated around 0 for the histogram. For the scatter plots, the distribution of the dots is more accumulated along the red line with coarser spatial resolutions (0.0125° > 0.0025° > 0.0005°), which is consistent with the statistical results from Tables 2–4. For the histogram, the red, blue, and yellow parts indicate the histogram of bias for retrieved ST at the spatial resolution of 0.0005°, 0.0025°, and 0.0125°, respectively. For ST at deeper depth, the distribution of the histograms is more accumulated around 0, and for coarse spatial resolutions, the bars within the histogram are clustered with less bias.

Figure 4. Comparisons between the observed and predicted ST. (a) Scatter plot between observed and predicted ST, and (b) histogram of bias for predicted ST.

4.3. Temporal performance of the Direct-DBN

Figure 5 is the time series of the Direct-DBN predicted ST profile at 0.0025°, in-situ ST, SM, and air temperature (Ta) measurements at four depths for one of the validation sites (site ID: MS3488), with the latitude is 31.843° N, and the longitude is 91.705° E, with an elevation of 4799 m. MS3488 is with grassland coverage, and the major soil type is sand (56.52%). The red, black, blue, and yellow lines represent the observed and predicted ST, SM, and Ta respectively. It can be noticed that the black line is highly overlapped with the red line, which illustrates that the proposed model has good accuracy. Furthermore, there is a typical annual cycle for ST at all four depths. ST in January is the lowest, and it gradually increases as time moves on and hits its peak at around July or August and then gradually decreases for all depths, which is identical to the physical facts. However, the proposed methods slightly underestimated ST with low values in January, especially at shallow layers (5 cm and 10 cm) compared with deep layers (20 cm and 40 cm). In this figure, it can be noticed that the ST of 5 cm < 10 cm < 20 cm < 40 cm, this result is consistent with the results discussed above. We also figured out that at the same depth, the peak of ST has slightly increased with time moved on from 2011 to 2020, especially for the shallow layers (5 cm and 10 cm).

Figure 5. Time series of predicted and observed ST at station MS3488.

5. Discussion

5.1. Evaluation of the Direct-DBN’s performance

In addition to the overall performances of the Direct-DBN models, the application for the specific month over the study period was also analyzed. Figure 6(a) ∼ 6(c) represents the R², RMSE, and MAE of the Direct-DBN model developed at the depths of 5, 10, 20, and 40 cm with a spatial resolution of 0.025° for each month, respectively. In Figure 6, blue, red, gray, and black indicate the statistical results of ST depth at 5, 10, 20, and 40 cm. It can be figured out that models developed at a depth of 40 cm have the best performance for most months of the year, with the highest R², lowest RMSE, and MAE. In general, 40 cm > 20 cm > 10 cm > 5 cm for the proposed approach's overall performance corresponds to the results discovered above. In addition, the proposed approach performs relatively better in summer (e.g. June, July, and August) at all four depths, and in winter (January and February) has comparatively worse results, especially for surface ST retrieval. And for specific months (e.g. April, November, and December), although the R² value is relatively high, RMSE and MAE values indicate that there were comparatively larger errors.

Figure 6. (a) R², (b) RMSE, and (c) MAE for different months.

5.2. Spatiotemporal distribution of ST at different depths and scales

Figure 7 illustrates the spatial distribution of annual averaged ST (0.0025°) retrieved via the Direct-DBN at four depths in 2011, 2014, 2017, and 2020. The yearly averaged ST ranges from −2.457 ${}^{\circ }\mathrm{C}$ to −0.320 ${}^{\circ }\mathrm{C}$ for the surface layer (5 cm), from −1.614 ${}^{\circ }\mathrm{C}$ to 0.335 ${}^{\circ }\mathrm{C}$ for the 10 cm depth, from 0.120 ${}^{\circ }\mathrm{C}$ to 1.939 ${}^{\circ }\mathrm{C}$ for the 20 cm depth, and from 1.387 ${}^{\circ }\mathrm{C}$ to 3.059 ${}^{\circ }\mathrm{C}$ for the 40 cm depth from 2011 to 2020. 2014 and 2020 have the lowest and highest annual averaged ST values for all four depths. For the four selected years, the annual mean ST of 2020 > 2017 > 2011 > 2014 of all four discussed depths. Besides that, for the same year, ST at 40 cm > 20 cm > 10 cm > 5 cm, which is consistent with the physical situation and in-situ measurements. As Figure 7 illustrates, the east part has relatively higher ST values than the west part, and the middle east region has the highest ST value at four depths for all the selected years. In addition, considering the DEM distribution, it can be figured out that regions with relatively high elevation have relatively low ST for all four depths.

Figure 7. Spatial distribution of predicted annual mean ST at a depth of 5, 10, 20, and 40 cm in 2011, 2014, 2017, and 2020. (unit: ${}^{\circ }\mathrm{C}$).

Figure 8 shows the spatial distribution of surface ST over the study area on July 24th, 2015, at three spatial resolutions (0.0125°, 0.0025°, and 0.0005°), with surface ST ranging from 3 ${}^{\circ }\mathrm{C}$ to 13.5 ${}^{\circ }\mathrm{C}$. The distributions of ST for three spatial resolutions are consistent over the study area; the mapping at 0.0025° and 0.0005° provide more detailed information, especially for regions with extremely high or low ST values. Moreover, ST retrieval at 0.0125° is slightly lower compared with the same regions at 0.0025° and 0.0005°; and ST at 0.0025° varies significantly compared with the two other scales.

Figure 8. Spatial distribution of estimated surface ST on July 24th, 2015, for three spatial resolutions (0.0125°, 0.0025°, and 0.0005°).

We chose two sub-regions, A over the west part and B over the east part of the study area, as marked by red and blue rectangles, to zoom in to compare the discrepancies of ST estimated at different spatial resolutions. For both sub-regions, the ST values of 0.0125° are slightly lower than the other two scales. The highest and lowest ST value of 0.0025° is higher or lower than that of 0.0005°, which is consistent with the spatial distribution of the study area. For 0.0125°, the zoomed-in spatial distribution of estimated ST is mostly a blur, and it is hard to recognize the spatial variation of ST. At the resolutions of 0.0025° and 0.0005°, it is clearer to obtain the spatial distribution of ST in detail. Although the spatial distribution of the three resolutions is generally consistent, the detailed distribution varies slightly. For the southeast part of sub-region A, the ST value at 0.0125° and 0.0025° is relatively low, while at 0.0005°, it is relatively higher compared with other regions. Since only the LST and EVI are downscaled to 0.0005° through the revised fusion method for blue, red, NIR, and TIR corresponding bands of MODIS and Landsat, other land surface properties that were applied to retrieve ST profile are at spatial resolutions over 250 m and were resampled to the same spatial resolutions as the downscaled observations. For the 0.0005° spatial resolution, the spatial distribution is not that continuous, especially at the boundaries of GLASS products pixels.

5.3. Increase in ST over the central of TP

Figure 9 illustrates ST’s warming rate at four depths discussed in this study from 2011 to 2020, respectively. The gray line indicates the time series of ST anomaly, and the red line represents the trend of the change of ST from 2011 to 2020. The time series of averaged ST over the study area at each depth for the specific day of the year was generated, and ST anomalies at each depth were calculated as illustrated by Eq. (19) ∼ (20). (19) $S{T}_D={\displaystyle \frac{S{T}_{2011D}+S{T}_{2012D}+\dots +S{T}_{2020D}}{10}}$(19) (20) $\begin{array}{c}ST{A}_D=ST-S{T}_D\end{array}$(20) where $S{T}_D$ is the averaged ST from 2011 to 2020 for the day of the year (DOY) $d$,$S{T}_{YearD}$ is the estimated ST through the methods proposed in this study on DOY $D$ of the specific year, $ST{A}_D$ is the ST anomaly on DOY $D$ of the year. The same goes for the four depths analyzed in this study.

Figure 9. ST’s rising trend at 5, 10, 20, and 40 cm.

Figure 9 has the same annual pattern as Figure 5, indicating the ST of deeper depths is higher, and the ST of the study area gradually increased with time moved on from 2011 to 2020. ST of the shallow depths more varies significantly than the deeper depths. Besides that, as the red line illustrates, there is a trend of increase in ST profile over long-term monitoring from 2011 to 2020, which is more evident at deep layers.

5.4. Limitations and future perspective

Studies focusing on ST variation over the TP typically utilized in-situ observations and reanalysis products, both of which cannot provide high-resolution ST (Hu et al. 2019; Li et al. 2022; Qin et al. 2020; Su et al., 2017; Yang et al. 2020; Zhao et al. 2021; Zhu et al. 2018). Temperature’s warming rate over TP is higher than the global rate and varies across different regions of the TP. Furthermore, the ST profile is an important indicator of permafrost degradation (Li et al. 2014; Li et al. 2022; Liu and Chen 2000; Liu et al. 2017; Ming et al. 2022; Qin et al. 2020; Tan et al. 2017; Wang et al. 2022). Given the importance of the TP in climate change and earth system science, monitoring high-resolution ST at different depths is of great significance. Compared to previous work, this study integrated remote sensing observations, data fusion techniques, and deep learning models to retrieve high spatiotemporal resolution ST estimates at four depths with considerable results. However, there are several limitations exist in this study. 1) The methodology proposed in this study of the same spatial resolutions performs relatively worse for months from January to April than the statistical results of other months. The uncertainty may owe to the ice coverage of the soil surface during winter, which impacts the reflectance/ radiance of the land surface and introduces uncertainties to the remotely sensed datasets applied. Besides that, the active layer thickness and permafrost condition affect the heat flow of soil and influence ST monitoring consequently. Considering the freeze–thaw cycle over the study area may improve the model's performance. 2) Besides that, in previous studies, the layer-wise method is the most widely used method for ST monitoring below the surface at different depths. However, in this study, the LW-DBN performs slightly worse than the Direct-DBN, the freeze–thaw status may accelerate the accumulation of errors at previous layers and contribute to the relatively poor performance of LW-DBN at a deeper depth, while the reason is still unclear yet. 3) There are other factors contributing to ST estimations (e.g. air temperature, humidity, wind speed, and solar radiation) owing to the impacts on the heat flow and energy balance, while these variables were not considered in the methodology proposed in this study due to the lack of high-resolution observations (Alizamir et al. 2020; Paul et al. 2004). 4) The fusing of MODIS and Landsat corresponding channels may also introduce bias, and the uncertainty of downscaled land surface observations impacts the models’ performances for ST retrieval.

The monitoring of ST at different depths plays an important role in Earth system science, helping to predict and understand ecosystem responses to climate change. Temperature’s rising rate over the TP is higher than the global average rate (Dai et al. 2020; Li et al. 2022). Monitoring the spatiotemporal variation of ST at different depths over the TP plays a significant role in earth surface processes and earth system science. ST helps predict and discuss the ecosystem variations to climate change and is an important variable to drive land surface models (Cao et al. 2020; Grünberg et al. 2020; Qin et al. 2020). Besides that, accurate ST estimation over TP is significant to numerous permafrost-related fields, such as thermal state, active-layer thickness, ground ice loss, and carbon release, methods proposed in this study may further be explored and utilized in permafrost science (Li et al. 2022; Yang et al. 2020). ST is considered one of the most essential factors for agricultural management and process, especially for ST at root depth (Alizamir et al. 2020; Lai, Farquharson, and Denton 2019). ST monitoring at different depths can be employed in plant growth and root conditions regulation as well to serve agricultural applications.

6. Conclusion

This study proposed an improved approach to estimating high-resolution ST at different depths with promising results over the central TP by fusing multi-remote sensing observation through the deep learning model. The MODIS and Landsat corresponding channels were fused to generate LST and EVI, which were integrated with ET, NPP, LAI, and FAPAR products and the interpolated SM profile to estimate the ST profile at four depths (5, 10, 20, and 40 cm) through deep learning technologies. The proposed methodologies were analyzed at three spatial scales, and four depths have performances of R² > 0.836 and MAE < 2.152 ${}^{\circ }\mathrm{C}$. The retrieved ST with the coarse resolution performs relatively better (0.0125° > 0.0025° > 0.0005°); the ST profile retrieved at 0.0005° is still with promising results achieving R² from 0.846 to 0.866 and MAE from 1.686 ${}^{\circ }\mathrm{C}$ to 2.152 ${}^{\circ }\mathrm{C}$ for all analyzed depths.

ST profile retrieving is still challenging, especially when the geographical environment is complicated. The methodologies proposed in this study retrieve the ST profile at four depths with high spatiotemporal resolutions. It has been observed from this study that the trend for ST’s rising is more obvious at deeper layers, which is crucial to global warming since the TP (the largest permafrost and high-altitude ecosystem) has been considered one of the most sensitive regions to climate change. ST is an important driving factor of numerous earth surface processes, accurate monitoring of ST over the TP is significant to evaluate and predict climate change. Methods proposed in this study retrieving spatiotemporal continuous ST at different depths over the central TP through remote sensing techniques are essential to permafrost fields and may be utilized in active-layer thickness, ground ice loss, and carbon release applications.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability

The data that support the findings of this study are openly available in NASA MODIS Portal (https://modis.gsfc.nasa.gov/), USGS Landsat Portal (https://earthexplorer.usgs.gov/), the Global land surface satellite products (http://www.glass.umd.edu/index.html), the Tibetan Plateau Science Data Center, National Earth System Science Data Sharing Infrastructure, National Science & Technology Infrastructure of China (http://tibet.geodata.cn).

Additional information

Funding

This work was supported by National Natural Science Foundation of China: [Grant Number 42171400]; National Natural Science Foundation of Guangdong: [Grant Number 2021A1515011324]; Henan Institute of Sun Yat-sen University: [Grant Number 2021-006]; Natural Resources of Guangdong: [Grant Number [2023]-25].