Recent Observations of Human-induced Asymmetric Effects on Climate in Very High-Altitude Area

Like urban heat islands (UHI), human-induced land degradation (HLD) is a phenomenon attributed to human activities, but this phenomenon occurs in non-urban areas. Although a large body of work has demonstrated that land-cover change influences local climate systems, little work has been done on separating the impact of HLD from naturally-occurring fluctuations in very high-altitude areas. We developed an innovative NDVI-difference method in order to evaluate HLD effects upon the climate system in the central Tibet Plateau. The results show that the minimum temperature increased at a significantly faster pace than the maximum temperature in the growing season at HLD meteorological stations, but this was reversed at stations with natural forces only. Further analysis revealed that abrupt changes of minimum temperature occurred five years earlier and amplitudes of these changes were 1.4 times larger than at stations with natural forces only. Therefore, our results complement other evidence that points to the fact that local effects from UHI contribute to climatic asymmetry observed between minimum and maximum temperature trends. Accordingly, we stress the need for consideration of non-urban factors from anthropogenic activities, such as human-induced land degradation, in understanding these asymmetric diurnal changes.


Introduction
Human activities, such as burning fossil fuels, releasing chemicals into the atmosphere from industrial activities, reducing the amount of forest cover, and rapid expansion of farming, are changing the balance of the climate system [1][2][3][4][5][6][7]. A prominent example is that that the rise of the minimum temperature has occurred at a higher rate than the rise of the maximum temperature since 1950, resulting in a broad decline in the diurnal temperature range (DTR) [8][9][10]. A large body of work has demonstrated that land-cover change provides an additional major forcing of the climate, through changes in the physical properties of the land surface [11][12][13][14]. Since land degradation is a humaninduced process leading to different amounts and rates of climate system variability, a pervasive issue in this topic has been how to separate human-induced climate change from natural climate fluctuations and henceforth help to quantify human influences on global warming [15][16][17][18][19][20]3]. Very high-altitude areas are ideal environments for monitoring human-induced climate change because they have low temperature conditions. However, human influences on the climatic system at high altitudes have a high degree of complexity and a high degree of uncertainty. Difficulties in taking direct measurements in very high altitude areas make this an even more challenging consideration.
Human-induced land degradation (HLD) is a man-made phenomenon that negatively affects the effective function of land within an ecosystem, to accept, store and recycle water, energy, and nutrients [21][22][23]. Changes in satellite-measured greenness are a good indicator of land degradation. The normalized different vegetation index (NDVI), which is a non-linear transformation of the visible (red) and near-infrared bands from remote sensing, has been proven to be capable of highlighting the area with comparatively reduced vegetation activity and thus monitoring the spatial-pattern and magnitude of land degradation [24][25][26][27][28][29]. In this study an innovative NDVI-difference method that utilizes spatially distributed information of time-series NDVI in the growing season, with verification of the actual vegetation conditions, was developed in order to evaluate HLD effects upon climate in the central Tibet Plateau (TP), where long-term pressures resulting from growing human populations and higher numbers of grazing livestock had been putting the actual and/or future capability of land in danger [30]. The Mann-Kendall (MK) test was used for detecting temperature trends and abrupt changes in time series.

Background and Methodology
There are clear suggestions from a number of high-altitude climate records from this century that the amplitude of temperature changes at high altitudes is greater than the observed global or hemispheric change; furthermore, a number of lines of evidence suggest that the warming signal in the tropics during the past few Table 1. (A) Temperature and (B) precipitation indicators for climatic regionalization of the Tibetan Plateau (TP), listing the number of days per year on which average temperature is above 10uC, the average temperature in the hottest month, the aridity index, and the annual precipitation [40]. decades increases with height. The IPCC Second Assessment Report [31], in its chapter on the impacts of climate change on mountainous regions [32], has recommended that ''future research needed to understand and predict effects of climatic change on mountain regions should represent balance and coordination between field studies, monitoring, experimental studies, and modeling.'' The Tibetan Plateau (TP), with the most prominent and complex terrain on the globe and an elevation averaging more than 4000 m above mean sea level [33], is often called the ''Third Pole'' because its geographic significance is akin to that of Antarctica and the Arctic [34]. The TP plays an important role in global atmospheric circulation through orographic and thermal forcing mechanisms. The TP, which is considered to be ''the driver and amplifier of global climate change,'' is the 'start-up region' for climate change in China, as well as for the world [35].
There is evidence that the degradation of grassland has become increasingly serious with growing population and livestock in recent years. Human-induced land degradation has become a widespread problem in the TP, with its large area and wideranging types of degradation. It was estimated that the total area of degraded land in the TP was approximately 4.0610 7 to 6.0610 7 ha in 1990s [36]. Serious degradation is particularly found around water sites and the camps where animals gather in the winter-spring. The causes of land degradation include inappropriate human activities like overgrazing, an imbalanced use of grasslands between winter-spring and summer -autumn and poor management of grasslands. To minimize the influence of non-vegetation land-cover, NDVI values below 0.1, representing snow, inland water bodies, desert and exposed soil, were removed. Furthermore NDVI was only calculated in growing season in order to reduce other factors not related to actual vegetation change. NDVI was sampled for 868 pixel windows centered on each meteorological station involved in this study. Such a pixel window, with a radius of about 32 km, was chosen because it is generally accepted that meteorological station observation data is influenced by mesoscale climate (1 km to 30 km) including proximity and size of vegetation cover types, urbanized areas and other factors [28,29]. It is assumed that each station NDVI is comprised of natural components and anthropogenic components due to averaging the corresponding 868 pixel window. For example, interannual changes in climate factors, especially temperature and precipitation, can profoundly influence vegetation growth and vegetation cover and are thus related to changes on NDVI. On the other hand, clearing of forest for agriculture also influences NDVI. The basic aim is to remove the climate signal in order to isolate the signal of human activities. We assume that NDVI in the 128 km radius area (32632 pixel window) around the station should not be sensitive to the HLD effect on the station. Hence natural components could be derived from the 32632 pixel window. Consequently NDVI difference, which is the 868 pixel window minus the 32632 pixel window, should be the best indicator available, in which a decline trend indicates intensity of human activities and occurrence of the HLD process. In particular, nightlights data from the Defense Meteorological Satellite Program Operational Linescan System (DMSP-OLS) (NOAA's National Geophysical Data Center) sensor was used to separate HLD stations from stations of urban expansion, which leads to a similar NDVI reduction. Finally HLD stations were identified by following the steps:

NDVI-difference method
Step 1. NDVI difference values around the station during growing season was aggregated and averaged for each year. Ordinary Least Squares (OLS) regression technique based on a linear regression model (Y = a+bX+e) was then applied on the dataset.
Step 2. Each 1 km grid in the satellite image of the DMSP-OLS data above the nightlights intensity threshold was classified into urban and non-urban categories. The decision tree to categorize the stations was based on analysis of the grid box group within 32 km of each station [15]. Step 3. Only stations with negative trend and no more than 25% urban grid boxes were classified as HLD.
Step 4. Results were validated through the auxiliary data of Landsat TM satellite image with 30 m spatial resolution from the United States Geological Survey (http://www.usgs.gov/).
No specific permissions were required for these locations/ activities, and the field studies did not involve endangered or protected species.

Identification of growing season
A number of methods, including NDVI threshold, lowest NDVI value from derivative, smoothing algorithms and model fit, have been devised for determining vegetation phenological stage [24,37]. In this study the model fit method, which was developed by Jönsson and Eklundh [37] for analyzing time-series NDVI data, was selected to indentify the growing season. It implements three processing methods based on least-squares fits to the upper envelope of the NDVI data: asymmetric Gaussian (AG), double logistic (DL) and adaptive Savitzky-Golay filtering (SG). The benefit of both asymmetric Gaussian and double logistic approaches is that they are less sensitive to incomplete time-series data with data gaps. In contrast, whilst adaptive Savitzky-Golay filtering is sensitive to data gaps, it can capture subtle and rapid changes in the time-series. Since the GIMMS data have received even better corrections for noise and artifacts [38,39] our method starts with SG to fit the time-series data. To smooth data and suppress disturbances, the SG method uses a filter, and replaces y i with a linear combination of nearby values in a window: where y i are each data value, c j are weights, i = 1, …, n. For each value of y i , the quadratic polynomial f (t) = c1+c2t+c3t2 was fit to all 2n+1 points in the moving window and the value y i was replaced with the value of the polynomial at position t i . The beginning of a season is defined from the SG functions as the point in time when the value has increased by 10% of the distance between the left minimum level and the maximum. The end of the season is defined in a similar way. As such the growing season is defined as the difference between the beginning and the end of a season.

Climatic Zones
The complex topography of the TP makes it difficult to analyze climate change at the scale of the whole region. Therefore, the TP is divided into different climatic zones in order to identify changes more clearly and accurately. The method, based on Lin and Wu [40], uses temperature and precipitation as the basic indicators ( Table 1). The TP is divided into different climatic zones Figure 2. Historical land degradation. Land degradation can be mapped and observed in physical terms, but can only be explained and understood when hidden social, political, and economic structures are analyzed [57]. Socioeconomic data from the 200 years from 1800 to 2000 were collected to explore the reasons for human-induced land degradation. It is clear that growing population and livestock have been putting high pressure on the study area. The flattening curve of population has changed since the 1950s and increased dramatically since the 1980s. Such growth is attributed by People's Republic of China officials to the improved quality of health and lifestyle of the average Tibetan since the beginning of reforms under the Chinese governance. In contrast, population livestock reached its climax in 1980 and has steadily declined since then, resulting in a downward trend in livestock per capita. Loss and degradation of grassland due to increasing human population and livestock pressures are the principal threats to this extremely fragile zone. doi:10.1371/journal.pone.0081535.g002 Figure 3. NDVI difference. Time-series NDVIs in 868 (32 km) and 32632 (128 km) pixel windows of the stations are shown at the top. Heavy lines are the result of smoothing with Savitzky-Golay filtering. Time-series NDVI differences and their trends are shown on the bottom. Zero values of trends from stations in areas with natural forces only (Group A) indicate that there are almost no differences between 32 km areas and 128 km areas, since no additional forces influence the vegetation except natural forces. Negative values of trends from stations in HLD areas (Group B) indicate that NDVI decreases more in the 32 km areas than in the 128 km areas, since pasture degradation from human activities occurred close to these stations. Note that there are two opposite initial states for NDVI difference: 32632 pixel window above 868 pixel window at Madoi station and 868 pixel windows above 32632 pixel window at other stations. doi:10.1371/journal.pone.0081535.g003 , which are centered on each station. Four ranks of degraded land were determined through remote sensing interpretation: slight, moderate, moderate to severe and severe. The major steps included (1) field survey for interpretation and mapping; (2) image preprocessing, including geometric rectification, image registration and atmospheric correction; (3) classification using visual interpretation of spectral features and texture types; and (4) final mapping after validation. It can be clearly seen that land degradation rate is high due to the pressure of human activities at HLD stations. For instance, total degraded area increased 5.6% in the period of 1994-2006 in Nagqu. The part that increased the most was the rank of severe degradation, about 51% of total degraded area. This result confirmed the declining trends of NDVI. doi:10.1371/journal.pone.0081535.g004 according to (1) the number of days per year on which the average temperature is above 10uC, with average temperature in the hottest month as an auxiliary indicator, and (2) an aridity index (Eq. 1), with annual precipitation as an auxiliary indicator:

K~PET=P ð2Þ
where K is the aridity index, PET is potential evapotranspiration and P is annual precipitation. Based on the temperature and precipitation indicators, the TP is divided into 11 climatic zones (Table 2, Fig. 1).

Mann-Kendall test
The Mann-Kendall test is a non-parametric statistical procedure that is well suited to analyzing trends in data over time [41]. One advantage of this test is that the data need not conform to any particular distribution. The second advantage of the test is its low sensitivity to abrupt breaks due to inhomogeneous time series [42]. As a result, the time series of DTR values at seven sites were analyzed for monotonous increasing or decreasing trends with the Mann-Kendall test.
The Mann-Kendall test can be viewed as a nonparametric test for zero slope of the first-order regression of time-ordered data versus time. Each data value is compared with all subsequent data values. If a data value from a later time period is higher than a data value from an earlier time period, the statistic S is incremented by 1. On the other hand, if the data value from a later time period is lower than a data value sampled earlier, S is decremented by 1. The net result of all such increments and decrements yields the final value of S [43]. The MK test used to be applied by considering the statistic S as [44]: where x j is the sequential data values, n is the length of the series and sgn and sgn The standard test statistic Z is computed as follows: The presence of a statistically significant trend is evaluated using the Z value. A positive (negative) value of Z indicates an upward (downward) trend.
The Sen slope [41] calculation is determined along with the Mann-Kendall test. Similar to Mann-Kendall test, the advantage of Sen's slope estimator is that it is insensitive to outliers or missing data. Therefore it is more rigorous than the usual regression slopes and thus provides a realistic measure of the trends in the data series.
Simply speaking, Sen's slope is the median of all differences between successive data values. The magnitude of the trend is computed as [45,46]: where T i is the slope of all data pairs, and x j and x k are data values at times j and k (j.k) correspondingly. The estimate of the slope of the trend is calculated as: A positive value of B means an increasing trend whereas a negative value means a decreasing trend.
The sequential version of the Mann-Kendall test was used to detect abrupt changes in climate against the long-term trend [47]. The null hypothesis H0 presumes that the sample under investigation shows no indication of a developing trend. This rank-based test considers the relations between all terms in the time series (61, 62, 63,…6n). The following steps are applied in order to accept or reject the null hypothesis [48]: Step 1. Definition of the test statistic. The test statistic t j variables are computed as follows: where nj denotes for each element x j (j.k) the number of cases where x j .x k , with j = 1,2,…..n and k = 1,2,…..j21. The distribution of t j is asymptotically normal with E(tj) = [j(j21)]/4 and Var(tj) = [j(j21)(2j+5)]/72.
Step 2. Calculation of the reduced variables. A reduced variable, called statistic u(t),is calculated for each of the test statistic variables t j as follows: Step

Meteorological data
Minimum and maximum temperatures from a total of 7 national meteorological stations for the period from 1 January 1982 to 31 December 2006 were used for this study. The observational data are derived from the ''Daily Surface Climate Variables of China'' catalog issued by the National Meteorological Figure 5. Growing season and related phenological measures for seven meteorological stations. The diagrams depict upper envelope Savitzky-Golay filtered data, with time in month steps. Dots represent original NDVI data. The thick solid lines show fitted functions of SG. The start point is defined as the date of the inflection point of the NDVI curve and end points as the date at which NDVI drops to the same level as the measured right minimum level. In between these dates, vegetation development stages are defined, covering the entire growing season (shaded area). doi:10.1371/journal.pone.0081535.g005 Information Center of the China Meteorological Administration (NMIC/CMA). The station locations span from 30.5uN to 34.9uN, and from 91.1uE to 102.1uE, whilst the altitude of the stations varies from 3471.4 m above sea level (a.s.l.) to 4507 m a.s.l. (Table 3).
Quality control (QC) of observational data is an extremely important factor in evaluating human-induced effects on climate in the study area. The QC must be undertaken prior to any data analysis to eliminate erroneous values. For example corrections must be performed on meteorological observed data if station metadata (coordinates, elevation, etc) change. In recent years researchers have paid more attention to atmospheric data assimilation techniques that could improve QC [49,50].
There were two stages of QC in this study. Firstly, the data were quality controlled by the provider of NMIC/CMA. They were checked by Feng et al. [51] for homogeneity and consistency in five tiers: high-low extreme checking for daily values, internal consistency checking, temporal outlier checking, spatial outlier checking and missing data checking. The data were then screened with the software package C3 extraQC (http://www.c3.urv.cat/ data.html): an expanded version of RClimDex (http://etccdi. pacificclimate.org/software.shtml, running with R 1.84 or a later version) http://cccma.seos.uvic.ca/ETCCDMIin order to (1) identify errors in the minimum and maximum temperature; (2) search for outliers in the minimum and maximum temperature; (3) use the generalized data plot to visually inspect the data and (4) assess data homogeneity.

Results and Analysis
Seven meteorological stations with an average altitude of around 4000 m in the same climatic zone in the central Tibet Plateau were classified into two groups: stations in areas of natural forces only (Group A) and stations in areas of HLD (Group B) ( Fig. 1) (Fig. 2). There is a transition zone between the area that has a high-pressure exerted from human activities (close to Group B) and the area with a low-pressure (close to Group A). Although NDVI differences trends in Group A remained unchanged, those in Group B showed clear declining trends (Fig. 3). The land degradation for stations were also visually interpreted and confirmed through the auxiliary data of Landsat TM satellite image with 30 m spatial resolution in the two periods of 1982-1994 and 1994-2006 (Fig. 4). Wessels [52] proposed the use of RESTREND (residual trend between actual NDVI and predicted NDVI) to distinguish land degradation from the effects of rainfall variability in South Africa. We use the RESTREND technique to validate the NDVI-difference method presented in this study. NDVI data and precipitation during the growing season of 1982-2006 are used in this technique and a comparison of results from NDVI-difference method and RESTREND is shown in Table 4. It can be seen that same results could be concluded for the 4 stations of Damxung, Qumarleb, Madoi and Qingshuihe, but the 3 stations of Nagqu, Maqu and Nangqian were not available using RESTREND. The reason may be that RESTREND is less sensitive to human-induced land degradation in alpine meadow in very high-altitude area of TP than NDVI-difference method, since RESTREND developed by Wessels was originally used in semiarid Karoo, Southern Africa.
A set of analyses were initially performed to compare the effects on trends of maximum temperature and minimum temperature during growing seasons for the period 1982-2006 for the two groups (Fig. 5). It was found that although maximum temperature and minimum temperature both increased, the minimum temperature increased at a significantly faster pace in Group B  Table 5). The confidence in the trend for the Mann-Kendall statistic is calculated using a Kendall probability table [53]. By assessing the S result along with the number of samples, n, the Kendall table provides the probability of rejecting the null hypothesis (H0 = no trend) for a given level of significance. We calculate a 'confidence level' percentage by subtracting the probability (p) from 1 (Confidence = 1-p%). Confidence of 90% represents a significance level of a = 0.1, and 95% confidence corresponds to a = 0.05. The resulting confidence in the trend is applied in the Mann Kendall trend analysis as outlined in Table 3.
The statistical probabilities of P-values were all above 95% except maximum temperature in Nangqian, and hence the null hypothesis (no trend) was rejected. Further analysis showed that there were abrupt changes for minimum temperature at all stations. However the year that abrupt changes happened was much earlier in Group B than in Group A, and the amplitude of this change was larger for Group B. The average abrupt year in Group B was about 5 years earlier: in 1996 vs. in 2001. The average amplitude change was 0.393% for Group B, about 1.4 times that (0.276%) for Group A (Fig. 6). These results reveal a significant and complex effect on the Tibet Plateau climate system from human activities, due to the multiple processes and feedbacks between degraded land surface from anthropogenic influences and the atmosphere. These results reveal differences between HLD and natural stations, which should be probed further in order to fully capture the role of human activities on the plateau climate system.

Conclusions
There are some interesting implications of the differences between HLD stations and the stations in areas with natural forces only. Firstly, annual trends in minimum temperature in the HLD group increase at a much faster rate than in the group from areas with natural forces only, indicating that HLD has a strong impact on the observed temperature. This is consistent with the conclusion that persistent land degradation can notably increase temperatures [54,55]. Secondly, the decreasing rate of annual DTR in the HLD group is particularly remarkable given that it is the reverse of the result from stations in areas with natural forces only. The result reinforces previous work suggesting that the reduction of vegetation cover and soil wetness may reduce the DTR by increasing nighttime surface air temperature [56]. This greater decrease may be partly attributed to (i) environmental conditions (e.g., vegetation) that reduce minimum temperature to a certain extent, thus minimizing the cooling mechanism in land degradation areas; (ii) weak evapotranspiration that prevents surface evaporative cooling effects in land degradation areas at very high-altitudes; and (iii) other factors that influence surface energy balance such as increased soil heating and long-wave surface forcing. These results demonstrate anthropogenic asymmetric changes to the climate, which reflect the complexity of the impact of human-caused land-cover change in the Earth's Third Pole. As such, more attention should be paid to issues such as HLD and more emphasis should be placed on characterizing and quantifying the effect of human-induced land-use change on asymmetric diurnal changes. Figure 6. Comparison of abrupt year and amplitude of minimum temperature, tested by the M-K method. The abrupt year is the year determined to have an abrupt change within the timespan of 1982 to 2006. In HLD stations Nagqu, Maqu, Nangqian and Damxung the abrupt year was 1995, 1998, 1997 and 1993, respectively. It is clear that the abrupt year occurred earlier than at the stations in areas of natural forces only: 2001 for Qumarleb, 1999 for Madoi and 2002 for Qingshuihe. We also compared the amplitude change for each station. HLD stations exhibited a larger amplitude jump between the two periods divided by the abrupt year. For instance Nagqu, Maqu, Nangqian and Damxung changed by 0.44%, 0.41%, 0.32% and 0.4% respectively, compared to 0.3%, 0.29% and 0.24% in Qumarleb, Madoi and Qingshuihe. doi:10.1371/journal.pone.0081535.g006