Dynamics of Land Surface Temperature in the Central Tien Shan Mountains

Permafrost conditions in mountain ranges are sensitive to regional land surface temperature (LST), among other factors. To explore that relationship, this study carried out 3 steps: (1) validated Moderate Resolution Imaging Spectroradiometer 1-km daily LST data using data measured in situ, (2) used the Harmonic Analysis of Time Series (HANTS) algorithm for fitting and removing the influence of clouds, and (3) analyzed the spatial and temporal characteristics of LST dynamics in the central Tien Shan mountain range based on remote-sensing data improved by covariance and empirical orthogonal function analysis. The results indicate that the in situ data present a basic reference for rebuilding invalid values in the retrieved data, and the data gap in daily LST products can be logically reconstructed with the HANTS algorithm. Major long-term and large-scale patterns can be well extracted with the reconstructed LST data. The most dynamic and sensitive LST areas occurred in the buffers around the periglacial areas. Areas above the periglacial line mainly exhibited a decrease in LST, while areas below it showed an increase. This suggests that the periglacial line of the central Tien Shan region has risen during the past decade. These findings can provide a reference for how periglacial areas respond to climate change and how this may affect hydrological and ecological processes.


Introduction
Earth's climate has warmed by approximately 0.6uC over the past 100 years, and the ecological responses to recent climate change attract increasing attention (Walther et al 2002). Climate warming has a direct impact on land surface temperature (LST) and speeds up the thaw of permafrost, which subsequently affects soil organic matter degradation, hydrology, and the carbon budget (Lawrence and Slater 2005;Zimov et al 2006). Continued or even accelerated future warming is likely to induce further retreat and degradation of high-elevation permafrost (Haeberli et al 1993).
As "water towers" fed by precipitation and "solid reservoirs," the high Central Eurasia mountains play important roles for the surrounding lowlands (Messerli et al 2004). LST is one of the key parameters in the physics of surface processes, combining surface-atmosphere interactions and energy fluxes (Wan et al 2002). The dynamics of permafrost in mountain ranges are sensitive to the regional LST and affected by multiple factors (eg elevation, slope, aspect, soil structure, snow, and vegetation cover); thus, they obviously often exhibit spatial heterogeneity.
Since little traditionally observed information exists for high-elevation mountain regions, scientific challenges in spatiotemporal surface temperature studies at different scales remain, and LST could serve as an indicator helping to determine the thermal regime of permafrost (Hachem et al 2009). Due to the remoteness of most permafrost areas, monitoring surface processes through remote sensing is desirable (Running et al 1994;Westermann et al 2011); it offers the most feasible, consistent, and accurate means of identifying land surface parameters (Sellers et al 1997;Gruber and Hoelzle 2001;Hall et al 2006). LST retrieval using infrared satellite measurements has been widely discussed (Prata 1993;Ulivieri et al 1994;Qin et al 2001;Trigo et al 2008). A number of studies of the retrieval of remotely sensed LST have been performed (Wan et al 2002;Mao et Trigo et al 2008), including multichannel approaches (Ulivieri et al 1994) and application of the mono-window algorithm to Landsat Thematic Mapper data (Qin et al 2001).
In recent years, Moderate Resolution Imaging Spectroradiometer (MODIS) LST data sets have been widely used for monitoring, mapping, and modeling the spatial distribution and active-layer dynamics of permafrost (Hachem et al 2009;Westermann and Langer 2010) as well as permafrost mass balance at meso to macro scales (Hall et al 2006;Sun and Kafatos 2007). After an initial evaluation we found the day-LST data to be the most applicable for bare and sparsely vegetated areas (Wan et al 2002;Wang et al 2007). The newly released MODIS LST Version 5 products provide an opportunity to improve this application. The accuracy of most algorithms of LST retrieval from thermal infrared is high (Qin et al 2001;Wan et al 2002). In mountainous permafrost regions, numerous climate and terrain features operate singly and in combination; in particular, the thermal infrared band is susceptible to the impact of clouds on the mountain peaks. Therefore, a field validation and continuous time series reconstruction is essential before application of MODIS LST products.
Using high-frequency automatic data recording, this study investigated the difference between satelliteretrieved and in situ-recorded data, which provided a reference for reconstructing data and analyzing spatiotemporal patterns and scale. This remote-sensing analysis enabled us to explore the complex surface temperature dynamics in terms of distribution and evolution of the collective features of heterogeneity in mountainous regions Xu et al 2013).
The aims of this study were (1) to validate daily MODIS LST data at 1-km resolution using data recorded in situ, (2) to rebuild cloud-free daily LST data using the Harmonic Analysis of Time Series (HANTS) algorithm, and (3) to analyze spatial and temporal LST patterns based on the reconstructed MODIS LST data. Statistical methods adopted for the analysis were covariance (CoV), slope of linear regression and empirical orthogonal function (EOF). This integration of statistical methods leads to a better understanding of LST regimes in mountainous regions and their implications for environmental change.

Study area
Located in the central region of the Eurasian continent, the alpine permafrost zone in the central Tien Shan ) belongs to the Asian high-mountain permafrost region (Marchenko et al 2007). Its elevation ranges mainly from 2500 to 7000 m in the Aksu River Basin, including the largest glaciers located on the highest mountain, Tomur Peak (7435 m). The field investigation area is on a tributary of the Aksu River near Tomur Peak (Figure 1). With a contribution of more than 70% to the total runoff of the Tarim River, the Aksu River has by far the largest impact on water resources and development of the region.
Stretching across the China-Kyrgyzstan border, the central Tien Shan is an important water source for both the north and south, and nearly 90% of the Tarim Basin population inhabits the outskirts of the mountains (Sun, Chang, and Opp 2010). There, sound water management has as much to do with water demands as it does with the effects of climate change on water availability. The identification of LST characteristics is fundamental for large-scale modeling of the state of the permafrost in the central Tien Shan. But the dynamics and amount of water contribution from the periglacial belt are largely unknown, and gauging data are scarce throughout the region. Given its extraordinary topographic heterogeneity and high-mountain ridges, the estimate of future freshwater availability in light of climate change is a key issue for agricultural planning, food security, and urbanization in western China.

Data In situ measurements
To perform a detailed LST study of the alpine permafrost, a monitoring network of 46 thermistor strings and 23 mini temperature data loggers (built-in PT1000 sensors with a resolution of 0.01uC) was installed on the southern regions of the central Tien Shan. The loggers were used to measure the temperature at multiple depths in the active layer at hourly intervals. To identify the factors with the largest influence on LST, the monitoring locations were carefully chosen in terms of elevation, slope, aspect, and other attributes. Because the in situ measurements are affected by snow and surface soil, for a better comparison with remotely sensed data we chose only time series of 5 data loggers that were less affected by these factors for this study.

MODIS LST data
Two MODIS sensors are aboard the Terra and Aqua platforms, which were launched in 1998 and 2002, respectively. These platforms have provided a new and improved capability for terrestrial satellite remote sensing aimed at meeting the needs of climate change research. In this study, MYD11A1 products were adopted. The MODIS LST products provide per-pixel temperature values in a sequence of swath-based to grid-based global image data. Due to cloud cover, there are always many invalid values in the daily MODIS LST products. This has been a major barrier for the application of remotely retrieved data. More than 7300 images of 1-km daily LST data were adopted as raw material for this analysis. Brief descriptions of the LST data products are available from Wan (1999Wan ( , 2009).

HANTS algorithm
The numerous data gaps in the MODIS daily LST product create a need for data interpolation. The HANTS algorithm has been applied to vegetation index and LST to rebuild cloud-free images by harmonic analysis (Roerink et al 2000;Julien et al 2006;Xu et al 2009). Clouds always have a negative influence on the Normalized Difference Vegetation Index, and therefore, taking the maximum value over a limited period tends to remove most cloud-contaminated observations. In HANTS, curve fitting is applied iteratively; first, a least squares curve is computed based on all data points, and next, the observations are compared to the curve.
Observations that are clearly below the curve are candidates for rejection due to cloud cover, and the points that have the greatest negative deviation from the curve are removed first. Next, a new curve is computed based on the remaining points, and the process is repeated. Pronounced negative outliers are removed by assigning them a weight of 0, and a new curve is computed. This sequence of iterations finally results in a smooth curve that approaches the upper boundary of the data points (Wen et al 2004). In this way, cloudy observations have been removed. In this study, the HANTS algorithm was used on daily LST data, and monthly mean data were then generated from the reconstructed daily LST data for long-term and largescale analysis.

CoV and slope
In statistics, CoV is simply a value calculated from the mean and the standard deviation of the LST time series in each pixel, using the following formula: where X i is the i th sample (value) of a variable (pixel), and n is the number of samples. This method has been widely used to determine the spatial variations of land surface variability (Weiss et al 2001). To finally pin down the LST dynamics, the interannual CoV derived from monthly LST values and the interannual CoV derived from LST values were calculated, respectively. The interannual CoV may reflect the overall change of the LST in the given period of time.
To detect and quantify the continuous changes of LST in each pixel, a linear regression analysis between LST time series and time from 2003 to 2012 was used. Such a linear regression makes it possible to reveal the long-term trend of CoV, which is indirectly but intimately related to LST dynamics. Mathematically, it determines the equation of the line that best fits the set of LST values for each pixel.

Y~a Xzb
where a is the slope of the best-fitting line. The slope a may act as a sustainability indicator for detecting changes in LST dynamics in the study area. If the slope of this regression exhibits a negative sign, showing a statistically decreasing trend over time, it may be concluded that the temporal variability of LST is declining; conversely, if the slope exhibits a positive sign, showing a statistically increasing trend, it may be concluded that the temporal variability of LST is increasing.

Empirical orthogonal function
The main purpose of the EOF is to carry out a linear transformation of the original data, producing a new set of orthogonal functions that exclude redundant information and extract the embedded patterns (Huete et al 2002). EOFs have been widely used in atmospheric science (Lorenz 1970;Dommenget and Latif 2002). Some of the most important oscillations in the climate system were identified using EOF analysis (Zhu 1996). For a spatiotemporal field, the mathematical form of the EOF can be defined as follows: where i 5 1, ..., m; j 5 1, ..., n; m is the number of sites (or grids); n is the time series length; Q ij are the i th components of the j th random vector for the centralized and normalized data; U ki are the weight coefficients representing the contribution of the k th component at the i th site (in other words, the components of the eigenvectors of the correlation matrix); z kj are the timedependent functions of the k th component of expansion (the so-called amplitude functions). Note that the weight coefficients U ki vary between the time series data sets (or between different sites) but are constant in time.
The EOFs are the eigenvectors. The relative importance of any individual EOF to the total variance in the field is measured by its associated eigenvalue. In practice, we often sort eigenvalues and corresponding eigenvectors in decreasing order using the first several EOFs to explain the principal variance. Each EOF is associated with a series of time coefficients that describe the time evolution of the particular EOF. The term "EOF" is also interchangeable with geographically weighted principal component analysis in geophysics.

Validation of the MODIS LST
We chose 5 in situ loggers and the MODIS LST imagery grids in which the loggers were situated for use in this analysis and then compared the LST values the in situ loggers recorded (around 14:00 h Beijing time) to valid MODIS-retrieved daily LST values from 1 September 2010 to 31 August 2012 in order to validate the MODIS data. Thus, the LST time series from the in situ loggers and MODIS retrieval are recorded at almost the same time of day. The relationship between MODIS-retrieved and in situ-measured data is presented in Figure 2; information about the in situ conditions and a comparison in CoV and correlation coefficient of the 2 series are given in Table 1.
As Figure 2 shows, the remote-sensing data series had a higher frequency of fluctuation than the logger data series. The values recorded by the MODIS LST and the loggers were within 4uC, and most of the correlation coefficients between the 2 data sources were greater than 0.9. To some extent, the in situ-measured data were consistent with the mean of the data retrieved via remote sensing. The deviation between the 2 data sets provides an important reference for interpolating missing data in the daily LST products. Daytime discrepancies are strongly impacted by the different heating rates of different elements within a pixel (eg, vegetation types, bare ground) (Trigo et al 2008). Major factors affecting this analysis are as follows: 1. The scale of the logger and MODIS LST grids differed.
Environmental variables can vary strongly within the grid cells. The MODIS LST data are an integration of a 1000 3 1000 m grid, and the logger data cover a much smaller area; thus, validating with a single measurement is difficult and susceptible to bias (Gubler et al 2011). 2. Snow cover occasionally affected the loggers (Bartlett et al 2004). Snow cover substantially increases the daily temperature for a logger sensor; this effect was very significant for the last 2 loggers during the winter.

HANTS simulation of the LST time series
Since the in situ-measured data are consistent to the mean of the data retrieved by remote sensing, rebuilding (interpolate) the missing data on this basis is reasonable. In this study, the HANTS algorithm was applied to the LST data on a per-pixel basis for each year between 2003 and 2010 for the entire study area. Examples of the HANTS-fitted and original LST time series are given in Figure 3. The results indicate that the data gap in the daily LST time series can be logically reconstructed.
Because of the small size of the filled gaps the data rebuilt by this algorithm did not greatly impact the long-term trends. So the reconstructed data are an ideal basis for monthly or long-term statistical analysis for the central Tien Shan region.

Spatiotemporal characteristics of LST
With the LST data reconstructed by HANTS fitting, the monthly and yearly mean LST can be calculated from the FIGURE 2 LST data as retrieved from MODIS and as measured by 5 in situ loggers (nonconsecutive days).   daily 1-km LST data. The statistical analysis indicates that, in general, the yearly mean LST shows a strong inverse relation with elevation in the study area ( Figure 4A). These LST/elevation regions exhibit a buffer with a width of about 1000 m at elevations above 3000 m; the range of the buffer becomes even wider at elevations below 3000 m. The 0uC isothermal region ranges from 3000 to 4200 m in elevation for the central Tien Shan. The distribution of CoV is a good way to determine spatial LST variability. In general, the spatial distribution of the CoV for yearly mean LST exhibits some high-value belts that are consistent with the geomorphology ( Figure 4B). The high-value regions show a clear buffer between the glacier cover and the grassland. A comparison with a digital elevation model for this region indicates that this dynamic buffer has an elevation ranging from 3000 to 4300 m, which is consistent with the periglacial areas in elevation. That means the periglacial areas are the most sensitive regions for LST dynamics. The thaw or advance of the permafrost and the changes in the active layer are much more obvious in these regions. Thus, the effort to understand the variability of permafrost under climate change should be focused on these buffer regions-the "hot lines" of permafrost dynamics.
As Figure 4C shows, the distributions of the CoV and elevation exhibit a bilaterally decreased pattern with a peak for areas with an elevation around 3800 m; the high-LST-dynamic regions are mainly distributed between 3000 and 4300 m (with CoV values higher than 0.5).
The slope of the intra-annual CoV reflects the dynamic trend of LST during the study period. As Figure 4D shows, the high value of slopes is mainly restricted to elevations of 3200-4200 m. This is somewhat consistent with the result of CoV analysis, with the difference that there are also some high-value regions for CoV on the north slopes of the central Tien Shan. Slope results further indicate that areas above the periglacial line exhibit a decreasing trend for LST change, while areas below the periglacial line show an increasing trend. These findings to some extent suggest that the periglacial line of the central Tien Shan region has risen, and as a consequence of atmospheric warming, the lower boundary of permafrost distribution in mountain ranges may have risen as well, causing local degradation of formerly frozen slopes. Synthetic analysis of CoV and its slope indicates that the most sensitive regions for LST dynamics are periglacial regions (3000-4300 m). This is consistent with earlier findings that climate change has led to an accelerated retreat of high mountain glaciers (Marchenko et al 2007;, which in turn has triggered more runoff and an increase in vegetation index .

EOF patterns and Principle Component evolution: spatial patterns and amplitudes
The EOF analysis efficiently extracted large-scale spatial patterns (EOFs) and the corresponding amplitudes, also called Principle Components (PCs), from the reconstructed LST monthly data set. Figure 5 presents the 4 leading EOFs. Spatial patterns for EOF1 and EOF2 exhibit a strong connection with the large-scale geomorphological features (elevation effects) of the central Tien Shan. Two patterns have obviously reversed characters in space (eg the high-value center located at the tops of the mountains for EOF1 and the reverse for EOF2). The spatial patterns of EOF3 and EOF4 present differences based on aspect (eg the high-value region of EOF3 has a southern aspect, while the high-value region of EOF4 has a northern aspect), which reflect the effect of local atmospheric circulation.
The amplitudes (PCs) for EOF1 and EOF2 have similar phases corresponding to the seasonal temperature cycles of the past decade, which indicate the fundamental impacts from climate processes ( Figure 6). PC2 is slightly later and weaker than PC1; this difference was especially pronounced in 2009. The frequency of PC3 and PC4 is higher than that of the first PC1 and PC2, and there are some clear differences in the reaction to the season rhythms. PC3 and PC4 phases also exhibit a somewhat opposite behavior. PC3 shows a decreasing trend from 2003 to 2006 and an increasing trend from 2007 to 2009, while PC4 exhibits almost the reverse. This reflects the amplitude shift of the impacts from the north and south slopes of the central Tien Shan.
For most pattern extraction, the leading patterns make a large contribution to the change as a whole, and these patterns are normally stable in evolutionary processes. From the amplitude of the EOFs, the patterns of fluctuation for EOF1 and EOF2 are relatively stable. This is consistent with earlier findings that the thermal state of permafrost in the Tien Shan region reflects the rise in average air temperature during the 20th century (Marchenko et al 2007). But the peripheries of these patterns are often sensitive regions as a result of patterns overlapping in a given period. This effect is highlighted in Figures 4B and 4D. At the same time, other patterns, which make a trivial contribution to the overall change, may contribute to important changes in a certain period and location, and proper attention should be paid to these phenomena.
Earlier studies indicated that climate warming is leading to accelerated ablation and retreat of high mountain glaciers (Lawrence and Slater 2005; Marchenko  (Woo et al 1994;, which may even cause a shift between runoff depth and precipitation . The findings of this study provide insight into how periglacial areas have responded to climate change. Based on long-term air temperature and snow cover data, there is a way to reconstruct the thermal state of permafrost and project these findings to evaluate the contribution of permafrost to water discharge in a given period.

Conclusions
Data analysis showed that the in situ LST data provided a useful reference for rebuilding invalid values in the retrieved LST data. The data gap in the 1-km daily LST time series can be logically reconstructed with the HANTS algorithm, which provides a good basis for studying the spatial and temporal patterns of LST dynamics. The long-term and large-scale patterns of LST in the high mountain region of central Tien Shan can be successfully identified with the reconstructed data and statistical methods discussed in this article. These largescale spatial patterns of LST are generally influenced by geomorphological features, such as elevation and aspect, and by local atmospheric circulation. The CoV of yearly LST and its intra-annual slope clearly show a sensitive dynamic buffer around the periglacial areas with elevation ranging from 3000 to 4300 m. LST mainly decreased above the periglacial line in the central Tien Shan and increased in areas below it. This indicates that the periglacial line region has risen during the past decade. As a consequence of atmospheric warming, the lower boundary of permafrost distribution in mountain ranges may have risen as well, causing local degradation of formerly frozen slopes.