The SRYR catchment area, situated in the upper stream of the Longyangxia Reservoir, covers approximately 145.3× 10³ km² in area. Geographically, it is located between 32° 10′ -36° 59′ N and 95° 54′ -103° 24′ E, and the average elevation is greater than 4000 m (Fig. 1). The study area is characterized by a typical continental plateau climate, with an annual average air temperature of -3.0° C to -4.1° C and an average annual precipitation of 300-700 mm. Three landform types exist in the SRYR: mountains, hills, and plains. Large areas in the eastern and central SRYR are covered with seasonally frozen soil interspersed with isolated patches of sporadically continuous permafrost. The permafrost occurs primarily on high mountains in the western part of the region (Wen et al., 2011).

Fig. 1

	Figure Option ViewDownloadNew Window
	Fig. 1 Topographic map of the source region of the Yellow River (SRYR)

We divided the study area into five regions according to the Qinghai-Tibetan Plateau natural zonation scheme (Zheng et al., 1979) for analyzing the spatial characteristics of the near-surface soil freeze-thaw cycle. These regions are the wide valley basin of the Yellow River source region (I), the Zoige plateau (II), the Golog and Yushu plateau gullies (III), the Huangnan Mountains (IV), and the eastern margin of the Qaidam Mountains (V) (Fig. 2).

We obtained the daily AMSR-E passive microwave remote sensing data from the National Snow and Ice Data Center (NSIDC) at the University of Colorado, USA. These data are provided in the Equal-Area Scalable Earth Grid projection (EASE-Grid) and have a spatial resolution of 25 km for all channels. The data are available twice daily from an ascending track (the satellite passed over the study region at 13:30 local time) and a descending track (the satellite passed over the study region at 01:30 local time) for 10 consecutive years from 19 June, 2002 to 3 October, 2011. The radiometer channel frequencies are 6.9, 10.7, 18.7, 23.8, 36.5 and 89.0 GHz, respectively, and each frequency is associated with dual vertically and horizontally polarized radiation measurements. In this study, the dual-indices algorithm was applied to the 18.5 and 36.5 GHz vertically polarized brightness temperature data (T_{B_18.5V} and T_{B_36.5V}) from both the ascending and descending tracks.

Some of the AMSR-E data were missing because of the unstable states of the radiometer and atmospheric conditions. The missing data included the 211^st-217^th days and the 256^th-262^nd days in 2002, the 304^th-308^th days in 2003, the 324^th day in 2004 and the 332^nd day in 2007. The daily AMSR-E data did not cover the entire study area at the middle and low latitudes; therefore, we had to splice the data from two adjacent days into one to cover the entire SRYR.

Fig. 2

Figure Option ViewDownloadNew Window

Fig. 2 Land cover and natural zonation of the SRYR. The purple solid lines represent the natural zonation boundaries. Natural zones: wide valley basin of the Yellow River source region (I), Zoige plateau (II), Golog and Yushu plateau gullies (III), Huangnan Mountains (IV), and eastern margin of the Qaidam Mountains (V). The red dotted line A represents the north-south profile line at 99.8° E for Figure 7, and the red dotted line B represents the east-west profile line at 34.5° N for Figure 6.

In this study, a land surface temperature (LST) dataset that consisted of data from eleven meteorological stations in the SRYR obtained from the Chinese Meteorological Science Data-sharing Service (Fig. 1) was used to calibrate and validate the thresholds of the dual-indices algorithm. The LST dataset included the daily average land surface temperature (LST_ave), the maximum land surface temperature (LST_max) and the minimum land surface temperature (LST_min) for the period 2002-2011. In addition to the LSTdata, the monthly average air temperature, the maximum air temperature, and the minimum air temperature were also collected for the period 2002-2011. Meteorological data from five meteorological stations, i.e., Qabqa, Guinan, Golog, Zoige and Madoi stations, were used to calibrate the dual-indices algorithm, and data from other stations (Xinghai, Henan, Maqu, Hongyuan, Jigzhi and Darlag stations) were used to validate the algorithm.

In this study, the Advanced Spaceborne Theemal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM), which has a spatial resolution of 30 m, was obtained from the Geospatial Data Cloud. A land cover map for 2009, which has a spatial resolution of 300 m and uses the UN Food and Agriculture Organization’ s land cover classification scheme, was obtained from the Cold and Arid Regions Sciences Data Center at Lanzhou (http://westdc.westgis.ac.cn). The above datasets were pre-processed using a system radiation correction and ground control point correction.

Profile lines can be used to analyze the spatial distribution of onset dates and durations of freeze-thaw cycles at regional scales (Bartsch et al., 2007). We created an east-west profile line at 34.5° N that crossed areas I, II, and III and a north-south profile line at 99.8° E that crossed areas IV and III. These lines intersected at the center of the SRYR and passed through areas with typical and varied landscapes based on the land use and natural zonation in the SRYR (Fig. 2).

We used the trend line method to analyze the changes in onset date and duration of the freeze-thaw cycles estimated using map algebra in ArcGIS. The computational formula is as follows:

${{\theta }_{slope}}=\frac{n\times \sum\limits_{j=1}^{n}{(j\times {{x}_{j}})-(\sum\limits_{j=1}^{n}{j)(\sum\limits_{j=1}^{n}{{{x}_{j}})}}}}{n\times (\sum\limits_{j=1}^{n}{{{j}^{2}})-(\sum\limits_{j=1}^{n}{j{{)}^{2}}}}}, $ (1)

where n is the number of years in the monitoring period (n=10 in this study); x_j, the quantity being studied in year j (onset dates or duration of the freeze-thaw cycle); and θ , the slope of the trend line (when θ > 0, the quantity being studied increases over the n years, and vice versa).

After the trend maps were created, the mean values were computed for each pixel, and the standard deviations were calculated as the average deviation of each pixel from the mean value. These calculations were implemented using the Spatial Analyst Toolbox in ArcGIS.

The dual-indices algorithm uses the 36.5 GHz vertical-polarization brightness temperature (T_{B_36.5V}), the spectral gradient (SG) between the 36.5 GHz (T_{B_36.5V}) and 18.5 GHz (T_{B_18.5V}) brightness temperatures and data from AMSR-E. This algorithm can be expressed as:

${{T}_{B\_36.5V}}< {{T}_{cutoff}}, $ (2)

$SG=\frac{{{T}_{B\_36.5V}}-{{T}_{B\_18.5V}}}{36.5-18.5}< 0, $ (3)

where T_cutoff is the threshold determined from the linear regressions between LST and T_{B_36.5V} for determining the near-surface soil freeze-thaw status.

LSTs tend to be close to the maximum values at the ascending local time (13:30) for AMSR-E and approach to the minimum values at the descending local time (01:30) (Xie et al., 2013; Chai et al., 2014). Therefore, we used LST_max and LST_min to calibrate and verify the classification algorithm for the ascending and descending tracks, respectively. Two linear regressions between the LSTs and T_{B_36.5V} were established for the period from 1 January, 2007 to 31 December, 2009 (Fig. 3), and robust correlations (R²=0.772 and R²=0.767, respectively) were identified between these quantities and the standard deviations were 3.87 and 3.32 K, respectively. According to the regression lines, the thresholds of the dual-indices algorithm were 257.59 K during the day and 261.28 K at night, which corresponded to the real soil freezing point (273.15 K). A pixel of the AMSR-E data was assigned to the thawed status during the day if the ascending T_{B_36.5V} and SGexceeded 257.59 and 0 K, respectively; otherwise, the pixel was assumed to be frozen. A pixel was classified into the frozen status at night if the descending T_{B_36.5V} and SG were less than 261.28 and 0 K, respectively; otherwise, the pixel was assumed to be thawed. Additionally, large bodies of water (lakes and rivers), areas of permanent snow and ice and deserts influence brightness temperatures and could be mistakenly classified. Therefore, areas up to 50 km with those land cover types (e.g., lake, river, snow, ice, and desert) were masked out in this study.

The following variables were used to evaluate the accuracy of the algorithm: (a) F_√ represents the number of correctly classified daily total frozen pixels; (b) F_a, the actual number of daily total frozen pixels; (c) T_√, the number of correctly classified daily total thawed pixels; (d) T_a, the actual number of daily total thawed pixels; (e) FTC_√, the number of correctly classified daily freeze-thaw cycle pixels; and (f) FTC_a, the actual total number of daily freeze-thaw cycle pixels.

The calculation method for the evaluation variables was as follows: (1) F_√/F_a ratio can be used for detecting the classification accuracy of the daily total frozen surface; (2) T_√/T_a ratio, the classification accuracy of the daily total thawed surface; (3) FTC_√/FTC_a ratio, the classification accuracy of the daily freeze-thaw cycle surface; and (4) (F_√+T_√+FTC_√)/(F_a+T_a+FTC_a) ratio, the overall classification accuracy.

Fig. 3

Figure Option ViewDownloadNew Window

Fig. 3 Relationship between land surface temperature (LST) and brightness temperature at 36.5 GHz (TB-36.5V) with vertical polarization based on Advanced Microwave Scanning Radiometer for the Earth Observing System (AMSR-E) data for daytime (a) and nighttime (b) passes. LSTmax, the maximum land surface temperature; LSTmin, the minimum land surface temperature; RMSE, root mean square error; F, degree of freedom.

We verified the algorithm with LST_max and LST_min measurements from six meteorological stations (Xinghai, Henan, Maqu, Hongyuan, Jigzhi and Darlag) during the period 2007-2009. Additionally, the same thresholds were used in the calibration algorithm. Based on 13, 152 validated data points, the classification accuracies of the daily total frozen surface, the daily total thawed surface, and the daily freeze-thaw cycle surface were 80.46%, 83.65%, and 78.07%, respectively, and the overall classification accuracy was 80.55%. These results indicate that the dual-indices algorithm can be used to detect the near-surface soil freeze-thaw cycle in the SRYR with sufficient accuracy.

Based on the aforementioned dual-indices algorithm and thresholds, we detected the daily near-surface freeze-thaw status derived from the AMSR-E data during the day and at night and then classified the status into three types. These types corresponded to the daily frozen surface (DFS) where the land surface was classified into the frozen status during both day and night; the daily thawed surface (DTS) where the land surface was classified into the thawed status during both day and night; and the daily surface freeze-thaw cycle (DFTS) where the land surface was classified into thawed during the day and frozen at night. Based on the time series information of the daily near-surface freeze-thaw status, four soil freeze-thaw stages in a year were defined: a spring freeze-thaw transition stage, a complete thawing stage, an autumn freeze-thaw transition stage and a complete freezing stage. If the surface belonged to one of three daily freeze-thaw statuses for three successive days, the periods of the surface status were determined as one of the four stages, and the first day of these three days was recorded as the onset date of the next freeze-thaw stage.

Figure 4 shows the monthly variations in near-surface soil freeze-thaw status from September 2010 to August 2011. As the temperature decreased after September, a DFS occurred near Ngroing Lake and Gyaring Lake in area I, while a DFTS occurred near Maduo County and Requ River. The DFS and DFTS subsequently expanded, moved southeastward, and gradually occupied area III. In January, which is the coldest month of the year, the DFS and DFTS were bounded by the A’ nyê maqê n Mountains, whereas the DFTS and DTS were bounded by the Xiqing Mountains. At that time, area I, which lies to the west of the SRYR, displayed a DFS status; areas III, IV and V, which lie in the middle of the SRYR, displayed a DFTS status; and area II, which lies to the southeast of the SRYR, displayed a DTS status. As the air temperature gradually increased, the western boundary of the DFS, the eastern boundary of the DTS, and the eastern and western boundaries of the DFTS moved westward. Later in March, the DFS disappeared, and a DFTS replaced the DFS in area I to the west of the A’ nyê maqê n Mountains. The DTS gradually increased and then occupied the entire SRYR in July and August. Therefore, the SRYR, except for the area with DTS throughout the year, underwent freeze-thaw cycles over long periods of time.

Fig. 4

	Figure Option ViewDownloadNew Window
	Fig. 4 Monthly variations in near-surface soil freeze-thaw status from September 2010 to August 2011 in the SRYR

The areas occupied by the three freeze-thaw statuses showed different trends from 2003 to 2010 (Fig. 5). The DFS gradually decreased from January to June and began to increase in September. Additionally, the peak of the curve for the month with the largest DFS moved backward from 2003 to 2010 (Fig. 5a), and the average DFS decreased by approximately 0.08%. In spring and autumn, the surface temperature dropped below zero at night, whereas the surface temperature exceeded zero during the day, resulting in freezing at night and thawing during the day. Being similar with the DFS, the DFTS occurred in autumn or winter, although its peak occurred during the early spring (Fig. 5b); and the surface gradually decreased and completely disappeared in summer. The DTS increased from January to July or August and then gradually decreased (Fig. 5c).

Fig. 5

	Figure Option ViewDownloadNew Window
	Fig. 5 Monthly variations in the daily frozen surface (DFS, a), daily freeze-thaw cycle surface (DFTS, b) and daily thawed surface (DTS, c) during the period 2003-2010

We analyzed the spatial distribution and inter-annual variation of the onset date of the four stages and the duration of the freeze-thaw cycles using the east-west profile line at 34.5° N (Fig. 6) and the north-south profile line at 99.8° E (Fig. 7). For the east-west profile line at 34.5° N, the dates of the onset of thawing (10-65^th days) (Fig. 6a) and complete thawing (28-115^th days) (Fig. 6b) moved forward from east to west as the elevation decreased, and the dates of the onset of freezing (260-320^th days) (Fig. 6c) and complete freezing (310-360^th days) (Fig. 6d) moved backward. The duration of the spring and autumn freeze-thaw transitions, which were determined by the onset of thawing/freezing and the onset of complete thawing/freezing, decreased from east to west. The duration was shorter in area II (spring: 45 days; autumn: 28 days) than in area III (spring: 52 days; autumn: 31 days); moreover, this duration was shorter than in area I (spring: 57 days; autumn: 35 days). The dates of the onset of thawing and complete thawing tended to occur earlier over the study period (2002-2011), and the dates of the onset of freezing and complete freezing tended to occur later over the study period. It should be noted that part of area II did not completely freeze in winter after 2006.

Fig. 6

	Figure Option ViewDownloadNew Window
	Fig. 6 Changes in the onset date and duration of the freeze-thaw cycles for the east-west profile line at 34.5° N during the period 2002-2011. Onset date of thawing (a), date of complete thawing (b), onset date of freezing (c), date of complete freezing (d), and elevations in areas I, II, and III (e).

For the north-south profile line at 99.8° E, the dates of the onset of thawing (10^th-55^th days) (Fig. 7a) and complete thawing (71^st-104^th days) (Fig. 7b) moved forward first and then backward from north to south. Similarly, the dates of the onset of freezing (280^th-318^th days) (Fig. 7c) and complete freezing (310^th-345^th days) (Fig. 7d) moved backward first and then forward. Turning point is at the transition zone between areas IV and III in the rolling mountain region of A’ nyê maqê n. The durations of the spring and autumn freeze-thaw transitions, which were disrupted by the elevated mountains, first increased and then decreased along this line from north to south. In area IV in the northern SRYR, the spring duration (approximately 52 days) and autumn duration (approximately 38 days) of freeze-thaw transitions increased from north to south, whereas an opposite trend was observed in area III. Similar to the east-west profile, the date of the onset of the freeze-thaw cycles along the north-south line showed a trend toward earlier thawing and later freezing from 2002 to 2011 against the background of climate warming in the SRYR.

Fig. 7

	Figure Option ViewDownloadNew Window
	Fig. 7 Changes in the onset time and duration of the freeze-thaw cycles for the north-south profile line at 99.8° E during the period 2002-2011. Onset date of thawing (a), date of complete thawing (b), onset date of freezing (c), date of complete freezing (d), and elevations in areas III and IV (e).

Figure 8 shows the spatial distribution of the onset date and duration of the freeze-thaw cycles. The dates of onset of thawing advanced or moved forward by 3.1 (± 1.4) days over the SRYR (Fig. 8a), and most of these advanced days appeared in area IV, followed by area III.The complete thawing date (Fig. 8c) was advanced by 2.46 (± 1.4) days and most of these advanced days appeared in area II, with a decreasing tendency in the number of advanced days from southeast to northwest. The dates of onset of freezing (Fig. 8e) and complete freezing (Fig. 8 g) were delayed by 0.9 (± 1.4) and 1.6 (± 1.1) days, respectively. Most of these delayed days appeared in parts of areas II and IV in the eastern SRYR, whereas no significant changes occurred in area I.

Fig. 8

Figure Option ViewDownloadNew Window

Fig. 8 Spatial distribution of trends in the onset date and duration of the freeze-thaw cycles (in days). Onset date of thawing (a), duration of the near-surface soil freeze-thaw transition during spring (b), date of complete thawing (c), duration of the near-surface soil thawing (d), Onset date of freezing (e), duration of the near-surface soil freeze-thaw transition during autumn (f), date of complete freezing (g), and duration of the near-surface soil freezing (h) in the SRYR during the period 2002-2011.

The duration of the spring freeze-thaw transition (Fig. 8b) in the middle of the SRYR increased by 2.1 (± 0.91) days, whereas that in western Ngoring Lake, Zoige Wetland and the northernmost SRYR the duration decreased by 1.21 (± 0.56) days. Similarly, the duration of the autumn freeze-thaw transition (Fig. 8f) in the central region increased by 1.1 (± 0.42) days, and the transition decreased by 1.21 (± 0.56) days in the other parts. Overall, the duration of complete thawing (Fig. 8d) increased by 0.72 (± 0.21) days/a with the maximum increase occurred in area IV and the minimum increase in western Ngoring Lake and northernmost area. Overall, the duration of complete freezing (Fig. 8h) decreased at a rate of 0.52 (± 0.26) days/a with the greatest decrease occurred in the northeastern SRYR at lower altitudes and the smallest increase in the southwest SRYR at higher altitude. The above results suggest that the changes in the freeze-thaw cycles were most severe in areas III and IV in the northeastern and southeastern parts of the SRYR in relatively low-lying regions.