A New Map of Permafrost Distribution on the Tibetan Plateau

The Tibetan Plateau (TP) has the largest areas of permafrost terrain in the midand low-latitude regions of the world. Some permafrost distribution maps have been compiled but, due to limited data sources, ambiguous criteria, 10 inadequate validation, and deficiency of high-quality spatial datasets, there is high uncertainty in the mapping of the permafrost distribution on the TP. We generated a new permafrost map based on freezing and thawing indices from modified Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperatures (LSTs) and validated this map using various ground-based datasets. The soil thermal properties of five soil types across the TP were estimated according to an empirical equation and soil properties (moisture content and bulk density). The Temperature at the Top 15 of Permafrost (TTOP) model was applied to simulate the permafrost distribution. Permafrost, seasonally frozen ground, and unfrozen ground covered areas of 1.0610 km (0.97-1.1510 km, 90% confidence interval) (40%), 1.4610 km (56%), and 0.0310 km (1%), respectively, excluding glaciers and lakes. Ground-based observations of the permafrost distribution across the five investigated regions (IRs, located in the transition zones of the permafrost and seasonally frozen ground) and three highway transects (across the entire permafrost regions from north to south) were used to 20 validate the model. Validation results showed that the kappa coefficient varied from 0.38 to 0.78 with a mean of 0.57 for the five IRs and 0.62 to 0.74 with a mean of 0.68 within the three transects. Compared with earlier studies, the TTOP modelling results show greater accuracy. The results provide more detailed information of the permafrost distribution and basic data for use in future research on the Tibetan Plateau permafrost.


Introduction
Permafrost is a major component of the cryosphere and is sensitive to climate changes (Wu et al., 2002b;Haeberli and Hohmann, 2008;Li et al., 2008;Gruber, 2012).Due to its unique and extremely high altitude (mean elevation over 4000 m) and low mean annual air temperatures (generally lower than −2 • C with an intra-annual amplitude over 20 • C in the permafrost region), the Tibetan Plateau (TP) possesses the largest areas of permafrost in the mid-and lowlatitude regions of the world (Zhou et al., 2000;Zhao et al., 2004Zhao et al., , 2010;;Yang et al., 2010).Permafrost and its dynamics complicate the water and energy exchange between soil and atmosphere and thereby introduce greater uncertainty into global climate models (GCMs) that predict climate change (Romanovsky et al., 2002;Smith and Riseborough, 2002;Cheng and Wu, 2007;Riseborough et al., 2008;Zhao et al., 2010).To generate better quantitative simulations, a more accurate TP permafrost distribution is needed.An accurate contemporary permafrost distribution map is also important as a baseline with which to estimate future permafrost degradation.
A significant amount of research has been conducted on the permafrost distribution of the TP, and many permafrost maps have been compiled to evaluate the distribution and thermal states of permafrost (Shi and Mi, 1988;Li and Cheng, 1996;Brown et al., 1997Brown et al., , 1998;;Qiu et al., 2000;D. Zou et al.: A new map of permafrost distribution on the Tibetan Plateau Wang et al., 2006).These maps have been used to study the responses of permafrost to climate change (Ran et al., 2012).However, considerable differences in the permafrost areas and boundaries occur in these maps, from 1.12 × 10 6 to 1.50 × 10 6 km 2 , due to different data collection periods, data sets, and methods (Yang et al., 2010;Ran et al., 2012).These maps represent different assessments of the permafrost distribution on the TP at different times.
Permafrost boundaries in the earlier maps have often been plotted on topographic maps by hand using conventional cartographic techniques (Tong and Li, 1983;Shi and Mi, 1988;Li and Cheng, 1996).The standard and most widely used map is the Map of Permafrost on the Qinghai-Tibetan Plateau (Li and Cheng, 1996).In this map, the permafrost boundaries were mainly determined using air temperature isotherms combined with field data, satellite images, and many relevant maps.After 2000, GIS techniques have been applied to permafrost mapping of the TP.Simple empirical models with minimal data requirements were established to consider the permafrost characteristics on the TP.These include the elevation model (Li and Cheng, 1999) and mean annual ground temperature (MAGT) (Nan et al., 2002).Meanwhile, some models with simplified physical processes applicable to high-latitude permafrost were transferred to simulate the permafrost distribution on the TP.Examples are the frost index (Nelson and Outcalt, 1983) and the temperature at the top of permafrost (TTOP) (Smith and Riseborough, 1996;Wu et al., 2002a).These models link permafrost temperature with surface temperature using seasonal surface transfer functions and subsurface thermal properties, which can provide reasonable assessments of permafrost distribution when the permafrost upper boundary conditions and regional soil thermal properties are satisfied.Recently, a global permafrost zonation index (PZI) was established based on the relationships between the air temperature and the occurrence of permafrost (Gruber, 2012).The PZI can represent broad spatial patterns but does not provide the actual extent of the permafrost.
Most temperature fields that have previously been used in these models were generated from spatially interpolated air temperature (Pang et al., 2011) or coarser-resolution (e.g.0.125 • × 0.125 • ) atmospheric reanalysis data (Gruber, 2012;Qin et al., 2015).Although air temperature produces inaccurate and low-resolution estimates of ground surface temperature (GST, defined as the surface or near-surface temperature of the ground (bedrock or surficial deposit) and measured in the uppermost centimetres of the ground), it was still widely used in practical applications because of limited GST observations.In these studies, the N factor has been the most effective way to transform the air temperature to the GST (Lunardini, 1978;Klene et al., 2001).With the recent development of infrared remote sensing technology, an increasing number of land surface temperature (LST, defined as the average temperature of an element of the exact surface of the Earth (e.g.surface of ground, vegetation canopy, or snow cover) calculated from measured radiance; Gillespie, 2014) products derived from different satellite images have been applied to global and regional permafrost distribution research (Kääb, 2008;Hachem et al., 2009;Nguyen et al., 2009;Langer et al., 2010;Westermann et al., 2012Westermann et al., , 2015)).These products are effective alternatives to GST, especially for regions with limited in situ observations, such as the TP.However, the LSTs observed by satellite sensors are instantaneous values at the passing times and must be transformed into mean daily temperature to serve as the thermal state of each day before use.Wang et al. (2011) averaged the twicedaily LST observations of the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors on board the Terra satellite to drive the TTOP model.Their results show a systematic bias with the ground observations because of different observation times (Wang et al., 2011).In addition, the limited availability of soil thermal property spatial data sets creates another problem when modelling permafrost distribution.Most soil surveys have been carried out in seasonally frozen ground or permafrost along the Qinghai-Tibet Highway (Li et al., 2014(Li et al., , 2015a) ) rather than in permafrost regions in the plateau hinterland due to the harsh climate and inconvenient access.Therefore, soil thermal properties have generally been estimated via soil types generated from a limited number of plateau geologic classification studies (Wang et al., 2011;Li et al., 2015b).Overall, there are an insufficient number of field investigations to use for modelling and validating the accuracy of the maps.
Field survey data sets have recently been obtained based on the project "Investigation of Permafrost and Its Environment over The Qinghai-Xizang (Tibet) Plateau" conducted by the Cryosphere Research Station on Qinghai-Xizang Plateau, Chinese Academy of Sciences (CAS).These data could provide validation of permafrost distribution maps.In addition, new information on remote sensing LST applications and spatial soil characteristics on the TP have been studied.An empirical model of daily mean LST was established and performed well in continuous permafrost regions of the central TP (Zou et al., 2014).Li et al. (2015b) studied the relationships between environmental factors and soil types in the permafrost region of the TP, and they used a decision-making tree to spatialize the soil types.The results exhibited good reliability and could be used to realize the spatialization of soil thermal properties.
This study aims to generate a new permafrost distribution map on the TP using remote sensing LST products and investigated soil properties.A multiple linear regression model based on MODIS LST was established, and ground-based LST observations were employed to calibrate the results.Soil thermal conductivities of each soil type on the TP were calculated via in situ observed soil moisture content and bulk density.The TTOP model was used to simulate the permafrost distribution, and the results were validated by the observed permafrost distributions of boreholes, five investigated regions (IRs, located in the transition zones of the permafrost and seasonally frozen ground), and three transects (across the entire permafrost regions from north to south).The TTOP modelling result was also compared with two recent benchmark maps (made in 1996 and 2006).
2 Materials and methods

Field survey data sets
The investigation of permafrost and its environments on the TP was conducted from 2009 to 2014.We studied five IRs -WenQuan (WQ) (Zhang et al., 2011(Zhang et al., , 2012) ) and Budongquan-Qingshuihe (B-Q) in the eastern TP, AErJin (AEJ) in the north-eastern TP, GaiZe (GZ) (Chen et al., 2016) in the southern TP, and XiKunLun (XKL) (Li et al., 2012) in the western TP (Fig. 1), which are located in the transition zones between permafrost and seasonally frozen ground with different climatic and geographic conditions.Groundbased observations, mechanical excavation, geophysical exploration (ground penetrating radar, GPR; time-domain electromagnetic, TEM), and borehole drilling were used, and comprehensive surveys of the permafrost distribution boundary, soil, vegetation, climate, and landform were carried out in all five IRs.The data sets of ground temperature profiles, spatial distribution of vegetation (Wang et al., 2016), and soil types (Li et al., 2014(Li et al., , 2015a) ) were obtained and a long-term permafrost monitoring network was established, including automatic weather station and borehole records.

Boreholes and soil pits
Field survey data sets including ground temperature, soil moisture content, and bulk density were obtained.The ground temperature, measured by temperature probes at different depths (generally set at 0.5 m intervals from 0 to 5 m, 1 m from 5 to 20 m, 2 m from 20 to 40 m, 5 m from 40 to 60 m, and 10 m greater than 60 m) in boreholes were used to determine the existence of permafrost.The soil samples were collected according to depth increments at each pit.The field bulk density (weight of the soil per unit volume) was measured by the clod method.Samples for moisture determination were stored in aluminum sampling boxes and carefully sealed to prevent moisture changes.The soil moisture content was expressed by weight as the ratio of the mass of water present to the dry weight of the soil sample (Wu et al., 2016).The sampling period was concentrated from July to October, so the weighted mean moisture content by depth was used to denote the mean state of each pit.The soil moisture content and bulk density were used to calculate the soil thermal conductivity.There was a total of 125 boreholes and 199 soil pits in five IRs of the field survey (Table 1).

Permafrost maps of five investigated regions
The permafrost maps of the five IRs were used as the validation data in this study.In a local region, the elevation and terrain factors have greater influence on permafrost occurrence than longitude and latitude, especially on mountainous permafrost (Riseborough et al., 2008).The lower limit of permafrost (LLP) was determined based on the linear regression relationship between MAGT and borehole elevation.The elevation where MAGT equals 0 • C was regarded as the LLP (Li et al., 2012;Zhang et al., 2012;Chen et al., 2016).In view of the influence of aspects on LLPs, the boreholes were classified into three types: north-facing, south-facing, and eastwest facing and the LLP of each type was determined.The permafrost distribution was generated based on the LLPs of different aspects and the digital elevation model (DEM) data, and a portion of the observed results of boreholes and geophysical methods (GPR and TEM) was reserved to validate the maps (Li et al., 2012;Zhang et al., 2012).For example, the results of GZ IR showed that the LLP was about 4950 m for north-facing, 5000 m for east-west-facing, and 5100 m for south-facing slopes (Chen et al., 2016).

Permafrost distribution of three highway transects
Three highway transects were established as follows (Fig. 1): National Highway 214 (Qinghai-Yunnan Highway, hereafter G214) from the northern Ela Mountain to Qingshuihe town, National Highway 109 (Qinghai-Xizang Highway, hereafter G109) from Xidatan to Nagqu, and National Highway 219 (Xinjiang-Xizang Highway, hereafter G219) from Kudi to Shiquanhe town.The overall transect lengths of G214, G109, and G219 were approximately 400, 750, and 900 km, respectively.Three transects across the entire permafrost regions from north to south in the eastern, central, and western TP were established.Many permafrost geological conditions were discovered in the process of the construction and renovation of the three highways, and many permafrost roadbed monitoring sections were subsequently set along the highways (Jin et al., 2008;Sheng et al., 2015).Based on these background data and our investigation results (geophysical and drilling exploration), the permafrost distribution limits and geothermal features of the three transects were generated and used as the validation data sets.

Existing two benchmark permafrost maps
One of the most widely used permafrost distribution benchmark maps is the Map of Permafrost on the Qinghai-Tibetan Plateau, compiled by the Lanzhou Institute of Glaciology and Geocryology, Chinese Academy of Sciences (hereafter TP-1996) to support basic research on cryospheric dynamics in China (Li and Cheng, 1996).TP-1996 synthesizes field data, literature, aerial photographs, satellite images, and many relevant maps and shows that the area of permafrost is 1.41 × 10 6 km 2 .The other is the Map of the Glaciers, Frozen Ground and Deserts in China, compiled by Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences (hereafter TP-2006) (Wang et al., 2006).In this map, the permafrost distribution was generated using a 0.5 • C MAGT isotherm as a threshold, which shows that the area of permafrost is 1.12 × 10 6 km 2 .
The MAGT was interpolated based on the relationship between elevation/latitude and the MAGT observation from all 76 boreholes along the Qinghai-Xizang Highway (Nan et al., 2002).

MODIS LST products
The MODIS LST data used in this study were the 1 km gridded clear-sky MOD11A2 (Terra MODIS) and MYD11A2 (Aqua MODIS) products (reprocessing version 5), which span from 2003 to 2012.The results of the radiance-based and temperature-based validation indicated that the accuracy of the global MODIS LST product is better than 1 • C in most cases, including lakes, homogeneous vegetation, and soils under clear-sky conditions (Wan et al., 2002(Wan et al., , 2004;;Coll et al., 2005;Wan and Li, 2008;Wan, 2008).Langer et al. (2010) and Westermann et al. (2011) focused on weekly averages and demonstrated an agreement generally better than 2 • C for MODIS LST in the summer season, at permafrost sites in Siberia and Svalbard, Norway, respectively.Both MOD11A2 and MYD11A2 provide two observations (daytime and night-time), which means that there are four LST observations per day for the same pixel.The temporal resolution of MOD11A2/MYD11A2 was 8 days.The LST values represent the 8-day average LST values (the missing values were ignored in the calculation) (Wan, 2009;Wan and Dozier, 1996), and there are theoretically 46 groups of LST values per year.While the 8-day MODIS LST products have more reliable data than the daily products, they still have numerous missing values when establishing the mean daily LST empirical models due to clouds or other factors (Prince et al., 1998).In this study, the Harmonic Analysis Time-Series (HANTS) algorithm, developed to deal with time series of irregularly spaced observations and to identify and fill in the missing values (Roerink et al., 2000), was applied to smoothen and reconstructing MODIS LST series on a per-pixel basis for the entire study area.The parameters set for the HANTS analysis were described in detail by Xu et al. (2013).
The full coverage of the entire TP requires a total of 13 swathes (h23v04, h23v05, h24v04, h24v05, h24v06, h25v04, h25v05, h25v06, h26v04, h26v05, h26v06, h27v05, and h27v06) of the MOD11A2/MYD11A2 products.The Terra overpass time is around 10:30(local time) in its descending mode and 22:30 in its ascending mode.The Aqua overpass time is around 13:30 in its ascending mode and 01:30 in its descending mode (https://modis.gsfc.nasa.gov/).The MODIS LSTs represent instantaneous observation values, and the overpass times of the satellites do not accurately correspond to standard meteorological observation times (Beijing time: 02:00, 08:00, 14:00, and 20:00) (China Meteorological Administration, 2003).Therefore, an arithmetic mean of the four LST observations with the same weights will produce a large deviation from the mean daily LST (Wang et al., 2011).We used a multiple linear regression to distribute different weights to each MODIS LST observation to establish the mean daily LST empirical model.The details of processing were described by Zou et al. (2014).The model validation at three permafrost sites showed that the determination, mean error, mean absolute error, and root mean squared error of mean daily LST were 0.91 to 0.93, −0.21 to 1, 2.28 to 2.42, and 2.96 to 3.05 • C, respectively.In this study, the empirical formula is as follows: LST daily = 0.18 × Terra day + 0.269 × Terra night + 0.143 × Aqua day + 0.435 × Aqua night + 0.896, (1) where LST daily is the mean daily LST, Terra day is the daytime LST observation of MOD11A2, Terra night is the nighttime LST observation of MOD11A2, Aqua day is the daytime LST observation of MYD11A2, and Aqua night is the nighttime LST observation of MYD11A2.
The calculations of the thawing indices (thawing degree days, TDD) and freezing indices (freezing degree days, FDD) were based on the 8-day average LST calculated from the previous processing.The procedures were realized using the IDL programming language, and the FDD and TDD from 2003 to 2012 were obtained and averaged as the model inputs.

Soil thermal properties
Soil thermal characteristics were modelled according to parameters measured from soil types encountered in the field.The classification of soil types was performed using the Decision Tree See 5.0 software and the Soil Land Inference Model (SoLIM) in conjunction with soil type and environment factor data (Li et al., 2014(Li et al., , 2015b)).According to the USDA soil taxonomy system, there are five soil orders on the TP as follows: Gelisols, Aridisols, Mollisols, Inceptisols, and Entisols.Considering the availability of soil sample parameters, the characteristics of sampling regions, and model applicability, the empirical model of soil thermal conductivity proposed by Kersten (1949) was adopted in this study.The equation of thawed soil thermal conductivity is (2) Furthermore, the equation of frozen soil thermal conductivity is where k t /k f is the thermal conductivity (W m −1 K −1 ) of thawed/frozen soil, ω is the soil moisture content (%), and γ d is the soil bulk density (kg m −3 ).Both ω and γ d were measured via soil samples collected in the field survey.The soil samples were classified according to soil orders; moisture content and bulk density values were averaged within soil orders to eliminate abnormal values (Table 2, the values show the mean with standard deviation of soil thermal parameters of each type).

Glacier and lake data
The spatial distribution and data on areas of glaciers and lakes on the TP were from the Second Glacier Inventory Dataset of China (Guo et al., 2014) and the Chinese Cryosphere Information System (Li, 1998) provided by the Cold and Arid Regions Science Data Center (http://westdc.westgis.ac.cn).

TTOP model
Considering the model's usefulness and sophistication, spatial scales, and available data sets (Riseborough et al., 2008), we selected the temperature at the top of permafrost (TTOP) model (Smith and Riseborough, 1996) to simulate the permafrost distribution on the TP.
The TTOP model can be expressed as follows: where P is the annual period (365 days).TDD (n t × I t ) is the ground surface thawing index, and FDD (n f × I f ) is the ground surface freezing index.n t and n f are n factors of the thawing and freezing seasons, and I t and I f are the air temperature thawing and freezing indices.  is thawing and freezing.We used the modified MODIS LST data processed in Sect.2.2.2 as the GST to derive the TDD and FDD, and the r k was calculated from the soil properties derived from processes in Sect.2.2.3.From Eq. ( 4), if the FDD is greater than , then TTOP will be below 0 • C and permafrost exists.This processing was realized in the ArcGIS software programme with the following expression: D = 1, TTOP ≤ 0 permafrost 0, TTOP > 0 seasonally frozen ground (5) The glacier and lake regions were excluded from the permafrost distribution modelling of the TTOP model.In addition to permafrost and seasonally frozen ground, unfrozen ground was also identified in this study.The unfrozen ground was defined as the region where the extreme minimum LST ≥ 0 • C. The night Aqua MODIS LST (observation time approximately 03:00) was employed as input data for the determination of unfrozen ground area.The uncertainty analysis of total permafrost area was conducted with R statistical software (version 3.3.1,www.r-project.org)using the percentile method, and we used a 90 % confidence interval to determine the range of the total permafrost area.The modelling scheme in this study is shown in Fig. 2.

Accuracy evaluation
The permafrost distribution of the borehole locations, five IRs, and three transects were used to estimate the accuracies of the three maps (TP-1996, TP-2006, and TP-2016).
The spatial distribution of borehole temperature data across a permafrost domain or seasonally frozen ground area was used as the criterion of advantages and disadvantages of results for the three time snapshots of 1996, 2006, and 2016.The permafrost distribution across the five IRs and three tran- The Cryosphere, 11, 2527-2542, 2017 sects were selected as the real values with which to validate the three maps.
To evaluate the agreement of the simulated permafrost distribution and the observed results, the kappa coefficient (K) (Cohen, 1960), which measures the degree of agreement, was selected for accuracy evaluation.
where the total number of pixels is n, and s is the number of pixels in which the simulation and investigated results agree.The number of investigated result pixels with permafrost is a 1 , and those without are a 0 , and the simulated map pixel numbers are b 1 and b 0 .Empirically and statistically arbitrary quality values for K have been proposed.Cohen (1960) suggested that K ≥ 0.8 signifies excellent agreement, 0.6 ≤ K< 0.8 represents substantial agreement, 0.4 ≤ K< 0.6 represents moderate agreement, 0.2 ≤ K< 0.4 represents fair agreement, and a lack of agreement corresponds to K< 0.2.

Permafrost distribution modelling of TTOP
Figure 3 shows the simulated permafrost distribution of the TTOP model on the TP (TP-2016).The distribution areas of permafrost and seasonally frozen ground were 1.06 × 10 6 km 2 with a 90 % confidence interval of 0.97-1.15× 10 6 km 2 , and 1.46 × 10 6 km 2 .This estimate excluded glaciers and lakes, which account for 40 and 56 % of the total TP area, respectively.The result shows that the permafrost distribution was centred in southern Qinghai and northern Tibet.The northern Qiangtang Plateau and the Kunlun Mountains were the regions with the most permafrost, which extends west and north-west to the Karakoram mountains.The permafrost continuity decreases gradually as the elevation decreases and the ground temperature increases with increasing distance from the central region.The geographic northern and southern boundaries of permafrost were Xidatan and Anduo from the mark sites of the Qinghai-Xizang Highway.There were a few areas of permafrost in the high mountains from Anduo to the southern Tibet Valley.Due to the existence of the Bayan Har Mountains and Anemaqen Mountain, the elevations of which are above 5000 m, there is permafrost occurrence in the eastern TP.Some unfrozen ground exists in the south-eastern margin of the TP, and the size of this area is approximately 0.03 × 10 6 km 2 (account for 1 % of the total TP area).

Validation with borehole observations
Boreholes can determine whether permafrost exists or not.
Figure 4 shows the spatial distribution of borehole locations at permafrost or seasonally frozen ground in five IRs of three maps.Different combinations were set up to analyse the differences between the three results.Columns a, b, and c show the results of TP-1996, TP-2006, and TP-2016, and rows 1, 2, 3, 4, and 5 show the results of XKL, GZ, AEJ, B-Q, and WQ IRs, respectively.The results show that TP-1996 is insensitive to the geographical boundaries across all five IRs, and there are many erroneous interpretations of both permafrost and seasonally frozen ground.TP-2006 had higher sensitivity to the boundaries, especially in WQ IR; however the recognition of the other four IRs is inadequate and the areas of permafrost distribution were overestimated.Compared to TP-1996 and TP-2006, TP-2016 performed better at identifying the geographic boundary of permafrost distribution, identifying almost all the boundaries of the five IRs correctly, especially for the seasonally frozen ground in the valley of the north-western XKL IR (Fig. 4c1) and that around the lakes of the eastern AEJ IR (Fig. 4c3).TP-2016 had some errors that were mainly affected by local terrain factors.These included the seasonally frozen ground distributed in valleys and a few permafrost spots at the margin, such as the two seasonally frozen ground boreholes in the northern AEJ IR (Fig. 4c3) and three permafrost boreholes at the south-western limit of GZ IR (Fig. 4c2).

Validation with five investigated regions (IRs)
The permafrost distributions of the five IRs were employed as true values to validate the modelling results of the three maps in order to analyse their performance in geographical boundary recognition ability.TP-1996 was the worst at recognizing the boundaries of permafrost in the five IRs.It misidentified all boundaries, with a low kappa coefficient (K< 0.2), due to greater misjudgment or overestimation of permafrost pixels.TP-2006 also performed poorly in the XKL, GZ, and AEJ IRs (K< 0.2) but performed better in the B-Q and WQ IRs, with a kappa coefficient reaching 0.63 and 0.77.TP-2016 had poor performance in the AEJ IR.The kappa coefficient was only 0.38, which is a slight improvement over estimates of the former two.In addition, it represents moderate agreement with the XKL and GZ IRs and substantial agreement with the B-Q and WQ IRs, which have kappa coefficients of 0.54, 0.48, 0.68, and 0.78, respectively.The average accuracies of TP-1996, TP-2006, and TP-2016 were 0.06, 0.35, and 0.57.TP-2016 performed best in the validation with the investigated permafrost distribution from both the individual and mean accuracies of the five IRs (Table 3).The TP-2016 performed better at identifying the permafrost boundary in the regions with complex terrain because of sharp changes in the LST within short distances, such as the WQ, B-Q, and XKL IRs.For GZ and AEJ IRs, where surface relief is much lower, the TP-2016 does not perform as well as the other three IRs.The worst performance in AEJ IR might also be due to no soil pits in the investigation and the soil thermal properties are inferred completely www.the-cryosphere.net/11/2527/2017/The Cryosphere, 11, 2527-2542, 2017  from the relationship between the environmental factors and the soil samples of the other four IRs.
The results of the AEJ IR and surrounding area were selected to compare the differences among the three maps (Fig. 5).In the AEJ IR, the investigated result shows that the seasonally frozen ground is mainly distributed at the northern valley and the eastern Ayakekumu Lake surrounding areas with permafrost.TP-2006 shows all judgements for permafrost in the AEJ IR, and the permafrost area is clearly overestimated.TP-1996 shows some seasonally frozen ground in the north-western AEJ IR, but the locations were misjudged.TP-2016 judged approximately 30 % seasonally frozen ground in the northern and eastern AEJ IR.Although there were few correct pixels, the locations in the eastern part were at the geographic boundary of permafrost.
The observed MAGT of the borehole closest to Ayakekumu Lake was 3 • C, which indicates the existence of seasonally frozen ground there.TP-2016 accurately modelled this phenomenon.In the regions around the AEJ IR, TP-2016 simulated the seasonally frozen ground around Aqikekule Lake (area approximately 350 km 2 ) and its source river, and this was not found in the other two maps.Most lakes on the TP are formed due to tectogenesis.The major axis basically remains consistent with the main structure directions and the secondary level fracture direction, and there generally exists penetrative or non-penetrative taliks under and around tectonic lakes (Zhou et al., 2000).TP-2016 also shows seasonally frozen ground in the mountainous region in proximity to the Pitileke River, while the other two maps did not identify this area.TP-2016 performed better at identifying the seasonally frozen ground formed by the surface water.
The permafrost distribution of TP-1996 and TP-2006 was modelled according to the relationship between temperature and 3-dimensional zonalities (longitude, latitude, and elevation) (Cheng, 1984).The higher weight given to elevation from the regression equation determined that it has a greater influence than that of longitude and latitude when interpolating temperature (air temperature or MAGT).The high continuity and low variability of the elevation difference in permafrost regions lead to results that appear more continuous.However, the temperature differences caused by local factors (e.g.lakes or rivers) are largely masked and this results in an excessive occurrence of the lower extrapolated temperature.This may explain the overestimated area of permafrost distribution in the previous TP-1996 andTP-2006.The use of remote sensing data can better reveal the spatial heterogeneity of LST.Relative to the two benchmark maps, the TP-2016 result driven by the processed MODIS LST in this paper was very sensitive to seasonally frozen ground formed by surface water, and the results show that there are many seasonally frozen ground areas surrounding lakes and major rivers that correspond to the previous studies (Lin et al., 2011;Niu et al., 2011).

Validation with three transects
The permafrost distribution of the three transects (G214, G109, and G219) of three maps were extracted for comparison with the investigated results to comprehensively evaluate their performance on the mainly permafrost developed regions of the TP.The accuracy statistics of the three maps for the three transects are listed in Table 4. TP-1996 had the worst accuracy of the three maps with an average kappa coefficient of 0.34.The accuracy of TP-2006 was higher than that of TP-1996 with an average kappa coefficient of 0.50.It performed well, especially for transect G109.TP-2016 has the highest accuracy.The kappa coefficients were 0.62, 0.69, and 0.74 for G214, G109, and G219, respectively, with an average of 0.68.TP-2016 performed best in the validation with the investigated permafrost distribution from both the individual and averaged accuracies of the three transects.In the three transects across all permafrost regions from north to south in the eastern, central, and western TP, which include most permafrost distribution characteristics in TP, the validation results should be a synthetic evaluation of the three maps.
Figure 6 shows the distributions of permafrost and seasonally frozen ground along the G109 transect of the three maps and the investigated result.The elevation and mark sites  were also added for analysis.For convenient comparison, the G109 transect was divided into five segments according to the investigated result as follows: two continuous permafrost regions (from XDT to southern FHSYK, and southern YSP to northern AD), one region of seasonally frozen ground only (from southern LDH to NQ) and two regions in which permafrost and seasonally frozen ground coexist (from WL to YSP, and AD to LDH).The comparison shows that the three maps performed well in two continuous permafrost regions.Almost all permafrost is identified correctly except for several seasonally frozen ground areas in the CMEH and BLH of TP-1996.In the region of seasonally frozen ground only, TP-1996 judged permafrost from AD to NQ, which is different from the investigated result and overestimated the permafrost area in this region.TP-2006 and TP-2016 identified that only seasonally frozen ground exists in this region, which is consistent with the investigated result.In two regions where permafrost and seasonally frozen ground coexist, a large difference was seen between the three maps and the investigated result.TP-2006 shows that continuous permafrost exists from XDT to northern AD, performed poorly in the recognition of the seasonally frozen ground, and overestimated the area of permafrost in the G109 transect.TP-1996 performed better than TP-2006 and recognized some of the seasonally frozen ground in TTH, TTH', YSP, and AD.However, the widths and locations reveal bias from that of the investigated result.TP-2016 identified almost all locations of seasonally frozen ground correctly with a smaller width difference, and was more consistent with the investigated result than the former two.Both TP-2006 and TP-2016 identified the sporadic permafrost in LDH, which was generally expected to be the southern limit of permafrost in previous studies.
In the G109 transect, seasonally frozen ground mainly exists due to the surface water effects, regional geologic structure/geothermal effects, and penetration/radiation effects, which cause a discontinuity in the plane and depth of the continuous distribution of permafrost (Zhou et al., 2000).Due to the large streamflow and high water temperature of TTH, TTH' and Buqu (flow through YSP) rivers, the penetrative taliks not only developed on the riverbed and high floodplain, but also expanded to the first or second terrace (the width generally reached 5-10 km).Additionally, a bare ground, gravel layer exists, and a higher mean annual air temperature was beneficial to precipitation infiltration, which created active thermal transfer conditions.Therefore, the seasonally frozen grounds existing in TTH and YSP were also affected by penetration/radiation effects.However, for the rivers with less streamflow and at higher latitudes, such as the CMEH and BHL rivers, the non-penetrative taliks are much smaller (generally < 100 m) and thus almost impossible to identify.The seasonally frozen ground in northern WL was mainly affected by regional geologic structure/geothermal effects, which has been validated by the results of engineer-ing geologic surveys of the Qinghai-Xizang Highway and Railway (Jin et al., 2008).From the kappa coefficients of the three maps and investigated result (Table 4) along the G109 transect, TP-2016 can better identify the seasonally frozen ground that is several kilometres wide and caused by local factors (surface water, geothermal, and permeation/radiation effects).

Spatial difference among the three maps
The kappa coefficients of each pair among the three maps were calculated (Table 5) to analyse the spatial difference.TP-1996 had low consistency with both TP-2006 and TP-2016.The kappa coefficients were 0.56 and 0.53, respectively, which indicates a large difference.TP-2006 had substantial agreement with TP-2016 and the kappa coefficient reached 0.71.The spatial differences between each pair among the three maps were compared (Fig. 7).Compared with TP-2006and TP-2016, TP-1996 overestimated the permafrost area, which was mainly distributed in the southeastern TP, south margin of continuous permafrost, and predominantly continuous and island permafrost in the Southern TP.In addition, TP-1996 misidentified some seasonally frozen ground on the continuous permafrost edge and the interior TP.The permafrost distribution area of TP-2006 was similar to that of TP-2016.Differences mainly existed in the interior TP regions, southern margin of continuous permafrost, and the regions surrounding the Bayan Har Mountains and eastern Nyainqêntanglha Mountains.

Discussion
TTOP was formulated with the modified MODIS LST, rather than ground surface temperature (GST) in this study.It is well known that MODIS LST observes a mixture of the vegetation canopy, snow cover, and ground surface, depending on the region and seasons.The snow cover and vegetation might have a significant influence on the relationship between the GST and MODIS LST, depending on the snow depth and duration (Zhang, 2005), and vegetation height and coverage.The snow cover distribution is spatially variable over the TP (Fig. 8a), with the most persistently snow-covered areas occurring in the south-eastern and western edge of the TP and in some alpine regions with elevations higher than 6000 m (Qin et al., 2006;Pu et al., 2007).Overall, the snow cover is rare, thin (< 3 cm) and has a short duration (mostly existing less than 1 day for a single snow event) due to the strong solar radiation and wind in the vast interior and the northern TP (Che et al., 2008;Huang et al., 2017), where the permafrost is most developed.Therefore, although thin snow cover might have a cooling effect on GST due to the high albedo of fresh snow and a rapid process of snowmelt (Zhang, 2005), the cooling effect may be of short duration and have very little effect on our study.The vegetation types www.the-cryosphere.net/11/2527/2017/The Cryosphere, 11, 2527-2542, 2017  in the alpine ecosystem of the permafrost region on the TP (Fig. 8b) are all composed of grassland and characterized by dwarf and sparsely distributed plants (Wang et al., 2016).The vegetation cover across most of the permafrost region was less than 30 % (Lehnert et al., 2015) and even less than 10 % in the middle and western TP.In view of the conditions of both snow cover and vegetation on the TP, there are only slight differences between the GST and MODIS LST on average, and even smaller differences in FDD and TDD in our study area.In addition, the HANTS algorithm might cause some bias under cloudy-sky conditions.Further evaluation of the algorithm was not performed in this study, because it has been proved to be an effective approach for filling in gaps in the MODIS LST data for the TP where clear-sky conditions dominated (Xu et al., 2013).
The data set used in the earliest maps (compiled in the 1980s and 1990s) included air temperature, field data, aerial photographs, satellite images, and many relevant maps (Tong and Li, 1983;Shi and Mi, 1988;Li and Cheng, 1996).The permafrost boundary was mainly based on a threshold of air temperature isotherms and modified in several regions (such as Qinghai-Xizang Highway, Qinghai-Yunnan Highway, and the Hengduan Mountains) using field data based on the authors' knowledge.The threshold was determined by the empirical statistical relationship between permafrost occurrence and meteorological observations in the eastern TP (Li and Cheng, 1996), while the universality of the threshold is questionable in the western TP due to insufficient data.In addition, high uncertainty exists in the use of air temperature interpolation because of the scarce, unevenly distributed monitoring sites.There were more sites in the eastern TP and fewer in the western TP, more sites at lower elevations and fewer at higher elevations, and very few sites in permafrost regions.This resulted in the low accuracy of extrapolated air temperature on the TP (Lin et al., 2002;Li et al., 2003), especially in the permafrost region.The permafrost maps were compiled with conventional cartographic techniques that plotted the permafrost boundaries on the topographic maps by hand (Tong and Li, 1983;Shi and Mi, 1988;Li and Cheng, 1996).The artefactual errors were very difficult to control and depended on the knowledge and skill of the mapper.These factors led to significant uncertainties in the maps.These maps place much emphasis on the broad concept of "possible" permafrost regions and this overestimated the actual permafrost areas (Wang et al., 2016).The permafrost mapping of TP-2006 was based on the MAGT that considered the characteristics of high-altitude permafrost.The regional MAGT was interpolated based on the relationship between elevation/latitude and the borehole observations along the Qinghai-Xizang Highway (Nan et al., 2002).The MAGT model performed more accurately in the central TP than in the eastern and western TP, which was demonstrated in the validation of the three transects.In view of the medium spatio-temporal resolution and sensitivity to spatial temperature heterogeneity of the MODIS LST data used in the mapping of TP-2016, it can accurately represent the spatial pattern of LST on the TP.In addition, the MODIS LST data were calibrated using ground-based LST observations obtained from automatic weather stations in typical permafrost regions (Zou et al., 2014), which correspond to actual climate conditions of the permafrost region.Moreover, the subsurface thermal properties derived from soil investigation data were also considered in the TTOP model.The improvement in upper boundary conditions of the permafrost model and the use of large quantities of reliable in situ observed data sets led to higher modelling accuracy.
In the earliest maps, only observational data from the field sites along Qinghai-Xizang Highway were used for map evaluation (Tong and Li, 1983;Shi and Mi, 1988;Li and Cheng, 1996).For TP-2006, the threshold of 0.5 • C MAGT was determined by a sensitivity analysis of comparison with the TP-1996, without independent validation (Nan et al., 2002).The validation in this study showed that the TP-2006 accuracy was higher than that of TP-1996. However, TP-2006 highlights the excessive elevation effects in the MAGT interpolation and masks the effects of local factors to some degree.The better performance of TP-2006 in the B-Q and WQ IRs might be explained by geomorphology similar to the Qinghai-Xizang Highway, because these two IRs were closer to the highway than the other three IRs.This suggests that the MAGT model could reflect the permafrost distribution when there are sufficient borehole ground temperature observations, and this is why we used it to model the permafrost distribution of five IRs.The validation results of the five IRs emphasized the performance when recognizing permafrost boundaries and that of the three transects emphasized the overall evaluation of the three maps.Overall, the validation results of both the five IRs and three transects suggested that TP-2016 performed the best and achieved the highest accuracy among the three maps.The results provide a standard permafrost distribution map on the TP based on current climate conditions.
The ground temperatures of permafrost on the TP have increased during the past several decades (Wu and Liu, 2004;Wu and Zhang, 2008;Zhao et al., 2010).This means a dis-equilibrium of permafrost under ongoing climate warming, and thereby any map based on a contemporary climate forcing is likely to underestimate the extent of permafrost.However, permafrost bodies have a long response time to atmospheric conditions (Riseborough, 2007;Romanovsky et al., 2010;Smith et al., 2010).The increasing rates of ground temperature were much lower in the TP than that in the circumpolar regions and much lower for the warm permafrost (Wu and Zhang, 2008;Smith et al., 2005;Zhao et al., 2010), which is mostly distributed near the permafrost boundaries.Moreover, the degradation of permafrost in these regions was characterized by a deepening of the active layer, rather than the disappearance of permafrost.The changes of the permafrost distribution on the TP might be very limited in the past several decades.Therefore, the spatial difference among the three maps might be mainly induced by the differences in methods and data sources.The TP-2016 could be used as the benchmark map for permafrost distribution on the TP, although more work is needed to improve the accuracy of surface forcing and the soil parameters.In addition, although the approach based on the relationship between current climate and permafrost occurrence is useful for mapping the distribution of the TP permafrost, it should be cautious when applying the transient responses of permafrost to climate change to modelling.

Conclusions
This study exploited the advantages of the medium spatiotemporal resolution of MODIS LST products to construct a database of mean daily LST of the TP.The permafrost distribution is simulated by the TTOP model combined with ground observations and soil investigated data sets.The model was validated against the permafrost distribution obtained from the borehole temperature data, five IRs and three transects and compared to two recent benchmark maps.From the validation with borehole temperature data, the suggested method of permafrost boundary identification shows a better result than the two maps, especially for the seasonally frozen ground in valleys and around lakes.The accuracy of the method validation shows that the TP-2016 case has the highest kappa coefficients for the five IRs and three transects.The average coefficients were 0.57 and 0.68.The modelling estimation shows that 1.06 × 10 6 km 2 of permafrost (0.97 × 10 6 -1.15 × 10 6 km 2 , 90 % confidence interval), 1.46 × 10 6 km 2 of seasonally frozen ground, and 0.03 × 10 6 km 2 of unfrozen ground could be on the TP.Compared with two recent benchmark maps, the TTOP model is superior in recognizing the boundary of permafrost, especially in the seasonally frozen ground areas caused by local factors.The new permafrost distribution map represents a basic data set for future permafrost research.

Figure 1 .
Figure 1.Spatial distribution of the field survey regions of the Tibetan Plateau (based on the permafrost distribution map made in 1996).

Figure 2 .
Figure 2. Flow diagram of the modelling scheme.

Figure 3 .
Figure 3. Spatial distribution of permafrost with the derived TTOP on the Tibetan Plateau.

Figure 4 .
Figure 4. Spatial distribution of boreholes in five IRs of three maps.

Figure 5 .
Figure 5.Comparison of the three maps in and around the AErJin investigated region (a TP-1996, b TP-2006, c TP-2016).

Figure 8 .
Figure 8. Annual average snow depth (a edited after Che et al., 2008) and vegetation types of the permafrost region (b edited after Wang et al., 2016) on the Tibetan Plateau.

Table 1 .
Field survey sample statistics in five investigated regions.

Table 2 .
Soil thermal parameters of each type on the Tibetan Plateau.

Table 3 .
Kappa coefficient statistics in five investigated regions of three maps.

Table 4 .
Kappa coefficient statistics in three transects of three maps.

Table 5 .
Kappa coefficient statistics for three maps.