Statistical evaluation of proxies for estimating the rainfall erosivity factor

Considering the high-temporal-resolution rainfall data requirements for calculating the Rainfall Erosivity factor (that is, the R-factor), studies have developed a large number of proxies for the R-factor (PR). This study aims to evaluate 15 widely used proxies, which were developed in various countries using daily, monthly, or yearly rainfall data, in terms of correlation and statistical equality with the R-factor by using the 6-min pluviographic data from 28 stations in Australia. Meng’s test was applied to rank the correlations. Although the Meng’s test indicated that the correlation between Rainfall Erosivity (R) and Rainfall Erosivity calculated by the proxy model (PR) generally increased with a finer time resolution of the rainfall data (in the order of year, month, and day), the 15 PRs under examination were all highly correlated with R (r > 0.62, p < 0.004), implying that all of them can be reasonably used as an R predictor. A direct estimation of the R-factor using PRs produced a mean relative error (MRE), root mean square error (RMSE), and Nash–Sutcliffe efficiency coefficient (NSE) with a mean of 50.0%, 1392 MJ mm ha−1 h−1 a−1, and 0.17, respectively. The linear calibrations improved the accuracy of the estimation and produced an MRE, RMSE, and NSE with a mean of 36.0%, 887 MJ mm ha−1 h−1 a−1, and 0.70, respectively. Finally, suitable proxies for instances where only daily, monthly, or yearly rainfall data are available were recommended.

Soil erosion can lead to land degradation and environmental deterioration [1][2][3] . Eroded sediments also have negative off-site effects, such as impaired water quality, altered aquatic habitats, limited channel navigability, reduced hydroelectric equipment lifespan, and increased flood risk [4][5][6] . Rainfall is a driver of water erosion. Unlike other natural factors affecting erosion, such as topography and vegetation, rainfall is less subject to human alteration or control. Rainfall thus represents an environmental contributor to soil loss and is of great importance in assessing or predicting soil erosion intensity 7 .
Rainfall Erosivity refers to the potential of a storm to erode the soil. The universal soil loss equation (USLE) 8,9 and the revised universal soil loss equation (RUSLE, RUSLE2) 10,11 , the most widely used models for predicting soil erosion, applied equations (Eqs. 1 and 2) to quantify the average annual Rainfall Erosivity, R (MJ mm ha −1 h −1 a −1 ). The discrepancy is reflected in the algorithm of the estimated unit rainfall kinetic energy (e r , MJ mm −1 ha −1 ), and Eqs. (3)(4)(5)  where n is the number of years of the record, m is the number of erosive events in each year, EI 30 is the Rainfall Erosivity index of a single event, E is the total rainfall energy (MJ ha −1 ), I 30 is the maximum 30 min rainfall intensity (mm h −1 ), and P r and i r are the rainfall amount (mm) and rainfall intensity of the rth time interval (mm h −1 ), respectively. Comparing the three equations for the R-factor, the results indicated that the RUSLE significantly underestimates the R value [12][13][14][15] . Nearing et al. 14 stated that the R value calculated by the RUSLE was 14% lower than that of the USLE and RUSLE2; however, the USLE and RUSLE2 results showed resemblances based on data from 56 stations from three countries: China, Italy, and the US. Studies show that the RUSLE used to estimate the R-factor is based on a conceptual mistake that results in underestimation of the R value; hence, Eq. (4) (RUSLE) should not be used for R calculation and Eq. (5) (RUSLE2) is recommended to be used instead 14 . The RUSLE2 is commonly utilized by the US government for conservation planning 16 and has been used in recent studies on Rainfall Erosivity 17,18 .
The EI 30 index is a standard algorithm for calculating Rainfall Erosivity. However, the algorithm requires pluviographic records at < 15 min intervals 9 ; moreover, the record period should be longer than 20 years for the purpose of containing dry and wet climatic regimes 19 . Such rainfall data are often unavailable, particularly in developing countries 20 . In addition, the algorithm is complicated and often laborious and time-consuming 21 .
Considering the limited data availability and tedious process, many proxies for the R-factor (PR) have been developed, using daily, monthly, and annual data instead of pluviographic data. Table 1 lists 15 proxies that have been extensively applied in the literature, three of which have been developed in the US, two in Spain, two in Malaysia, five in China, two in Australia, and one in Honduras.
To assess candidate proxies, studies have often used the Nash-Sutcliffe efficiency coefficient 22 and the relative error to evaluate the discrepancies between the R-factor and the proxies 7,23-27 . However, the most suitable proxy was determined simply through a direct comparison of the obtained Nash-Sutcliffe coefficient and the relative error, which vary with the collected sample data and are random variables. Because the comparison among random variables should be done in the framework of statistics, a direct comparison of the obtained Nash-Sutcliffe coefficient or relative error is questionable.
To the best of our knowledge, no studies have conducted statistical evaluations of the established PR. The objective of this study was to statistically evaluate the 15 R-proxies listed in Table 1 using rainfall data from Australia. We ranked the correlations between the proxies and R-factor using Meng's test. Based on the statistical evaluations, we recommended proxies for instances where only daily, monthly, or yearly rainfall data are available.

Methodology
Data source. We used 6-min pluviographic rainfall data from 28 stations in Australia (Table 2), which were provided by the Australian Government Bureau of Meteorology (http:// www. bom. gov. au). Australia has a diverse climate, including tropical monsoon climate in the northern part, Mediterranean climate in the southwestern corner of the country, oceanic and humid subtropical climate in the southeastern part, and arid to semiarid climate in the interior mainland. Rainfall has high spatial variation, ranging from < 50 mm to > 3000 mm, with a mean of 465 mm ( Table 2). Detailed information on the selected stations is listed in Table 2, including location information, length of record, mean annual precipitation (MAP), and R value. The selected stations spanned from tropical to subtropical and temperate, and from humid to semi-arid and arid regions, with a MAP varying from 195 to 1236.2 mm. Most of the selected stations (22 out of the 28 stations) had a MAP ranging between 400 and 800 mm because soil erosion is rather limited in very humid or very dry areas. The missing rainfall data for each station varied, and for the purpose of the quality of the rainfall data, the maximum missing rate of each year should not be greater than 90% for each station. Hence, the data for the selected period were unequal. Rainfall data from the 28 stations were for the period 1960-2011, and the data periods ranged from 20 to 31 years, with a mean of 23 years for the 28 selected stations.
Calculation of the R-factor using the RUSLE2 and its proxies. The R value is generally underestimated using the RUSLE2, as discussed in the Introduction section. Hence, the R-factor was calculated using Eqs.
(1), (2), and (5) (RUSLE2), and the 15 PRs (the proxies for the R-factor; Established proxies for the R-factor. PR i is the ith proxy for the R-factor. AR i , MR j , and HMR k (MJ mm ha −1 h −1 a −1 ) are the Rainfall Erosivity in the ith year, jth month, and kth half-month, respectively. p j , p a , p i , and p are the rainfall in month j, annual rainfall, average rainfall in month i, and highest monthly rainfall, respectively. rain 10 (mm) represents monthly rainfall exceeding 10 mm for days, and day 10 represents the monthly rainfall exceeding 10 mm for days. (P d > 12.7(12, 0) ) m (mm) represents daily rainfall of ≥ 12.7 (12, 0) mm on the mth day. P dmax (mm) represents the maximum daily rainfall in an average year. H (mm) and L are the altitude and latitude in decimal degrees of the stations, respectively. S (mm) is the average rainfall for the summer half of the year (November-April, Bureau of Meteorology, 1989). MAP (mm) is the mean annual precipitation. α, β, and η are regional constants. F and MFI represent the Fournier index 48  Daily 6.97(rain 10 ) j − 11.23 day 10 j Malaysia 32 1.7265 m α = 0.3937 warm season (10-12months, 1-4months) 0.3101coldseason(5−9months)  29 designed a Z-test to compare two correlations (r 1 and r 2 ) that have a common dependent variable. The equations are as follows: where z ri is the Fisher z-transformed value for r i (i = 1 or 2), z ri = 1 2 ln 1+r i 1−r i z ri = 1 2 ln 1+r i 1−r i , N is the sample size, r x is the correlation between the two independent variables, r 2 i = (r 1 2 + r 2 2 ) / 2, and f should be set to 1 if it is greater than 1.
To compare k correlations with a common dependent variable (r 1 , r 2 ,…, r i ,…, r k ; k > 2), Meng's test becomes www.nature.com/scientificreports/ where r z r represents the correlation coefficient between z ri and λ i , z ri is the Fisher z-transformed value for r i , and λ i is the contrast weight assigned to each z ri and sums to zero; for example, λ i should be − 3, 1, 1, and 1 if r 1 is under examination and k = 4. χ 2 (k − 1) represents χ 2 distributed with k-1 degrees of freedom.
where N is the sample size, z r is the mean of z ri , r x is the median intercorrelation among the k independent variables, and r 2 i is the mean of r i 2 . By comparing r i and the average of the other k-1 correlation coefficients, Meng's test determines whether r i is a contrast in the two-tailed test. In the case of the one-tailed test, Meng's test determines whether r i is a less-or better-correlated contrast. This study used a stepwise Meng's test procedure to rank the correlations between R and PR i . We first labeled each of the correlations between R and PR as better-correlated contrast, less-correlated contrast, or non-contrast by means of the one-tailed Meng's Z-test, whereby the correlations under examination were grouped into two classes: the less-correlated class and the other class. Then, we repeated the procedure for each of the two classes, and the procedure did not terminate until no contrast was detected among any subclass. The significance level was set at 0.05. Bonferroni correction applies when several bivariate tests are performed simultaneously. We do not use Eq. (9) multiple times simultaneously, so that the Bonferroni correction is not necessary. In this way, we ranked the correlation between R and PR.
In addition to the correlation coefficient, we also used the MRE (%), the RMSE (MJ mm ha −1 h −1 a −1 ), and NSE to assess the performance of the regression equation. The MRE, RMSE, and NSE are defined as follows: where R i and PR i are the R and PR values of the ith station, respectively, n is the number of stations, and R(−−) is the average value of R i .

Results and discussion
Calculation results and R and PR deviations. The calculated R-factor of the 28 stations varied between 290.8 MJ mm ha −1 h −1 a −1 and 7153.7 MJ mm ha −1 h −1 a −1 with a mean of 2424.1 MJ mm ha −1 h −1 a −1 . As shown in Fig. 1, all examined proxies were surprisingly found to be highly correlated with the R-factor (r > 0.62, p < 0.004), implying that all of them can be reasonably used as a predictor of the R-factor in Australia. This finding is at odds with the common notion that regional differences in physiographic settings prevent a model of the R-factor from application to other areas. The reason might be that the PR models selected were established in different regions and developed by the R value calculated by the USLE or RUSLE method. The PR value is in fact an approximation of the R-factor calculation of the USLE or RUSLE. Therefore, the PR value is significantly related to the R value, and the difference between PR and R reflects the correction coefficient, which is the regional difference, and therefore all PR values were highly correlated with the R values.
The linear regression equation in Fig. 1 can be applied to modify the deviation between the PR and R. Table 3 lists the RMSE, MRE, and NSE for the estimation of R values that directly used PR and modified PR. The direct estimation showed high deviation, the MRE ranged from 18.3 to 153.0%, the RMSE ranged from 561 to 2279 MJ mm ha −1 h −1 a −1 , and the NSE ranged from − 2.1 to 0.89. The revised calculations, as shown in Table 3, indicated the MRE ranged from 13.1 to 77.4%, the RMSE ranged from 302 to 1482 MJ mm ha −1 h −1 a −1 , and the NSE ranged from 0.23 0.97. The RMSE (p = 0.002) and the MRE (p = 0.01) were notably reduced, and the NSE (p = 0.009) significantly increased under the correction, which further proved that the error was string-reduced when the modified PR was used. However, not all PRs showed similar expression, such as the RMSE, MRE, and NSE of PR 11 , PR 12 , PR 13 , and PR 14 , which changed slightly after calibration.

Correlation comparison between R and PR S using Meng's test.
According to the step-by-step test method designed in this study, the flow chart is shown in Fig. 2, Eq. (9) divided the correlation coefficients of 15 PRs and R into two classes: L 1 (r 1 , r 2 , r 6 , r 9 ), and O 1 (r 3 , r 4 , r 5 , r 7 , r 8 , r 10 , r 11 , r 12 , r 13 , r 14 , r 15 ), and no significant difference was observed in the L 1 class (p > 0.06). The correlation coefficients in the O 1 class were separated into L 2 (r 3 , r 4 , r 8 , r 10 ) and O 2 (r 5 , r 7 , r 11 , r 12 , r 13 , r 14 , r 15 ) classes, and no contrast was detected in the L 2 class (p > 0.35). The correlation coefficients in the O 2 class were divided into L 3 (r 5 , r 7 , r 12 ) and O 3 (r 11 , r 13 , r 14 , r 15 ) classes, and the coefficients in the L 3 class showed statistical equality under Meng's test. The correlation coefficients of PR 14 and PR 15 were notably greater than those of r 11 and r 13, and further statistical tests indicated no significant differences (9) Z = r z r χ 2 (k − 1) www.nature.com/scientificreports/ between r 11 and r 13, r 14 , and r 15 , respectively. After four rounds (Fig. 2), the stepwise Meng's test procedure ranked the correlation between PR i and R, r i , as follows: As expected, the correlation generally became stronger with finer time resolution of the rainfall data, and the daily proxies were stronger correlated with R than the monthly ones, which in turn were stronger correlated (14) r 1 = r 2 = r 6 = r 9 < r 3 = r 4 = r 8 = r 10 < r 5 = r 7 = r 12 < r 11 = r 13 < r 14 = r 15  (Table 3) showed a very similar rank to the correlation in Eq. (14).
Recommending the PR for rainfall data with different resolutions. Two yearly proxies, PR 1 and PR 2 showed strong correlation with R (Fig. 1a,b). The t-test indicated that the coefficient of PR 1 was equal to 1, and the constant term of PR 1 was, significantly, equal to 0. It showed that the R value in Australia can be directly estimated by PR 1 , and Table 3 shows that no striking variation was observed among the MRE, RMSE and NSE when estimating the R value by revised PR 1 . The coefficient and constant term of PR 2 were significantly unequal to 1 and 0, respectively. It showed that the estimation of the R value in Australia directly by PR 2 would result in a large discrepancy. Based on Eq. (14), no statistical difference was observed between r 1 and r 2 . Hence, we suggest the use of PR 1 and revised PR 2 as predictors of the R-factor when only yearly rainfall data are available. Among the monthly proxies, PR 4 , PR 5 , and PR 7 are based on the modified Fournier index (MFI), whereas PR 6 is based on the Fournier index (F) ( Table 1). Both indices have been widely applied to estimate the R-factor 31-34 . Renard et al. 10 and Yue et al. (2014) found that the MFI was stronger correlated with MAP than the F. Arnoldus 35 and Hernando and Romana 36 reported a stronger correlation between the MFI and R than between the F and R. The Australian data we used agree well with the observations above. Although the MFI and F were both highly correlated with MAP and R (Fig. 3), Meng's test indicated that the MFI was stronger correlated with both MAP and R (p < 0.048). The three MFI-based monthly proxies, PR 4 , PR 5 , and PR 7 , even showed a stronger correlation with R than PR 6 (Eq. (14)), a daily proxy for the R-factor. Hence, we concluded that the MFI is superior to the F in developing a proxy for the R-factor. The t-test of the regression equation of PR 3 -PR 8 indicated that the coefficients of the equations were significantly unequal to 1. Table 3 shows that directly using the monthly PRs to estimate the R-factor produced large errors. Therefore, when only monthly rainfall data are available, it is better to use the revised PR 5 and PR 7 to predict the R-factor.
Among the daily proxies, PR 14 and PR 15 were developed in Australia, and their correlation with R ranked first among the 15 proxies. PR 11 and PR 13 , two Chinese proxies, ranked second in terms of their correlation with R. The reason for the superior performance of PR 14 and PR 15 may not simply be due to the fact that they were established using Australian data. Capturing the periodic variation within a year, the model developed by Yu and Rosewell 37 (the original version of PR 14 and PR 15 ) performed best among the 8 R-proxies in northeastern Spain 7 .
The t-test of the regression equation of PR 10 -PR 13 indicated that the coefficient and the constant were equal to 1 (p > 0.06) and 0 (p > 0.33), respectively. Table 3 illustrates that no remarkable improvement was observed in the deviation of the R value when estimated by the revised equations. Consequently, PR 10 -PR 13 can be used to estimate the R value. The coefficient and constant term of PR 9 were significantly unequal to 1 and 0, respectively. The revised PR 9 was recommended when estimating the R-factor in Australia. The coefficient of PR 14 was significantly equal to 1 (p = 0.87) and the constant was notably greater than 0 (p = 0.02), which indicated that the equation for PR 14 in Fig. 1k probably leads to a large deviation in estimating R values for dry areas. The constant term of PR 15 was equal to 0 (p = 0.45) and the coefficient was significantly greater than 1 (p = 0.00), which indicated  Fig. 1l. We suggest the use of revised PR 14 and PR 15 as surrogates of the R-factor when daily rainfall data are available. However, with MREs of 26.1% and 23.6%, respectively (Fig. 1i,j), the accuracy remains satisfactory while using PR 11 and PR 13 as a surrogate of the R-factor.

Conclusion
This study statistically evaluated 15 developed PR using rainfall data from Australia. All proxies were found to be significantly correlated with R (p < 0.004), and thus can be used to predict the R-factor by means of linear calibration models. Meng's test ranked the correlations between R and PR (see Eq. (14)), indicating that the daily proxies were stronger correlated with R than the monthly ones, which in turn were stronger correlated than the yearly ones. Furthermore, the t-test of the regression equations showed that five proxies (PR 1 , PR 10 , PR 11 , PR 12 , and PR 13 ) were statistically equal to the R-factor, implying that they can be directly used as an R surrogate without calibration. A direct estimation of the R-factor using PRs produced an MRE between 18.3 and 153.0%, with a mean of 50.0%. The linear calibrations produced an MRE ranging from 13.1 to 77.4%, with a mean of 36%. Based on the statistical evaluations above, we recommended proxies for instances where only daily, monthly, or yearly rainfall data are available. Our rainfall data were collected from varying climate zones (from tropical to temperate, and from humid to semi-arid and arid). Hence, although future work in other regions is needed, our results may be applicable to other regions around the world.

Figure 2.
Results of comparing correlations between R and PR by the stepwise Meng's test procedure. The numbers within parentheses are the p values resulting from the Meng's tests, the L class represents the lesscorrelated class, and the O class represents the classes other than those in the L class, including the non-contrast and the better-correlated contrast. Bold and red texts denote the better-than-average correlated.