Spatiotemporal variation and tendency analysis on rainfall erosivity in the Loess Plateau of China

Soil erosion is one of the most important ecological problems in the world. Nearly a third of arable land has been exposed to soil erosion, which has increased at a rate of 1.0 × 10⁷ hm² per year (Pimentel et al. 1995). Soil erosion removes fertile soil and seriously threatens agricultural production and food security. The eroded soil usually enters rivers, causing siltation and water eutrophication, and thus endangers water safety and affects normal production and life. Predicting and protecting against soil erosion are of great significance to the construction of an ecological civilization and to the sustainable development of the economy and society. Consequently, mechanisms of soil erosion and measures to control it have been hot research topics (Lal 1989; Pimentel 2006).

Rainfall is the driving force and prerequisite for soil erosion (Zhang et al. 2002). The Universal Soil Loss Equation (USLE) devised by Wischmeier & Smith (1978) and the Revised Universal Soil Loss Equation (RUSLE) proposed by Renard & Freimund (1994) and Renard et al. (1997) for the United States are widely used to predict soil erosion (Anees et al. 2017). The rainfall erosivity (R-factor) in these models indicates the ability of rainfall to cause soil erosion and is considered the best factor for studying the response of soil erosion to rainfall changes (Nearing et al. 2004). Meusburger et al. (2012) collected rainfall data for 1989–2010 and used the product of the total kinetic energy of rainfall (E) and the 30-min maximum rainfall intensity (I₃₀) as the measurement index of rainfall erosivity to study the temporal and spatial variation of rainfall erosivity in Switzerland. Da Silva (2004) used rainfall data in conjunction with geographic information system (GIS) for spatial interpolation to map the spatial distribution of rainfall erosivity in Brazil. Some other researchers have conducted spatiotemporal analyses of rainfall erosivity in different regions (Angulo-Martínez & Beguería 2009; Panagos et al. 2015; Xie et al. 2016). A high-resolution distribution map of rainfall erosivity has been established using global data, and the global rainfall erosivity was estimated to be about 2,190.0 MJ mm hm⁻² h⁻¹ (Panagos et al. 2017). Due to the difficulty in obtaining rainfall data with a high temporal resolution, Vrieling et al. (2010) used satellite data and combined rainfall intensity and rainfall erosivity to establish a method for forecasting rainfall erosivity, thereby demonstrating a potential tool for soil erosion prediction in data-poor areas. Previous studies have enhanced our understanding of spatiotemporal variations in rainfall erosivity, but these studies were constrained by the limited availability of daily rainfall data and by the length of rainfall data. Studying available long-term rainfall data and estimating rainfall erosivity could provide information on the potential trends of soil erosion, especially in areas that are lacking data.

In the Loess Plateau of China, intensive rainfall and sparse vegetation can lead to relatively serious soil erosion. The soil erosion modulus can reach 1.0 × 10⁷ kg km⁻² a⁻¹. A large amount of sediment enters the Yellow River and causes many problems, including river channel siltation and the deterioration of aquatic ecosystems (Liu & Liu 2010). Many researchers have analyzed the spatial and temporal distribution of rainfall erosion on the Loess Plateau in different periods (Xin et al. 2011; Abd Elbasit et al. 2013; Yang & Lu 2015). Fu et al. (2005) calculated all the RUSLE parameters using GIS based on regional data from the Yanhe watershed and concluded that the annual soil loss was 0.5–2.0 × 10⁷ kg km⁻². Sun et al. (2014) analyzed the effects of topography and land-use patterns on soil erosion in the Loess Plateau and suggested that ‘Grain-to-Green Program’ is effective for preventing soil erosion. Xu (2005) carried out a quantitative analysis of the correlations between vegetation coverage (C_f), annual rainfall erosivity (R_e), and annual precipitation (P_m). Xu's results indicated that R_e will increase rapidly when P_m > 300 mm, and that when P_m > 530 mm, the rate at which R_e increases with P_m becomes higher. Current research generally focuses on specific watersheds where daily rainfall data are easily available. Consequently, there are relatively few papers addressing the long-term spatiotemporal variation of rainfall erosivity in the Loess Plateau. More analysis is needed for the region as a whole based on available long-term daily rainfall data. The results of the analysis could form the basis for soil loss prediction in the region. Moreover, the impact of climate change on soil erosion needs to be understood (Almagro et al. 2017). Few studies, however, have been carried out that have made predictions and analyses of variations in rainfall erosivity under projected climate change in the Loess Plateau of China. Consequently, a climate change perspective and decision basis cannot be provided for decision makers.

The overall aim of this study was to address these knowledge gaps. Precipitation data from 1971 to 2010 were collected from 39 representative stations in the Loess Plateau of China to analyze the spatiotemporal variability of rainfall erosivity. The main objectives of this study are:

1. to determine the spatial and temporal variability of annual rainfall erosivity in the Loess Plateau;

2. to predict the change tendency of future rainfall erosivity using a typical global circulation model (GCM).

The Loess Plateau is located between 100°52′–114°31′E and 33°37′–41°25′N, including most or part of Shanxi, Shaanxi, Gansu, Qinghai, and Henan Provinces and the Inner Mongolia and Ningxia Regions. The total area is about 6.35 × 10⁵ km². Most areas of the plateau have semi-humid or semi-arid climates, and the climates in different areas are quite different. The annual mean temperature is 3.6–14.3 °C, and the annual evaporation is 1,400–2,000 mm. The main crops grown are wheat, corn, soybeans, and sorghum. The study area is characterized by complex landforms and deep soil layers (mainly loose and easily eroded dark loessial soil and loessal soil). The annual precipitation ranges from 200 to 750 mm, with large inter-annual and seasonal variations. There are frequent storms, which are the main driving force of soil erosion. The annual mean soil erosion modulus of the Loess Plateau is 0.5–1.0 × 10⁷ kg km⁻² a⁻¹, and the thickness of the soil layer lost every year is about 1 cm.

We used data from the 39 representative meteorological stations located in the Loess Plateau. The daily rainfall data from 1971 to 2010 were selected for investigating the spatiotemporal variability of precipitation. The location of each meteorological station is shown in Figure 1: SRTM (Shuttle Radar Topography Mission) data from http://srtm.csi.cgiar.org/srtmdata/, with spatial resolution approximately 30 meters on the line of the equator. Daily rainfall data were obtained from the China Meteorological Data Service Center (http://data.cma.cn/). The Linfen meteorological station does not have data for 2000, and the Yongji station does not have data for 1981–1990.

In addition to the data from the above meteorological stations, data were also obtained using a GCMs (Hadley Center Global Environment Model version 2 (HadGEM2-ES)) forced by two Representative Concentration Pathway (RCP) scenarios – RCP4.5 and RCP8.5 (Chou et al. 2014; Yan et al. 2015). The HadGEM2-ES is an earth system category GCM that was developed by the Hadley Center, and the resolution is about 1.875° for longitude and 1.25° for latitude (Collins et al. 2011). The RCP scenarios were developed by the research community by considering emissions, concentrations, and land-use trajectories and are labeled according to the expected values of global radiative forcing in 2100. RCP4.5 is considered an intermediate scenario that assumes emissions reduction during the 21st century by employing clean technologies and stringent climate policies. This scenario predicts a global forcing radiation of ∼4.5W m² and concentrations of ∼650 ppm CO2-eq at stabilization after 2100. In contrast, RCP8.5 is a kind of business-as-usual scenario characterized by no implementation of climate policies, low rates of technology development, high energy intensities, and a reliance on fossil fuels that leads to higher emissions than RCP4.5 over time. Global forcing radiation and atmospheric concentrations of CO2-eq are projected to reach >8.5 W m⁻² and >1,370 ppm in 2100, respectively (Almagro et al. 2017). Daily precipitations for 2020–2100 are available for the RCP scenarios. The boundaries of the Loess Plateau were taken to 33.75°–41.25°N and 101.25°–114.375°E. Daily rainfall data for 48 sample spots were obtained (Figure 1) from https://esgf-data.dkrz.de/search/cmip5-dkrz/.

In the present study, the average rainfall erosivity refers to the average of data from different meteorological stations over a given time, while the annual mean rainfall erosivity refers to the average of data over a period of time for a specific location.

MATLAB R2018b (The MathWorks, Inc., USA) and ORIGIN 2020 (OriginLab Corporation, USA) were used for statistical analysis and mapping of data. The ordinary Kriging interpolation method and the geostatistical modules of ArcGIS 10.2 (Esri, Inc., USA) were used for the spatial interpolation of precipitation and rainfall erosivity. In addition, the non-parametric Mann–Kendall trend analysis method was used to test the significance of changes in meteorological elements (Hamed & Ramachandra Rao 1998; Ahmad et al. 2015).

Data from the Linfen and Yongji meteorological stations were excluded from the analysis because data from some years were missing. Data of the remaining 37 meteorological stations were averaged to obtain the precipitation variation for the Loess Plateau from 1971 to 1990 and from 1991 to 2010 (Figure 2(a) and 2(b)). The mean precipitation from 1971 to 1990 and from 1991 to 2010 was 425.26 and 399.87 mm, respectively. The annual precipitation ranges were 205.39 mm (1971–1990) and 251.99 mm (1991–2010), and the coefficients of variation were 0.15 and 0.14, respectively. In general, years with a mean precipitation of >500 mm in the Loess Plateau from 1971 to 2010 accounted for only 7.5% of the years, and years with less than 400 mm accounted for 45.0% (1971–1990 accounted for 15.0%). The annual precipitation showed an overall decreasing trend from 1971 to 2010.

The mean erosive rainfall on the Loess Plateau from 1971 to 1990 and from 1991 to 2010 was 242.59 and 229.49 mm, respectively, which accounted for about 57.0% of the mean precipitation, and showed a decreasing trend from 1971 to 2010. There is a strong linear regression relationship between mean erosive rainfall and precipitation (R² = 0.93, P < 0.05). The results of the Mann–Kendall test showed that the erosive rainfall trend from 1971 to 2010 has increased from 1973 to 1981 and decreased from 1982 to 2010 (except for 1985 and 1996). Both trends are insignificant and changeable, which indicates that erosive rainfall is not only affected by the precipitation amount but also reflects the annual changes in precipitation.

The intra-annual distribution of rainfall erosivity for the years from 1971 to 2010 is shown in Figure 3(a) and 3(b). The rainfall erosivity during the year formed a single-peak distribution, first increasing and then decreasing. The annual precipitation distribution and the rainfall erosivity followed a similar pattern: in January and February, the rainfall erosivity was almost zero, it then gradually increased, reaching a peak in July and August, and then gradually decreased to zero in November. The total rainfall erosivity in July to August accounted for 61.5% (1971–1990) and 59.4% (1991–2010) of the year's rainfall erosivity, indicating that summer rainfall, especially in July and August, was the main cause of soil erosion. In 1991–2010, the precipitation decreased slightly from that of 1971 to 1990, but the maximum still occurred in July and August, and was about 200.0 MJ mm hm⁻² h⁻¹ for both series.

The average rainfall erosivity for 1971–1990 and 1991–2010 was 1,278.68 and 1,200.60 MJ mm hm⁻² h⁻¹, respectively, and the coefficients of variation were 0.24 and 0.19, respectively. Rainfall erosivity was essentially stable at around 1,200 MJ mm hm⁻² h⁻¹ from 2004 to 2010 (Figure 4(b)). Overall, rainfall erosivity fluctuated with the increases and decreases in annual precipitation, but there were some years where the correlation between rainfall erosivity and precipitation was not strong. Taking 1989 and 1990 as an example, the rainfall erosivity was basically the same, but the difference in precipitation was 17.5%. This shows that rainfall erosivity is also affected by rainfall distribution and other factors. The Mann–Kendall tests showed that, in general, the inter-annual changes in rainfall erosivity increased from 1975 to 1981 and decreased from 1982 to 2010. This is similar to the trends for erosive rainfall, but, as shown in Figure 4(a) and 4(b), the annual decrease rate in rainfall erosivity during 1991–2010 was slightly slower than that during 1971–1990, and eventually stabilized.

As shown in Figure 5, the annual mean rainfall erosivity decreased from southeast to northwest. For the area as a whole, the annual mean rainfall erosivity was 1,286.02 and 1,201.46 MJ mm hm⁻² h⁻¹ a⁻¹ in 1971–1990 and 1991–2010, respectively. The maximum values obtained (for the Wutai Mountain and Anze station in the province of Shanxi) were 2,433.14 MJ mm hm⁻² h⁻¹ a⁻¹ (1971–1990) and 2,144.56 MJ mm hm⁻² h⁻¹ a⁻¹ (1991–2010). The minimum values obtained for the two periods (for the Jingyuan station, Gansu province) were 277.57 and 298.94 MJ mm hm⁻² h⁻¹ a⁻¹, respectively. The annual mean rainfall erosivity from 1991 to 2010 was slightly lower than that from 1971 to 1990 and showed a very high decrease in the northern part of Shanxi Province. The areas with the greatest rainfall erosivity were, however, concentrated in the southeast of the Loess Plateau. Consequently, most of Shanxi Province, central and southern Shaanxi were susceptible to rainfall erosion from 1971 to 2010. In these areas, soil conservation measures should be strengthened to reduce soil erosion.

Precipitation data for representative areas of the Loess Plateau are shown in Table 1. (The data were obtained using ordinary Kriging interpolation. Because there are very few representative meteorological stations that could be selected in Henan and Qinghai Provinces, and the areas of these Provinces are relatively small, these Provinces are not included in Tables 1 and 2). The annual mean precipitation from 1971 to 1990 was 398.50 ± 114.80 mm a⁻¹. The annual mean precipitation was about 500.0 mm a⁻¹ for Shanxi and Shaanxi Provinces, and below 300.0 mm a⁻¹ for the Ningxia and Inner Mongolia Regions, while the annual mean precipitation value distribution from 1971 to 1990 was relatively discrete in Gansu Province, which was about 170.60–598.97 mm a⁻¹. The annual mean precipitation from 1991 to 2010 was 143.45–585.32 mm a⁻¹, which was 5.5% lower than that from 1971 to 1990. The province with the highest mean precipitation changed from Shanxi Province in 1971–1990 to Shaanxi Province in 1991–2010, while Ningxia still had the lowest. The standard deviation (STD) of the average precipitation in Shanxi Province decreased by 32.8% from 1971–1990 to 1991–2010. The reduction in the other regions was, however, much smaller.

Under the RCP4.5 scenario, the annual mean precipitation from 2020 to 2100 would be 505.12 ± 130.02 mm a⁻¹, which is 26.8 and 34.1% higher than the means for 1971–1990 and 1991–2010, respectively. The maximum precipitation in RCP4.5 would reach 1,050.60 mm a⁻¹. The spatial distribution of precipitation similar to the spatial distribution for 1971–1990: Shanxi Province had the largest annual mean precipitation of 680.54 mm a⁻¹, while the provinces with the lowest precipitation changed from Ningxia (1971–2010) to Inner Mongolia, with a precipitation of only 341.16 mm a⁻¹. The range of precipitation under the RCP8.5 scenario was 219.93–1,012.40 mm a⁻¹, which was an increase of about 35.9% compared with 1991–2010. The regional variations in annual mean precipitation were essentially the same as those in the RCP4.5 scenario, with a maximum of 1,000.0 mm a⁻¹ (Gansu and Shaanxi) and a minimum of only 219.93 mm a⁻¹ (Inner Mongolia).

From historical data (Table 2), the annual mean rainfall erosivity in the Loess Plateau from 1971 to 1990 was calculated as 251.49–2,405.43 MJ mm hm⁻² h⁻¹ a⁻¹, with an average of 1,137.75 ± 487.57 MJ mm hm⁻² h⁻¹ a⁻¹. Shanxi Province had the highest rainfall erosivity, which reached 2,405.43 MJ mm hm⁻² h⁻¹ a⁻¹ from 1971 to 1990, followed by Shaanxi Province with 1,557.78 MJ mm hm⁻² h⁻¹ a⁻¹. The rainfall erosivity in Ningxia was low, and the STD was also low. The STD for Gansu was the largest; this indicates a large spatial difference, which is similar to the results obtained from the analysis of data from meteorological stations for 1991–2010. The annual mean rainfall erosivity from 1991 to 2010 was 1,082.32 ± 418.69 MJ mm hm⁻² h⁻¹ a⁻¹. Compared with 1971–1990, except for Ningxia, which showed an erosivity increase of 16.6% for 1991–2010, other provinces and regions all showed decreases. For example, erosivity in Shanxi and Shaanxi decreased by 10.1 and 3.0%, respectively, which represents small spatial differences.

The calculated rainfall erosivity in future climate scenarios was lower than that obtained from the historical data. Under the RCP4.5 and RCP8.5 scenarios, the rainfall erosivity in 2020–2100 would be 769.64 ± 433.32 and 868.49 ± 507.97 MJ mm hm⁻² h⁻¹ a⁻¹, which is a 28.9 and 19.8% decrease, respectively, compared with 1991–2010. The areas with the highest rainfall erosivity are still in Shanxi Province in the RCP scenarios, while the lowest areas change from being in Ningxia (according to the historical data) to Inner Mongolia. These results are similar to the results for the distribution of precipitation. The minimum values for erosivity from 2020 to 2100 would be 161.32 MJ mm hm⁻² h⁻¹ a⁻¹ (RCP4.5) and 135.53 MJ mm hm⁻² h⁻¹ a⁻¹ (RCP8.5). In the RCP8.5 scenario, the rainfall erosivity ranges from 135.53 to 3,494.06 MJ mm hm⁻² h⁻¹ a⁻¹ and is 12.8% higher than the RCP4.5 values. The value distribution under RCP8.5 is, however, more discrete than under RCP4.5. The maximum erosivity in Shanxi Province increases by 35.8% from RCP4.5 to RCP8.5, while the increases in the other provinces and regions are small. This arises because heavy rainfall becomes more concentrated in Shanxi Province in the RCP8.5 scenario.

The rainfall erosivity in different periods was divided into two parts according to the length of the series, after the annual mean rainfall erosivity was interpolated by ordinary Kriging in ArcGIS in every part, and the rainfall erosivity variation figure (Figure 6) was obtained using raster calculation. The variation in the rainfall erosivity from 1971–1980 to 1981–1990 was −604.89 to 394.74 MJ mm hm⁻² h⁻¹ a⁻¹. In the southern Gansu Province, the central part of Shaanxi Province, and the northern Shanxi Province, rainfall erosivity from 1981 to 1990 was less than that in 1971–1980, while the rainfall erosivity in Ningxia and some parts of Gansu Province and the southern part of Shanxi Province increased the most. The spatial variation of rainfall erosivity from 1991 to 2010 was not the same as that from 1971 to 1990. After the 2000s, rainfall erosivity increased significantly in southern Shaanxi Province and the central part of Shanxi Province, while the other areas showed an overall decreasing trend. For example, rainfall erosivity in Qinghai Province and Inner Mongolia decreased by about 429.88 MJ mm hm⁻² h⁻¹ a⁻¹.

Under the RCP4.5 scenario, the variation of rainfall erosivity from 2020–2060 to 2061–2100 was −66.89 to 595.59 MJ mm hm⁻² h⁻¹ a⁻¹. The areas with increased rainfall erosivity are concentrated in the southern Loess Plateau, particularly in southern Shanxi and in the central and southern part of Shaanxi Province, while the rainfall erosivity in the other regions in 2061–2100 is slightly less than that in 2020–2060. In contrast, under the RCP8.5 scenario, the rainfall erosivity increased greatly and was concentrated in the southeast of the Loess Plateau (southern Shanxi and northern Henan Provinces). In the other regions, rainfall erosivity in RCP8.5 also increased, and the overall increase was around 99.19–2,009.03 MJ mm hm⁻² h⁻¹ a⁻¹. This indicates that the rainfall erosivity increased during the second half of the 21st century under the RCP8.5 scenario, and hence, the degree of soil erosion in the southeast of the Loess Plateau would also increase.

The mean precipitation on the Loess Plateau was 412.57 mm from 1971 to 2010 and showed an insignificant decreasing trend, with an annual decrease of about 0.88 mm. These results are similar to those in the study by Wu et al. (2016). The mean precipitation in 2020–2100, however, increased by 30.3% (RCP4.5) and 32.1% (RCP8.5) compared with that in 1971–2010. The spatial distribution of future precipitation indicates that precipitation will increase in the southeast of the study area and decrease in the northwest. Consequently, Shanxi Province may have an increased risk of soil erosion.

Precipitation distribution causes annual rainfall erosivity to appear as a unimodal distribution (Figure 3). The maximum half-month rainfall erosivity between 1971–1990 and 1991–2010 was about 200.0 MJ mm hm⁻² h⁻¹ and occurs mainly in July to August, which indicates that heavy rainfall during the summer is the main cause of soil erosion in the Loess Plateau. The comprehensive utilization of engineering measures (e.g. slope managements), vegetation measures (e.g. afforestation), and agricultural measures (e.g. no-tillage) are the key elements of ecological environmental development that are required to prevent soil erosion in the Loess Plateau.

In terms of inter-annual changes, the annual mean rainfall erosivity from 1971 to 1990 was 1,139.75 ± 487.57 MJ mm hm⁻² h⁻¹ a⁻¹, with an insignificant decreasing rate of 1.83 MJ mm hm⁻² h⁻¹ a⁻¹ (Figure 4). When compared with the erosivity from 1971 to 1990, the annual mean rainfall erosivity from 1991 to 2010 decreased by 5.0%. Similarly, the rainfall erosivity from 2020 to 2100 was 28.9% (RCP4.5) and 19.8% (RCP8.5) less than that in 1991–2010. Almagro et al. (2017) also obtained similar results for Brazil. The similarity in results may have been caused by the rainfall in the future climate scenarios being relatively dispersed, with few short-duration heavy rainfall events. An alternative explanation is that the HadGEM2-ES model outputs are based on meteorological factors and ignore the effect of topography on rainfall intensity (Djebou et al. 2014). Consequently, the model outputs may lead to underestimations of the erosivity of future rainfall. For further research, comprehensive analysis should be combined with factors such as topography and land-use data. However, despite the aforementioned limitations, the results of the present study can still be used as a reference for the study of soil erosion under climate change in the future.

To improve the analysis of the relationship between rainfall erosivity and precipitation over various time scales, the present study made use of the Standardized Precipitation Index (SPI) and used the SPI PROGRAM to calculate the SPI value (obtained from https://drought.unl.edu/droughtmonitoring/SPI/SPIProgram.aspx). The time scale of the historical data was 36 and 60 months, and 60 and 96 months for the estimated data. The results are shown in Supplementary material, Appendix B. Based on the SPI value and the classification criteria of McKee et al. (1993), cases were selected for a moderate drought year (−1.49 ≤ SPI ≤ −1.00), a normal year (−0.99 ≤ SPI ≤ 0.99), and a moderate humid year (1.00 ≤ SPI ≤ 1.49). The selected hydrological years and SPI values are shown in Table 3.

It can be seen from Figure 7 that the monthly rainfall erosivity has a good corresponding relationship with precipitation. However, the precipitation amount is not the determining condition for the variation in rainfall erosivity. For example, the precipitation in June 2005 decreased by 11.1% compared with May 2005, but the rainfall erosivity increased by 74.9%. The increase was mainly the result of the greater rainfall intensity in summer. The threshold of erosive rainfall adopted in this paper was a daily rainfall of P ≥ 12 mm. Consequently, a lot of the precipitation was excluded and could not be used as erosive rainfall when calculating rainfall erosivity. The use of the erosive rainfall threshold also leads to deviations between rainfall erosivity and precipitation (Xie et al. 2016).

The percentage of precipitation from April to June in moderate drought years was relatively large and accounted for 34.0% of the annual precipitation (1971–1990), while it only accounted for 20.0% of the annual precipitation in normal and humid years. In moderate drought years, the precipitation from July to August accounted for about 25.3% of the annual precipitation, while in normal and moderate humid years, it reached more than 35.0% (1991–2010). The historical data show that the monthly distribution of rainfall erosivity during different hydrological years is essentially the same. In drought years, rainfall erosivity increased from April and reached a maximum in May, while in normal and humid years, the maxima were often reached in August and September, and then gradually decreased and reached zero in November. The rainfall erosivity from July to September in drought years accounted for about 50.0% for the whole year, while in normal and humid years, this period accounted for 75.7–90.0% of the year's total.

A discussion of the data obtained from the RCP scenarios follows. When compared with historical data, it can be seen that precipitation derived from the RCP scenarios is more dispersed and is mainly concentrated in April to September, with the peak generally occurring in July and August (Figure 7(c) and 7(d)). Although the precipitation increased from April to June, most of the precipitation events were of long-duration and low-rainfall intensity, with daily precipitation less than 12 mm. This precipitation, therefore, had little effect on the distribution of annual rainfall erosivity. The long-duration, low-intensity rainfall events were the results of climate change and other factors. Precipitation from April to June accounted for about 30.0% of the annual precipitation on the whole, while rainfall erosivity during this period only accounted for about 25.0% of the annual rainfall erosivity. Except for the heavy rainfall erosivity that occurred in July in humid years, the monthly rainfall erosivity values for different hydrological years were essentially below 500.0 MJ mm hm⁻² h⁻¹. In a similar fashion to the results from the historical data, 83.3% of the months had rainfall erosivity below 200.0 MJ mm hm⁻² h⁻¹ in the RCP4.5 and RCP8.5 scenarios.

In the 1980s, rainfall erosivity decreased in southern Gansu, central Shaanxi, and northern Shanxi and increased slightly in Ningxia and some parts of Gansu and southern Shanxi. After 2000, it increased in central and southern Shaanxi and central Shanxi but decreased in Inner Mongolia. However, on the whole, the spatial distribution of rainfall erosivity in the Loess Plateau shows a decreasing trend from the southeast to the northwest. Most parts of Shanxi Province, and the central and southern parts of Shaanxi remained the areas most threatened by rainfall erosion in the Loess Plateau, and thus requires ecological management interventions to prevent soil erosion.

In the future climate scenarios, Shanxi continues to be the province with the greatest rainfall erosivity, while the area with the lowest rainfall erosivity changes Ningxia (as determined by historical data) to Inner Mongolia. This represents a relatively regular distribution of highs in the southeast and lows in the northwest. Areas with increased rainfall erosion after 2060 are also concentrated in the south (RCP4.5) and southeast (RCP8.5). In the RCP4.5 and RCP8.5 scenarios, the rainfall erosivity in 2020–2100 decreased by 28.9 and 19.8%, respectively, compared with that in 1991–2010, respectively, but the range of variations is larger (Table 2). The maximum annual mean rainfall erosivity reaches 3,494.06 MJ mm hm⁻² h⁻¹ a⁻¹, while the minimum is only 135.53 MJ mm hm⁻² h⁻¹ a⁻¹.

To analyze the causes of increased precipitation and reduced rainfall erosivity in future climate modes and to further study the quantitative relationship between rainfall erosivity and precipitation, the present study compared annual rainfall erosivity and annual precipitation. (Missing data were recorded as 0, and the historical data comprised 39 × 40 = 1,560 pairs; 26 sample spots in the Loess Plateau were selected for the predicted data, and 26 × 81 = 2,106 pairs of data were generated in the RCP4.5 and RCP8.5 scenarios, respectively.) The fitting results are shown in Figure 8.

When compared with the historical data, it was evident that the annual precipitation in the RCP4.5 scenario was mainly concentrated in the 300–900 mm range, and the annual rainfall erosivity was generally below 2,500 MJ mm hm⁻² h⁻¹. In contrast to the RCP4.5 scenario, the annual rainfall erosivity in the RCP8.5 scenario was mainly concentrated below 3,000 MJ mm hm⁻² h⁻¹. The logistic function fits the RCP4.5 data well (R² = 0.49, P < 0.05), but the function does not converge with the RCP8.5 data. The Allometric1 function in ORIGIN 2020 was used to do the fitting, and a good fitting effect was also obtained for the RCP8.5 data (R² = 0.55 (P < 0.05), Figure 8(c)). For the same rainfall amount, when P < 442.51 mm, the rainfall erosivity of the RCP4.5 scenario is greater than that of RCP8.5, and vice versa. Under the future climate modes, the annual rainfall erosivity gradually increases and becomes more discrete with an increase in precipitation, which indicates that the rainfall data generated by the HadGEM2-ES model are complex. However, as Figure 8 shows, when the annual precipitation is less than 1,000.0 mm, the 95% confidence band is narrow, which indicates that two functions are applicable for the prediction of rainfall erosivity using data from HadGEM2-ES. It is clear, therefore, that the detailed changes of future rainfall erosivity require further quantitative research.

Rainfall data from 1971 to 2010 for 39 typical meteorological stations in the Loess Plateau were collected to calculate the spatiotemporal variation of rainfall erosivity. Additionally, future trends in rainfall erosivity were predicted using the RCP4.5 and RCP8.5 scenarios. The main conclusions are as follows:

1. From 1971 to 2010, the mean rainfall erosivity on the Loess Plateau was 1,239.64 MJ mm hm⁻² h⁻¹. Overall, annual rainfall erosivity exhibited a slightly decreasing trend. On a monthly basis, rainfall erosivity appeared as a unimodal distribution throughout the year and was mainly concentrated in July to August. There were some differences, however, in the distribution during different hydrological years.

2. The change in annual mean rainfall erosivity was different for different time series, while the spatial distribution of rainfall erosivity showed an overall decrease from the southeast to the northwest. The province with the highest rainfall erosivity was Shanxi, while Ningxia had the lowest. Southern Shanxi, and central and southern Shaanxi were areas susceptible to rainfall erosion in the Loess Plateau.

3. Rainfall erosivity obtained from future climate scenarios is lower than that obtained from historical data. While the precipitation amount increased in future scenarios, the relationship between precipitation and rainfall erosivity was relatively discrete. Future rainfall erosivity presents a relatively regular spatial distribution pattern: high in the southeast and low in the northwest. Areas with intense rainfall erosion on the Loess Plateau were concentrated in the south (RCP4.5) and southeast (RCP8.5).

This research was financially supported by the National Natural Science Foundation of China (Grants 41530858 and 41771305) and the National Key Research and Development Program of China (2017YFC0504702). We thank Paul Seward, PhD, from Liwen Bianji, Edanz Group China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.

All relevant data are included in the paper or its Supplementary Information.