The Richardson daily rainfall erosivity estimation model modified by Zhang et al. (2002) was applied to calculate RE, which reads R h a l f   a   m o n t h = α ∑ k = 1 m ( P k ) β . (1)

In the abovementioned equation, R h a l f   a   m o n t h (ΜJ·mm·hm⁻²·h⁻¹·a⁻¹) is the RE value of half a month. The threshold value of 12 mm is consistent with the Chinese standard for rainfall erosivity (Xie et al., 2016). K is the number of days of RE within half a month. P k is the daily precipitation greater than or equal to 12 mm on the K day within half a month. The seasonal distribution of erosive forces reflected by semimonthly periods is consistent with the USLE and RUSLE models. α and β are model parameters and can be calculated as follows: α = 21.586 β − 7.1892 , (2) β = 0.8363 + 18.144 P d 12 + 24.455 P y 12 , (3) where P_d12 (mm) refers to the daily average precipitation greater than 12 mm and P_y12 (mm) refers to the average annual precipitation above or equal to 12 mm per day.

The RE of each half a month of each year is calculated using Eqs 1–3. The monthly rainfall, annual rainfall erosivity (ARE), and average annual rainfall erosivity (AARE) can be obtained by summarizing.

Trend analysis, mutation analysis, and periodic analysis were employed to analyze the changes in RE on the timescale. The Mann-Kendall method (M-K) was applied to calculate the Z-value of statistics, and the changing trend in RE was analyzed and discovered through the positive and negative Z-values (Kendall, 1948; Alashan 2020; Wang 2020). The M-K test has been widely applied to detect sudden changes in values in long time series of climate data. The climate data are large and span long periods of time. Therefore, we normally regard years as the basic unit for research and analysis (Sa’adi et al., 2019; Machiwal et al., 2019). If z is above 0, the sequence will show an upward trend, and if z is below 0, it will display a downward trend. For example, if |z| is greater than or equal to 1.96, the mutation points of the time series of RE will be identified by the standardized cumulative anomaly, and the rank test method will be applied to further accurately test the mutation points. When the statistic |U| is below 1.96, the mutation points identified by standardized cumulative anomaly will be excluded. Wavelet analysis was employed to analyze the periodic change of RE. We spatially interpolated the precipitation and RE using a semi-covariance function and co-kriging interpolation. As to the co-kriging interpolation, we introduced elevation data as an auxiliary variable and the resolution was 2 km (Xiao et al., 2018; Belkhiri et al., 2020). The cross-verification method is the use of points around each measured point to predict the real value, and it compares the predicted value with the real value (Sharma et al., 2022). The accuracy of the co-kriging interpolation results was assessed by mean relative error (MRE), mean absolute error (MAE), and root mean square error (RMSE), and the smaller the value is, the higher the accuracy becomes (Qi et al., 2020).

Topographic relief (m) refers to the height difference between the highest and lowest points within a particular area (Shiraishi et al., 2020). In this study, the mean–variable point methods were applied to calculate the optimal statistical unit (18 × 18) of topographic relief in the Qinba Mountains, which is the optimal statistical window for neighborhood analysis to obtain the topographic relief. T R = A L T / 1000 + { [ Max ( H ) − Min ( H ) ] × [ 1 − P ( A ) / A ] } / 500 , (4) where T R is the topographic relief; A L T is the average altitude in the region; Max ( H ) and Min ( H ) are the highest and lowest values of altitude in the region, respectively; P(A) is the area of flat land in the region, and areas with slopes of less than 2° are defined as flat in this article; A is the total area.

The geographical detector method developed by Wang and Xu (2017) was used to conduct factor detection and analysis of RE and calculate the degree of influence of various environmental factors on RE in the region. The calculation equation is as follows: q = 1 − 1 N δ 2 ∑ h = 1 L N h δ h 2 . (5) In this equation, q is the detection factor (rainfall, temperature, NDVI, DEM, longitude, and latitude) detection force value, and q ∈ [ 0,1 ] . The larger the q value becomes, more is the influence that factor A has on RE in the study area. N and N h are, respectively, the sample numbers of the Qinba Mountain region. δ h 2 is the discrete variance of A in sample h, and L is the type of the factors.

The AARE of the Qinba Mountains from 1970 to 2017 was 4,197.85 MJ mm hm⁻² h⁻¹ a⁻¹, showing a slight variation with a mutation coefficient of 0.16. As illustrated in Figure 2A, the overall trend in RE is fluctuating, displaying a three-stage change characteristic of increasing first, then decreasing, and finally increasing. In the first stage (1970–1984), the RE showed a rising trend of fluctuation, and it then gradually decreased in the second stage (1984–1995). During the third stage of 1995–2017, the RE fluctuated from 3,157.31 to 3,954.75 MJ mm hm⁻² h⁻¹ a⁻¹. According to Figure 3, the wavelet analysis was introduced to conduct periodic analysis of RE. It was found that RE presented a periodic change of 5–10 years.

The variation trend of annual precipitation (Figure 2B) is consistent with annual RE (Figure 2A). During the past 50 years, the maximum precipitation was 1,046.98 mm in 1983, and the maximum RE was 5,264.11 MJ mm hm⁻² h⁻¹a⁻¹ in 1980. The lowest precipitation was 560.76 mm in 1997, and the minimum RE was 3,157.31 MJ mm hm⁻² h⁻¹a⁻¹ in 1995. The results demonstrated that the maximum precipitation did not coincide with the occurrence time of the maximum RE, and the maximum value of RE was 3 years earlier, and the minimum value was 2 years earlier, compared with precipitation. The maximum and minimum precipitation values differed by 1.86 times, but the maximum and minimum values of RE differed by 1.66 times, meaning that the variation range of precipitation is more extensive than RE. The variation coefficient of annual RE (0.16) was higher than that of the annual precipitation (0.12), illustrating that the interannual change of RE was more significant than the interannual change of precipitation. Figure 3 indicates that the amplitude changes of UF and UB are similar, suggesting a fluctuation trend of first increasing, then decreasing, and again increasing. The positive series (UF) represented a fluctuation trend from 1970 to 1983, a slow decline from 1983 to 2003, and a gradual rise from 2004 to 2017. The reverse series (UB) showed a declining trend from 1970 to 1979, a rising trend from 1980 to 1991, and an ascending trend from 1992 to 2017. The change in UB was more intense. UF and UB, respectively, had an intersection point between 1988 and 1989, 2006 and 2007, and 2007 and 2008. The positive sequence (UF) exceeded the critical value of 1.96 in 1984, and the inverse sequence (UB) went beyond the critical value of −1.96 in 2009 because the mutation point in 2007 and 2008 was the closest to 2010 and displayed a trend of significant decline, which further indicates that RE mutations in the Qinba mountains had occurred in 2008.

Considering the seasonal differences of precipitation, spring, summer, autumn, and winter are correspondingly divided into March–May, June–August, September–November, and December–February. As shown in Figure 4, RE slumps in spring and summer and soars in autumn, and the changes in winter are small and can be ignored. From 1970 to 2017, in spring, the average RE was 243.20 MJ mm hm⁻² h⁻¹ a⁻¹, accounting for 5.79% of the annual average RE, and the average RE decreased by 9.16 MJ mm⁻² h⁻¹ a⁻¹ per 10a, with a Z-value of −2.60, showing a significant decreasing trend. In summer, the average RE was 1,660.73 MJ mm hm⁻² h⁻¹ a⁻¹, accounting for 39.56%, and the changing trend was first decreasing, then increasing, and again decreasing. In 1999, the maximum values of RE appeared in spring. The variation law of autumn annual RE was similar to that of interannual RE, which was 2,044.84 MJ mm hm⁻² h⁻¹ a⁻¹, accounting for 48.71% of the annual average RE. The average RE increased by 217.65 MJ mm hm⁻² h⁻¹ a⁻¹ per 10a, exhibiting a significant increasing trend (z = 1.36). In winter, since snowfall did not directly affect soil erosion, snow is not considered in the calculation of winter RE. As a result, the proportion of average RE was relatively small (5.94%).

Based on the ArcGIS 10.3 platform, semi-covariance functions and kriging interpolation were examined to interpolate the average annual rainfall erosivity (AARE) and average annual precipitation (AAP) in the Qinba Mountains from 1970 to 2017. The semi-covariance function shows that the AARE is a spherical model, and AAP is a Gaussian model. The accuracy of the kriging interpolation was verified through the cross-verification method, and the results showed that MAE and RMSE were less than 1, demonstrating the credibility of the AARE and AAP interpolation results. As illustrated in Figure 5A, the spatial distribution characteristics of AARE were high in the southeast and low in the northwest, which was consistent with the distribution rule of AAP (Figure 5B). The AAP gradually dropped from above the average precipitation line of 800 mm in the southeast to below the average precipitation line of 200 mm in the northwest, and the AARE dropped from above the average RE line of 9,000 MJ mm hm⁻² h⁻¹ a⁻¹ in the southeast to below the average RE line of 3,000 MJ mm hm⁻² h⁻¹ a⁻¹ in the northwest. The distribution center and the high value of the AARE were consistent and distributed in the southeast. The high-value distribution center of AAP was consistent with the AARE, which was distributed in the southeastern part. The low value was distributed in the northwestern part, but the distribution area of AARE was more widespread.

To further analyze the trend of the annual RE changes, the statistical Mann-Kendall Z-value and coefficient of variation (Cv) of 118 weather stations from 1970 to 2017 were calculated. Meanwhile, semi-covariance functions and kriging interpolation were applied for interpolation. The Mann-Kendall Z-value is more suitable for the linear model, and Cv fits the spherical model more. The cross-verification results display that the MAE is 0.001 and RMSE is 1.02 in the Mann-Kendall Z-value and that the MAE is 0.08 and RMSE is 0.14 in the Cv (Figure 6).

As Figure 6A demonstrates, there are pronounced regional differences in the distribution of the RE coefficient. The spatial distribution characteristics are high in the middle and low in the surrounding region, showing that in the central region, the RE was increased and that water conservation measures should be implemented to prevent further soil erosion in the central region (Muhire et al., 2015). RE was declining around the study area boundary, but most of them failed to pass the significance test (|Z| was below 1.96). Figure 6B also illustrates that the variation coefficient of RE was relatively stable (Cv was below 0.2), showing a distribution pattern of being high in the west and low in the east. Located in a high-altitude area, the terrain in the western region is complex, and the coefficient of variation was relatively large. Also, the annual rainfall was unstable, and the interannual RE varied greatly. The variation coefficient decreased from west to east in the Qinba Mountains, indicating that the precipitation tended to be stable from west to east, and the interannual variation of RE gradually decreased, mainly caused by the dual influence of topography and climate.

Four lines of latitude and longitude (32.67°N, 33.96°N, 107.05°E, and 108°E) were selected to analyze the trends in the AARE of the Qinba Mountains in the longitudinal and latitudinal directions. Figure 7 indicates a trend of increasing and then decreasing, with increasing longitude at different latitudes, and the trend line displays a sharper increase at 33.96°N than at 32.67°N. The highest value of AARE (9,085.31 MJ mm hm⁻² h⁻¹ a⁻¹) at 32.67°N appears at 108.41°E, while the highest value (5,750.45 MJ mm hm⁻² h⁻¹ a⁻¹) at 33.96°N appears at 107.96°E, with a difference of 1.57 times and a deviation of 0.45°E. Meanwhile, both 32.67°N and 33.96°N indicate that RE is more noticeable in areas with flatter topography, mainly due to the flatter topography being influenced by the monsoonal climate with abundant rainfall (northern hemisphere). The AARE tends to drop with increasing latitude at 107.05°E and 110.52°E (Figure 8). The highest value is 9,593.04 MJ mm hm⁻² h⁻¹ a⁻¹ and the lowest is 4,245.61 MJ mm hm⁻² h⁻¹ a⁻¹ at 107.05°E, with a decrease of 2.26 times. At 110.52°E, the highest (8,022.83 MJ mm hm⁻² h⁻¹ a⁻¹) and the lowest (4090.49 MJ mm hm⁻² h⁻¹ a⁻¹) values of RE are decreased by 1.96 times. Compared with 110.52°E, the decreasing trend of AARE is more apparent at 107.05°E.

Based on the research by Shiraishi et al. (2020) and considering the complex geomorphic features of the Qinba Mountains, the topographic relief (m) was divided into four categories: micro-topography (below or equal to 140), small topography (140–250), medium topography (250–370), and large topography (above or equal to 370). The topographic relief degree is mainly small topography (39.08%), distributed alternately with the middle relief area. The second is the medium topography (28.88%), distributed in the southwest and central parts. The third is the large topography, accounting for 10.37%, mainly distributed in the southwest with a high elevation and along the ridgeline from west to east like a belt. Micro-topography (21.67%) is mainly on the border of the central Hanzhong Basin and eastern border of the Qinba Mountains (Figure 9; Table 1).

Figures 10A–D display the corresponding AAP changes that exist in different topographic reliefs. The AAP in the micro-topography showed a distribution pattern of being high in the southeast and low in the northwest, with fewer areas of high AAP and low AAP concentrated in the northwest of the study area. The APP distribution pattern in the small topography is similar to that in the medium topography, decreasing from southeast to northwest. The main reason is that the Qinba Mountains pass through the Qinling–Huaihe Line, and the study area is a subtropical monsoon climate, resulting in high precipitation. Due to the small area of large topography in the topographic relief, the AAP was scattered with a decreasing trend from south to north. The AAP low-value area was concentrated in the southwest, while the AAP high-value area was not evident. The distribution of the AARE over different topographic reliefs is shown in Figures 11A–D. The RE in the micro-topography gradually decreases from the southeast to northwest, which is consistent with the variation of the AAP in the topographic relief. The sum of the RE distribution in the small and medium topographies accounts for 67.96%, with a decreasing trend from the east to west and from the south to north, mainly because of the influence of the monsoon climate on the eastern and southern areas and the superior hydrothermal conditions there. However, the northern and western areas are less affected by the monsoon climate, and the topography is complex, leading to a lower value of RE.