Estimating Surface Soil Heat Flux in Permafrost Regions Using Remote Sensing-Based Models on the Northern Qinghai-Tibetan Plateau under Clear-Sky Conditions

The ground surface soil heat flux (G0) quantifies the energy transfer between the atmosphere and the ground through the land surface. However; it is difficult to obtain the spatial distribution of G0 in permafrost regions because of the limitation of in situ observation and complication of ground surface conditions. This study aims at developing an improved G0 parameterization scheme applicable to permafrost regions of the Qinghai-Tibet Plateau under clear-sky conditions. We validated several existing remote sensing-based models to estimate G0 by analyzing in situ measurement data. Based on the validation of previous models on G0; we added the solar time angle to the G0 parameterization scheme; which considered the phase difference problem. The maximum values of RMSE and MAE between “measured G0” and simulated G0 using the improved parameterization scheme and in situ data were calculated to be 6.102 W/m2 and 5.382 W/m2; respectively. When the error of the remotely sensed land surface temperature is less than 1 K and the surface albedo measured is less than 0.02; the accuracy of estimates based on remote sensing data for G0 will be less than 5%. MODIS data (surface reflectance; land surface temperature; and emissivity) were used to calculate G0 in a 10 x 10 km region around Tanggula site; which is located in the continuous permafrost region with long-term records of meteorological and permafrost parameters. The results obtained by the improved scheme and MODIS data were consistent with the observation. This study enhances our understanding of the impacts of climate change on the ground thermal regime of permafrost and the land surface processes between atmosphere and ground surface in cold regions.


Introduction
The total area of the Qinghai-Tibet Plateau (QTP) is approximately 2.57 million km 2 and it has an average elevation of over 4000 m [1].This explains why it is called 'the roof of the world' or the Earth's third pole.It is estimated that 54.3% of the QTP is located in permafrost regions, which is vulnerable to thawing due to ongoing climate change [2].The permafrost changes under the scenarios of climatic warming would have significant impacts on the land surface processes, hydrological cycles, and ecosystem from regional to global scale [3].It is of great significance to improve the ground surface soil heat flux (G 0 ) parameterization scheme to understand the impacts of climate change on the ground thermal regime in cold regions.
Ground surface soil heat flux (G 0 ) is an important parameter characterizing the energy balance between the ground surface and the atmosphere.The soil's state and thermal performance have significant impact on G 0 [4].Freeze-thaw cycles are particularly significant in determining the energy exchange profile between ground surface and atmosphere in the region and through the year [5][6][7][8].Phase change in the active layer during thawing periods consumes a large amount of heat [9].Soil moisture content also increases during thawing periods, which can regulate the energy flux and temperature gradient between the permafrost soils and the atmosphere to some extent [10].However, during freezing periods of ground surface, the temperature gradient between the ground and air across the freezing zone increases relatively due to the phase change, which enhances the heat release intensity of the system to the atmosphere.G 0 manifests the energy exchange between the active layer and the atmosphere and it affects the energy budget of the active layer and permafrost, so it is important to understand the freeze-thaw process that takes place within the ground surface in the QTP region, as well as degradation or development of permafrost through time [11].This type of analysis requires the effective parameterization of ground surface soil heat flux [12].
The single-point method of obtaining a value for ground surface soil heat flux includes observation and calculation [13].The observation methods mainly involve PlateCal (Combination of Heat Flux Plate Measurements and Calorimetry) [13][14][15][16][17] and TCAV (Averaging Soil Thermocouple Probe) [14,18].Studies by Yu et al. [19] and Xu et al. [14], however, demonstrated that it is still challenging to accurately measure the surface soil heat flux using the heat plates.Surface soil heat flux can also be calculated using the soil temperature and soil moisture measurements [14].Calculation methods include the Harmonic analysis Method (HM) [20], the Lambda Method (Incorporation of Thermal Conduction and Convection) [21], and the TDEC (the correction of thermal diffusion equation) [22].
Soil temperature and moisture data cannot be obtained directly through optical remote sensing imaging [23][24][25][26], but the measured net radiation flux (R n ) can be divided into three parts: latent heat flux, sensible heat flux, and surface soil heat flux (G 0 ), so if the ratio of G 0 to R n is known, then it can be used to estimate G 0 from remote sensing data.It is important to use an accurate ratio to accurately estimate G 0 .Previous studies suggested that the G 0 /R n ratio at the surface of bare soil is approximately 0.3 [27][28][29], and may reach to 0.5 for dry soil surfaces [29].For vegetated surfaces, studies have shown that G 0 /R n is related to the local vegetation indices, which can be obtained from remote sensing data products.For vegetation surfaces under full coverage, G 0 /R n ranges between 0.05 and 0.1 [30].Reginato et al. calculated G 0 /R n according to crop height using an empirical formula [31].Choudhury et al. obtained G 0 /R n as an exponential function of the leaf area index (LAI), yielding a correlation coefficient of 0.90 [18].Jaskson et al. calculated G 0 /R n using an exponential function based on the normalized difference vegetation index (NDVI) [32].Su proposed that G 0 /R n can be calculated for different ground surfaces, including for regions with full vegetation coverage, for completely bare land, and for those regions with varied vegetation coverage [33].The NDVI is sensitive to solar position and to the variability in soil reflectance factors [34][35][36], so there are some uncertainties in the estimates of NDVI itself.Menenti et al. [37] and Bastiaanssen [23] established an improved scheme for obtaining the G 0 /R n ratio, related to parameters such as land surface temperature and surface albedo.Ma et al. [38,39] proposed the use of MSAVI (Modified Soil Adjusted Vegetation) in arid and semi-arid and high altitudes regions, instead of NDVI.
There are also case studies that focus on the issue of phase differences between R n and G 0 [40,41].The phase difference is related to the soil type, soil moisture content, and vegetation condition [40].Due to the existence of permafrost under the active layer (i.e., the effect of a water-proof layer of permafrost), the soil moisture content in the active layer in summer is higher than that of shallow soil on seasonally frozen ground.As a result, vegetation coverage can be higher over permafrost than over seasonally frozen ground.The phase differences over permafrost are therefore different to those over Remote Sens. 2019, 11, 416 3 of 26 seasonally frozen ground.This difference can be significant for the permafrost regions of the Tibetan Plateau and is, therefore, explored in this study.When combined with the widely used MODIS data, accurate estimates of the ratios to be used can provide estimates of G 0 close to "true value" in areas with high spatial resolution (1 km).
The performance of several remote sensing models used to estimate G 0 were compared in this study, and the model with the least error was selected as the G 0 scheme for the QTP.Furthermore, the solar time angle was added into the G 0 parameterization scheme to solve the phase difference problem.The improved G 0 parameterization scheme was then validated by testing it against data directly measured across a 10 × 10 km area of permafrost at the Tanggula (TGL) station.Since the MODIS data integrity is susceptible to the cloud, and considering the conditional limitation of the solar time angle (R n > 0), the improved G 0 parameterization scheme is only applicable to the daytime of the clear-sky days.

Study Region and Data Collection
Observation data were collected from seven automatic weather stations: (Xidatan (XDT), Tanggula (TGL), Ayake (AYK), Zuonaihu (ZNH), and MS3478 (NPAM) in the permafrost region and Amdo, MS3637 (SPAM) in the seasonally frozen ground regions) on the Qinghai-Tibet Plateau (QTP) (Table 1).A 10 × 10 km permafrost region (33 • 0.8 N~33 • 6.6 N, 91 • 51 E~91 • 56.7 E) at Tanggula (TGL) was selected for the validation case study (Figure 1).The background maps on the left and right of Figure 1 are the spatial distribution map of frozen ground on Qinghai-Xizang Plateau [42] and ASTER Global Digital Elevation Model respectively.Tibetan Plateau and is, therefore, explored in this study.When combined with the widely used MODIS data, accurate estimates of the ratios to be used can provide estimates of G0 close to "true value" in areas with high spatial resolution (1 km).The performance of several remote sensing models used to estimate G0 were compared in this study, and the model with the least error was selected as the G0 scheme for the QTP.Furthermore, the solar time angle was added into the G0 parameterization scheme to solve the phase difference problem.The improved G0 parameterization scheme was then validated by testing it against data directly measured across a 10 × 10 km area of permafrost at the Tanggula (TGL) station.Since the MODIS data integrity is susceptible to the cloud, and considering the conditional limitation of the solar time angle (Rn > 0), the improved G0 parameterization scheme is only applicable to the daytime of the clear-sky days.

Study Region and Data Collection
Observation data were collected from seven automatic weather stations: (Xidatan (XDT), Tanggula (TGL), Ayake (AYK), Zuonaihu (ZNH), and MS3478 (NPAM) in the permafrost region and Amdo, MS3637 (SPAM) in the seasonally frozen ground regions) on the Qinghai-Tibet Plateau (QTP) (Table 1).A 10 × 10 km permafrost region (33° 0.8′ N~ 33° 6.6′ N, 91° 51′ E~ 91° 56.7′ E) at Tanggula (TGL) was selected for the validation case study (Figure 1).The background maps on the left and right of Figure . 1 are the spatial distribution map of frozen ground on Qinghai -Xizang Plateau [42] and ASTER Global Digital Elevation Model respectively.AYK is located at the boundary region between the permafrost and the seasonally frozen ground on the QTP, which is covered by sparse desert grasslands.ZNH is located in the continuous permafrost region and is covered by desert grasslands and sparse vegetation.TGL is located in the continuous permafrost region and is covered by alpine steppe.XDT is located in the sporadic permafrost region near the northern lower limit of the permafrost on the QTP.NPAM is located in the permafrost region [42], where the surface consists of sandy clay containing fine stones.This area is covered by alpine meadows about 15 cm high.Both Amdo and SPAM are located in the seasonally frozen ground regions, and both sites are covered with alpine meadows [43].AYK is located at the boundary region between the permafrost and the seasonally frozen ground on the QTP, which is covered by sparse desert grasslands.ZNH is located in the continuous permafrost region and is covered by desert grasslands and sparse vegetation.TGL is located in the continuous permafrost region and is covered by alpine steppe.XDT is located in the sporadic permafrost region near the northern lower limit of the permafrost on the QTP.NPAM is located in the permafrost region [42], where the surface consists of sandy clay containing fine stones.This area is covered by alpine meadows about 15 cm high.Both Amdo and SPAM are located in the seasonally frozen ground regions, and both sites are covered with alpine meadows [43].The in situ measurement data and remotely sensed data were compared in this study.The in situ observation data included surface radiation values, land surface temperatures, soil heat flux measured at the AWSs, soil temperature profiles, and soil moisture profile (Table 1).The GAME-Tibet (Asian Monsoon Experiment on the Tibetan Plateau) was carried out from May to September in 1998, so measured data for the three sites (Amdo, NPAM and SPAM) were obtained from the GAME-Tibet project, collected during this period.The temporal resolution of the data collected from the SPAM site is every 1 h, while the temporal resolution of the data obtained from other six stations is 30 min.
Remote sensing data included the MODIS satellite data products and NOAA-14/AVHRR satellite data.The attributes of the remotely sensed datasets are shown in Table 2. NOAA-14/AVHRR data provides spectral reflectance at red and short-wave infrared bands.These are used to estimate the MSAVI (Modified Soil-Adjusted Vegetation Index).MOD/MYD09 data products are used to estimate the vegetation indices (NDVI, MSAVI and vegetation coverage f c ) and the albedo.MOD/MYD11_L2 data products are used to provide the land surface temperature and estimates of the land surface emissivity.The MOD/MYD03 geolocation data product provides the latitude and longitude information allowing geocalibration of the MOD/MYD11_L2 and MOD/MYD09 data products.Algorithms for the retrieval of land surface parameters from the remote sensing datasets are described in detail in Sections 3.2 and 5.2.

Data Processing Methods
G 0 cannot be directly observed at a site and Section 3.1 describes the acquisition of G 0 at sites in detail.The main G 0 parameterization schemes were tested by comparing the synthetic results with the automatic weather station (AWS) data measured in situ.By comparing the AWS data with the remote sensing products produced vegetation indexes (NDVI, MSAVI and f c ), the most effective G 0 parameterization scheme was identified (Section 3.2).The solar hour angle was also added in to the selected G 0 parameterization scheme to improve the performance of the model (Section 3.3).Figure 2 demonstrates all the processes.

Calculation of G0 at AWS Sites
Two methods were used to calculate the surface soil heat flux at the sites (Figure 2).The first method was to calculate G0 from the soil heat flux measured at a depth of 5 cm using a heat flux plate.The soil temperature at the 5-cm depth was used to carry out the calculation as follows [15]: where G5 (W•m -2 ) is the measured soil heat flux at a depth of 5 cm; Cs is the volumetric soil heat capacity (a value of 1.18 × 10 6 J•m -3 •K -1 was used in this study [44]); T is the soil temperature at a depth of 5 cm; and t stands for the time of each measurement.We can call this the Plates method.The other method used the Thermal Diffusion Equation and Correction (TDEC) algorithm, as demonstrated by Yang and Wang [22].

Calculation of G 0 at AWS Sites
Two methods were used to calculate the surface soil heat flux at the sites (Figure 2).The first method was to calculate G 0 from the soil heat flux measured at a depth of 5 cm using a heat flux plate.The soil temperature at the 5-cm depth was used to carry out the calculation as follows [15]: where G 5 (W•m −2 ) is the measured soil heat flux at a depth of 5 cm; C s is the volumetric soil heat capacity (a value of 1.18 × 10 6 J•m −3 •K −1 was used in this study [44]); T is the soil temperature at a depth of 5 cm; and t stands for the time of each measurement.We can call this the Plates method.The other method used the Thermal Diffusion Equation and Correction (TDEC) algorithm, as demonstrated by Yang and Wang [22].The G 0 calculated by the Plates method somewhat underestimates the actual value (Figure 3a).Plates (Equation ( 1)) assumes that the volumetric soil heat capacity c s is the same between the surface and underground depths of 5 cm (1.18 × 10 6 J•m −3 •K −1 ) and that the soil temperatures at different depths have the same rates of change at any time.The rate of change of soil temperature with time over a set period, however, increases with the decrease of soil depth (Figure 3b).The TDEC method has been shown to be insensitive to the soil thermal conductivity coefficient, and its reliability has been verified by the twin experiment [22].Some studies on the plateau have shown that the TDEC method has better accuracy [12,13,45].In this study, the values of G 0 at the Amdo, NPAM, SPAM, and TGL sites were acquired using the TDEC method.G 0 at the other sites was calculated by the Plates method because of the limited availability of measured data (Table 1).The TGL site data could be used for both the Plates and the TDEC.Data from the four sites (XDT, AYK, ZNH, and TGL in 2014), which could be used with the Plates method, was used to calibrate the four existing models, to evaluate their applicability on the plateau, and to select the best applicable scheme (Sections 3.2 and 4.1).Data from four sites (Amdo, NPAM, SPAM and TGL) satisfying TDEC were then used to obtain an improved scheme based on the applicable scheme (Sections 3.3 and 4.2) and to verify the reliability of the improved scheme at those sites (Section 5.1).
where G5 (W•m -2 ) is the measured soil heat flux at a depth of 5 cm; Cs is the volumetric soil heat capacity (a value of 1.18 × 10 6 J•m -3 •K -1 was used in this study [44]); T is the soil temperature at a depth of 5 cm; and t stands for the time of each measurement.We can call this the Plates method.The other method used the Thermal Diffusion Equation and Correction (TDEC) algorithm, as demonstrated by Yang and Wang [22].

Model Name
Mathematical Form * T s is the land surface temperature, α and α are the instantaneous value and daily mean value of the land surface albedo, respectively, NDVI is the normalized difference vegetation index, MSAVI is the modified soil-adjusted vegetation index, f c is the vegetation coverage.Γ s and Γ c are the ratio of soil heat flux to net radiation for bare soil and a fully covered vegetation canopy, respectively.
The stations (TGL, XDT, ZNH, AYK, NPAM and Amdo) do not provide measured land surface temperature data (T s ), but T s can be obtained using long wave radiation data as follows [47]: where T s is the land surface temperature ( • C), R u and R d are the incoming and outgoing long wave radiation powers, respectively (W/m 2 ).σ is the Stefan-Boltzmann constant (5.67 ε is the land surface emissivity, which is approximately 0.95 on the plateau [44].
The net radiation R n can be obtained from short-wave and long-wave radiation and the land surface albedo α is the ratio of downward to upward shortwave radiation: Both the NDVI [48] and the MSAVI [49] use combined values of red (RED) and near red (NIR) bands of remote sensing data: RED and NIR are land surface reflectance values of band 1 and band 2 of MYD09 or AVHRR-1 and AVHRR-2 of NOAA-14/AVHRR, respectively.Vegetation coverage f c can be obtained via the NDVI [50]: where NDVI max and NDVI min are the NDVI values for bare soil and full vegetation coverage, respectively.The coefficients in these original models should be validated with in situ measurement data for the different cases, so that there will be eight models (four prototypes and four calibrated).The accuracies of the values produced by the models were evaluated using the mean absolute error (MAE) and the root-mean-square error (RMSE).
In the eight specific forms of the G 0 parameterization scheme, the best one will be selected by MAE and RMSE as the basic scheme to be improved, and it can be expressed as: where, Γ best (G 0 /R n ) is the best of the eight models.

Improvement of the G 0 Parameterization Scheme
Observations (Figure 4) showed that both R n and G 0 undergo obvious daily fluctuations, and the change of surface soil heat flux G 0 lags behind the change of net radiation R n [51], which means that there is a phase difference between R n and G 0 .Previous studies used the solar hour angle to solve the issue of phase difference between R n and G 0 during the daytime [40]: where R n is the net radiation, "R n > 0" represents daytime; a is the amplitude parameter, which must be determined based on measured data; b is the period of 1 day, and b = 86,400 s (i.e., 24 h); c is the phase shift.Observed values for R n and G 0 must be used to determine the value of c. Colaizzi et al. used 10,800 s as the value of c in their paper [41]; t is the solar time angle(s) [52].We then used the solar time angle to solve the problem of the phase shift between R n and G 0 on the QTP under clear-sky conditions, and can then derive the preliminary form of the improved scheme based on Equations ( 8) and ( 9): There are only two unknowns; the phase shift c and the amplitude parameter a in Equation (10).Both c and a can be determined using the data for the four TDEC sites.The coefficients of determination (R 2 ) of the linear relationship between f (Γ best , R n , c, t) and G 0 in the case of different phase shifts are obtained from the in situ measured data.The overlapping value of the phase shifts that their corresponding R 2 are greater than R 2 without considering the phase shift at each site is taken as the optimal phase offset c between G 0 and f (Γ best , R n , c, t) in permafrost regions of the QTP.The same amount of measured data from each site located at permafrost regions of the QTP were selected to form an observed dataset.The value of the amplitude parameter a can then be obtained by introducing the measured dataset and the known c value into Equation (10).
MAE and RMSE as the basic scheme to be improved, and it can be expressed as: where, Γbest (G0/Rn) is the best of the eight models.4) and the modified versions of the four schemes are shown in Appendix A. The prototypes of these models are shown in Table 3. From the mathematical forms of the SEBAL adj and Ma adj parameterization schemes (Table 4), it can be seen that the simulated G 0 of the two schemes will always be less than 0 in the daytime when the vegetation index (NDVI or MSAVI) is high and land surface temperature Ts is higher than 0, which is inconsistent with the actual situation., which is inconsistent with the actual situation.In the case, it is unfortunately not possible to use the two adjusted schemes.

Model Name
Mathematical Form According to the observation data, the corresponding "measured" G 0 and T s have negative values when net radiation R n is greater than 0 (daytime).However, from the mathematical forms of these schemes (Tables 3 and 4), we can know that the schemes Moran, Moran adj , SEBS and SEBS adj cannot simulate the negative G 0 and the schemes SEBAL and Ma can simulate the negative G 0 (Whether the simulated G 0 by the schemes SEBAL and Ma is negative or not depends entirely on the corresponding "measured" T s ).It can be seen that both schemes SEBAL and Ma performed better than the schemes Moran and SEBS.In addition, some studies proposed using MSAVI (Modified Soil Adjusted Vegetation) instead of NDVI in arid and semi-arid and high altitudes regions [38,39].Therefore, The Ma scheme is suitable for the QTP.
Both the MAE and RMSE (Table 5) obtained using the validation data in Table A1 of Appendix B clearly showed that the modified Moran adj and SEBS adj parameterization schemes performed better than the original Moran and SEBS schemes; however, the Moran and SEBS schemes were significantly worse than the Ma and SEBAL schemes for producing MAE and RMSE values, and the estimated G 0 values using the Ma scheme were slightly better than those using SEBAL, overall.In addition, Ma et al. [38] determined the original form of the Ma scheme by using data from 6 stations on the northern QTP including Amdo, NPAM and SPAM.The Ma scheme is the most suitable for the QTP and was verified by 7 sites (TGL, XDT, ZNH, AYK, Amdo, NPAM, and SPAM).
According to the mathematical form and RMSE/MAE of each scheme, we choose the scheme Ma as the best form Γ best of Γ(G 0 /R n )s on the QTP.So the improved scheme can be called the Ma impr scheme and can be expressed as: The specific values of the phase shifts c and the coefficient a in the Ma impr scheme (Equation ( 11)) can be obtained from observation data from each site.

Determination of the Optimal Phase Shifts c on the QTP
The coefficients of determination (R 2 ) of the linear relationship between f (Γ best , R n , c, t) and G 0 in the case of different phase shifts are obtained from the in situ measured data.The phase shift corresponding to maximum R 2 is taken as the optimal phase offset c between G 0 and f (Γ best , R n , c, t) at a site.Taking TGL as an example, the relationship between f (Γ Ma , R n , c, t) and G 0 in the case of different phase shifts c (Figure 5) is obtained from the observation data from TGL (Table A2 in Appendix B), and the best phase shift c at TGL can be identified.In the same way, the optimal phase shifts (Table 6) are obtained for the three sites of NPAM, Amdo, and SPAM.The best phase shifts for TGL, NPAM, Amdo, and SPAM are -3h, -5h, and no cos ("no cos" means that the scheme does not include a solar time angle that is to say that it is by default the Ma scheme).The best phase shift is in fact "no cos", which does not mean that there is no phase shift between G 0 and R n , but only that the phase shift between G 0 and R n in the seasonally frozen soil areas is small (Figure 6).Considering the surface temperature, surface albedo, and MSAVI, the phase shift effect is removed.There is a certain difference in the optimal phase shift c at each site, which may be related to the type of frozen ground at each location.Different types of frozen ground may cause different time differences in surface energy conversion.After considering the solar time angle, the coefficient of determination R 2 at TGL and NPAM in the permafrost region does increase, while R 2 reduces at SPAM and Amdo, which are located in the seasonally frozen ground region.The Ma scheme still produces the best estimates for the seasonally frozen ground region.The Ma scheme is, therefore, used in the seasonally frozen region, while the model needs further improvement in the permafrost region.
between G0 and Rn in the seasonally frozen soil areas is small (Figure 6).Considering the surface temperature, surface albedo, and MSAVI, the phase shift effect is removed.There is a certain difference in the optimal phase shift c at each site, which may be related to the type of frozen ground at each location.Different types of frozen ground may cause different time differences in surface energy conversion.After considering the solar time angle, the coefficient of determination R 2 at TGL and NPAM in the permafrost region does increase, while R 2 reduces at SPAM and Amdo, which are located in the seasonally frozen ground region.The Ma scheme still produces the best estimates for the seasonally frozen ground region.The Ma scheme is, therefore, used in the seasonally frozen region, while the model needs further improvement in the permafrost region.The sign of the phase shift c represents the chronological direction of translation of the net radiation Rn.Since the change of surface soil heat flux G0 lags behind the change of net radiation Rn, Rn needs to shift to the right, and the cosine of the rightward translation should be the negative value of the amount of translation, so the sign of the phase shift c should be negative.The sign of the phase shift c represents the chronological direction of translation of the net radiation R n .Since the change of surface soil heat flux G 0 lags behind the change of net radiation R n , R n needs to shift to the right, and the cosine of the rightward translation should be the negative value of the amount of translation, so the sign of the phase shift c should be negative.

AWSs f(no cos) f(0h) f(-1h) f(-2h) f(-3h) f(-4h) f(-5h) f(-6h)
In Table 6, the best phase shift at TGL is also −3 h, while the best phase shift c at NPAM is −5 h.When the phase shift c of TGL is greater than −3 h, the R 2 will be less than 0.587, which is the coefficient of determination without cosine function.When the phase shifts of TGL and NPAM are −3 h, their coefficients of determination are higher than the coefficients of determination without cosine function.In this case, it was decided to use −3 h as the best phase shift for the permafrost, i.e., the permafrost phase shift c is equal to −10,800 s (−3 h).
Regardless of whether or not phase offset or how much phase offset is considered, the simulation results are not good in the negative part of G 0 (Figure 5).It should be related to the actual relationship between the "measured" G 0 and measured R n at sites.Most remote sensing-based models to estimate G 0 are established on the basis of the relationship between G 0 and net radiation R n .These models consider other land surface parameters, such as land surface temperature, albedo and vegetation index, in order to better simulate G 0 /R n values.The purpose of this manuscript is to establish a G 0 parameterization scheme in sunny days on the QTP.The standard for daytime selection is that net radiation R n is greater than 0. It can be seen that the relationship between the "measured" negative G 0 and the simulated results f (Γ Ma , R n , c, t) (It can be approximated as G 0 analog value, Figure 5) is very similar to that between the "measured" negative G 0 and observed R n (the red box area in Figure 6).Therefore, it is caused by the relationship between the "measured" negative G 0 and measured R n that the simulation results of remote sensing-based models in the negative part of G 0 are not satisfactory during daytime of the clear-sky days.According to in situ measurement data of sites, the majority of the "observed" G 0 are larger than 0. So the poor effect of negative simulation will affect the application of these models in daytime.
phase shift c is equal to -10,800 s (-3 h).
Regardless of whether or not phase offset or how much phase offset is considered, the simulation results are not good in the negative part of G0 (Figure .5).It should be related to the actual relationship between the "measured" G0 and measured Rn at sites.Most remote sensing-based models to estimate G0 are established on the basis of the relationship between G0 and net radiation Rn.These models consider other land surface parameters, such as land surface temperature, albedo and vegetation index, in order to better simulate G0/Rn values.The purpose of this manuscript is to establish a G0 parameterization scheme in sunny days on the QTP.The standard for daytime selection is that net radiation Rn is greater than 0. It can be seen that the relationship between the "measured" negative G0 and the simulated results f (ΓMa, Rn, c, t) (It can be approximated as G0 analog value, Figure .5) is very similar to that between the "measured" negative G0 and observed Rn (the red box area in Figure .6).Therefore, it is caused by the relationship between the "measured" negative G0 and measured Rn that the simulation results of remote sensing-based models in the negative part of G0 are not satisfactory during daytime of the clear-sky days.According to in situ measurement data of sites, the majority of the "observed" G0 are larger than 0. So the poor effect of negative simulation will affect the application of these models in daytime.

Lookup for the Coefficient a
When the phase shift was equal to -10,800s in Equation (11), only the coefficient a is unknown and its value can be calculated based on the fitted data for clear-sky conditions (Table B3 in Appendix B).The coefficient a at TGL (a=1.3232) was close to the NPAM coefficient a (a=1.2173).It is reasonable to apply the TGL and NPAM values across all of the sites, and 1.2686 can be used as the average permafrost coefficient a.
From Figure 7, the instantaneous value of the surface albedo α in the time range from 10:00 to 16:00 is approximately equal to the daily average value of the surface albedo α ̅ on the same day.The Maimpr scheme can then be simplified as Equation ( 12):

Lookup for the Coefficient a
When the phase shift was equal to -10,800s in Equation (11), only the coefficient a is unknown and its value can be calculated based on the fitted data for clear-sky conditions (Table A3 in Appendix B).The coefficient a at TGL (a = 1.3232) was close to the NPAM coefficient a (a = 1.2173).It is reasonable to apply the TGL and NPAM values across all of the sites, and 1.2686 can be used as the average permafrost coefficient a.
From Figure 7, the instantaneous value of the surface albedo α in the time range from 10:00 to 16:00 is approximately equal to the daily average value of the surface albedo α on the same day.The Ma impr scheme can then be simplified as Equation ( 12): 16:00 is approximately equal to the daily average value of the surface albedo α ̅ on the same day.The Maimpr scheme can then be simplified as Equation ( 12):

The Uncertainty Analysis of the Improved G0 Parameterization Scheme Maimpr
The uncertainty of the scheme was analyzed using data from the TGL and NPAM sites and using the R software (R x64 3.5.0,https://www.r-project.org/), as shown in Figure 8 and Table 8.

The Uncertainty Analysis of the Improved G 0 Parameterization Scheme Ma impr
The uncertainty of the scheme was analyzed using data from the TGL and NPAM sites and using the R software (R x64 3.5.0,https://www.r-project.org/), as shown in Figure 8 and Table 8.

The Uncertainty Analysis of the Improved G0 Parameterization Scheme Maimpr
The uncertainty of the scheme was analyzed using data from the TGL and NPAM sites and using the R software (R x64 3.5.0,https://www.r-project.org/), as shown in Figure 8 and Table 8.At the TGL and NPAM sites (Figure 8 and Table 7), the slope and worst intercept absolute value from the linear relationship between G 0 (TDEC) and the simulated values of G 0 using the Ma impr scheme are around 1 and around 5 respectively, within the 95% confidence interval.The p values are Remote Sens. 2019, 11, 416 14 of 26 all less than 0.01.The RMSE values of G 0 estimated by the Ma impr scheme at the TGL and NPAM sites are 6.102 W/m 2 and 5.756 W/m 2 , respectively, the MAE values are 5.382 W/m 2 and 4.878 W/m 2 , respectively, and the values of RMSE and MAE are less than 7 W/m 2 at both sites.In summary, the simulated results using the G 0 parameterization Ma impr scheme are reliable when the input parameters are measured data.As some parameters for the Ma impr scheme are derived from remote sensing data, the uncertainty of the overall Ma impr scheme depends on its sensitivity to uncertainty in those parameters, in addition to the uncertainty of the scheme itself.For example, the vegetation index MSAVI is obtained through remote sensing data.G 0 estimates made using the Ma impr scheme in permafrost regions are affected by land surface temperature T s , land surface albedo α, and net radiation flux R n (Equation ( 12)), while R n is affected by land surface temperature T s and the land surface albedo α (Equation ( 15)).The errors in MODIS land surface temperature measurements and land surface albedo estimated using Equation ( 14) and MOD/MYD09, are 1K [53] and 0.02 [54], respectively.The sensitivity of the Ma impr scheme to land surface temperature T s and land surface albedo α, based on remotely observed data, can be derived from the percentage change of simulated G 0 values when increasing or reducing the "observed" T s by 1K or the measured α by 0.02 [55]: where, VR +accuracy and VR −accuracy indicate the percentage change of simulated G 0 values by increasing and reducing the "observed" T s by 1K or the measured α by 0.02, respectively; VRaccuracy is the ultimate sensitivity coefficient; G base is the simulated G 0 value from observed data corresponding to the unchanged land surface temperature Ts and unchanged land surface albedo α; G +accuracy and G -accuracy represent the simulated G 0 value when the "observed" T s is increased or decreased by 1K or the measured α is increased or decreased by 0.02, respectively.It is evident from Table 8 that there are differences in the sensitivity of the G 0 parameterization scheme Ma impr to land surface temperature T s or land surface albedo α at various stations, but the percentage change of the simulated G 0 values caused by the error of satellite-derived T s or α will not exceed 5%.

The Performance of The Improved Parameterization Scheme Ma impr in the Study Area
The purpose of the G 0 parameterization scheme is to estimate the G 0 values across the study region using satellite-derived land surface parameters.Owing to the impact of cloud coverage, which resulted in an incomplete MODIS dataset, the study area was reduced to 100 km 2 at TGL (Figure 1).The detailed process is shown in Figure 9.
measured α is increased or decreased by 0.02, respectively.8 that there are differences in the sensitivity of the G0 parameterization scheme Maimpr to land surface temperature Ts or land surface albedo α at various stations, but the percentage change of the simulated G0 values caused by the error of satellite-derived Ts or α will not exceed 5%.

The Performance of The Improved Parameterization Scheme Maimpr in the Study Area
The purpose of the G0 parameterization scheme is to estimate the G0 values across the study region using satellite-derived land surface parameters.Owing to the impact of cloud coverage, which resulted in an incomplete MODIS dataset, the study area was reduced to 100 km 2 at TGL (Figure . 1).The detailed process is shown in Figure 9. G 0 in the study area was estimated by Equation (12).The solar time angle t was obtained from the satellite transmitting time and the longitude of the pixel in question [52].The satellite-derived vegetation indices are obtained as described in Section 3.2.The land surface albedo α is obtained using an empirical formula and band 1~5 and band 7 data in the MODIS product MYD09 [54][55][56]: (14) where R 1 , R 2 , R 3 , R 4 , R 5 , and R 7 are the reflectances of bands 1, 2, 3, 4, 5, and 7 of the MODIS data product MYD09.The residual standard error (RSE) of the empirical formula is not more than 0.02, as validated by Shunlin Liang et al. [54].
The net radiation R n can be obtained through the surface radiation balance equation [57]: where α is the land surface albedo (Equation ( 14)), σ is Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 •K −4 ).DSR and DLR are the incoming short and long wave radiation fluxes, respectively (W/m 2 ).DSR and DLR in the study area were only a little different, and they could be replaced by actual values of DSR and DLR at the TGL station measured when the satellite transited [58].The land surface temperature T s was directly obtained from the MODIS land surface temperature data product MOD11_L2.The land surface emissivity ε s was estimated using the MYD11_L2 data product and a nonlinear formula [59]: where ε 31 and ε 32 are the emissivity of band 31 and band 32 in MYD11_L2, respectively.The residual standard error (RSE) of the empirical formula was approximately 0.027, as validated by Shunlin Liang [59].
The following image (taken at Beijing time: 2014.07.24 14:40, Figure 10) was used as a case to discuss the performance of the improved scheme Ma impr in the study area.
the satellite transmitting time and the longitude of the pixel in question [52].The satellite-derived vegetation indices are obtained as described in Section 3.2.The land surface albedo α is obtained using an empirical formula and band 1~5 and band 7 data in the MODIS product MYD09 [54][55][56] where R1, R2, R3, R4, R5, and R7 are the reflectances of bands 1, 2, 3, 4, 5, and 7 of the MODIS data product MYD09.The residual standard error (RSE) of the empirical formula is not more than 0.02, as validated by Shunlin Liang et al. [54].
The net radiation Rn can be obtained through the surface radiation balance equation [57]: where α is the land surface albedo (Equation ( 14)), σ is Stefan-Boltzmann constant (5.67 × 10 - 8 W• m -2 •K -4 ).DSR and DLR are the incoming short and long wave radiation fluxes, respectively (W/m 2 ).DSR and DLR in the study area were only a little different, and they could be replaced by actual values of DSR and DLR at the TGL station measured when the satellite transited [58].The land surface temperature Ts was directly obtained from the MODIS land surface temperature data product MOD11_L2.The land surface emissivity εs was estimated using the MYD11_L2 data product and a nonlinear formula [59]: where ε31 and ε32 are the emissivity of band 31 and band 32 in MYD11_L2, respectively.The residual standard error (RSE) of the empirical formula was approximately 0.027, as validated by Shunlin Liang [59].
The following image (taken at Beijing time: 2014.07.24 14:40, Figure 10) was used as a case to discuss the performance of the improved scheme Maimpr in the study area.In the study region depicted, G0 is low in the southeast and high in the northwest (Figure 10a), and was ranged between 140-180 W/m 2 , which accounted for 82% of the total data (Figure 10b).In the study region depicted, G 0 is low in the southeast and high in the northwest (Figure 10a), and was ranged between 140-180 W/m 2 , which accounted for 82% of the total data (Figure 10b).
The land surface emissivity LSE gap was small (Figure 11d) in the study area, and the net radiation R n was mainly affected by the land surface temperature T s and the land surface albedo α (Equation ( 15)).In addition, the solar time angle was not significantly different, therefore, the spatial distribution of the ground surface soil heat flux G 0 was mainly determined by T s , α, and MSAVI (Equation ( 12)).
There was a significant positive correlation between the "observed" G 0 and the measured T s in the data measured at the TGL station (Figure 12a), while the "observed" G 0 decreased with the increasing "measured" land surface albedo α (Figure 12b).This is consistent with previous findings [60,61].The relationship between the spatial distribution patterns for T s (Figure 11a) and G 0 (Figure 10a) in the study area was also consistent with the above analysis.The values in the spatial distribution pattern were low in the southeast and high in the northwest.However, the spatial distributions of the surface albedo α (Figure 11b) and the surface soil heat flux G 0 (Figure 10a) were negatively correlated.
From south to north, α gradually decreases, while G 0 gradually increases.These results demonstrate the reliability of the improved G 0 Ma impr parameterization scheme based on remotely sensed data since they are in line with the actual situation in the study region.The land surface emissivity LSE gap was small (Figure 11d) in the study area, and the net radiation Rn was mainly affected by the land surface temperature Ts and the land surface albedo α (Equation ( 15)).In addition, the solar time angle was not significantly different, therefore, the spatial distribution of the ground surface soil heat flux G0 was mainly determined by Ts, α, and MSAVI (Equation ( 12)).The land surface emissivity LSE gap was small (Figure 11d) in the study area, and the net radiation Rn was mainly affected by the land surface temperature Ts and the land surface albedo α (Equation ( 15)).In addition, the solar time angle was not significantly different, therefore, the spatial distribution of the ground surface soil heat flux G0 was mainly determined by Ts, α, and MSAVI (Equation ( 12)).

Comparison of the G 0 Parameterization Schemes Ma impr and Ma
The two schemes Ma impr and Ma were tested against the validation data (Table A3 in Appendix B), and the performance of the two schemes was compared.
The correlation coefficient a between "observed" G 0 (TDEC) data and G 0 values simulated by the Ma scheme for each station using remotely measured data was less than 1, while the correlation coefficient a between "observed" G 0 (TDEC) and G 0 estimated by the Ma impr scheme for each station was greater than 1, indicating that the Ma model overestimated values in permafrost regions, while Ma impr slightly underestimated values (Figure 13).The coefficient of determination R 2 for the Ma impr scheme was greater than the R 2 of the Ma scheme, and the coefficient α of the Ma impr scheme was closer to 1.This is another demonstration that the improved Ma impr scheme performed better than the Ma scheme in the permafrost region.[60,61].The relationship between the spatial distribution patterns for Ts (Figure 11a) and G0 (Figure 10a) in the study area was also consistent with the above analysis.The values in the spatial distribution pattern were low in the southeast and high in the northwest.However, the spatial distributions of the surface albedo α (Figure 11b) and the surface soil heat flux G0 (Figure 10a) were negatively correlated.From south to north, α gradually decreases, while G0 gradually increases.These results demonstrate the reliability of the improved G0 Maimpr parameterization scheme based on remotely sensed data since they are in line with the actual situation in the study region.

Comparison of the G0 Parameterization Schemes Maimpr and Ma
The two schemes Maimpr and Ma were tested against the validation data (Table B3 in Appendix B), and the performance of the two schemes was compared.
The correlation coefficient a between "observed" G0 (TDEC) data and G0 values simulated by the Ma scheme for each station using remotely measured data was less than 1, while the correlation coefficient a between "observed" G0 (TDEC) and G0 estimated by the Maimpr scheme for each station was greater than 1, indicating that the Ma model overestimated values in permafrost regions, while Maimpr slightly underestimated values (Figure 13).The coefficient of determination R 2 for the Maimpr scheme was greater than the R 2 of the Ma scheme, and the coefficient α of the Maimpr scheme was closer to 1.This is another demonstration that the improved Maimpr scheme performed better than the Ma scheme in the permafrost region.The same conclusion was reached from the remotely sensed data.Data measured at the TGL site in summer 2014 were also used because of the lack of remote sensing data in the study area, and the images produced at three moments (Beijing time: 2014.6.3015:25, 2014.7.24 14:40, and 2014.9.18 15:25) were selected for study.
From Table 9, it can be concluded that the result of the Ma impr scheme were closer to the "measured" G 0 (TDEC) data than the results of the Ma scheme, when using either directly measured data or AQUA/MODIS data.
The negative G 0 are not well simulated by the schemes Ma impr and Ma, while the performance of the scheme Ma impr is obviously better than that of the scheme Ma in the respect (Figure 14, which is a line chart converted from the data of Figure 13).Therefore, from the effect of simulating negative G 0 , the scheme Ma impr is also better than the scheme Ma in permafrost regions on the QTP.Ts is land surface temperature, α represents surface albedo, DSR and DLR are the down shortwave radiation and the long wave radiation respectively, Rn is the net radiation at the land surface, MSAVI is the modified soil Adjusted Vegetation index, and G0 is the surface soil heat flux.Obs is the measured value of the corresponding parameter, RS represents the parameter value inverted by MODIS data.G0_TDEC is the "measured" G0 by the TDEC algorithm using the observed data.G0_Maimpr and G0_Ma are the simulated G0 calculated by the Ma and Maimpr schemes using insuit data, respectively.G0 MODIS _Maimpr and G0 MODIS _Ma are the simulated G0 calculated by the Ma and Maimpr schemes using MODIS data, respectively.
From Table 9, it can be concluded that the result of the Maimpr scheme were closer to the "measured" G0 (TDEC) data than the results of the Ma scheme, when using either directly measured data or AQUA/MODIS data.

Conclusions
Four existing G 0 parameterization schemes used on the QTP were evaluated against observation data.To introduce the phase difference between G 0 and R n , an improved G 0 parameterization scheme Ma impr , which factors in the solar time angle, was established.Data from sites across a 100 km 2 permafrost region were used to verify the reliability of the Ma impr scheme.The following conclusions were drawn: The phase shifts between G 0 and R n in permafrost regions and seasonally frozen ground regions differ considerably.In permafrost regions, the shift is larger than it is in seasonally frozen ground regions.The improved G 0 Ma impr parameterization scheme was established based on the Ma scheme, by representing the characteristic impact of phase shifts in permafrost regions of the QTP.
The Ma impr scheme produced satisfactory reliability.At both sites in the study area, the maximum values of RMSE and MAE for the Ma impr scheme were 6.102 W/m 2 and 5.382 W/m 2 , respectively.The G 0 obtained using MODIS data by the Ma impr scheme, satellite-derived land surface temperature T s and the satellite-derived land surface albedo α in the study area were also consistent with the actual situation.The negative G 0 are not well simulated by the scheme Ma impr , while the scheme Ma impr is obviously better than others in the respect.
A parameterization scheme fixed to the actual situation contributes to the related research into G 0 in permafrost regions and leads to a deeper understanding these areas.The Ma impr scheme nevertheless needs further improvement because of the limited number of observation sites used in this study.A longer time series and multi-site observations would improve the understanding of energy balance processes, especially in the permafrost regions.* TGL(77) represents the 77 samples of data from the TGL site included in the overall dataset when TGL and NPAM are taken as a whole.

Figure B1
. The coefficient a of each site when the phase shift c is equal to -10800 s (in order to ensure that the result will not be biased towards a certain site), N on graph (c) represents a portion of the sample data, including the same amount of data from TGL and NPAM (Table B3), not a simple combination of all of the data from the two sites).The coefficient a of each site when the phase shift c is equal to -10800 s (in order to ensure that the result will not be biased towards a certain site), N on graph (c) represents a portion of the sample data, including the same amount of data from TGL and NPAM (Table A3), not a simple combination of all of the data from the two sites).

Figure 1 .
Figure 1.Map of the distribution of observation sites.

Figure 1 .
Figure 1.Map of the distribution of observation sites.
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 27 selected G0 parameterization scheme to improve the performance of the model (Section 3.3).Figure2demonstrates all the processes.

Figure 2 .
Figure 2. Flow chart of the of optimizing process for G0 parameterization scheme

Figure 3 .
Figure 3. Linear relationship between G0_TDEC and G0_Plates (a) and the temporal variation characteristics of soil temperature at the TGL site, at depths of 5 cm, 10 cm, and 20 cm (b)

Figure 2 .
Figure 2. Flow chart of the of optimizing process for G 0 parameterization scheme.

Figure 3 .Figure 3 .
Figure 3. Linear relationship between G0_TDEC and G0_Plates (a) and the temporal variation characteristics of soil temperature at the TGL site, at depths of 5 cm, 10 cm, and 20 cm (b) Figure 3. Linear relationship between G 0 _TDEC and G 0 _Plates (a) and the temporal variation characteristics of soil temperature at the TGL site, at depths of 5 cm, 10 cm, and 20 cm (b).

Figure 4 .
Figure 4. Daily fluctuations of net radiation Rn and ground surface soil heat flux G0 at TGL, NPAM, Amdo and SPAM sites on a sunny day

Figure 4 .
Figure 4. Daily fluctuations of net radiation R n and ground surface soil heat flux G 0 at TGL, NPAM, Amdo and SPAM sites on a sunny day.

1 .
The Best Form Γ best of Γ(G 0 /R n )s on the QTP SEBAL adj , Ma adj , Moran adj , and SEBS adj are the modified versions of the SEBAL, Ma, Moran, and SEBS models respectively (Table

Figure 5 . 6 .
Figure 5.The relationship between G0 and f (ΓMa, Rn, c, t) under the different phase shifts in TGL

Figure 5 .Table 6 .
Figure 5.The relationship between G 0 and f (Γ Ma , R n , c, t) under the different phase shifts in TGL.Table 6.The coefficient of determination R 2 of f (Γ Ma , R n , c, t) and G 0 at the four stations with different phase shifts c.

Figure 6 .
Figure 6.The relationship between G0 and Rn(t) under the different phase shifts at SPAM.

Figure 6 .
Figure 6.The relationship between G 0 and R n (t) under the different phase shifts at SPAM.

Figure 7 .
Figure 7. Relationship between the land surface albedo instantaneous value α and the daily mean value α ̅ in permafrost and seasonally frozen ground

Figure 7 .
Figure 7. Relationship between the land surface albedo instantaneous value α and the daily mean value α in permafrost and seasonally frozen ground.

Figure 7 .
Figure 7. Relationship between the land surface albedo instantaneous value α and the daily mean value α ̅ in permafrost and seasonally frozen ground

Figure 8 .
Figure 8.The linear relationship between the land surface soil heat flux G 0 calculated by the TDEC method (G 0 _TDEC) and the parameterization scheme Ma impr (G 0 _Ma impr ) using measured data from the TGL and NPAM sites (the red line represents a straight line with a slope of 1 and an intercept of 0).The blue line represents the 95% confidence interval line.Land surface temperature T s is obtained by Equation (2) using observed long wave radiation.

Figure 9 . 9 .
Figure 9.The inversion process of ground surface soil heat flux G0 in the study area Figure 9.The inversion process of ground surface soil heat flux G 0 in the study area.

Figure 10 .
Figure 10.The spatial (a) and numerical (b) distribution of ground surface soil heat flux G0 in the study area

Figure 10 .
Figure 10.The spatial (a) and numerical (b) distribution of ground surface soil heat flux G 0 in the study area.

Figure 11 .
Figure 11.The spatial distribution characteristics of the land surface temperature LST (a), land surface albedo Albedo (b), MSAVI (c), and land surface emissivity LSE (d) in the study area

Figure 11 . 27 Figure 11 .
Figure 11.The spatial distribution characteristics of the land surface temperature LST (a), land surface albedo Albedo (b), MSAVI (c), and land surface emissivity LSE (d) in the study area.

Figure 12 .
Figure 12.The relationship between the "observed" G 0 and the measured land surface temperature T s (a) and the relationship between the "observed" G 0 and the measured land surface albedo α (b) at the TGL station on sunny summer days of 2010-2012.

Figure 13 .
Figure 13.Comparison between the Maimpr and Ma schemes at each site

Figure 13 .
Figure 13.Comparison between the Ma impr and Ma schemes at each site.
T s is land surface temperature, α represents surface albedo, DSR and DLR are the down shortwave radiation and the long wave radiation respectively, R n is the net radiation at the land surface, MSAVI is the modified soil Adjusted Vegetation index, and G 0 is the surface soil heat flux.Obs is the measured value of the corresponding parameter, RS represents the parameter value inverted by MODIS data.G 0 _TDEC is the "measured" G 0 by the TDEC algorithm using the observed data.G 0 _Ma impr and G 0 _Ma are the simulated G 0 calculated by the Ma and Ma impr schemes using in-suit data, respectively.G 0 MODIS _Ma impr and G 0 MODIS _Ma are the simulated G 0 calculated by the Ma and Ma impr schemes using MODIS data, respectively.

Figure 14 . 14 .
Figure 14.Comparison between the Maimpr and Ma schemes at each site Figure 14.Comparison between the Ma impr and Ma schemes at each site.

Figure A4.
Figure A4.The coefficient a of each site when the phase shift c is equal to -10800 s (in order to ensure that the result will not be biased towards a certain site), N on graph (c) represents a portion of the sample data, including the same amount of data from TGL and NPAM (TableA3), not a simple combination of all of the data from the two sites).

Table 1 .
Observation instruments and observation parameters.

Table 2 .
Specifications of remote sensing datasets.

Table 3 .
Common remote sensing-based models for estimating Γ

Table 4 .
Adjusted Remote Sensing Model of Estimating Γ

Table 7 .
The uncertainties in the simulated G 0 by the Ma impr scheme using in situ data from TGL and NPAM stations.

Table 8 .
Average sensitivity of the G 0 parameterization scheme Ma impr to land surface temperature T s and land surface albedo α.

Table 9 .
Land surface parameters and ground surface soil heat flux G 0 inverted in the images taken at three selected times (hh:mn dd.mm) in 2014.