Spatial and Temporal Variations of Land Surface Temperature Over the Tibetan Plateau Based on Harmonic Analysis

Abstract Land surface temperature (LST) is an essential parameter in the physics of land surface processes. The spatiotemporal variations of LST on the Tibetan Plateau were studied using AQUA Moderate Resolution Imaging Spectroradiometer LST data. Considering the data gaps in remotely sensed LST products caused by cloud contamination, the harmonic analysis of time series (HANTS) algorithm was used to eliminate the influence of cloud cover and to describe the periodical signals of LST. Observed air temperature data from 79 weather stations were employed to evaluate the fitting performance of the HANTS algorithm. Results indicate that HANTS can effectively fit the LST time series and remove the influence of cloud cover. Based on the HANTS-derived mean term and annual harmonics, annual mean LST, seasonal fluctuation, and peak time of the LST annual cycle are discussed. The spatial distribution of annual mean LST generally exhibits consistency with altitude in the study area, and the spatial distribution of seasonal oscillation is closely related to precipitation. However, the timing of the peak LST does not exhibit an obvious regular pattern. The LST characteristics of different land cover types were also studied. Bare land has the highest mean LST and exhibits remarkable seasonal fluctuation. Snow, ice, or both show the lowest mean temperature, and forest shows the weakest seasonality. Different land cover types also reflect different peak occurrences of the LST annual cycle, with grassland showing the lowest annual phase value. This paper provides detailed information on the LST variations on the Tibetan Plateau, with the cloud contamination removed. The HANTS algorithm is demonstrated to be effective for understanding spatiotemporal variations of remotely sensed LST, especially for regions over which dense clouds cause large gaps in the LST data.


Introduction
Land surface temperature (LST) is an essential parameter in the physics of land surface processes: it plays a key role in the energy and water transfers between the ground and the atmosphere (Jin et al 1997). LST is controlled by solar radiation and the land-atmosphere heat exchange (Manzo-Delgado et al 2004). Therefore, its spatial and temporal distributions reflect not only the variations of climate factors but also the land surface characteristics. An in-depth understanding of the spatial and temporal variations of LST is important to a range of research fields, including climate, vegetation, and hydrology.
Traditional LST data are collected by meteorological observation and then interpolated to gridded data. However, in some areas, especially in remote areas, spatial interpolation cannot provide satisfactory results from point observations due to the sparse density and uneven distribution of meteorological stations. Satellite remote sensing provides a straightforward and consistent way to observe LST over large scales with more spatially detailed information. During the past decades, LST retrieval has become a hot topic in remote sensing studies. Some classical approaches, such as split-window algorithm, monowindow algorithm, and single-channel algorithm, have been developed and have achieved satisfactory accuracy (Price 1984;Becker and Li 1990;Wan and Li 1997;Qin et al 2001;Jimenez-Munoz and Sobrino 2003). Based on LST data derived from remote sensing, several studies have been conducted to study the spatial and temporal variations of LST. Manzo-Delgado et al (2004) analyzed LST changes in central Mexico in the dry seasons from 1996 to 2000 and the relationships between LST and forest fires. Julien et al (2006) studied LST dynamics over the European continent from 1982 to 1999 using Advanced Very High Resolution Radiometer data.

Systems knowledge
Mountain Research and Development (MRD) An international, peer-reviewed open access journal published by the International Mountain Society (IMS) www.mrd-journal.org Westermann et al (2011) analyzed the spatial and temporal variations of the surface temperatures of higharctic tundra in the summer and assessed the long-term impact of snow cover, soil moisture, and surface cover on surface temperature. Van De Kerchove et al (2011) assessed the spatiotemporal variability of LST in the Russian Altai Mountains and its relationship with physiographic variables.
As the highest and largest highland in the world, the Tibetan Plateau plays an important role in regional and global climate through thermal forcing mechanisms (Yanai et al 1992). The plateau's surface absorbs a large amount of incoming solar radiation energy and undergoes seasonal variations of surface heat fluxes (Ye and Gao 1979). The quantitative analysis of LST is important for understanding the land-atmosphere energy balance of the plateau and its effects on climate. In recent decades, many studies have been performed to analyze the temperature changes of the Tibetan Plateau using meteorological observation data (Lin and Zhao 1996;Liu and Chen 2000;Du 2001;Xu et al 2003;Li et al 2005;Duan et al 2006;Liu et al 2006Liu et al , 2009Ding and Zhang 2008;Wang et al 2008a;You et al 2008;Cuo et al 2012). Due to the sparse density and uneven distribution of meteorological stations on the plateau, the spatial pattern of temperature cannot be properly depicted. Some efforts have been devoted to analyzing LST variations by remote sensing (Jiang and Wang 2001;Ma et al 2003;Oku et al 2006;Salama et al 2012;Tseng et al 2011;Zhong et al 2011). Most of these studies used coarse resolution remotely sensed LST data, which cannot represent detailed spatial patterns. In addition, the incompleteness of the remotely sensed LST caused by cloud cover has not been considered in these studies.
Considering the Tibetan Plateau's unique topography and its thermal forcing mechanism on atmospheric circulation, a study of the LST variations in this region is essential for understanding the energy and water cycle; therefore, it has important scientific significance in regional and global climate research. The main objective of this paper is to analyze the spatiotemporal variations of LST on the Tibetan Plateau using harmonic analysis, which can effectively eliminate cloud contaminations and depict fundamental periodic variations. The 1-km resolution Moderate Resolution Imaging Spectroradiometer (MODIS) LST product from 2003 to 2010 was employed in this study.

Study area
The Tibetan Plateau is located in western China and ranges from 73.44u to 103.98uE and 26.87u to 39.81uN. The plateau has a total area of approximately 2,530,000 km 2 and covers more than one fourth of the land area of China ( Figure 1). Known as the ''third pole of the earth'' (Qiu 2008), the Tibetan Plateau is the highest and largest plateau in the world, with an average elevation of more than 4000 m above sea level. As a result of the high altitude, the climate of the plateau is quite different from that of other regions at nearly the same latitudes. The mean annual air temperature varies between 22.8 and 11.9uC, depending on the topographic conditions. The annual precipitation decreases from more than 800 mm in the southeast to lower than 200 mm in the northwest, with most precipitation occurring from May to September. Influenced by the climate and altitude, land cover on the Tibetan Plateau also shows a clear geographic distribution, following the order of forest, grassland, and desert from the southeast to the northwest (Wang et al 2005).

Data set
The AQUA MODIS LST product (MYD11A2) spanning 2003 to 2010 was obtained from the Warehouse Inventory Search Tool (WIST 2011). MYD11A2 is an 8-day composite daytime and nighttime LST product with a 1-km spatial resolution that is derived from AQUA MODIS data. The LST product is first computed using the generalized splitwindow algorithm proposed by Wan and Dozier (1996) and then produced by averaging the LST values under clear-sky conditions over an 8-day period (Wan et al 2002) for the purpose of removing cloudy pixels. In this study, the daytime LST was employed with an overpass time of approximately 1:30 PM (local solar time).
First, the MODIS daily LST product was computed using the generalized split-window algorithm (Wan and Dozier 1996): T 31 zT 32 2 where T s is LST; T 31 and T 32 are the brightness temperatures in MODIS band 31 and 32, respectively; e 31 and e 32 are the surface emissivities in these 2 bands estimated using the classification-based emissivity method (Snyder et al 1998), according to land cover types determined by the MODIS land cover product (MOD12Q1) and snow cover product (MOD10_L2); and C, A 1 , A 2 , A 3 , B 1 , B 2 , and B 3 are regression coefficients given by multidimensional lookup tables. Then, MYD11A2 was generated by averaging the daily LST values under clear-sky conditions over an 8-day period (Wan et al 2002) for the purpose of removing cloudy pixels. Each year of the MYD11A2 LST data contained 46 eight-day images.
The MODIS LST product (V005) were evaluated using in situ observations from the Asia-Australia Monsoon Project on Tibetan Plateau (CAMP/Tibet), with an average error of 2.36 K . Compared with the validation results of MODIS LST in other regions (Wan 2008;Wang et al 2008b;Crosman and Horel 2009), the accuracy error in the Tibetan Plateau is somewhat higher, which may be due to the complicated topography and atmosphere conditions, but it is still acceptable for practical applications.
The corresponding MODIS land cover product (MCD12Q1) was also used in this study. The MCD12Q1 product incorporates 5 land cover classification schemes derived from supervised decision-tree classification methods (Friedl et al 2002). We used the International Geosphere-Biosphere Program (IGBP) classification scheme, which contains 17 land cover classes (water, evergreen needleleaf forest, evergreen broadleaf forest, deciduous needleleaf forest, deciduous broadleaf forest, mixed forests, closed shrubland, open shrubland, woody savanna, savanna, grassland, permanent wetland, cropland, urban and built up, cropland/natural vegetation mosaic, snow and ice, and barren or sparsely vegetated). Using MODIS Reprojection Tool software, all MODIS data were mosaicked and reprojected. In addition, elevation data of the plateau were acquired from the Shuttle Radar Topography Mission 90-m digital elevation model (SRTM3 DEM, version 4.1) and were resampled to a 1-km resolution to coincide with the MODIS data.
Daily maximum air temperature (T a ) observation data of 79 weather stations from 2003 to 2010 in the research area ( Figure 1) were employed to evaluate the fitting performance of harmonic analysis. These T a data were averaged over each 8-day period for comparison with the MODIS LST data.

Methodology
Although MYD11A2 is an 8-day composite LST product, it includes numerous invalid values (identified as fill values) caused by long-term cloud cover and anomalous values caused by unidentified cloud cover. These problematic values need to be handled carefully because they do not represent actual surface conditions. To eliminate the influence of these cloud-affected values in the LST data and extract the characteristics of the LST dynamics, a harmonic analysis algorithm, named the harmonic analysis of time series (HANTS), was performed on the LST time series. The HANTS algorithm (Menenti et al 1993;Verhoef et al 1996;Roerink et al 2000) was developed to remove cloud-contaminated values in normalized difference vegetation index (NDVI) time series. The HANTS algorithm allows for greater flexibility in the choice of frequencies and the length of the time series than the fast Fourier transform (FFT) algorithm. In addition, it is relatively easy to exclude certain points from the time series because the samples are not required to be equidistant in time. By removing the obvious outliers in the time series, the extracted harmonics are far more reliable than the straightforward FFT algorithm (Roerink et al 2000).
The HANTS algorithm applies a least squares curvefitting procedure based on harmonic components by considering the most important frequencies in the time profiles. The fitted curve in the HANTS transform is described as the sum of its mean value and several cosine functions with different frequencies, where y(t) is the fitted curve value at time t, a 0 is the average value of the time series, N is the number of harmonics, a i is the amplitude of harmonic i, v i is the frequency of harmonic i, and h i refers to the phase of harmonic i. For each harmonic, the amplitude and phase of the cosine function is determined during an iterative fitting procedure. The HANTS iterative process can be summarized as follows. First, a least squares curve is computed based on all data points within the valid range. Then, the observation points are compared with the fitted curve. The points that significantly deviate from the current curve (eg cloudy and missing pixels) are identified and removed. Afterward, the procedure is repeated based on the remaining points until no points exceed the predefined error threshold or the number of remaining points is too small. Finally, the iteration leads to a smooth temporal curve. In this way, the cloud contamination is removed, and the amplitudes and phases of the harmonics are extracted (Roerink et al 2000;Julien et al 2006).
The HANTS algorithm was presented for the application of NDVI time series and has been successfully applied for reconstructing cloud-free NDVI data sets and for studying vegetation phenology (Verhoef et al 1996;Roerink et al 2000;Wen et al 2004;Julien et al 2006;Julien and Sobrino 2010;de Jong et al 2011). Moreover, the HANTS transform can also be applied to LST time series because LST fluctuations exhibit periodic patterns and can theoretically be fitted by several harmonics.

HANTS simulation of the LST time series
To perform the HANTS analysis, 5 parameters should be set: the number of frequencies, a high/low suppression flag, a valid data range, the fit error tolerance, and the degree of overdeterminedness. A detailed description of these 5 HANTS parameters can be found in Roerink et al (2000). In this study, the number of frequencies was set to 2 (frequency 1 5 annual; frequency 2 5 semiannual). The high/low suppression flag was set to low, considering that cloud contamination leads to abnormal low LST values. The valid data range was from 250 to 70uC. The fit error tolerance was set to 6uC, and the degree of overdeterminedness was set to 7.
The HANTS algorithm was applied to the LST data on a per pixel basis in each year between 2003 and 2010 for the entire study area. Two examples of the HANTS-fitted and the original LST time series are given in Figure 2, whose locations are also labeled in Figure 1 as Pt 1 and Pt 2, respectively. In this figure, the light gray bars are the original LST temporal profiles, and the solid line is the HANTS-fitted LST temporal profile. The outliers in the original LST temporal profile that were rejected according to the fit error tolerance in the HANTS iterative procedures are shown as dark gray bars. Figure 2A represents an LST temporal curve that contains only 3 data gaps, and Figure 2B represents an LST temporal curve that contains 9 data gaps. The simulation results show that HANTS can properly and robustly fit the LST temporal curves even if there are conspicuous data gaps.
Observed maximum air temperature from 79 weather stations was compared with the original and HANTSfitted LST to evaluate the fitting performance. Figure 3A shows the relationship between valid values in the original LST and corresponding T a values. There is an obvious linear correlation between these 2 data, with a correlation coefficient of 0.7885 (P , 0.01). The LST values are generally higher than T a , but some LST values are lower than T a . Considering the overpass time of AQUA MODIS, these pixels most likely contain undetected residual clouds, which result in abnormally low surface temperature. Figure 3B shows the relationship between HANTS-fitted LST and T a . The light gray triangles are the fitted LST values for the invalid values in the original LST caused by cloud cover. These triangles show similar relationship between LST and T a as the dark gray diamonds in the figure, indicating that no significant inconsistencies were introduced in these fitted values. The correlation coefficient between HANTS-fitted LST and T a is 0.7751 (P , 0.01), which is slightly lower than that between original LST and T a and suggests the good performance of HANTS fitting. In addition, the abnormally low values in the original LST do not show very low surface temperatures relative to air temperature, indicating that the influence of residual clouds on the original LST are effectively eliminated. The success of HANTS fitting can also be testified by LST spatial distributions. Figure 4 shows the original and HANTS-fitted LST map in the period of 361 days of year (DOY) 2010. In the original LST map, there are several large data gaps in the northeastern, northwestern, and southeastern part of the plateau, which are caused by cloud cover. The HANTS-fitted LST map shows good consistency with the original LST map for the clear-sky pixels, and the missing values in the original LST are effectively filled.
The fitting results were saved as 5 harmonic components: the mean LST and the amplitude and the phase of the first (frequency 5 1) and the second (frequency 5 2) harmonics. The mean LST and amplitudes were represented in units of degrees Celsius, and the phases were represented in units of day of year. An example of the fitted LST curve and its harmonic components is given in Figure 5. The HANTS-fitted LST curve and the mean LST term are shown on the left axis, while the first and second LST harmonic curves are shown on the right axis. The HANTS-fitted curve is achieved according to the summation of the mean LST and the 2 harmonics. The first and second harmonics represent the annual and biannual cycles of the LST time series, respectively, according to their frequencies. Their amplitudes denote the periodic fluctuation range of the LST time series based on the mean value, and their phases illustrate the peak time of the LST periodic cycles, which indicate the temporal character of the LST fluctuation. Figure 5 proves that the first harmonic shows a much larger amplitude than the second harmonic. Moreover, the mean amplitude of the first harmonic of all pixels is 12.14uC, which is also much higher than that of the second harmonic (2.86uC). This result suggests that the first harmonic (annual term with frequency 5 1) contributes far more to the overall LST time series than the second harmonic (semiannual term with frequency 5 2), indicating a strong unimodal periodic pattern of the LST temporal series on the Tibetan Plateau. The second harmonic contributes less to the overall fit but helps improve the HANTS fit performance.

Spatiotemporal variations of LST
The HANTS simulation result of the LST time series is composed of the mean value of the LST temporal profile, amplitudes, and phases of the 2 harmonics. According to the preceding analysis, the mean and annual harmonic contribute more to the overall LST time series than the semiannual harmonic. Therefore, the 8-year (2003-2010) averaged mean LST and annual harmonic were used to represent the spatiotemporal variations of LST in the Tibetan Plateau. The mean term represents the annual mean LST value. The annual harmonic depicts the annual cycle of LST, its amplitude illustrates the seasonal fluctuation magnitude, and its phase illustrates the peak time of the LST annual cycle. Figure 6 provides the 8-year (2003-2010) averaged mean LST on the Tibetan Plateau. The study area exhibits a highly contrasted spatial pattern of annual mean LST values, suggesting significant regional variabilities. The annual mean LSTs in most areas were generally within the range of 0 to 30uC, but some pixels were even larger than 30uC or lower than 0uC. Consulting the DEM map of the plateau (Figure 1), the spatial distribution of annual mean LST is influenced by the altitude. Areas with lower altitude generally have higher annual mean LST, indicating that higher altitude leads to lower temperature. The Qaidam Basin located in the northeast shows higher annual mean temperatures above 28uC due to its lower altitude and desert landscape. However, the annual mean LSTs of the South Tibetan Valley in the southwest, with relatively high altitudes ranging from 4500 to 5000 m, are also among the highest in the study area. This can be attributed to the area's lower latitude and valley terrain and to the influence of the warm Southwest Monsoon. The Kunlun Mountains, located in the northwest, and the Nyainqentanglha Mountains, located in the south, show lower annual mean temperatures below 3uC. Large lakes in the plateau, such as Qinghai Lake and Nam Lake, also show relatively lower annual temperatures. Figure 7A and B represents the 8-year (2003-2010) averaged amplitude and phase, respectively, of the annual cycle of the LST series on the Tibetan Plateau. Figure 7A clearly reveals that the annual amplitude shows a spatial tendency that increases from the southeast to the northwest, which shows a similar spatial pattern with precipitation ( Figure 7C). According to the surface heat balance, LST is controlled by solar radiation and surface thermal properties. High precipitation caused by the monsoon indicates heavy clouds, mainly in summer. The occurrence of clouds reduces solar radiation and in turn reduces the surface temperature, especially in summer. Lower summer temperatures decrease lower seasonal  fluctuations. Moreover, in this arid and semiarid area, higher precipitation usually results in higher vegetation cover and higher soil moisture, which are more resistant to seasonal oscillations of LST. Despite this general southeast-northwest difference, the most striking seasonal oscillation occurs at the Qaidam Basin, which is ascribed to its sparse vegetation cover, low soil moisture, and low altitude. The annual amplitude in this basin is usually higher than 15uC. In the Yarlung Zangbo Grand Canyon located in the southeast, the annual amplitude is mostly lower than 5uC, indicating slight seasonality. As the main water vapor channel of the Tibetan Plateau (Lin and Wu 1990), this area is highly influenced by the warm-humid Southwest Monsoon and shows a warm-temperate monsoon climate. Figure 7B demonstrates the spatial distributions of peak time for the LST annual cycle. Most of the study area has annual phase values ranging from 148 to 204 DOY, suggesting they achieve their peak value of the LST annual cycle at June and July.

MountainResearch
In contrast to the amplitude, the spatial distribution of the annual phase does not present an obvious regular pattern, so it may be controlled by complex interactions of environmental variables, including solar radiance, underlying surface type, topography, vegetation abundance, climate, and soil moisture conditions. On the whole, the annual phase exhibits an inconspicuous southeast-northwest spatial tendency. The largest annual phases were found over large water bodies, suggesting that these areas have the latest peak time of LST due to their high heat capacities. High mountains, such as the Karakoram Mountains and Nyainqentanglha Mountains, also exhibit delayed peak times of LST for their permanent snow cover. Compared with the DEM map, there is an inconspicuous analogy between the spatial distributions of altitude and the annual harmonic, suggesting that altitude is not the decisive factor for the annual variations of the LST peak time.

LST variations for different land cover types
To analyze the LST dynamic characters of the different land cover types in the study area, the original MODIS land cover map was reclassified by aggregating the 17 IGBP categories into 5 major land cover classes: forest, grassland, water, snow/ice, and bare land ( Figure 7D). Based on the land cover map, the average values of the annual mean LST and the amplitude and phase of the first harmonic for each land cover type were calculated (Table 1). Table 1 shows that bare land exhibits the highest annual mean LST (18.00uC), followed by grassland (17.24uC), and forest (16.62uC).
Vegetation cover has a dampening effect on LST due to the cooling effect of shading and evapotranspiration. For example, bare land displays a high temperature that can be related to its high insolation (caused by low cloud cover), low heat capacity, and lack of evaporative cooling (caused by low vegetation cover). In contrast, water displays a relatively low temperature because of its rather high thermal capacity and evapotranspiration. Thus, water shows a much lower mean temperature of 6.75uC, and snow/ice shows the lowest mean temperature of 0.53uC.
For the amplitude of the first harmonic, bare land also shows the most distinct annual fluctuation of LST, with an amplitude of 15.63uC, due to its low heat capacity and low evapotranspiration. Grassland, water, and snow/ice show lower seasonalities, with amplitudes of 10.40, 10.02, and 9.81uC, respectively. However, forest exhibits the lowest seasonal oscillation, with an amplitude of 9.48uC. The lowest amplitude of forest is probably due to the evaporation and shelter effect of the vegetation canopy, which efficiently reduce temperature in growing season, especially in summer. For the phase of the first harmonic, grassland and forest show relatively lower respective values of 170.88 and 174.40 DOY, indicating that their maximum LSTs occur on June 19 and June 23, respectively. Forest and grassland show earlier peak times than other land types, which can be explained by the high cloud cover and high vegetation abundance in the growing season (from July to September) that reduce the summer temperature. Bare land reaches its peak LST on 181.12 DOY, approximately 10 days later than grassland. Snow/ice reaches its highest LST on 191.74 DOY, approximately 10 days later than bare land. Water shows the largest amplitude of 200.95 DOY, suggesting that its peak LST occurs on July 19. The lag effect of water surface temperature results from its high heat capacity.

Conclusions
In this paper, the spatiotemporal variations of daytime LST in the Tibetan Plateau were studied using AQUA MODIS 1-km LST data from 2003 to 2010. The harmonic analysis method HANTS was introduced to suppress cloud contamination and discriminate significant periodic variations of LST. Evaluation results indicate that the HANTS algorithm can provide satisfactory fitting for the LST time series and eliminate the influence of cloud cover on the remotely sensed LST data.
The results of the HANTS analysis clearly show that the LST seasonal variation in the study area has a strongly unimodal periodic pattern. Therefore, most of the temporal variance was constrained to the mean term and the annual harmonic. The mean term of the harmonic analysis represents the annual mean LST value, which exhibits a significant regional variability in the study area. The spatial pattern of mean LST is generally influenced by the distributions of altitude and underlying surface. The amplitude of the annual harmonic describes seasonal fluctuations of LST and displays a remarkable southeastnorthwest decreasing gradient, which corresponds to the distribution of precipitation. The phase of the annual harmonic describes the peak time of the LST annual cycle, which does not present such an obvious regular pattern as amplitude. Most of the area achieves its peak value of LST in June and July, and water bodies display the latest peak time. The geographic distributions of both seasonal fluctuation and peak time show little correspondence with altitude. The characteristics of the LST variations for different land cover types were also studied. Bare land shows the highest annual mean LST, followed by grassland, forest, water, and snow and ice. Bare land also exhibits the most prominent seasonal oscillations. Other land covers, such as grassland, water, and snow and ice, show lower seasonalities relative to bare land, and forest exhibits the lowest seasonal fluctuation. For the annual phase, vegetation covers including grassland and forest show their early peak of the LST annual cycle. The bare land and the snow/ice classes reach their peak LST approximately 10 and 20 days later than grassland, respectively. Water shows the latest LST peak, representing an approximately 1-month lag compared with vegetation cover.
This study gives spatially detailed information on the spatiotemporal variability of LST in the Tibetan Plateau, which is meaningful to understand the surface thermal condition and its influence on vegetation, frozen ground, and climate on the plateau. In addition, the HANTS algorithm is shown to be an effective approach for understanding the spatiotemporal variations of LST, especially for the regions with dense clouds that cause a large number of gaps in the LST data.