An Assessment of Using Remote Sensing-based Models to Estimate Ground Surface Soil Heat Flux on the Tibetan Plateau during the Freeze-thaw Process

The ground surface soil heat flux (G0) is very important to simulate the changes of frozen ground and the active layer thickness; in addition, the freeze-thaw cycle will also affect G0 on the Tibetan Plateau (TP). As G0 could not be measured directly and soil heat flux is difficult to be observed on the TP in situ due to its high altitude and cold environment, most of previous studies have directly applied existing remote sensing-based models to estimate G0 without assessing whether the selected model is the best one of those models for those study regions. We use in-situ observation data collected at 12 sites combined with Moderate Resolution Imaging Spectroradiometer (MODIS) data (MOD13Q1, MODLT1D, MOD09CMG, and MCD15A2H) and the China meteorological forcing dataset (CMFD-SRad and CMFD-LRad) to validate the main models during the freeze-thaw process. The results show that during the three stages (complete freezing (CF), daily freeze-thaw cycle (DFT), and complete thawing (CT)) of the freeze-thaw cycle, the root mean square error (RMSE) between the models' G0 simulated value and the corresponding G0 "measured value" is the largest in the CT phase and smallest in the CF phase. The simulated results of the second group schemes (SEBAL, Ma, SEBALadj, and Maadj) were slightly underestimated, more stable, and closer to the measured values than the first group schemes (Choudhury, Clawson, SEBS, Choudhuryadj, Clawsonadj, and SEBSadj). The Maadj scheme is the one with the smallest RMSE among all the schemes and could be directly applied across the entire TP. Then, four possible reasons leading to the errors of the main schemes were analyzed. The soil moisture affecting the ratio G0/Rn and the phase shift between G0 and net radiation Rn are not considered in the schemes directly; the scheme cannot completely and correctly capture the direction of G0; and the input data of the schemes to estimate the regional G0 maybe bring some errors into the simulated results. The results are expected to provide a basis for selecting remote sensing-based models to simulate G0 in frozen ground dynamics and to calculate evapotranspiration on the TP during the freeze-thaw process. The scheme Maadj suitable for the TP was also offered in the study. We proposed several improvement directions of remote sensing-based models in order to enhance understanding of the energy exchange between the ground surface and the atmosphere.


Introduction
The Tibetan Plateau (TP) faces the subtropical zone in the south and the middle latitudes in the north. It spans over 25 • longitude from east to west, accounting for about a quarter of China's land area, with an average elevation of over 4000 m [1,2]. The TP, known as the third pole of the earth, affects climate change in important ways. This occurs because of the dynamic action of the high terrain and the thermal influence of the strong heat source of the earth's surface and its upper atmospheric column [3,4]. The energy and water cycle in the TP is one of the main parts of the global energy and water cycles [5,6]. A large area of frozen ground covers some of the TP. With global warming, researchers are very concerned about the effects of climate warming on frozen ground's changes. Understanding the ground surface soil heat flux is one of the starting points needed to clarify these effects.
Measurement of the ground surface soil heat flux (G 0 ) is an important component of calculating the surface energy balance and represents the energy released into or from the soil. Estimating G 0 is a concern for almost all the energy balance analyses done by atmospheric boundary layer observatories and is part of the study of the effects of climate change on frozen ground [7]. Due to the particularity of the plateau climate, a process of daily and seasonal freeze-thaw occurs in the surface soil. This will not only alter the thermal conditions in soil but also cause major differences in soil heat conduction and in the energy exchange between the soil and the nearby upper atmosphere during the state of freezing and ablation. The thermal conductivity and heat capacity of frozen ground are quite different from those of the ablation soil [8][9][10]. The seasonal freeze-thaw cycles of near-surface soil significantly affect the energy balance of the surface. The freezing and thawing of the soil on the TP reflect the changes in the water-heat exchange between the plateau land and atmosphere. The existence of a thin frozen layer will fundamentally affect the direct water-heat exchange between the atmosphere and deep soil [11]. As the energy exchange between soil and atmosphere mainly occurs at the interface layer, the freezing and thawing of shallow soil is particularly important when studying G 0 on the TP.
Measurements of G 0 can be obtained using the actual observed soil temperature, soil moisture, and shallow soil heat flux data at experimental sites [12], but soil temperature and moisture data cannot be obtained directly through optical remote sensing imagery [13,14]. The basic idea of estimating G 0 using remote sensing-based models is based on the fact that if the ratio of G 0 to net radiation (R n ) can be acquired, then we can employ the ratio to estimate G 0 using remote sensing data. Different forms of remote sensing-based models used for estimating G 0 come from the estimation of G 0 /R n using different land surface parameters. Some use only vegetation data to estimate this ratio, such as the NDVI-based (normalized difference vegetation index) scheme of Clawson [15,16], the LAI-based (leaf area index) scheme of Choudhury [17], and f c -based (fractional canopy coverage) scheme known as the surface energy Balance System (SEBS) [18]. Others also consider land surface temperature and surface albedo in addition to a vegetation index, such as the schemes of the Surface Energy Balance for Land (SEBAL) [13] and Ma [19]. Ma et al. [19,20] proposed the use of modified soil adjusted vegetation index (MSAVI) in arid and semi-arid as well as in high altitudes regions, instead of NDVI. Therefore, the scheme of Ma is similar to the scheme of SEBAL. It is important to use an accurate ratio (G 0 /R n ) to accurately estimate G 0 ; most methods directly introduce a certain scheme to estimate G 0 and lack verification of the model, which may lead to some estimated errors caused by inappropriate model selection. Meanwhile, a few studies [21,22] have tested and verified remote sensing-based models; the evaluations of these models have mainly been concentrated in the complete melting stage (in summer) during the freeze-thaw cycle, which also has a great influence on simulating G 0 [23]. Therefore, it can be seen that analyzing the regional adaptability of remote sensing-based models and determining the appropriate model in a specific region is a very important step in G 0 estimation and research using remote sensing data during the freeze-thaw cycle across the TP [24].
In-situ measurement data of 12 observed sites, Moderate Resolution Imaging Spectroradiometer (MODIS) data (specifically MOD13Q1, MODLT1D, MOD09CMG, and MCD15A2H), and the China meteorological forcing dataset (CMFD-SRad and CMFD-LRad) were applied in this study. We adjusted Remote Sens. 2020, 12, 501 3 of 24 and evaluated the applicability of five remote sensing-based models estimating G 0 (Choudhury, Clawson, SEBS, SEBAL, and Ma) during the freeze-thaw cycle on the TP. Additionally, the sensitivity of these schemes to land surface parameters (land surface temperature, land surface albedo, and vegetation indices) was analyzed. We also discussed the possible sources of error in the models. It can provide some substantial evidence for improving remote sensing-based models to estimate G 0 on the TP.

Study Regions
Observational data were collected from 12 automatic weather stations listed in Figure 1. This background map (Figure 1) shows the spatial distribution of frozen ground on the TP [25]. The Binggou site is located in the alpine meadow valley of the Binggou sub-basin in the upper reaches of the Heihe River. The east and west sides are continuous hills and mountains. The surrounding area is relatively flat and covered by a mixture of sparse grassland and riverbed gravel. The AYKMS site is located at the sparse desert grasslands of the boundary region between the permafrost and the seasonally frozen ground of the TP. The XDTMS sparse desert grasslands lie in the alpine steppe of the permafrost region near the northern lower limit of the permafrost on the TP. The TGLMS site is located in the alpine meadow of the predominantly continuous permafrost region. The observed area around the BGGRASS site is relatively flat and covered with alpine grassland, which is 2-4 cm tall. Alpine grassland is the main land surface type of Biru site. The Jiali site falls within the sub-frigid zone and semi-humid monsoon climate regions of the Plateau, with sufficient sunshine and strong radiation. The Naqu site is located in the alpine grassland of Naqu Valley and lies within the sub-frigid zone and semi-humid climate regions of Plateau [26]. The Nierong Site is affected by the terrain; the cold and dry climate falls within the sub-arid and semi-arid monsoon climate zone of Plateau. The terrain around the site is broad, expansive, and relatively flat and supports alpine meadow [27]. The Linzhou site is located in a valley surrounded by relatively flat of the alpine desert. Shiquanhe is located in the western part of the TP and has sparse and short grassland that is generally only a few centimeters tall. This part of the plateau experiences a cold zone monsoon arid climate. The solar radiation is strong with frequently strong surface wind speeds; the air is dry with less rain than other areas and experiences large daily temperature differences [28]. Amdo station has a relatively smooth surface with a very open view where sparse plateau meadow covers the ground.    Figure 1. Distribution of study sites on the Tibetan Plateau (source of permafrost distribution data: [25]) For abbreviations of site names and further details see Table 1.

In-Situ Measurement Data
The above 12 sites used in this study provided data for net radiation (Rn, downward/upward shortwave and longwave radiation), precipitation, soil temperature and humidity profiles, land surface temperature, and soil heat flux ( Table 1). The data acquisition interval was 30 min. The observed data of the Binggou site were acquired from the Cold and Arid Regions Sciences Data Center (http://westdc.westgis.ac.cn/data/30617ca4-f46e-4254-a40a-5d0fbea58ed2). Meanwhile, the data for the AYKMS, TGLMS, and XDTMS sites were provided by the Cryosphere Research Station on the Qinghai-Tibet Plateau, Chinese Academy of Sciences (http://www.crs.ac.cn/), while the measured data of the other eight sites were downloaded from the National Meteorological Information Center (http://tipex.data.cma.cn/tipex). Precipitation data were used to eliminate data  Table 1.

In-Situ Measurement Data
The above 12 sites used in this study provided data for net radiation (R n , downward/upward shortwave and longwave radiation), precipitation, soil temperature and humidity profiles, land surface temperature, and soil heat flux ( Table 1). The data acquisition interval was 30 min. The observed data of the Binggou site were acquired from the Cold and Arid Regions Sciences Data Center Remote Sens. 2020, 12, 501 6 of 24 (http://westdc.westgis.ac.cn/data/30617ca4-f46e-4254-a40a-5d0fbea58ed2). Meanwhile, the data for the AYKMS, TGLMS, and XDTMS sites were provided by the Cryosphere Research Station on the Qinghai-Tibet Plateau, Chinese Academy of Sciences (http://www.crs.ac.cn/), while the measured data of the other eight sites were downloaded from the National Meteorological Information Center (http: //tipex.data.cma.cn/tipex). Precipitation data were used to eliminate data based on unusual weather. Soil heat flux, soil temperature, and soil moisture data were used to calculate G 0 . Downward/upward shortwave radiation was used to calculate the land surface albedo in addition to R n . In addition, R n , land surface temperature, and land surface albedo acted as input parameters for the remote sensing-based models in estimating G 0 .
Measured land surface temperature data (T s ) were unavailable for seven stations (Binggou, Amdo, BGGRASS, Biru, Jiali, Naqu, and Nierong), but the in-situ LST measurements were estimated from the upwelling and downwelling longwave radiation data using the Stefan-Boltzmann law [29,30]: where T s is the land surface temperature (unit: • C), R L ↓ and R L ↑ are the incoming and outgoing longwave radiation powers, respectively (unit: W/m 2 ), σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W·m −2 ·K −4 ), and ε s is the land surface emissivity, which is approximately 0.987 at Binggou [31] and 0.98 at the other six sites [32].

Remote Sensing Data (MODIS) and China Meteorological Forcing Dataset
Vegetation indices used in the study include NDVI, modified soil adjusted vegetation index (MSAVI), f c , and LAI. These data were not observation-based but retrieved from remote sensing data (Table 2). MODIS product MCD15A2H including LAI is an 8-day composite dataset with 500 m pixel size. The algorithm chooses the "best" pixel available from all the acquisitions of both MODIS sensors located on NASA's Terra and Aqua satellites from within the 8-day period [33]. The NDVI data were obtained from the MODIS vegetation index product MOD13Q1; these data were produced on 16-day intervals and at 250 m spatial resolution. The modified soil adjusted vegetation index (MSAVI) was also calculated from the product MOD13Q1 by Equation (2) [34]: where RED and NIR are land surface reflectance values of RED and NIR bands of MOD13Q1, respectively.
The fractional canopy coverage f c can be obtained by NDVI [35]: where NDVI max and NDVI min are the NDVI values for bare soil and full vegetation, respectively. The product MODLT1D was made using MOD11A1 through splicing, cutting, and projection conversion (http://www.gscloud.cn) with 1 km spatial resolution, provided along with the daytime and nighttime surface temperature bands associated using quality control assessments. The MOD09CMG product provides an estimate of the surface spectral reflectance of MODIS Bands 1 through 7, resampled to 5600 m pixel resolution and corrected for atmospheric conditions such as gasses, aerosols, and Rayleigh scattering. China meteorological forcing dataset (CMFD, https://data.tpdc.ac.cn) developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences, was used as the model input in this study. The dataset has merged the observations from 740 operational stations of the China Meteorological Administration with the corresponding Princeton Global Meteorological Forcing Dataset [36][37][38]. The parameters used in this study included downward shortwave (SRad) and downward longwave (LRad) radiation. These two parameters have a 10 km spatial resolution and a temporal resolution of 3 h ( Table 2). All parameters provided or retrieved by MODIS data were converted to the same projection coordinate as CMFD datasets and resampled to 10 × 10 km. The transit time of the Terra satellite was local time 10:30 a.m., and TP was in the range of 75 • E to 105 • E (Beijing standard time 11:30-13:30); we chose the average values of CMFD corresponding to Beijing standard time 11:00 and 14:00. Figure 2 shows the general steps used for deriving the G 0 in this paper. The green box demonstrates the process of calculating G 0 at sites based on the measured data (details in Section 3.1). The content in the purple irregular wireframe shows the process of adjusting and selecting the G 0 parameterization schemes (Section 3.2), and that in the red irregular wireframe presents the steps for estimating G 0 over the TP by the selected scheme using MODIS data and forcing datasets (Section 3.3). Figure 2 shows the general steps used for deriving the G0 in this paper. The green box demonstrates the process of calculating G0 at sites based on the measured data (details in Section 3.1). The content in the purple irregular wireframe shows the process of adjusting and selecting the G0 parameterization schemes (Sections 3.2), and that in the red irregular wireframe presents the steps for estimating G0 over the TP by the selected scheme using MODIS data and forcing datasets (Section 3.3).

Estimate G0 at a Site
The G0 estimation started with the soil 1-D heat conduction equation [7]: where ρs is the soil density (unit: kg/m 3 ), cs is the soil specific heat (unit: J•kg −1 •K −1 ), T is the soil temperature (unit: K), t is time (unit: s), z is the soil depth (unit: m, the downward direction is taken as positive), ks is the soil heat conduction coefficient (unit: W•K −1 •m −1 ), and G is the soil heat flux (unit: W/m 2 , the direction down here is positive).
By integrating Equation (4), the G0 is obtained: where zref is the reference depth (it was 0.10 m in this study), and G(zref) is the soil heat flux at the depth of the reference location.
For sites on the TP, there is a freeze-thaw process in the soil near the surface, so the ice thermal capacity should also be considered [7,39], so the discrete form of Equation (6) can be expressed as

Estimate G 0 at a Site
The G 0 estimation started with the soil 1-D heat conduction equation [7]: where ρ s is the soil density (unit: kg/m 3 ), c s is the soil specific heat (unit: J·kg −1 ·K −1 ), T is the soil temperature (unit: K), t is time (unit: s), z is the soil depth (unit: m, the downward direction is taken as positive), k s is the soil heat conduction coefficient (unit: W·K −1 ·m −1 ), and G is the soil heat flux (unit: W/m 2 , the direction down here is positive). By integrating Equation (4), the G 0 is obtained: where z ref is the reference depth (it was 0.10 m in this study), and G(z ref ) is the soil heat flux at the depth of the reference location. For sites on the TP, there is a freeze-thaw process in the soil near the surface, so the ice thermal capacity should also be considered [7,39], so the discrete form of Equation (6) can be expressed as where θ 0.05m is the unfrozen water content at 0.05 m depth (unit: m 3 ·m −3 ), ρ dry c dry is the thermal capacity of dry soil (about 0.90 × 10 6 J·m −3 ·K −1 ) [40], ρ w c w is the thermal capacity of liquid water Remote Sens. 2020, 12, 501 9 of 24 (4.2 × 10 6 J·m −3 ·K −1 ), and ρ i c i is the ice thermal capacity (1.89 × 10 6 J·m −3 ·K −1 ). There were three different stages for the calculation of the soil ice content at 5-cm depth ϑ i,0.05m : (1) When the soil layer above 0.05 m was completely thawed, no ice existed in the layer (ϑ i,0.05m = 0); (2) When the daily freeze-thaw stage of the soil layer was above 0.05 m, the soil ice content ϑ i,0.05m was approximately obtained by the difference of the steady unfrozen water content before and after the daily freeze-thaw stage of this soil layer. ϑ i,0.05m can be calculated using Equation (8) [41]: (3) When the soil layer above 5 cm was completely frozen, ϑ i,0.05m equals the soil ice content of the last day of the daily freeze-thaw stage. It should be noted that the ice content was variable even if the soil was frozen, but such ice content variation was small in the frozen stages, so the ice content was approximately regarded as a constant during the frozen stages in the calculations [39]. The daily freeze-thaw stages were estimated based on the soil temperatures T 0.05m and unfrozen water contents θ 0.05m [39].

Introduction and Adjustments to G 0 Parameterization Schemes
Several empirical remote sensing models can be used to estimate Γ(G 0 /R n ), such as Choudhury [17], Clawson [15], SEBAL [13], Ma [19], and SEBS [18]. Choudhury et al. [17] expressed the relationship (the Choudhury scheme) between G 0 and R n as an exponential function of the LAI with a correlation coefficient of 0.93: A detailed description of the other four schemes is presented in Table 3 of Yang et al. [23]. It should be noted that the Moran scheme is named the Clawson scheme in this study, which uses an exponential function of NDVI to the ratio of G 0 /R n [15,16]. In addition, Yang et al. [23] found that the ratio Γ s of soil heat flux to net radiation for bare soil in the SEBS scheme should be less than 0.315 on the TP, which needs to be adjusted later. Table 3. The original and adjusted versions of remote sensing-based models estimating G 0 /R n .

Model Name Mathematical Form References
This study Choudhury adj G 0 /R n = 0.267 · e 0.27·LAI This study 13·NDVI [42] Clawson adj G 0 /R n = 0.238 · e 0.78·NDVI This study This study evaluates the applicability of each scheme over the TP from their errors and their sensitivity to land surface parameters during the freeze-thaw process.
The accuracies of the values produced by the models were evaluated using the root mean square error (RMSE) and mean bias error (MBE), defined as follows: where P j is the modeled value, G j is the measured value, and m is the number of samples. We can choose the best one from the 10 models based on the RMSE and MBE. These remote sensing-based schemes estimate G 0 based on the proportional relationship between G 0 and R n . The difference between the schemes is that some models only use a vegetation index to estimate G 0 /R n , while others use land surface temperature and albedo as well as a vegetation index to estimate G 0 /R n . In addition, the net radiation R n was obtained through Equation (13). The mathematical formulas of the schemes (Table 3) and Equation (13) were merged. It can be seen that their sensitivities to land surface temperature T s , land surface albedo α, and the vegetation index VI (LAI, NDVI, MSAVI, and f c ) can be derived from the percentage change of simulated G 0 values when increasing or reducing the observed T s by 1 K, the measured α by 0.02, or VI by 0.1 [43]: |G−accuracy−Gbase| N VR accuracy = max VR +accuracy , VR −accuracy (12) where VR +accuracy and VR −accuracy indicate the average change of simulated G 0 values by increasing and reducing the "observed" T s by 1 K or the measured α by 0.02, respectively; VR accuracy is the ultimate sensitivity coefficient; G base is the simulated G 0 value from observed data corresponding to the unchanged land surface temperature T s , unchanged land surface albedo α, and unchanged vegetation index. Meanwhile, G +accuracy and G −accuracy represent the simulated G 0 value when the "observed" T s is increased or decreased by 1 K or the measured α is increased or decreased by 0.02 or remotely sensed VI is increased or decreased by 0.1, respectively. N is the number of samples.

Estimating Regional G 0 on the TP
We used the evaluation indices (RMSE and MBE) in Section 3.2 to select the optimal scheme on the TP and then obtain all the parameters by remote sensing data that were needed for the optimal scheme; these were input into the scheme to obtain the spatial distribution characteristics of G 0 on the TP. All schemes required the use of net radiation, land surface albedo, surface temperature, and vegetation index as inputs; the remote sensing inversion algorithm for the vegetation index was introduced in Section 2.2.2. Here, we only introduce the methods for obtaining the remaining land surface parameters by remote sensing and reanalyzed data.
The net radiation R n can be obtained through Equation (13) [44]: where R S ↑ and R S ↓ are the upwelling and downwelling shortwave radiation (W·m −2 ), respectively, and R L ↓ and R L ↑ are the downwelling and upwelling longwave radiation (W·m −2 ), respectively. The product MODLT1D can provide T s (the scale-factor is 0.02 and the unit is Kelvin). The instantaneous value of land surface albedo α in the time range from 10:00 to 16:00 is approximately equal to the daily average value of land surface albedo α on the same day [23]. In addition, the land surface instantaneous albedo α can be derived from bands 1-5 and band 7 data in the MODIS product MOD09CMG and an empirical Equation (14) [45]: The errors in MODIS land surface temperature measurements and land surface albedo estimated using Equation (14) and MOD/MYD09, are 1 K [46] and 0.02 [47], respectively. ε s is the land surface emissivity, which was calculated using Equation (15) [48]:

Errors of the Schemes during the Freeze-Thaw Process
The errors of the models were evaluated by RMSE and MBE using the data from four sites (TGL, Naqu, Biru, and Nierong) of the verification group ( Figure 3). The RMSE values of the 10 schemes (five original schemes and five corresponding adjusted schemes) were recorded in the three stages of a freeze-thaw cycle. We can see that, overall, the RMSE values of the adjusted schemes based on the measured data of the eight sites are smaller than those of the original schemes. This is particularly obvious for the Choudhury scheme. The RMSE of the adjusted Choudhury adj scheme is significantly smaller than the original Choudhury scheme. This phenomenon once again shows that if measured data are used, it is best to re-determine the empirical coefficients in a remote sensing-based model based on the measured data and then use it to conduct related research.
The RMSE of each scheme was significantly different throughout the freeze-thaw cycle ( Figure 3a); the RMSE values of the three stages are roughly presented as a whole: CF < DFT < CT. The order may depend on whether a particular scheme is better able to correctly capture the direction of G 0 when G 0 and R n are different in the three phases. When the direction of G 0 and R n is different (Table 6), the data volume of the three stages is CT > DFT ≈ CF, and the amount of analog data that can correctly capture the G 0 direction is CF > DFT > CT, which directly leads to the smallest RMSE in the CF stage and the largest RMSE in the CT stage. The MBE of the CT stage, whether overestimated or underestimated, is less than that of the other two stages, while the MBE values of the other two stages are near to each other (Figure 3b). when G0 and Rn are different in the three phases. When the direction of G0 and Rn is different (Table  6), the data volume of the three stages is CT > DFT ≈ CF, and the amount of analog data that can correctly capture the G0 direction is CF > DFT > CT, which directly leads to the smallest RMSE in the CF stage and the largest RMSE in the CT stage. The MBE of the CT stage, whether overestimated or underestimated, is less than that of the other two stages, while the MBE values of the other two stages are near to each other (Figure 3b). For convenience, we divided the schemes into two groups according to the factors considered in them. The first group of schemes simulates G0/Rn value using a vegetation index alone (Choudhury, Clawson, SEBS, Choudhuryadj, Clawsonadj, and SEBSadj). Meanwhile, the second group schemes are based on Ts, α, and a vegetation index (NDVI or MSAVI) to estimate the G0/Rn value (SEBAL, Ma, For convenience, we divided the schemes into two groups according to the factors considered in them. The first group of schemes simulates G 0 /R n value using a vegetation index alone (Choudhury, Clawson, SEBS, Choudhury adj , Clawson adj , and SEBS adj ). Meanwhile, the second group schemes are based on T s , α, and a vegetation index (NDVI or MSAVI) to estimate the G 0 /R n value (SEBAL, Ma, SEBAL adj , and Ma adj ). An individual one of the first group of schemes was slightly better in a certain stage of the freeze-thaw cycle, but overall, the RMSE values of the second group schemes was smaller ( Figure 3a). Overall, the first group schemes simulations were obviously overestimated, and those of another group were weakly underestimated (Figure 3b). This study considers that the Ma adj scheme is more suitable for estimating G 0 on the TP. In addition, MSAVI is more suitable for arid and semi-arid or high-elevation areas than NDVI [19,20,49]. The eight sites that determined the scheme not only are very evenly distributed on the TP but also include the main land surface and frozen ground types, which can represent the overall situation of the TP; the RMSEs of the schemes in any of the three stages (CF, DFT, and CM) were also relatively close with a minimum of difference. If no measured data are available, the Ma adj scheme can be considered for direct use on the TP.

Regional G 0 on the TP
This study mainly focused on the applicability of the G 0 parameterization scheme on the TP during the freeze-thaw cycle. Due to data limitations, this study selected the day with the largest proportion of valid data per month during a four-month period (July and October 2014, January and April 2015) to show the simulated G 0 by the scheme Ma adj during the freeze-thaw cycle. We selected the four days (July 24, 2014, October 25, 2014, January 21, 2015, and April 26, 2015) according to the quality control documents of the remote sensing land surface temperature product MODLT1D (Figure 4). The CMFD downward shortwave (CMFD-SRad) and longwave (CMFD-LRad) radiation data have proven to have reasonable accuracy with the in-situ measurements over the TP, which means that the two CMFD datasets can be used as model inputs [38]. The method to estimate net radiation R n over the TP by regional data (MODIS products and CMFD datasets) is described in detail in Section 3.3. Although some data are missing, it can roughly reflect the seasonal differences of land surface temperature ( Figure 5), net radiation (Figure 6), and ground surface soil heat flux (Figure 7) on the TP, which are higher in July, are lower in October and January, and become high again in April. means that the two CMFD datasets can be used as model inputs [38]. The method to estimate net radiation Rn over the TP by regional data (MODIS products and CMFD datasets) is described in detail in Section 3.3. Although some data are missing, it can roughly reflect the seasonal differences of land surface temperature (Figure 5), net radiation (Figure 6), and ground surface soil heat flux (Figure 7) on the TP, which are higher in July, are lower in October and January, and become high again in April.        Due to the absence of both remote sensing and measured data and consideration of a spatial homogeneity (10 km), only 11 measured G0(Rn) values corresponding to the simulated G0(Rn) values were available. As shown in Figure 9a,b, the estimates of net radiation and ground surface soil heat flux show reasonable agreement with the in-situ measurements. The RMSEs for the net radiation flux and soil heat flux are 114.43 and 45.99 W/m 2 . The RMSE (45.99 W/m 2 ) obtained from the regional data (MODIS and CMFD) as the input of the scheme Maadj is close to that (Figure 3a) obtained from the measured data as the input of the scheme Maadj. The slope of the fitting line between "observed" and simulated G0 values was slightly less than but very close to 1 (0.83; Figure 8b), This tells us there was a slight underestimation of the simulated results; The coefficient of determination R 2 is 0.87. These indicators (RMSE, MBE, and R 2 ) indicate that satisfactory results can be obtained by using the Maadj scheme. It should be noted that some bias exists between the estimated G0 and ground measurements. Due to the absence of both remote sensing and measured data and consideration of a spatial homogeneity (10 km), only 11 measured G 0 (R n ) values corresponding to the simulated G 0 (R n ) values were available. As shown in Figure 9a,b, the estimates of net radiation and ground surface soil heat flux show reasonable agreement with the in-situ measurements. The RMSEs for the net radiation flux and soil heat flux are 114.43 and 45.99 W/m 2 . The RMSE (45.99 W/m 2 ) obtained from the regional data (MODIS and CMFD) as the input of the scheme Ma adj is close to that (Figure 3a) obtained from the measured data as the input of the scheme Ma adj . The slope of the fitting line between "observed" and simulated G 0 values was slightly less than but very close to 1 (0.83; Figure 8b), This tells us there was a slight underestimation of the simulated results; The coefficient of determination R 2 is 0.87. These indicators (RMSE, MBE, and R 2 ) indicate that satisfactory results can be obtained by using the Ma adj scheme. It should be noted that some bias exists between the estimated G 0 and ground measurements. The causes of the scheme errors are analyzed in detail in Section 5.2.

The Sensitivity of Remote Sensing-Based Schemes to Land Surface Parameters
According to the sensitivity analysis indexes in Section 3.2.2, the sensitivity values of three stages of freeze-thaw process of each scheme can be determined. The Choudhury scheme is used as an example to illustrate the sensitivity of this scheme to three surface parameters (Table 4). There are three cases for each land surface parameter: base + accuracy, base, and base − accuracy, so each scheme has 26 variations (3 3 -1 = 26, three to the third power indicates that there are three changes in each land surface parameter (+ accuracy, base, and − accuracy), and three land surface parameters (Ts, α, or VI) are considered in common; subtracting 1 is the case when three land surface parameters are all original values), where the largest VR value is used as VRaccuracy of the scheme.   Figure 8a,b represent the sites and dates corresponding to the red circles, shown in the table (c)). Note: RMSE, root mean square error; MBE, mean bias error; R 2 , the coefficient of determination.

The Sensitivity of Remote Sensing-Based Schemes to Land Surface Parameters
According to the sensitivity analysis indexes in Section 3.2.2, the sensitivity values of three stages of freeze-thaw process of each scheme can be determined. The Choudhury scheme is used as an example to illustrate the sensitivity of this scheme to three surface parameters (Table 4). There are three cases for each land surface parameter: base + accuracy, base, and base − accuracy, so each scheme has 26 variations (3 3 -1 = 26, three to the third power indicates that there are three changes in each land surface parameter (+ accuracy, base, and − accuracy), and three land surface parameters (T s , α, or VI) are considered in common; subtracting 1 is the case when three land surface parameters are all original values), where the largest VR value is used as VR accuracy of the scheme.
The same idea is to get VR accuracy of each scheme (Table 5); we know that, overall, these schemes that are based on T s , α, and VI are more stable than others based on only VI. That is, when remotely sensed products are used as input parameters of remote sensing-based schemes, the results simulated by the schemes based on T s , α, and VI are more stable. In addition, VR accuracy values of the adjusted schemes (Choudhury adj , Clawson adj , and SEBS adj ) based on the vegetation index are significantly reduced (Table 5), which tells us their simulated results are more stable than those of their original schemes (Choudhury, Clawson, and SEBS). The VR accuracy values of the two adjusted schemes (SEBAL adj and Ma adj ) are larger than that of the two original schemes (SEBAL and Ma), although the difference is still no more than 2 W/m 2 , which shows that the two adjusted schemes are more sensitive to land surface parameters within acceptable limits. Note: The sign "+" ("−") before T s , α, or LAI represents that the parameter is increased (decreased) by 1 K, 0.02, or 0.1, respectively. N represents the number of samples. CF, DFT, and CT represent complete freezing stage, daily freeze-thaw cycle stage, and complete thawing stage of a freeze-thaw cycle, respectively.

Analysis of the Causes of the Scheme Errors
Through a review of previous research [13,15,[17][18][19]23,38] and the above findings of the present study, we knew that remote sensing-based models can provide satisfactory simulation results, and their application is convenient and simple. However, some errors do occur in the modeled results. Here, we discuss the sources of error from the following four aspects: Soil moisture has a certain effect on the G 0 /R n ratio. The existing studies confirmed that the variability in G 0 /R n on hourly timescales is partly caused by variations in soil moisture [50][51][52][53]. This dependence arises from the competing effects of changing evaporation rates and soil thermal conductivity as the soil dries [54]. The drier the soil, the larger the G 0 /R n [24]. In addition, a significant difference exists in soil moisture during a freeze-thaw cycle, with the smallest soil moisture in the CF stage and the largest in the CT stage ( Figure 9). However, none of these schemes considered the effects of soil moisture directly, which will result in the production of some errors. Although T s is sensitive to soil moisture [24], it cannot completely replace soil moisture. From this point of view, the schemes considering T s also have more advantages.
Remote Sens. 2020, 12, 501 18 of 25 Soil moisture has a certain effect on the G0/Rn ratio. The existing studies confirmed that the variability in G0/Rn on hourly timescales is partly caused by variations in soil moisture [50][51][52][53]. This dependence arises from the competing effects of changing evaporation rates and soil thermal conductivity as the soil dries [54]. The drier the soil, the larger the G0/Rn [24]. In addition, a significant difference exists in soil moisture during a freeze-thaw cycle, with the smallest soil moisture in the CF stage and the largest in the CT stage ( Figure 9). However, none of these schemes considered the effects of soil moisture directly, which will result in the production of some errors. Although Ts is sensitive to soil moisture [24], it cannot completely replace soil moisture. From this point of view, the schemes considering Ts also have more advantages. The phase shift between G0 and Rn is also one of the sources of errors in all schemes. Some studies have shown that a phase shift exists between G0 and Rn [23,24,[54][55][56]. Ground surface soil heat flux and net radiation flux do not have the same diurnal variation shape. The soil heat flux peak values are usually later than the net radiation flux peak values, which was not taken into account in these parameterizations. In the present study, taking the measured data of TGL station as an example, no matter which stage of the freeze-thaw process was analyzed, an obvious phase shift occurred between G0 and Rn. That is, the peak value of diurnal variation in G0 always arrived later than the peak value of diurnal variation in Rn ( Figure 10). Relevant studies have also been conducted on the TP, which The phase shift between G 0 and R n is also one of the sources of errors in all schemes. Some studies have shown that a phase shift exists between G 0 and R n [23,24,[54][55][56]. Ground surface soil heat flux and net radiation flux do not have the same diurnal variation shape. The soil heat flux peak values are usually later than the net radiation flux peak values, which was not taken into account in these parameterizations. In the present study, taking the measured data of TGL station as an example, no matter which stage of the freeze-thaw process was analyzed, an obvious phase shift occurred between G 0 and R n . That is, the peak value of diurnal variation in G 0 always arrived later than the peak value of diurnal variation in R n ( Figure 10). Relevant studies have also been conducted on the TP, which concluded that the scheme considering the phase shift between G 0 and R n does have a lower estimated error than before in the CT stage [23]. concluded that the scheme considering the phase shift between G0 and Rn does have a lower estimated error than before in the CT stage [23]. These remote sensing-based models are all based on G0/Rn, so when the G0 and Rn directions are opposite, the simulation results of each scheme are not very good. Table 6 shows that when G0 has an opposite trend to that of Rn, the vegetation index-based scheme basically cannot correctly capture the direction of G0. In contrast, the schemes that are based on Ts, α, and VI can reflect the direction of part of G0 correctly and were used in the three stages of the freeze-thaw cycle. Here, G0 manifests the energy exchange between the active layer and the atmosphere while it affects the energy budget of the active layer and permafrost, so its direction is crucial for the active layer thickness and permafrost [23]. When the direction of G0 is positive (the downward direction of G0 is positive), this finding indicates that the active layer absorbs heat from the atmosphere [56]. Over time, soil temperature will rise, which will increase the thickness of active layer and promote the degradation of permafrost. Vice versa, the active layer releases heat into the atmosphere [56], and the decrease of soil temperature is beneficial to the preservation and development of permafrost. Considering the correctness of the G0 direction, although all schemes cannot fully capture the direction of G0, a scheme based on Ts, α, and VI is better than a scheme based on vegetation index alone. Table 6. Simulation effect of each scheme on the direction of ground surface soil heat flux (G0) when G0 is the opposite direction from net radiation (Rn).

Scheme CT DFT CF
G0 is the opposite direction from Rn 5102 2117 2270 Under the premise that G0 is opposite to Rn, the simulated G0 direction is consistent with the "measured" G0  These remote sensing-based models are all based on G 0 /R n , so when the G 0 and R n directions are opposite, the simulation results of each scheme are not very good. Table 6 shows that when G 0 has an opposite trend to that of R n , the vegetation index-based scheme basically cannot correctly capture the direction of G 0 . In contrast, the schemes that are based on T s , α, and VI can reflect the direction of part of G 0 correctly and were used in the three stages of the freeze-thaw cycle. Here, G 0 manifests the energy exchange between the active layer and the atmosphere while it affects the energy budget of the active layer and permafrost, so its direction is crucial for the active layer thickness and permafrost [23]. When the direction of G 0 is positive (the downward direction of G 0 is positive), this finding indicates that the active layer absorbs heat from the atmosphere [56]. Over time, soil temperature will rise, which will increase the thickness of active layer and promote the degradation of permafrost. Vice versa, the active layer releases heat into the atmosphere [56], and the decrease of soil temperature is beneficial to the preservation and development of permafrost. Considering the correctness of the G 0 direction, although all schemes cannot fully capture the direction of G 0 , a scheme based on T s , α, and VI is better than a scheme based on vegetation index alone. Table 6. Simulation effect of each scheme on the direction of ground surface soil heat flux (G 0 ) when G 0 is the opposite direction from net radiation (R n ).

Scheme CT DFT CF
G 0 is the opposite direction from R n 5102 2117 2270 Under the premise that G 0 is opposite to R n , the simulated G 0 direction is consistent with the "measured" G 0 Remote sensing data will result in certain errors as presented in Table 5. Although the table shows the uncertainty of each scheme, it also shows the influence of an increase or decrease in the accuracy of remotely sensed land surface temperature, surface albedo, and vegetation indices on the simulated G 0 ; therefore, the accuracy of remotely sensed products has a certain effect on the simulation results of the main schemes.

Conclusions
Based on the observed TP and vegetation index data of 12 sites obtained from remote sensing products, we modified five remote sensing-based models to estimate the ground surface soil heat flux on the Tibetan Plateau during the freeze-thaw process and obtained 10 schemes (five original schemes and five modified schemes). We further analyzed RMSE values of 10 schemes, their sensitivity to land surface parameters, and the sources of their errors. The following conclusions are drawn.
(1) The RMSE of each scheme was significantly different during the freeze-thaw cycle on the TP; meanwhile, the RMSE values of the three stages can be roughly presented as a whole: CF < DFT < CT. The MBE of the CT stage, whether overestimated or underestimated, was less than that of the other two stages, while the MBE values of the other two stages were close. Overall, RMSE values of the second group of schemes (SEBAL, Ma, SEBAL adj , and Ma adj ) were smaller than that of the first group of schemes (Choudhury, Clawson, SEBS, Choudhury adj , Clawson adj , and SEBS adj ). The MBE values tell us the G 0 values simulated by the first group of schemes were obviously overestimated, and those of the second group of schemes are weakly underestimated.
(2) In the absence of measured data, the Ma adj scheme can be considered for estimating G 0 on the TP; its RMSE and MBE were the least among all the schemes. (3) AV accuracy of the second group of schemes was less and their simulated results were relatively more stable than those of the first group of schemes during the freeze-thaw process on the TP. (4) We discussed four possible reasons for the errors of the main schemes. Soil moisture, the phase difference between G 0 and R n , whether each main scheme can correctly capture the propagation direction of G 0 , and the accuracy of remote sensing data resulted in errors for the scheme analog results.  (

3) SEBSadj
As there is no site with full vegetation coverage, only the ratio Γs of G0 to Rn in the SEBS scheme can be adjusted in this study. It can be 0.20 in bare soil on the TP ( Figure A2). Additionally, this ratio (0.20) was close to 0.25 when obtained by Yang et al. [23]; therefore, the adjusted form of SEBSadj of the SEBS scheme can be expressed as Equation (A7): (A7) Figure A2. Adjustment scatter plot for the ratio Γs of ground surface soil heat flux G0 to net radiation Rn in bare soil on the TP (R, the correlation coefficient; N, the number of samples.). (

3) SEBS adj
As there is no site with full vegetation coverage, only the ratio Γ s of G 0 to R n in the SEBS scheme can be adjusted in this study. It can be 0.20 in bare soil on the TP ( Figure A2). Additionally, this ratio (0.20) was close to 0.25 when obtained by Yang et al. [23]; therefore, the adjusted form of SEBS adj of the SEBS scheme can be expressed as Equation (A7): (

3) SEBSadj
As there is no site with full vegetation coverage, only the ratio Γs of G0 to Rn in the SEBS scheme can be adjusted in this study. It can be 0.20 in bare soil on the TP ( Figure A2). Additionally, this ratio (0.20) was close to 0.25 when obtained by Yang et al. [23]; therefore, the adjusted form of SEBSadj of the SEBS scheme can be expressed as Equation (A7): (A7) Figure A2. Adjustment scatter plot for the ratio Γs of ground surface soil heat flux G0 to net radiation Rn in bare soil on the TP (R, the correlation coefficient; N, the number of samples.). Figure A2. Adjustment scatter plot for the ratio Γ s of ground surface soil heat flux G 0 to net radiation R n in bare soil on the TP (R, the correlation coefficient; N, the number of samples.).