Estimating High Resolution Daily Air Temperature Based on Remote Sensing Products and Climate Reanalysis Datasets over Glacierized Basins: A Case Study in the Langtang Valley, Nepal

Abstract

Near surface air temperature (Ta) is one of the key input parameters in land surface models and hydrological models as it affects most biogeophysical and biogeochemical processes of the earth surface system. For distributed hydrological modeling over glacierized basins, obtaining high resolution Ta forcing is one of the major challenges. In this study, we proposed a new high resolution daily Ta estimation scheme under both clear and cloudy sky conditions through integrating the moderate-resolution imaging spectroradiometer (MODIS) land surface temperature (LST) and China Meteorological Administration (CMA) land data assimilation system (CLDAS) reanalyzed daily Ta. Spatio-temporal continuous MODIS LST was reconstructed through the data interpolating empirical orthogonal functions (DINEOF) method. Multi-variable regression models were developed at CLDAS scale and then used to estimate Ta at MODIS scale. The new Ta estimation scheme was tested over the Langtang Valley, Nepal as a demonstrating case study. Observations from two automatic weather stations at Kyanging and Yala located in the Langtang Valley from 2012 to 2014 were used to validate the accuracy of Ta estimation. The RMSEs are 2.05, 1.88, and 3.63 K, and the biases are 0.42, −0.68 and −2.86 K for daily maximum, mean and minimum Ta, respectively, at the Kyanging station. At the Yala station, the RMSE values are 4.53, 2.68 and 2.36 K, and biases are 4.03, 1.96 and −0.35 K for the estimated daily maximum, mean and minimum Ta, respectively. Moreover, the proposed scheme can produce reasonable spatial distribution pattern of Ta at the Langtang Valley. Our results show the proposed Ta estimation scheme is promising for integration with distributed hydrological model for glacier melting simulation over glacierized basins.

1. Introduction

Near surface air temperature (Ta), usually referring to the sheltered thermometer temperature at 2 m above ground, is one of the key input parameters in land surface models and hydrological models [1,2] and has been widely used in various fields, such as glacier melting modeling, crop growth simulation, evapotranspiration estimation, drought monitoring, urban heat island and global climate change [3,4,5,6,7,8,9]. Ta information can be obtained through three approaches: the station observation, climate reanalysis datasets and remote sensing. Station observation is the most direct and accurate method and is always used as the ground truth to validate the Ta estimated from other approaches. Advantages of station observation include long time records and high observation frequencies. The longest records of surface air temperature at some station can be over one hundred years [10]. At manually operated stations, observations are recorded four times per day (00:00, 06:00, 12:00 and 18:00 in UTC time) while the observation frequency at automatic weather stations can be up to one minute. However, the density of the ground stations is still very low and the spatial distribution of current stations is quite non-uniform. The situation is even worse at glacierized basins due to their harsh natural condition. The climate reanalysis datasets are generally produced by the atmospheric data assimilation system in which the numerical model ingests all available observations from the ground stations, sounders and satellites. The climate reanalysis datasets have been widely used to feed land surface models and large scale hydrological models. Limited by the computational resources and input data sources, the spatial resolution of climate reanalysis datasets is usually relatively coarse, ranging from tens to hundreds of kilometers. For example, the ERA-Interim [11] produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) has a resolution about 80 km, and the MERRA (Modern-Era Retrospective analysis for Research and Applications) [12] from the National Aeronautics and Space Administration (NASA) has a spatial resolution about 0.5°. Therefore, dynamical or statistical downscaling is needed when these climate reanalysis datasets are applied to small basins [13] which again introduces large uncertainties. In the era of satellite, remote sensing technique provides another promising way to obtain near surface temperature information with the strengths of large spatial coverage, repeat visiting and relatively high resolution.

The remote sensing based Ta estimation methods can be classified into four categories: direct inversion, physical, semi-empirical and statistical regression methods. In the direct inversion method, atmospheric temperatures at different height levels are directly retrieved from the brightness temperature of satellite channels and Ta can then be obtained using the temperature lapse rate. For example, the atmospheric temperature and moisture profiles products (MOD07_L2/MYD07_L2) from the moderate-resolution imaging spectroradiometer (MODIS) are produced at 20 vertical levels by using synthetic regression retrieval algorithm under clear sky, in which the regression coefficients are derived using a fast radiative transfer model with atmospheric characteristics taken from a global profile dataset [14]. Bisht et al. [15] and Zhu et al. [16] obtained the Ta by using the temperature lapse rate method from the MODIS temperature profile products. The main problem of the direct inversion method is that the spatial resolution and inversion accuracy are relatively low and the lapse rate is assumed to be vertically constant which may be inconsistent with the actual situation [17]. The physical approach is based on the surface energy balance equation in which the temperature gradient, i.e., the difference between Ta and land surface temperature (LST), is related to the surface sensible heat flux with some resistance terms [18,19]. The physical approach usually uses the remotely sensed LST and net radiation (R_n) products as input parameters and could obtain a reasonable Ta estimation result. Sun et al. [19] proposed a Ta estimation method with LST, Rn and land surface properties, such as normalized difference vegetation index (NDVI), surface resistance, etc. However, some parameters in this method (such as the ground heat flux to net radiation ratio and crop water stress index) highly depend on NDVI, which limits its application in the glacierized basins. Semi-empirical methods establish the relationship between Ta and LST using regression formula with some theoretical basis [20,21,22]. One of the widely used semi-empirical methods is the temperature–vegetation index (TVX) method, which assumes that the Ta equals to the canopy temperature when vegetation coverage reaches its maximum [22,23,24]. In vegetation covered area, the MODIS LST and NDVI is always negatively correlated due to the vegetation cooling effect. After Ta is obtained from station observation at the maximum NDVI region, the sliding window method is adopted to obtain the Ta in each pixel of the satellite images. The TVX method also needs ground observation and could not be used in ungauged basins and vegetation sparse watersheds. Statistical regression methods usually try to build the relationship between station observed Ta and remotely sensed LST, elevation, NDVI, surface albedo and other explaining variables (such as longitude/latitude) using multi-variable linear regression model or machine learning approaches [25,26,27,28,29,30,31,32,33,34,35]. Using MODIS LST products and observed daily mean temperature from 95 meteorological stations in the Tibetan Plateau, Zhang et al. [36] compared the performances of several regression methods, including MLR (multiple linear regression), PLS (partial least squares regression), BPNN (back propagation neural network), SVR (support vector regression), RF (random forests) and CR (Cubist regression), and 15 different combination schemes of the four MODIS observations per day over the Tibet Plateau. All the statistical regression methods need station observations to train the models which limits their potential applications in ungauged basins.

The Hindu-Kush-Himalayan (HKH) region is the water towers of Asia as many water resources are stored in the high altitude snow and glaciers [37,38]. The HKH region is extremely sensitive to climate change and evident impacts from global warming, including glacier retreat [39,40,41] and permafrost degradation [42,43], can be observed in this region. Distributed hydrological modeling is a valuable tool to assess the climate change impact on the glacier/snow melting and the resulting runoff. However, the application of high resolution distributed hydrological modeling is limited in the HKH region mainly due to lack of high resolution meteorological forcing. For most applications in the glacierized basins of HKH region, Ta is spatially interpolated from several, or even single, observation stations using the temperature lapse rate, which can hardly be calibrated for ungauged glacierized basins. Here we propose a new framework to estimating high-resolution Ta over a glacierized basin using remote sensing products and climate reanalysis dataset which can be easily extended to other ungauged glacierized basins of the HKH region.

Following this introduction, we present information about the study area, datasets used in this study as well as details of proposed Ta estimation framework in Section 2. Experimental results and analysis are presented in Section 3. We discuss the effects of terrain and cloud on Ta estimation as well as possible improvement in Section 4 and conclude our study in Section 5.

2. Materials and Methods

2.1. Study Area

Our main study area is the Langtang Valley located in the boundary region between Nepal and China (Figure 1), which is in the heart region of HKH. The Langtang Valley is a glacierized basin bounded by longitude from 85°35′E to 85°48′E and latitude from 28°08′N to 28°23′N with elevation ranging from 3800 to 7234 m and annual mean Ta of about 276 K. The total area of the Langtang Valley is 360 km², 166 km² of which is glacierized area. The elevation and land cover maps around the Langtang Valley are presented in Figure 1. The digital elevation model (DEM) data is from the shuttle radar topography mission (SRTM) provided by NASA. The land cover dataset is the GlobCover-2009 provided by European Space Agency (ESA). The Nepalese side of broader study area (whole area in the bounding box) is in the transition zone from plain to mountainous area with high altitude gradient and complicated topography. The southern part is low-elevated plain area, while the northern part is high-elevated mountainous area. The main land cover types in our study area include permanent snow and ice, grassland, closed needle leaved evergreen forest, bare areas and rainfed croplands, among which the permanent snow and ice accounts for 42.3% of the total area.

2.2. Datasets

2.2.1. MODIS LST Product

The MODIS with 36 spectral bands ranging from visible to thermal infrared wavelength is a key instrument aboard on the Terra and Aqua satellites launched in 1999 and 2002, respectively. For most area of the earth, the MODIS sensors could provide four times per day observations at local time around 1:30, 10:30, 13:30 and 22:30. Due to the high observation frequency and multi-level products in different temporal/spatial resolutions, the MODIS products have been widely used in various applications. We use the MOD11A1/MYD11A1 LST products provided by the Level-1 Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/), at the Goddard Space Flight Center. The current version of those LST products is version 6 with a 1 km spatial resolution. The algorithm adopted by MOD11A1/MYD11A1 LST retrieval is to use the classification-based emissivity method to estimate the emissivity in MODIS band 31 (10.78–11.28 μm) and band 32 (11.77–12.27 μm), and then use the generalized split-window LST retrieval algorithm to generate the LST products under clear sky, while under cloudy sky, the LST is not retrieved in the product [44,45]. According to the ground validation results of these LST products, the accuracy is about 1 K under clear sky, and could meet the accuracy requirement of most modeling applications on LST [46,47]. The LST products used in this study range from 1 January 2012 to 31 December 2014, and cover the Langtang basin (Figure 1). LST products from both Aqua and Terra in both daytime and nighttime are used. To facilitate product processing, all LST data is resampled to 0.01° using the nearest neighbor method with the help of the MRT (MODIS Reprojection Tool) provided by NASA.

2.2.2. The CLDAS Climate Reanalysis Dataset

The climate reanalysis dataset used in this study is from the China Meteorological Administration (CMA) land data assimilation system (CLDAS) version 2.0 provided by the China metrological data service center (http://data.cma.cn/) [48,49]. CLDAS is developed by the National Meteorological Information Center (NMIC) of CMA using data fusion and assimilation technology. The spatial coverage of CLDAS is the Eastern Asia bounded by longitude from 60°E to 160°E and latitude from 0° to 65°N with a spatial resolution of 0.0625° and CLDAS has an hourly temporal frequency. The main input datasets in CLDAS are the quality controlled ground observation data from more than thirty thousand automatic weather stations in China, the background forecasting from the Global Forecast System (GFS) of the National Centers for Environmental Prediction (NCEP) (NCEP/GFS) and the Fengyun-2 (FY-2) satellite data. Through fusing multi-source data from ground measurements, satellite observations and numerical model products, the CLDAS provides high-quality gridded hourly Ta, air pressure, humidity, wind speed, precipitation and surface shortwave radiation. The fusion of temperature, air pressure, humidity and wind speed is mainly achieved by the Space-Time Multiscale Analysis System (STMAS), which is an upgraded version of Local Analysis and Prediction System (LAPS) of the National Oceanic and Atmospheric Administration (NOAA) [50,51] using the multi-grid sequential variational method. Validation study of this product showed the bias and root mean square error (RMSE) of Ta are less than 1 K and 1.5 K, respectively, benchmarked with national automatic ground station observed air temperature in China [52].

2.2.3. Ground Observation at the Langtang Valley

Observations collected at two automatic weather stations located at the Langtang Valley (Figure 1) are used to validate the estimated Ta from remote sensing products and climate reanalysis datasets. The Yala base camp station (Yala) locates at 85.61208°E, 28.23252°N with an elevation of 5213 m and the Kyanging station locates at 85.56169°E, 28.21081°N with an elevation of 3869 m [53]. Hourly data from these two stations are provided by the International Centre for Integrated Mountain Development (ICIMOD) regional database system (http://rds.icimod.org/). Data records from the Yala station cover the time period from 22 March 2012 to 31 December 2014 while those from the Kyanging station are from 8 May 2012 to 31 December 2014. As we aim at estimating daily maximum, minimum and mean temperatures in this study, the hourly data from these two stations are averaged to obtain the ground-based daily maximum, minimum and mean temperatures which will later be used to validate the estimated Ta from remote sensing and reanalysis datasets.

2.3. Methods

The main technical framework for Ta estimation is as follows: (a) reconstructing MODIS LST products using the data interpolating empirical orthogonal functions (DINEOF) method to obtain spatio-temporal continuous LST; (b) upscaling the reconstructed LST, DEM and duration of daylight (DDL) to the scale of CLDAS climate reanalysis dataset and establishing the statistical relationship between Ta and explaining variables at the CLDAS scale; and (c) applying the resulted relationship at the original MODIS LST scale to obtain high resolution Ta estimation. More details about the technical framework and accuracy assessment method are provided in the following subsections.

2.3.1. The LST Gap-Filling Method

Due to the influence of cloud contamination and gaps in orbit converge, missing data problem in MODIS LST products is a common issue, which limits the application of LST. To overcome this problem, the DINEOF method was adopted to reconstruct LST in this study. The DINEOF method is proposed by Beckers et al. in 2003 [54], and has been successfully used in sea surface temperature reconstruction [55,56,57]. The DINEOF is developed based on the Empirical Orthogonal Functions to solve the missing data problem appearing in many geoscience fields. Compared with the traditional interpolation methods, the DINEOF needs fewer input parameters, has high computational efficiency and does not need to compute the spatial/temporal correlation. The DINEOF also has no strict assumption on the distribution of the original dataset. Therefore, it could be applied even under high missing data rate condition [58]. The basic principle of the DINEOF method for MODIS LST reconstruction is as following.

Assuming the overpassing time number of each LST product is T in our case study and each overpass has S pixels, a LST matrix with a size of T × S can be organized as following:

$$\mathit{LSTM}=\left(\begin{array}{c}\begin{array}{cc}\begin{array}{ccc}LS{T}_{1,1}& LS{T}_{1,2}& \cdots \end{array}& LS{T}_{1,T}\end{array}\\ \begin{array}{cc}\begin{array}{ccc}LS{T}_{2,1}& LS{T}_{2,2}& \cdots \end{array}& LS{T}_{2,T}\end{array}\\ ⋮\\ \begin{array}{cc}\begin{array}{ccc}LS{T}_{S,1}& LS{T}_{S,2}& \cdots \end{array}& LS{T}_{S,T}\end{array}\end{array}\right),$$

(1)

The matrix LSTM could be decomposed using the singular value decomposition (SVD) method as following:

$$\mathit{LSTM}=\mathit{U}\sum \mathit{V},$$

(2)

where ∑ contains the singular values of the matrix LSTM, U contains the EOFs which represent the spatial distribution mode of LSTM, and V contains the principal components (PCs) which represent the temporal distribution of LSTM.

The original LSTM could be accurately reconstructed by using Equation (2) with decomposed U, V and ∑ from LSTM. However, in the actual reconstruction process, to filter the outliers and fill the gaps, only the first k components of U and V are used depending on the predefined reconstruction accuracy. We have tested the suitability of the DINEOF method in LST reconstruction on the Tibetan Plateau and obtained good reconstruction accuracy [58].

2.3.2. The Air Temperature Estimation Method

Ta was estimated using the reconstructed MODIS LST products, DEM and DDL as input parameters. The DDL was used to reflect the annual variation of Ta, and DEM was used to reflect the lapse of temperature with elevation, while the LST was used to represent the spatial distribution pattern of Ta. By using the DINEOF method, we obtained spatial and temporal continuous MODIS LST products four times per day. The LST, DEM and DDL were then aggregated to the CLDAS scale. A multi-variable linear regression model was then trained for daily maximum, minimum and mean temperatures at CLDAS scale using the CLDAS Ta and up-scaled input parameters:

$$Ta={\mathrm{a}}_1*LS{T}_{MODDAY}+{\mathrm{a}}_2*LS{T}_{MODNIG}+{\mathrm{a}}_3*LS{T}_{MYDDAY}+{\mathrm{a}}_4*LS{T}_{MYDNIG}+\mathrm{b}*DEM+\mathrm{c}*DDL+\mathrm{d},$$

(3)

where $LS{T}_{MODDAY}$, $LS{T}_{MODNIG}$, $LS{T}_{MYDDAY}$, and $LS{T}_{MYDNIG}$ represent the reconstructed LST for Terra and Aqua satellites in day and nighttime, respectively. The DDL is calculated from local latitude ($\phi$) and day of year (DOY) following the method of Chen et al. [59]:

$$DDL=\frac{24}{\pi }\arccos (\tan (\frac{\phi }{180}\pi )\tan (\frac{23.45\pi }{180}\sin (\frac{2\pi (284 + DOY)}{365}))),$$

(4)

We further assumed the regression coefficients are scale-independent and directly applied the model to the original MODIS LST scale to obtain high-resolution Ta. The technical workflow for Ta estimation is summarized in Figure 2.

We used both simple multi-variable linear regression (SML) and piecewise multi-variable linear regression (PML) methods to fit the model in Equation (3). The PML method was used separately for summer (from 1 June to 30 September) and winter (from 1 October to 31 May) seasons to consider potential seasonal change in the relationship between Ta and other variables. Due to the high spatial variability of Ta in the mountainous area, a local regression method instead of a global method was adopted to build the model. To find the appropriate window size, 3 × 3, 5 × 5, and 7 × 7 window sizes have been tested using the daily maximum, minimum and mean Ta from CLDAS data, and the statistical results of the models showed that, the best window size is 5 × 5, which will be adopted in both SML and PML to obtain the high resolution Ta at the Langtang Valley in later study.

2.3.3. Statistical Metrics for Accuracy Assessment

The metrics used to evaluate the Ta estimation performance are RMSE and bias, which are defined as follows:

$$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{(T{a}_{pi}-T{a}_{oi})}^2},$$

(5)

$$bias=\frac{1}{n}\sum _{i=1}^n(T{a}_{pi}-T{a}_{oi}),$$

(6)

where $T{a}_{pi}\mathrm{and}T{a}_{oi}$ are the estimated and station observed Ta, respectively.

3. Results and Analysis

3.1. The MODIS LST Reconstrution Result

The spatial distribution patterns of annual mean LST before and after reconstruction are presented in Figure 3 and the temporal variations of mean LST across the study region before and after reconstruction are presented in Figure 4. The spatial patterns of annual mean LST before and after reconstruction are quite similar, and they are highly consistent with the spatial distribution of topography and land cover data. The seasonal variations of LST before and after reconstruction agree with each other very well, and the reconstructed data are more stable in the time domain. Larger fluctuations are observed in the time series of mean LST across the study region before reconstruction due to the impact from cloud. We excluded the cloud-contaminated area when calculating the area mean LST in the study region before reconstruction. Excluding data from these unsampled areas directly leads to the fluctuations in the area mean LST time series before reconstruction, as shown in Figure 4.

In Figure 5, the spatial distribution patterns of LST before and after reconstruction on 19 February 2012 are presented. As illustrated, the DINEOF method has the capability to recover the LST very well. In the daytime, the distribution of original LST products is fragmentized, but after the reconstruction, the spatially continuous LST products are generated. In the nighttime, the spatial distribution pattern of the original MOD11A1 LST is relatively better due to a lower missing data rate, but the western part of nighttime MYD11A1 is missing since that area is not covered by the satellite. The filled nighttime MYD11A1 LST is consistent well with the corresponding nighttime MOD11A1, which means that the DINEOF method has the ability to deal with large missing data problem such as the gaps observed in the MODIS LST products. The spatial patterns of the recovered LST in both daytime and nighttime agree well with those of the DEM and land cover map.

3.2. Spatial Distribution Pattern of Estimated Ta

Figure 6 and Figure 7 present the spatial distribution patterns of the estimated daily mean/maximum/minimum Ta from CLDAS, SML and PML on 1 March 2013 and 8 August 2013, which are in the winter and summer, respectively. From both Figure 6 and Figure 7, it could be discovered that CLDAS Ta could not represent the detailed distribution of Ta at the Langtang Valley due to its relatively coarse spatial resolution, which limits its application in hydrological modeling over this basin. Results obtained by the SML and PML models could recover more details about the spatial distribution of Ta. In the wintertime, both the SML and PML models could recover the spatial distribution of Ta very well, and the results from these two models are very similar. The spatial distribution patterns of Ta obtained by SML and PML models are generally consistent with the spatial pattern of DEM and land cover map. In the summer, the SML model estimated Ta has smaller spatial variability than that from the PML model. The Ta estimated by PML model could provide more details of spatial distribution of Ta than SML and its general pattern is more consistent with the original CLDAS Ta, which indicates the PML method has better performance than SML method and could be used to estimate high resolution Ta from MODIS LST and climate reanalysis dataset over this glacierized basin.

3.3. Accuracy Assessment of the Estimated Ta

We validated the accuracy of the estimated Ta at the Langtang Valley benchmarked with observations at two automatic weather stations, i.e., the Kyanging and Yala stations. The validation result has been summarized in

Table 1. The validation result of CLDAS, SML and PML estimated Ta benchmarked with station observation. Bolded values represent the best performances among the three Ta estimates. CLDAS represents China Meteorological Administration (CMA) land data assimilation system. SML is simple multi-variable linear regression method and PML is piecewise multi-variable linear regression method.

			Maximum		Mean		Minimum
			RMSE (K)	Bias (K)	RMSE (K)	Bias (K)	RMSE (K)	Bias (K)
Kyanging	Whole Year	CLDAS	2.07	−0.90	2.87	−2.14	5.13	−4.19
		SML	2.21	0.17	2.55	−1.22	4.77	−3.83
		PML	2.05	0.42	1.88	−0.68	3.63	−2.86
	summer	CLDAS	1.27	0.46	0.94	−0.74	2.23	−2.09
		SML	2.23	−0.75	3.20	−2.46	5.72	−5.09
		PML	2.01	1.16	1.40	−0.40	2.84	−2.39
	winter	CLDAS	2.40	−1.66	3.51	−2.92	6.18	−5.36
		SML	2.19	0.70	2.09	−0.51	4.13	−3.11
		PML	2.07	−0.01	2.11	−0.85	4.01	−3.13
Yala	Whole Year	CLDAS	7.26	7.03	5.23	4.95	3.65	2.64
		SML	4.81	4.30	2.93	2.12	2.84	−0.37
		PML	4.53	4.03	2.68	1.96	2.36	−0.35
	summer	CLDAS	7.28	7.15	5.29	5.24	3.71	3.59
		SML	5.07	4.84	2.83	2.45	1.85	0.28
		PML	4.26	4.05	2.39	2.11	1.36	0.43
	winter	CLDAS	6.09	4.67	4.82	2.60	5.00	−0.08
		SML	4.61	3.90	3.01	1.87	3.39	−0.85
		PML	4.72	4.02	2.87	1.85	2.89	−0.93

. The RMSEs of CLDAS Ta are larger than those for the SML and PML estimated Ta, only with an exception for the estimated Ta at the Kyanging station in the summer. For the daily maximum Ta, the PML estimated Ta at the Kyanging and Yala stations has the smallest RMSE and bias, except in the summer at the Kyanging station and in the winter at the Yala station. For the daily mean Ta, Ta estimated by the PML at the Kyanging and Yala stations has the smallest RMSEs and biases, except in the summer and winter at the Kyanging station. For the daily minimum Ta, Ta estimated by the PML in both summer and winter has the smallest RMSEs, except that the CLDAS has the smallest RMSE at the Kyanging station in the summer. Overall, Ta estimated by the PML has the best validation performance. The RMSE and bias are 2.05 K and 0.42 K at the Kyanging station and 4.53 K and 4.03 K at the Yala station for the daily maximum Ta, are 1.88 K and −0.68 K at the Kyanging station and 2.68 K and 1.96 K at the Yaya station for the daily mean Ta, and are 3.63 K and −2.86 K at the Kyanging station and 2.36 K and −0.35 K at the Yala station for the daily minimum Ta.

The time series of Ta obtained from CLDAS, SML and PML methods are compared with the station observation in Figure 8. The daily maximum, mean and minimum Ta from CLDAS at the Yala station are higher than the station observation, while at the Kyanging station, the daily minimum Ta from CLDAS is lower than the station observation and the maximum and mean Ta from CLDAS are more consistent with the station observation. The SML estimated Ta matches well with station observation at both Kyanging and Yala stations in winter, but in summer, the performance of SML is not as good as in winter. The PML estimated Ta at the Kyanging station is consistent well with the station observation only with a slight underestimation of daily minimum Ta in summer. The PML estimated Ta at the Yala station also matches well with the station observation except an overestimation of maximum Ta in summer. These validation results show that our proposed PML Ta estimating scheme has good performance in estimating high spatial resolution daily Ta at the Langtang Valley, except about a 5 K overestimation of maximum Ta at the Yala station in summer.

4. Discussion

4.1. Impact of Weather Condition on the Ta Estimation Accuracy

Weather condition is an important factor influencing the application of remote sensing products, especially those from visible and thermal bands. As MODIS LST is only available under clear sky in the land products (MOD11 and MYD11), many previous studies [25,27,35,36,59] estimated Ta by integrating MODIS LST and station observations under clear sky. To estimate Ta under both clear and cloudy sky conditions, Zhu et al. [16] used the clear sky LST from the MODIS cloud product (MOD06) and Ta from the MODIS atmosphere profile products (MOD07) to build the relationship between LST and Ta at 5-km scale, and then applied this relationship under cloudy sky. The LST provided by the MOD06 cloud product is a combination of data from various sources including MOD11 product, National Centers for Environmental Prediction (NCEP) gridded reanalysis and Data Assimilation Office (DAO) data [60]. The validation RMSE of estimated daytime Ta by Zhu et al. under clear sky is 2.5 K and under both clear and cloudy sky conditions is 2.85 K at eastern part of the Qaidam Basin in China. In this paper, we proposed a novel scheme to estimate high resolution daily mean, maximum and minimum Ta by integrating the reconstructed MOD11A1/MYD11A1 LST and CLDAS Ta reanalysis. We tested the scheme over the Langtang Valley, which is a glacierized basin in the HKH region. The performance of our method under clear and cloudy sky at the Kyanging and Yala stations is presented in Figure 9. At the Kanging station, the PML estimated Ta has larger RMSE (2.19 K) and bias (−1.19 K) under clear sky than under cloudy sky, but at the Yala station, the PML estimated Ta has larger RMSE (2.78 K) and bias (2.13 K) under cloudy sky than under clear sky. Overall, the estimated RMSEs and biases under clear and cloudy sky have no large differences, which means that our proposed method has the ability to obtaining high resolution Ta under both clear and cloudy sky conditions over the Langtang Valley.

4.2. Lapse Rate Derived from the Station Observed, CLDAS and Remote Sensing Based Ta

Elevation is one of the key factors that determine the spatial distribution of Ta, especially in mountainous area. Traditionally, distributed Ta over a glacierzed basin is obtained through using a constant temperature lapse rate for hydrological simulation. It has been recognized that the temperature lapse rate has significant temporal and spatial variabilities [17,61,62,63]. Figure 10 shows the temperature lapse rate derived from the Kyangjing and Yala station observed Ta, compared with those derived from the gridded CLDAS and remote sensing based Ta over the Langtang Valley. The three temperature lapse rate estimates are generally consistent with each other showing obvious seasonal variability (larger absolute lapse rates in winter and smaller ones in summer). The temperature lapse rate derived from the remote sensing based Ta ranges from −4 to −7 K/km, which is generally consistent with the reported observation values from independent station observations [61]. Due to the influence of terrain on air flow and surface radiation, the lapse rate varies with elevation and terrain aspect [64,65]. The lapse rates at different elevation bands and aspects derived from the remote sensing based daily Ta are presented in

Table 2. Temperature lapse rate at different elevation bands derived from the remote sensing based daily Ta.

Elevation (km)	Lapse Rate (K/km)
<2.5	−3.35
2.5–3.5	−4.50
3.5–4.5	−5.12
4.5–5.5	−5.25
5.5–6.5	−5.27
>6.5	−5.67

and

Table 3. Temperature lapse rate at different terrain aspects derived from the remote sensing based daily Ta.

Aspect	Lapse Rate (K/km)
North	−4.7
Northeast	−5.5
East	−5.4
Southeast	−5.3
South	−5.1
Southwest	−5.3
West	−5.3
Northwest	−5.0

. The results show that the absolute lapse rate increases with the elevation and also increases from northern aspect to northeast aspect, decreases from northeast aspect to southern aspect, increases from southern to western aspect and decreases from western to northern in a clockwise direction. The variation of lapse rate with aspect may be the result of the combined effect of the southwest monsoon and larger solar radiation intercepted by the southern slope.

4.3. Relationship between Ta and Runoff at Glacierized Basins

At glacierized basins, the runoff is quite sensitive to the variation of Ta [39]. The temperature-index model is one of the simplest and most widely used methods to simulate the ice/snow melting process [66,67,68,69,70]. According to the temperature-index model, the melting of ice/snow per day is equal to ${\mathrm{DDF}}_{\mathrm{snow}/\mathrm{ice}}*\mathrm{T}$, where DDF_ice/snow is the degree day factor for ice or snow, T is equal to 0 when the daily mean air temperature is less then 0 °C, otherwise T is equal to the daily mean air temperature. The DDF for ice ranges from 5.4 to 20 mm·d⁻¹·°C⁻¹, while the DDF for snow is generally smaller than that for ice, ranging from 2.7 to 11.6 mm·d⁻¹·°C⁻¹, according to the review by Hock [9]. Zhang et al. [71] reported a mean DDF of 7.1 mm·d⁻¹·°C⁻¹ for glacier ice in the western China. For the Nepal Himalaya region, the DDF for glacier ice ranges from 6 to 8 mm·d⁻¹·°C⁻¹ [72,73]. This means when daily mean air temperature is positive, 1 °C daily mean air temperature rise will lead to a runoff increase about 6 to 8 mm per day from ice melting water. Therefore, glacier melting is very sensitive to climate warming and uncertainties in Ta estimation would directly lead to uncertainties in runoff simulation at glacierized basins. Specific for this study, uncertainties in estimated Ta may be attributed into three aspects: (1) the error introduced by the referencing CLDAS dataset; (2) the error in MODIS land surface temperature retrieval; and (3) the error introduced by the downscaling algorithm. How these uncertainties in Ta contribute to uncertainty in glacier melting runoff simulation deserves future exploration.

Figure 11 shows the relationship between daily mean Ta and discharge measured at the Kyanging hydrological station. Limited by the data in hand, only the kyanging station observed Ta and discharge in 2004 and the PML estimated Ta at the Kyanging station and station observed discharge in 2004 and from 2008 to 2010 are presented in this figure. Both the station observed Ta and the PML based Ta show exponential relationships with the observed discharge, which is generally consistent with the findings by Chen et al. [74].

5. Summary and Conclusions

Remote sensing based Ta estimation is an important approach to obtain distributed Ta information at glacierized basins. Traditionally, people use the temperature lapse rate to obtain distributed Ta for hydrological simulation at glacierized basins from several or even one station. At least three main problems exist in this method. The first one is that this method assumes the lapse rate is a constant over the whole watershed, and only considers the influence of elevation on Ta and ignores the impacts from other factors, such as the land cover and terrain aspect, which limits the application of this method over very heterogeneous surface. The second problem is that, in most applications, the temperature lapse rate is assumed to be unchanged throughout the year, which may be not the reality. Due to the impacts from annual variation of vegetation phenology and incoming solar radiation, the lapse rate may be not constant across the whole year. The third problem is that this method can be hardly calibrated for many ungauged glacierized basins. To overcome these issues faced by the traditional Ta interpolation method for hydrological simulation at glacierized basins, a new scheme for estimating Ta using MODIS LST, DEM, DDL and reanalysis climate dataset as input parameters is proposed in this study. The main input parameters of our method are from remote sensing products and climate reanalysis dataset, and no additional station observations are needed for calibration, which makes the proposed method applicable in ungauged basins.

Cloud contamination is an inevitable problem in all terrestrial remote sensing products from visible to thermal infrared channels. The appearance of cloud will cause the missing data problem in MODIS LST products and there are also uncertainties in the cloud identify algorithm. Cloud shadow will also change the surface energy balance and may introduce large error in the parameter retrievals. In this study, to obtain the spatio-temporal continuous LST, multi-year average LST and the DINEOF method have been used to reconstruct the MODIS LST. Although the reconstructed LST could satisfy the requirements for Ta estimation, some errors are also introduced in the reconstruction process, and may finally influence the estimation accuracy of Ta.

To validate the new proposed method, we tested our scheme over the Langtang Valley and validated the estimated Ta using observations from two automatic weather stations at the Kyanging and Yala stations. The validation results for the SML regression method show that, at the Kyanging station, the RMSE and bias are 2.21 and 0.17 K for daily maximum Ta, are 2.55 and −1.22 K for daily mean Ta, and are 4.77 and −3.83 K for daily minimum Ta, respectively. At the Yala station, the RMSE and bias are 4.81 and 4.30 K for daily maximum Ta, are 2.93 and 2.12 K for daily mean Ta, and are 2.84 and −0.37 K for daily minimum Ta, respectively. To consider the seasonal variability of the relationship between Ta and LST, PML method is further adopted to build models for summer and winter, separately. The validation results for the PML regression method show that, at the Kyanging station, the RMSE and bias for daily maximum Ta are 2.05 and 0.42 K, are 1.88 and −0.68 K for daily mean Ta, and are 3.36 and −2.86 K for minimum Ta, respectively. At the Yala station, the RMSE and bias for daily maximum Ta are 4.53 and 4.03 K, are 2.68 and 1.96 K for daily mean Ta, and are 2.36 and −0.35 K for daily minimum Ta, respectively. The Ta obtained by the PML method has better performance than the SML method except for the mean Ta at the Kyanging station and the maximum Ta at the Yala station in winter. Our results show that the proposed method is applicable for Ta estimation for hydrological simulation at ungauged glacierized basins.