Spatial and Temporal Patterns of Rainfall Erosivity in the Tibetan Plateau

: The Tibetan Plateau is inﬂuenced by global climate change which results in frequent melting of glaciers and snow, and in heavy rainfalls. These conditions may increase the risk of soil erosion, but prediction is not feasible due to scarcity of rainfall data in the high altitudes of the region. In this study, daily precipitation data from 1 January 1981 to 31 December 2015 were selected for 38 meteorological stations in the Tibetan Plateau, and annual and seasonal rainfall erosivity were calculated for each station. Additionally, we used the Mann–Kendall trend test, Sen’s slope, trend coe ﬃ cient, and climate tendency rate indicators to detect the temporal variation trend of rainfall erosivity. The results showed that the spatial distribution of rainfall erosivity in the Tibetan Plateau exhibited a signiﬁcant decreasing trend from southeast to northwest. The average annual rainfall erosivity is 714 MJ · mm · ha − 1 · h − 1 , and varies from 61 to 1776 MJ · mm · ha − 1 · h − 1 . Rainfall erosivity was mainly concentrated in summer and autumn, accounting for 67.5% and 18.5%, respectively. In addition, annual, spring, and summer rainfall erosivity were increasing, with spring rainfall erosivity highly signiﬁcant. Temporal and spatial patterns of rainfall erosivity indicated that the risk of soil erosion was relatively high in the Hengduan mountains in the eastern Tibetan Plateau, as well as in the Yarlung Zangbo River Valley and its vicinity.


Introduction
Soil erosion is a global environmental problem, which leads to land degradation, siltation of reservoirs, and eutrophication of water bodies, among others [1][2][3]. Formation mechanisms and succession processes of soil erosion are affected by rainfall erosivity, which is highly correlated with the product (EI) of the total storm energy (E) and the maximum 30 min intensity (I 30 ), both derived from data by Wischmeier and Smith (1958) [4] and Wischmeier (1959) [5]. The concept of rainfall erosivity developed further by Hudson (1971) [6] and Wischmeier and Smith (1978) [7] describes erosivity as the average of annual summations of storm EI 30 . Rainfall erosivity is the basic factor in the Universal Soil Loss Equation [7] and the Revised Universal Soil Loss Equation [8]. Additionally, many empirical soil erosion predictions also use rainfall erosivity [9][10][11][12]. Although the EI 30 index has been accepted worldwide, its calculation requires hyetograph data for a storm. Therefore, the use of the EI 30 index has been limited by a lack of high temporal resolution rainfall data. Many statistical models have been established based on annual rainfall [13][14][15][16][17], monthly rainfall [18][19][20], daily rainfall [21][22][23][24][25], and even hourly rainfall [26,27]. Hourly rainfall data cannot be obtained from many of the national meteorological observation stations, which limits their application. By contrast, daily rainfall data can guarantee the accuracy of rainfall erosivity calculations due to the level of detail on rainfall variability [21]. Daily rainfall data is widely used in the calculation of rainfall erosivity [21][22][23][24][25].
The spatial heterogeneity of soil erosion is very high and has multiple scale characteristics. It depends on the diversity and complexity of factors affecting soil erosion. Rainfall erosivity is one of the basic factors responsible for soil erosion [28]. However, it is not distributed uniformly throughout the year, and at regional scale, knowledge of seasonal and even semi-monthly distribution of rainfall erosivity is critical to the accuracy of soil erosion calculations [8,11]. Therefore, the spatial and temporal distribution of rainfall erosivity concerns various countries and regions [12,15,[29][30][31][32]. Spatial variability in rainfall erosivity in China is relatively large, with rainfall erosivity decreasing from the southeast to the northwest [33,34]. Therefore, compared with the northwest region, the risk of erosion in the southeastern region was significantly greater, and the area is also receiving the attention of scientists. The global climate system is undergoing a change characterized by significant warming [35]. The Tibetan Plateau is sensitive to climate change and ecologically fragile [36]. Climate in the Tibetan Plateau has changed dramatically in the early 21st century, with the continuing severe warming and increasing precipitation [37]; changes include increasing snowmelt and more frequent heavy precipitation events. Vegetation changes in the Tibetan Plateau and its response to climate change [38], precipitation and climate change [39], runoff and soil erosion [40][41][42], as well as watershed non-point source pollution in the main farming area [43], have attracted the attention of researchers. Protecting the grasslands of the Tibetan Plateau is of great importance in limiting global climate change [44]. However, grassland degradation phenomena have been found due to the unreasonable use of grassland by human beings and changing global climate [45]. Besides, aeolian desertification is impeding sustainable socio-economic development [46]. Furthermore, soil erosion in this area is relatively weak compared to the eastern part of China, and recovery from soil erosion is lengthy and difficult once it occurs.
The service span of water conservancy facilities and flood control capabilities were reduced by severe soil erosion that led to the siltation of channel beds and thus a reduction in reservoir capacity in the Tibetan Plateau [47]. The Tibetan Plateau, especially the gorge area in the eastern part, has a large elevation difference where severe soil erosion restricts the development of agriculture and animal husbandry [42]. Sparse distribution of meteorological stations in this region prevented data collection on rain characteristics; however, surface runoff generated by snowmelt can also cause soil erosion. Although researchers have done some studies on rainfall erosivity in the Tibetan Plateau, it is far from enough [48][49][50][51]. Dynamic monitoring of soil erosion in China still requires detailed data on rainfall erosivity in the area.
To address this, daily precipitation data were used from 1 January 1981 to 31 December 2015 for a total of 38 stations in Tibet, and the annual and seasonal rainfall erosivity was calculated. The objectives of the study were to (a) identify long-term trends in rainfall erosivity, and (b) map the spatiotemporal patterns in rainfall erosivity. The results of this paper are intended to optimize the quantitative prediction of soil erosion and soil and water conservation planning services in the Tibetan Plateau.

Study Area
The Tibetan Plateau is part of and located in the southwestern part of the Qinghai-Tibet Plateau. It has a total area of about 1.2 million km 2 and an average elevation of 4000 m, which gave it its designation as the "roof of the world". Most of the Tibetan Plateau is arid and semi-arid. Precipitation has an uneven distribution and large spatial differentiation. The southeast area has annual precipitation of 600-800 mm, while the western area suffers from drought with annual precipitation of less than 200 mm. Precipitation from May through October accounts for more than 80% of the annual total. The spatial distribution of precipitation indicates a decreasing trend from the southeast to the northwest [52]. Mean annual temperature varies between −3.0 and 11.8 • C, and daily temperature fluctuations are often as much as about 15 • C [53]. The region experiences soil erosion, desertification, and mountain hazards. Grasslands make up the largest proportion of land cover, accounting for 56.72% of the total land area in the region. The proportion of cultivated land is 0.42% and it is mainly found in southern Tibet, with a small amount distributed in the east and southeast. Due to the high altitude in the northwestern region, the snow cover and frozen soil make this area unsuitable for human habitation, and the level of human activities is thus low. Thirty-eight meteorological stations used in this study are primarily found in the eastern and southern parts of the Tibetan Plateau, with very few in the northwest (Figure 1). has an uneven distribution and large spatial differentiation. The southeast area has annual precipitation of 600-800 mm, while the western area suffers from drought with annual precipitation of less than 200 mm. Precipitation from May through October accounts for more than 80% of the annual total. The spatial distribution of precipitation indicates a decreasing trend from the southeast to the northwest [52]. Mean annual temperature varies between −3.0 and 11.8 °C, and daily temperature fluctuations are often as much as about 15 °C [53]. The region experiences soil erosion, desertification, and mountain hazards. Grasslands make up the largest proportion of land cover, accounting for 56.72% of the total land area in the region. The proportion of cultivated land is 0.42% and it is mainly found in southern Tibet, with a small amount distributed in the east and southeast. Due to the high altitude in the northwestern region, the snow cover and frozen soil make this area unsuitable for human habitation, and the level of human activities is thus low. Thirty-eight meteorological stations used in this study are primarily found in the eastern and southern parts of the Tibetan Plateau, with very few in the northwest ( Figure 1).

Collection of Rainfall Data
Daily precipitation data from 1 January 1981 to 31 December 2015 for a total of 38 stations from the National Meteorological Information Center of the China Meteorological Administration were used in the present study [54]. The altitude of these stations varies from 2327 to 4800 m, with the lowest elevation at Zayu station and the highest at Amdo station. The northernmost and westernmost station was Shiquanhe, the southernmost was Pagri, and the easternmost was Markam station ( Figure  1).

Calculation of Rainfall Erosivity
Rainfall data were checked for quality. The cold and warm season rainfall estimation model was selected to calculate rainfall erosivity [34,55]. The model takes into full account the characteristics of rainfall in China; that is, the warm season (May-September) with multi-convective, short-duration heavy rain, and the cold season (October-December, January-April) with many frontal, longduration and heavy rains. The model has been widely used, especially in the first national water census in China [11]. The rainfall erosivity model for estimating the daily rainfall during the warm and cold seasons has the following form:

Collection of Rainfall Data
Daily precipitation data from 1 January 1981 to 31 December 2015 for a total of 38 stations from the National Meteorological Information Center of the China Meteorological Administration were used in the present study [54]. The altitude of these stations varies from 2327 to 4800 m, with the lowest elevation at Zayu station and the highest at Amdo station. The northernmost and westernmost station was Shiquanhe, the southernmost was Pagri, and the easternmost was Markam station (Figure 1).

Calculation of Rainfall Erosivity
Rainfall data were checked for quality. The cold and warm season rainfall estimation model was selected to calculate rainfall erosivity [34,55]. The model takes into full account the characteristics of rainfall in China; that is, the warm season (May-September) with multi-convective, short-duration heavy rain, and the cold season (October-December, January-April) with many frontal, long-duration and heavy rains. The model has been widely used, especially in the first national water census in China [11]. The rainfall erosivity model for estimating the daily rainfall during the warm and cold seasons has the following form: where R is the average annual rainfall erosivity in MJ·mm·ha −1 ·h −1 ; k = 1, 2, . . . , 24, which divide the year into 24 half months; R hal f month k is rain erosivity of the k th half-month in MJ·mm·ha −1 ·h −1 ; I = 1, 2, . . . , n, and n indicates the rainfall data year; j = 0,1, . . . , m, m is the number of the erosion-causing rainy days in the k th half month of the i th year; p i,j,k is the j amount of erosive daily rainfall (rainfall ≥12 mm) in the k th month of the i th year in mm, with j = 0 indicating no erosive rainfall resulting in p i, j,k = 0; in the dry season α is 0.3101, and in the rainy season α is 0.3937; and WR hal f month k is the ratio of the average rainfall erosivity of the k th month to the average annual rainfall erosivity.

Trend Coefficient and Climate Tendency Rate
To determine temporal variability in rainfall erosivity, we applied the trend coefficient (TC) and climate tendency rate (CTR) indicators, which are widely used to detect and assess, respectively, the direction and extent of long-term change in climate factors [56,57]. The formula for the trend coefficient is as follows: where r xt is the trend coefficient and n is the number of years. x i is rainfall erosivity of the i-th year.
x is the average annual rainfall erosivity for the n years. When the trend coefficient r xt is positive, it indicates a linearly increasing trend of rainfall erosivity, and a negative value indicates a linearly decreasing trend. A linear equation is usually used to indicate trends in meteorological elements [58,59]. The climate tendency rate was calculated as follows: where a 1 × 10 is the climate tendency rate over a 10-year period (MJ·mm·ha −1 ·h −1 ·10a −1 ), indicating the variation of rainfall erosivity per 10 years. According to the linear regression theory, where σ x is the mean square error of the rainfall erosivity and σ t is the mean square error of sequence 1, 2, . . . , n.

Mann-Kendall Trend Test
The Mann-Kendall test [60,61] was used to detect the significance level and the abrupt-change point of the long-term variability in R factors; the specific calculation method was as follows: Construct a matrix for the time series x which has a sample of n: Define the statistic UF k given the assumption that the time series is random and independent: where UF 1 = 0, and E(s k ) and Var(s k ) are the average and variance of the cumulative number s k , respectively. When x1, x2, . . . , xn are mutually independent and of the same continuous distribution, they can be calculated from the equations below: Var(s k ) = n(n − 1)(2n + 5) 72 .
UF i is the standard normal distribution, which is calculated according to the time sequence x1, x2, . . . , xn. Given the significance level α, when|UF i | > Uα, there is a significant trend in the sequence.
Repeat the calculation to the inverted sequence xn, xn−1, . . . , x1, and make UB k = −UF k (k = n, n − 1, . . . 1), UB = 0. Intersection of UF and UB is the abrupt change point. Sen (1968) [62] developed the non-parametric procedure for estimating the slope of the trend in the sample of N pairs of data:

Sen's Slope Estimator
where x j and x k are the data values at times j and k (j > k), respectively.
If there is only one datum in each time period, then N = n(n − 1)/2, where n is the number of time periods. If there are multiple observations in one or more time periods, then N < n(n − 1)/2, where n is the total number of observations. The N values of Q i are ranked from smallest to largest and the median of slope or Sen's slope estimator is computed as shown: The sign of Q med reflects the direction of data trend, while its value indicates the steepness of the trend.
The confidence interval about the time slope [62,63] can be computed as follows: where Var (S) is defined in Equation (13) and Z 1−α/2 is obtained from the standard normal distribution table. In this study, the confidence interval was computed at two significance levels (α = 0.01 and α = 0.05). Then, M 1 = (N − C α )/2 and M 2 = (N + C α )/2 are calculated. The lower and upper limits of the confidence interval, Q min and Q max , are the M 1 largest and the (M 2+1 ) th largest of the N ordered slope estimates [63].
The slope Q med is statistically different than zero if the two limits (Q min and Q max ) have the same sign.

Temporal Variability in Rainfall Erosivity
Sen's slope estimator analysis indicated that rainfall erosivity in the Tibetan Plateau exhibited an increasing trend from 1980 to 2015, but the trend was not significant. Sen's slope of annual rainfall erosivity was 3.05 (p = 0.32). However, variability at each station was high, and 13 stations exhibited a decreasing trend (34% of the stations), and 25 stations an increasing trend (66% of the stations) ( Figure 2a, Table 1). The climate tendency rate had a significantly increasing trend at the Gerze station in the northern part of the Tibetan Plateau, Shenza station in the middle, and Markam station in the eastern part of the Minjiang River basin (at the 0.05 significance level), and the rate increased by 20 MJ·mm·ha −1 ·h −1 ·10a −1 , 32 MJ·mm·ha −1 ·h −1 ·10a −1 , and 50 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively, at the three stations. Increasing rainfall erosivity was mainly found in the southern Tibetan valley, Yarlung Zangbo River basin, and Hengduan Mountains in the Tibetan Plateau, especially in the Lancang River Basin. An annual rainfall erosivity anomaly showed that it fluctuated. The rainfall erosivity was small before 1987, and the value increased after 2010. In the study period, two-thirds of the annual rainfall erosivity was higher than the average; the minimum appeared in 1992, and the maximum appeared in 1998 ( Figure 2b). The results of the M-K test showed a non-significant increasing trend in rainfall erosivity in the Tibetan Plateau since 1984 ( Figure 2c). The confidence interval about the time slope [62,63] can be computed as follows: where Var (S) is defined in Equation (13) and Z1−α/2 is obtained from the standard normal distribution table.
In this study, the confidence interval was computed at two significance levels (α = 0.01 and α = 0.05). Then, are calculated. The lower and upper limits of the confidence interval, Qmin and Qmax, are the M1 largest and the (M2+1)th largest of the N ordered slope estimates [63]. The slope Qmed is statistically different than zero if the two limits (Qmin and Qmax) have the same sign. Both the Mann-Kendall statistical test and Sen's slope have been frequently used to detect the significance of trends in hydro-meteorological time series [64][65][66][67][68][69].

Temporal Variability in Rainfall Erosivity
Sen's slope estimator analysis indicated that rainfall erosivity in the Tibetan Plateau exhibited an increasing trend from 1980 to 2015, but the trend was not significant. Sen's slope of annual rainfall erosivity was 3.05 (p = 0.32). However, variability at each station was high, and 13 stations exhibited a decreasing trend (34% of the stations), and 25 stations an increasing trend (66% of the stations) ( Figure 2a, Table 1). The climate tendency rate had a significantly increasing trend at the Gerze station in the northern part of the Tibetan Plateau, Shenza station in the middle, and Markam station in the eastern part of the Minjiang River basin (at the 0.05 significance level), and the rate increased by 20 MJ·mm·ha −1 ·h −1 ·10a −1 , 32 MJ·mm·ha −1 ·h −1 ·10a −1 , and 50 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively, at the three stations. Increasing rainfall erosivity was mainly found in the southern Tibetan valley, Yarlung Zangbo River basin, and Hengduan Mountains in the Tibetan Plateau, especially in the Lancang River Basin. An annual rainfall erosivity anomaly showed that it fluctuated. The rainfall erosivity was small before 1987, and the value increased after 2010. In the study period, two-thirds of the annual rainfall erosivity was higher than the average; the minimum appeared in 1992, and the maximum appeared in 1998 ( Figure 2b). The results of the M-K test showed a non-significant increasing trend in rainfall erosivity in the Tibetan Plateau since 1984 ( Figure 2c).   The distribution of rainfall erosivity over the 24 and a half months is an input factor for many soil erosion model calculations, including the Chinese Soil Loss Equation (CSLE) [70] and for calculating vegetation cover and biological measure factors. Data from 1981 to 2010 indicated that rainfall erosivity in the Tibetan Plateau was mainly concentrated in June-September, accounting for 81% of the year ( Figure 3). This was closely related to the monsoon climate of the Tibetan Plateau, as monsoons account for 58.5% of the annual precipitation [52]. The proportion of rainfall in May-October accounted for 90% of the rain in the eastern part of Tibet [49]. Therefore, the rainy season was also a frequent period of soil water erosion.  The distribution of rainfall erosivity over the 24 and a half months is an input factor for many soil erosion model calculations, including the Chinese Soil Loss Equation (CSLE) [70] and for calculating vegetation cover and biological measure factors. Data from 1981 to 2010 indicated that rainfall erosivity in the Tibetan Plateau was mainly concentrated in June-September, accounting for 81% of the year (Figure 3). This was closely related to the monsoon climate of the Tibetan Plateau, as monsoons account for 58.5% of the annual precipitation [52]. The proportion of rainfall in May-October accounted for 90% of the rain in the eastern part of Tibet [49]. Therefore, the rainy season was also a frequent period of soil water erosion.  Table 2). The climate tendency rate indicated a significant increasing trend (at the 0.05 significance level) at Cona and Nyingchi stations in the southern part of the Tibetan Plateau, and a highly significant increase (at 0.01 significance level) at the Sog station in the upper reaches of the Nu River Basin; the climate tendency rate increased by 10, 20, and 11 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively, for the three stations. Sen's slope analysis indicated that rainfall erosivity in the Tibetan Plateau increased during the study period. Sen's slope value for spring rainfall erosivity was 1.85 (p = 0.005 < 0.01). Rainfall erosivity with a decreasing trend was found in the western part of the Tibetan Plateau where the Zuogong station had an insignificant decreasing trend; its erosivity decreased by 40 MJ·mm·ha −1 ·h −1 ·10a −1 ( Table 2). A spring rainfall erosivity anomaly showed that rainfall erosivity was smaller before 2000, and larger after 2000. During the study period, spring rainfall erosivity was below average 67% of the time; the minimum appeared in 1987 and the maximum in 2010 (Figure 4b Table 2). The climate tendency rate indicated a significant increasing trend (at the 0.05 significance level) at Cona and Nyingchi stations in the southern part of the Tibetan Plateau, and a highly significant increase (at 0.01 significance level) at the Sog station in the upper reaches of the Nu River Basin; the climate tendency rate increased by 10, 20, and 11 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively, for the three stations. Sen's slope analysis indicated that rainfall erosivity in the Tibetan Plateau increased during the study period. Sen's slope value for spring rainfall erosivity was 1.85 (p = 0.005 < 0.01). Rainfall erosivity with a decreasing trend was found in the western part of the Tibetan Plateau where the Zuogong station had an insignificant decreasing trend; its erosivity decreased by 40 MJ·mm·ha −1 ·h −1 ·10a −1 ( Table 2). A spring rainfall erosivity anomaly showed that rainfall erosivity was smaller before 2000, and larger after 2000. During the study period, spring rainfall erosivity was below average 67% of the time; the minimum appeared in 1987 and the maximum in 2010 (Figure 4b). The M-K test results showed an increasing trend in spring erosivity in the Tibetan Plateau since 1987, a significant increasing trend after 2010, and the abrupt year occurred in 2005 (Figure 4c).  Increasing rainfall erosivity in summer was observed at more stations than did decreasing, with 21 of 38 meteorological stations (55%) showing increasing trends in rainfall erosivity (Figure 5a, Table 3). Based on the climate tendency rate, Shenza station in the middle of the Tibetan Plateau, and the Qamdo station of the Lancang River Basin in the eastern Hengduan Mountains showed a significant trend (at 0.05 significance level) with an increase of 31 and 41 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively; the trend at Gerze station in the west and Markam station in the Lancang River Basin in the eastern Hengduan Mountains was highly significant (at the 0.01 significance level) at 24 and 69 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively. The results of Sen's slope estimator indicated an increasing but not significant trend during the study period. The Sen's slope value of summer rainfall erosivity was 1.41 (p = 0.41). Rainfall erosivity showing a decreasing trend was observed mainly in the western part of the Tibetan Plateau, the southern part, and the western part of Hengduan Mountains, with the largest decrease at Gyangze and Lhorong stations where the trend coefficient decreased by 34 and 24 MJ·mm·ha −1 ·h −1 ·10a -1 , respectively (Table 3). A summer rainfall erosivity anomaly showed that summer rainfall erosivity was above average in most years (57%), with the minimum in 1983 and the maximum in 1998 (Figure 5b). The M-K test results indicated an insignificant increasing trend of summer rainfall erosivity in the Tibetan Plateau since 1983 (Figure 5c).  Increasing rainfall erosivity in summer was observed at more stations than did decreasing, with 21 of 38 meteorological stations (55%) showing increasing trends in rainfall erosivity (Figure 5a, Table 3). Based on the climate tendency rate, Shenza station in the middle of the Tibetan Plateau, and the Qamdo station of the Lancang River Basin in the eastern Hengduan Mountains showed a significant trend (at 0.05 significance level) with an increase of 31 and 41 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively; the trend at Gerze station in the west and Markam station in the Lancang River Basin in the eastern Hengduan Mountains was highly significant (at the 0.01 significance level) at 24 and 69 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively. The results of Sen's slope estimator indicated an increasing but not significant trend during the study period. The Sen's slope value of summer rainfall erosivity was 1.41 (p = 0.41). Rainfall erosivity showing a decreasing trend was observed mainly in the western part of the Tibetan Plateau, the southern part, and the western part of Hengduan Mountains, with the largest decrease at Gyangze and Lhorong stations where the trend coefficient decreased by 34 and 24 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively (Table 3). A summer rainfall erosivity anomaly showed that summer rainfall erosivity was above average in most years (57%), with the minimum in 1983 and the maximum in 1998 (Figure 5b). The M-K test results indicated an insignificant increasing trend of summer rainfall erosivity in the Tibetan Plateau since 1983 (Figure 5c).  The number of stations with increasing rainfall erosivity in autumn was slightly lower than that with a decreasing trend; 17 of the 38 meteorological stations (45%) showed an increasing trend for rainfall erosivity (Figure 6a, Table 4). Bases on the climate tendency rate, the Konggar station in the Yarlung Zangbo River valley in the southern part of the Tibetan Plateau (near Zetang station) showed a significant increasing trend (at the 0.05 significance level), and the Lhunze station showed a highly significant increasing trend (at the 0.01 significance level), with an increase of 16 and 13 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively. The results of autumn rainfall erosivity indicated an insignificant decreasing trend during the study period. The Sen's slope value of autumn rainfall erosivity was −0.02 (p = 0.98). In addition to the Yarlung Zangbo River Basin, autumn rainfall erosivity in other areas was mainly decreasing (Table 3). Autumn rainfall erosivity at the Bome and Zuogong stations in the Hengduan Mountains in the eastern part of the Tibetan Plateau decreased by 46 and 9 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively. An autumn rainfall erosivity anomaly shows that autumn rainfall erosivity was above average in most years (50%), with the minimum in 1992 and the maximum in 1985 (Figure 6b). The M-K test results indicated an insignificant increasing trend in the Tibetan Plateau from 1984 to 2003, and a slowly increasing trend after 2008 (Figure 5c).  The number of stations with increasing rainfall erosivity in autumn was slightly lower than that with a decreasing trend; 17 of the 38 meteorological stations (45%) showed an increasing trend for rainfall erosivity (Figure 6a, Table 4). Bases on the climate tendency rate, the Konggar station in the Yarlung Zangbo River valley in the southern part of the Tibetan Plateau (near Zetang station) showed a significant increasing trend (at the 0.05 significance level), and the Lhunze station showed a highly significant increasing trend (at the 0.01 significance level), with an increase of 16 and 13 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively. The results of autumn rainfall erosivity indicated an insignificant decreasing trend during the study period. The Sen's slope value of autumn rainfall erosivity was −0.02 (p = 0.98). In addition to the Yarlung Zangbo River Basin, autumn rainfall erosivity in other areas was mainly decreasing (Table 3). Autumn rainfall erosivity at the Bome and Zuogong stations in the Hengduan Mountains in the eastern part of the Tibetan Plateau decreased by 46 and 9 MJ·mm·ha −1 ·h −1 ·10a −1 , respectively. An autumn rainfall erosivity anomaly shows that autumn rainfall erosivity was above average in most years (50%), with the minimum in 1992 and the maximum in 1985 (Figure 6b). The M-K test results indicated an insignificant increasing trend in the Tibetan Plateau from 1984 to 2003, and a slowly increasing trend after 2008 (Figure 5c).  Winter rainfall erosivity in the Tibetan Plateau was very low, and about 66% of the stations had a value of 0 ( Figure 7a, Table 5). The results indicated an insignificant decreasing trend in winter rainfall erosivity during the study period. Meanwhile, the Sen's slope value of winter rainfall erosivity was −0.03 (p = 0.89). Based on the climate tendency rate, the largest decreasing trend was observed at the Lhunze station and Qamdo in the Hengduan Mountains in the eastern part of the Tibetan Plateau, as well as at the Cona station in the southern region, where it decreased almost 100%. The anomaly showed that winter rainfall erosivity was below average in most years (66%), with the minimum in 1992 and the maximum in 1989 (Figure 7b). The M-K test results demonstrated that the winter rainfall erosivity in the Tibetan Plateau had an insignificant decreasing trend since 1991 (Figure 7c).  Winter rainfall erosivity in the Tibetan Plateau was very low, and about 66% of the stations had a value of 0 ( Figure 7a, Table 5). The results indicated an insignificant decreasing trend in winter rainfall erosivity during the study period. Meanwhile, the Sen's slope value of winter rainfall erosivity was −0.03 (p = 0.89). Based on the climate tendency rate, the largest decreasing trend was observed at the Lhunze station and Qamdo in the Hengduan Mountains in the eastern part of the Tibetan Plateau, as well as at the Cona station in the southern region, where it decreased almost 100%. The anomaly showed that winter rainfall erosivity was below average in most years (66%), with the minimum in 1992 and the maximum in 1989 (Figure 7b). The M-K test results demonstrated that the winter rainfall erosivity in the Tibetan Plateau had an insignificant decreasing trend since 1991 (Figure 7c).

Spatial Distribution of Rainfall Erosivity in the Tibetan Plateau
Spatial distribution of rainfall erosivity in the Tibetan Plateau showed a decreasing trend from the southeast to the northwest (Figure 8). The average annual rainfall erosivity at 38 meteorological stations varied from 61 to 1776 MJ·mm·ha -1 ·h -1 during the study period, with an average of 714 MJ·mm·ha −1 ·h −1 . The largest annual rainfall erosivity was observed at the Bome station in the eastern part of the Tibetan Plateau, and the smallest at the Shiquanhe station in the western region ( Figure  8). Rainfall erosivity was <500 MJ·mm·ha −1 ·h −1 at 24% of the stations, 500-1000 MJ·mm·ha −1 ·h −1 at 55%, and >1000 MJ·mm·ha −1 ·h −1 at 21% of the stations. Relatively high rainfall erosivity was mainly distributed in the Hengduan Mountains in the eastern part of the Tibetan Plateau and in the lowelevation areas between the Yarlung Zangbo and Nu rivers.

Spatial Distribution of Rainfall Erosivity in the Tibetan Plateau
Spatial distribution of rainfall erosivity in the Tibetan Plateau showed a decreasing trend from the southeast to the northwest (Figure 8). The average annual rainfall erosivity at 38 meteorological stations varied from 61 to 1776 MJ·mm·ha −1 ·h −1 during the study period, with an average of 714 MJ·mm·ha −1 ·h −1 . The largest annual rainfall erosivity was observed at the Bome station in the eastern part of the Tibetan Plateau, and the smallest at the Shiquanhe station in the western region ( Figure 8). Rainfall erosivity was <500 MJ·mm·ha −1 ·h −1 at 24% of the stations, 500-1000 MJ·mm·ha −1 ·h −1 at 55%, and >1000 MJ·mm·ha −1 ·h −1 at 21% of the stations. Relatively high rainfall erosivity was mainly distributed in the Hengduan Mountains in the eastern part of the Tibetan Plateau and in the low-elevation areas between the Yarlung Zangbo and Nu rivers. The seasonal distribution of rainfall erosivity varied greatly across the Tibetan Plateau; the average summer rainfall erosivity was 482 MJ·mm·ha −1 ·h −1 , accounting for a maximum of 67.5% of the annual rainfall erosivity, followed by autumn and spring, accounting for 18.5% and 11.5%, respectively. The proportion of winter rainfall erosivity was the smallest at only 2.5%. Seasonal variability in rainfall erosivity varied among meteorological stations, but generally followed rain distribution in summer > autumn > spring > winter (Figure 9). The proportions of spring to annual total rainfall erosivity differed among meteorological stations. The Shiquanhe station had the smallest proportion at 0 and the Zayu station in the Hengduan Mountains had a proportion of 46.5%. Meanwhile, the proportion in the south was higher than in the north, and higher in the east than in the west, resulting in a slight decreasing trend from the southeast to the northwest (Figure 9a). The proportion of summer rainfall erosivity varied from the smallest of 16.4% at Nyalam, to the largest of 94.7% at Shiquanhe station. In addition, 71% of the stations accounted for more than 70% of summer rainfall erosivity, and summer rainfall erosivity in most meteorological stations in the southern Tibet Valley contributed more than 80% to the total (Figure 9b). The proportion of autumn rainfall erosivity at each meteorological station varied between 5.3% and 39.4% with the highest at Purang station, followed by the Nyalan and Cona stations; these stations were located on the southernmost edge of the Tibetan Plateau (Figure 9c). The proportion of winter rainfall erosivity to the total was small, with less than 4% of the stations contributing >90% of the winter erosivity. In that, the Nyalam station contributed the most at 33.9%, followed by Purang, both located on the southwestern edge of the Tibetan Plateau (Figure 9d). The seasonal distribution of rainfall erosivity varied greatly across the Tibetan Plateau; the average summer rainfall erosivity was 482 MJ·mm·ha −1 ·h −1 , accounting for a maximum of 67.5% of the annual rainfall erosivity, followed by autumn and spring, accounting for 18.5% and 11.5%, respectively. The proportion of winter rainfall erosivity was the smallest at only 2.5%. Seasonal variability in rainfall erosivity varied among meteorological stations, but generally followed rain distribution in summer > autumn > spring > winter (Figure 9). The proportions of spring to annual total rainfall erosivity differed among meteorological stations. The Shiquanhe station had the smallest proportion at 0 and the Zayu station in the Hengduan Mountains had a proportion of 46.5%. Meanwhile, the proportion in the south was higher than in the north, and higher in the east than in the west, resulting in a slight decreasing trend from the southeast to the northwest (Figure 9a). The proportion of summer rainfall erosivity varied from the smallest of 16.4% at Nyalam, to the largest of 94.7% at Shiquanhe station. In addition, 71% of the stations accounted for more than 70% of summer rainfall erosivity, and summer rainfall erosivity in most meteorological stations in the southern Tibet Valley contributed more than 80% to the total (Figure 9b). The proportion of autumn rainfall erosivity at each meteorological station varied between 5.3% and 39.4% with the highest at Purang station, followed by the Nyalan and Cona stations; these stations were located on the southernmost edge of the Tibetan Plateau (Figure 9c). The proportion of winter rainfall erosivity to the total was small, with less than 4% of the stations contributing >90% of the winter erosivity. In that, the Nyalam station contributed the most at 33.9%, followed by Purang, both located on the southwestern edge of the Tibetan Plateau (Figure 9d).

Discussion
Rainfall erosivity in the Tibetan Plateau from 1981 to 2015 varied from 61 to 1776 MJ·mm·ha −1 ·h −1 , with an average of 714 MJ·mm·ha −1 ·h −1 . This value is higher than the inland areas of northwest China, but smaller than the areas of east and south China [34]. Yan et al. [48] assessed temporal and spatial distribution of rainfall erosivity in the Tibetan Plateau using TRMM 3B42 data from 2000 to 2008 and showed a decreasing trend during the study period. However, our results showed that from 1980 to 2015, the annual rainfall erosivity indicated an insignificant increase. If we put our research period in 2000-2008, the annual rainfall erosivity also showed a slight decrease. Fan et al. [53] found that rainfall erosivity varied from 32 to 12,189 MJ·mm·ha −1 ·h −1 , with an average of 768 MJ·mm·ha −1 ·h −1 during 2000-2010, which is closed to our findings. Gu et al. [49] also estimated rainfall erosivity in the Hengduan Mountainous Region of Eastern Tibet. Their results showed that rainfall erosivity of Qamdo station based on monthly rainfall was about 300 MJ·mm·ha −1 ·h −1 , which was half of our value. Data precision and models cause the difference. Daily rainfall data is better than monthly rainfall data. Spatial distribution is characterized by a decreasing trend from the southeast to the northwest, which is in accord with previous studies [33,34,50,51].
Topography and monsoon affect rainfall distribution in the study area. The higher the altitude, the less the rainfall (Figure 10a). High altitude and atmospheric circulation patterns have created a cold, dry, and fragile environment, which is very sensitive to soil erosion. Weak rainfall erosivity can cause very severe soil erosion, and terrain has an important influence on the distribution of rainfall erosivity in the Tibetan Plateau. Previous studies have shown that altitude, slope, and aspect affect transport of water vapor and the spatial distribution of precipitation in the study area [71]. The relationship between the average annual rainfall erosivity and altitude of meteorological stations revealed a high correlation coefficient (−0.63), with rainfall erosivity decreasing with increasing altitude (Table 6, Figure 10b). The main factors affecting rainfall erosivity were rainfall and rainfall

Discussion
Rainfall erosivity in the Tibetan Plateau from 1981 to 2015 varied from 61 to 1776 MJ·mm·ha −1 ·h −1 , with an average of 714 MJ·mm·ha −1 ·h −1 . This value is higher than the inland areas of northwest China, but smaller than the areas of east and south China [34]. Yan et al. [48] assessed temporal and spatial distribution of rainfall erosivity in the Tibetan Plateau using TRMM 3B42 data from 2000 to 2008 and showed a decreasing trend during the study period. However, our results showed that from 1980 to 2015, the annual rainfall erosivity indicated an insignificant increase. If we put our research period in 2000-2008, the annual rainfall erosivity also showed a slight decrease. Fan et al. [53] found that rainfall erosivity varied from 32 to 12,189 MJ·mm·ha −1 ·h −1 , with an average of 768 MJ·mm·ha −1 ·h −1 during 2000-2010, which is closed to our findings. Gu et al. [49] also estimated rainfall erosivity in the Hengduan Mountainous Region of Eastern Tibet. Their results showed that rainfall erosivity of Qamdo station based on monthly rainfall was about 300 MJ·mm·ha −1 ·h −1 , which was half of our value. Data precision and models cause the difference. Daily rainfall data is better than monthly rainfall data. Spatial distribution is characterized by a decreasing trend from the southeast to the northwest, which is in accord with previous studies [33,34,50,51].
Topography and monsoon affect rainfall distribution in the study area. The higher the altitude, the less the rainfall (Figure 10a). High altitude and atmospheric circulation patterns have created a cold, dry, and fragile environment, which is very sensitive to soil erosion. Weak rainfall erosivity can cause very severe soil erosion, and terrain has an important influence on the distribution of rainfall erosivity in the Tibetan Plateau. Previous studies have shown that altitude, slope, and aspect affect transport of water vapor and the spatial distribution of precipitation in the study area [71]. The relationship between the average annual rainfall erosivity and altitude of meteorological stations revealed a high correlation coefficient (−0.63), with rainfall erosivity decreasing with increasing altitude (Table 6, Figure 10b). The main factors affecting rainfall erosivity were rainfall and rainfall intensity. There is a strong correlation between rainfall and rainfall erosivity in the study area (Figure 10c). The correlation coefficient is as high as 0.92.  The north-south pattern of the mountains and river valleys in Hengduan Mountains in eastern Tibet was conducive to the entry of water vapor from the Indian Ocean, and the rainfall was relatively abundant in that area. Therefore, the annual rainfall erosivity was also relatively high in eastern Tibet, with annual, spring, and summer rainfall erosivity increasing during the study period. In northwestern Tibet, precipitation was extremely small, and the corresponding rainfall erosivity was low; this was due to the obstruction by the mountain system of the Himalayas coupled with high altitude and low temperature. In southern Tibet and in the Himalayas, annual rainfall erosivity was small, mainly due to the local topography of the leeward slope of the mountain. However, in the Yarlung Zangbo river valley, rainfall erosivity was higher due to a lower terrain and sufficient water and heat. Additionally, precipitation in the rainy season has a pronounced vertical differentiation in the plateau, and with the increase in altitude, precipitation significantly decreases, and then increases [71], leading to complexity in the spatial distribution of rainfall erosivity.
Rainfall erosivity in the Tibetan Plateau exhibited a decreasing trend from the southeast to the northwest, due to mainly to the monsoon. Specifically, the trend was affected by the southeast monsoon from the Pacific Ocean, and by the southwest monsoon from the Indian Ocean [72]. The climate of the Tibetan Plateau is diverse, from southeast to northwest, and includes tropical, subtropical, temperate plateau, plateau sub-frigid, and plateau cold zones. The southern Tibetan valley belongs to the subtropical monsoon climate, with relatively abundant precipitation, and relatively high rainfall erosivity. The annual, spring, and summer rainfall erosivity showed an increasing trend in the Tibetan Plateau, with spring erosivity showing an especially significant increase, which was closely related to climate change in this area. The region experienced a warming and more humid trend [73]. A previous study has shown that there was a slight increase in annual rainfall erosivity in most parts of China. Besides, the central and eastern Qinghai-Tibet Plateau is the region with the most significant increase in rainfall erosivity [33]. Our results also showed that the  Besides, the central and eastern Qinghai-Tibet Plateau is the region with the most significant increase in rainfall erosivity [33]. Our results also showed that the rainfall erosivity of most meteorological stations on the Tibetan Plateau was increasing. This variation may increase the risk of soil erosion in the region. The Tibetan Plateau is located in the most active geological and structural tectonic belt with the most varied geological history and the strongest tectonic movement in China; physical weathering of the ground is strong, the freezing and thawing effects are widely distributed, and the material on slopes is unstable, leading to an extremely fragile environment. As shown in this study, the risk of soil erosion in the Hengduan Mountains in the eastern Tibetan Plateau and the Yarlung Zangbo River Valley and its vicinity were high. Rainfall erosivity in this area was not only high, but also showed an increased trend, especially in spring. As the temperature rises, the snow melts and erosion caused by rainfall makes the risk of soil erosion in spring also high. Previous studies also showed that the region with the largest annual erosion modulus in Tibet was the "Three Rivers Basin" in the east (Jinsha River, Lancang River, and Nu River), followed by the Yarlung River valley [74]. The terrain of the Tibetan Plateau is complex; heavy rainfall is likely to cause soil erosion, as well as floods, landslides, and mudslides. Annual rainfall erosivity in the Hengduan Mountains in eastern Tibet was relatively high, and mainly concentrated in summer. The area was high, and with sleep slopes and deep river gorges, and the mountains and valleys coexist. Coupled with flood disasters, it is a region with high risk of soil erosion. Due to the influence of the southwest monsoon in the Indian Ocean, the southern Tibetan valley has excellent hydrothermal conditions and can grow many subtropical crops, but human activities undermine environmental protection, which results in creation of important sources of soil erosion from low-coverage grasslands and sloping farmland in the valley.

Conclusions
The Tibetan Plateau is an ecologically fragile area. Climate in the Tibetan Plateau has changed significantly in the early 21st century, especially due to significant warming. The risk of soil erosion brought about by climate change may be significant. In this study, we collected daily precipitation data from 1 January 1981 to 31 December 2015 at 38 meteorological stations in the Tibetan Plateau, and calculated annual and seasonal rainfall erosivity. Temporal and spatial variability patterns of rainfall erosivity were analyzed. We concluded that: (1) Rainfall erosivity in the Tibetan Plateau from 1981 to 2015 varied from 61 to 1776 MJ·mm·ha −1 ·h −1 , with an average of 714 MJ·mm·ha −1 ·h −1 . The spatial distribution of rainfall erosivity differed significantly across the area, with the value decreasing from the southeast to the northwest.
(2) Summer rainfall erosivity accounted for 67.5% of the total annual rain erosivity, followed by autumn and spring, accounting for 18.5% and 11.5%, respectively. Between 1981 to 2015, and especially since the mid-1980s, annual and summer rainfall erosivity in the Tibetan Plateau exhibited an increasing trend (increased by 50% and 60%); spring erosivity had a highly significant increasing (increased by 94%), while autumn and winter a slightly decreasing trend (less than 10%).
(3) The spatial distribution and temporal variability in rainfall erosivity were affected by the special topography and monsoon of the Tibetan Plateau, as well as by increasing warming and humidification in this area. In contrast, soil erosion risk was high in the Hengduan Mountains in the eastern Tibetan Plateau and in the Yarlung Zangbo River Valley and its vicinity.