Rainfall erosivity estimation based on rainfall data collected over a range of temporal resolutions

Introduction Conclusions References


Introduction
Soil erosion prediction models are effective tools for helping to guide and inform soil conservation planning and practice. The most widely used soil erosion models used for conserva-tion planning are derived from the Universal Soil Loss Equation (USLE) Smith, 1965, 1978). These models include the USLE, the Revised USLE (RUSLE) (Renard et al., 1997), and RUSLE2 (Foster, 2004). Adaptations of the USLE have also been developed for use in other parts of the world, including, for example, Germany (Schwertmann et al., 1990), Russia (Larionov, 1993), and China . For example, the Chinese Soil Loss Equation (CSLE) was used in the first national water erosion sample survey in China (Liu et al., 2013).
These models have in common a rainfall erosivity factor (R), which reflects the potential capability of rainfall to cause soil loss from hillslopes, and which is one of the most important basic factors for estimating soil erosion. In its simplest form, the R factor is an average annual value, calculated as a summation of event-based energy-intensity values (EI 30 ) for a location divided by the number of years over which the data were collected. EI 30 is defined as the product of kinetic energy of rainfall and the maximum contiguous 30 min rainfall intensity during the rainfall event. It is the basic rainfall erosivity index that was developed by Wischmeier (1958) originally for the USLE, and is still widely used in other erosion prediction models (e.g., RUSLE, RUSLE2), with some modifications and improvements. Wischmeier (1976) suggested that more than 20 years' rainfall data are needed to calculate average annual erosivity to include relatively dry and wet periods.
Determination of the maximum contiguous 30 min rainfall intensity during the rainfall event is a relatively straightforward process, although it requires a temporally detailed rainfall record (e.g., 5 min) for a storm. Determination of the kinetic energy of a storm is more complex.
Kinetic energy (KE) is generally suggested to indicate the ability of a raindrop to detach soil particles from a soil mass (e.g., Nearing and Bradford, 1985). Since the direct measurement of KE requires sophisticated and costly instruments, several different estimating methods have been developed that estimate KE based on rainfall intensity (I) using logarithmic, exponential, or power functions. The original 1978 release of the USLE utilized a logarithmic function (Wischmeier and Smith, 1978) that was based on rainfall energy data published by Laws and Parsons (1943). Brown and Foster (1987) re-evaluated this relationship and recommended the use of an exponential relationship, which was subsequently used in RUSLE (Renard et al., 1997). For RUSLE2 Foster (2004) used the exponent value of −0.082, rather than the −0.05 value used in RUSLE, as follows: where e r is the estimated unit rainfall kinetic energy (MJ ha −1 mm −1 ) and i r is the rainfall intensity (mm h −1 ) at any given time within a rainfall event (usually taken as 1 min for computational purposes, with average-intensity representative of the time increment). This was based largely on work of McGregor and Mutchler (1976) and McGregor et al. (1995), who found that the RUSLE equation gave values that were too low. The energy term used in RUSLE2 gives results on the order of those from the original USLE method. The temporal resolution of rainfall data available across the world does not always allow for a direct computation of rainfall kinetic energy (Sadeghi et al., 2011;Sadeghi and Tavangar, 2015;Oliveira et al., 2012;Panagos et al., 2015;Zhang and Fu, 2003), even within countries with extensive rainfall monitoring programs. In the United States, for example, intra-storm, temporally detailed data (historically taken on pen recording charts, now taken as 1 min digital data) are only available at limited stations, whereas daily data are common (Nicks and Lane, 1995;Flanagan et al., 2001). Yet there is a need for developing models for application in all areas of the world in order to produce erosivity maps that can be used for evaluating soil erosion rates (e.g., Sadeghi et al., 2011, Sadeghi andTavangar, 2015;Oliveira et al., 2012;Panagos et al., 2015;Zhang and Fu, 2003). For that reason many efforts have been undertaken to estimate rainfall erosivity by using daily (Richardson et al., 1983;Yu, 1998;Capolongo et al., 2008;Yin et al., 2007;Zhang et al., 2002a, b;Xie et al., 2001Xie et al., , 2015, monthly (Arnoldus, 1977;Renard and Freimund, 1994;Yu and Rosewell, 1996;Ferro et al., 1999;Wu, 1994;Zhou et al., 1995), or annual rainfall data (Lo et al., 1985;Renard and Freimund, 1994;Yu and Rosewell, 1996;Bonilla and Vidal, 2011;Zhang and Fu, 2003;Wang, 1987;Sun, 1990). Generally the technique has been to develop a simple empirical relationship between erosivity and coarse resolution rainfall based on limited finer resolution data and then to extend the analyses to wider areas and longer periods with coarser temporal resolution rainfall data (Angulo-Martinez and Begueria, 2012;Ma et al., 2014;Ramos and Duran, 2014;Sanchez-Moreno et al., 2014).
Several studies evaluated different timescales of erosivity using different temporal resolutions of rainfall data. In Europe, Panagos et al. (2015) undertook the task to develop an erosivity map for Europe based on data from 1541 precipitation stations with temporal resolutions of 5 to 60 min. To use data that had been reported at the different time resolutions they had to apply adjustment factors to the data, which they reported to have introduced some uncertainty into the estimations. Sadeghi and Tavangar (2015) evaluated various erosivity estimation indices, including Fournier (Fournier, 1960), modified Fournier (Arnoldus, 1977), Roose (1977) and Lo et al. (1985), using data from 14 stations in Iran. They evaluated annual, seasonal, and monthly information. Similarly, the work in Brazil summarized by Oliveira et al. (2012) highlighted several studies that used various estimations of erosivity based on various types of data and interpolations. Other innovative ways have been advanced to produce better mappings of erosivity, including the use of daily (Fan et al., 2013) or 3 h (Vrieling et al., 2010(Vrieling et al., , 2014 data from the Tropical Rainfall Measuring Mission (TRMM) Multi-satellite Precipitation Analysis (TMPA) precipitation data.
In China the specifications for surface meteorological observations by the China Meteorological Administration (China Meteorological Administration, 2003) have required since the 1950s that the maximum 60 and 10 min rainfall amounts, (P 60 ) day and (P 10 ) day , be compiled; hence, these data are readily available in China. The measurements were made using siphon-method, self-recording rain gauges. Because of this, there is an interest in China to utilize the maximum daily 10 and 60 min rainfall intensities, (I 10 ) day and (I 60 ) day , to calculate erosivity.
The objectives of this study were threefold: (1) calibrate methods of estimating erosivity for timescales ranging from daily to average annual based on different temporal resolutions of rainfall data from 11 calibration stations with 1 min resolution data; (2) compare models in this study with those published in previous research, based on seven independent validation stations using the same data types; and (3) determine the most accurate methods, based on these data, for calculating different timescales of erosivity when different temporal resolutions of rainfall data are available. Note that, in this paper, we use the term timescales when discussing the erosivity values (equation outputs) and resolution (equation inputs) when referring to the rainfall input data, for clarity. Although several studies have been conducted on this topic in the past, no study used as comprehensive a data set collected over this wide geographic area of China to evaluate the wide range of erosivity timescales needed for erosion work, and utilizing such a wide range of temporal resolution rainfall data as the independent variable.

Data
Data collected at 18 stations by the Meteorological Bureaus of Heilongjiang, Shanxi, Shaanxi, Sichuan, Hubei, Fujian, and Yunnan provinces and the municipality of Beijing were used (Fig. 1, Table 1). These stations were distributed over the eastern half of China; 1 min resolution rainfall data (Data M) were obtained by using a siphon, self-recording rain gauge. The data collection period began in 1971 for Wuzhai (53663) and Yangcheng (53975) in Shanxi Province and from 1961 for the remaining 16 stations; the data records ended in 2000 for all stations. Quality control of Data M was done to select the best observation years using the more complete data sets of daily rainfall totals, Data D, which were observed by simple rain gauges at the same stations. Data M was compared with Data D on a day-by-day basis, and those days with deviation exceeding a certain criterion were marked as questionable and were not used in this analysis (Wang et al., 2004). The criterion used was that the data were considered good when the absolute deviation between Data M and Data D was less than 0.5 mm when the daily rainfall amount was less than 5 mm and no more than 10 % when the daily rainfall amount was greater than or equal to 5 mm. Data M in the earlier years of record tended to have more days with missing or suspicious observations. These totals of Data M and Data D were compared year-by-year to determine which years could be designated as common years for use in this study, with an effective year having a relative deviation for yearly rainfall amount of no more than 15 %. There were at least 29 common years for all 18 stations, and seven stations had common years of at least 38 years ( Table 1). Note that though there were missing data in the information used, Data D was only used for quality control purposes and the data used in the analysis, Data M, were internally consistent in that only the data from common years were used in all comparisons and evaluations reported. Data M was used to calculate the event-based EI 30 values as a function of the calculated kinetic energy and maximum 30 min rainfall intensity (Foster, 2004). These were treated as observed values and summed to obtain the erosivity factors, R, for daily, month (individual month totals), year (individual year totals), average monthly (one value for each month at each station), and average annual (one value for each station) timescales. Total rainfall event depth values were also compiled into the other temporal resolutions of rainfall data, including correspondent daily, month, year, average monthly, and average annual resolutions. For the eight stations in the northern part of China (including stations in Heilongjiang, Shanxi, Shaanxi provinces and Beijing municipality), only the periods from May through September were used because the siphon, self-recording rain gauges were not utilized in the winter to avoid freeze damage. Percentages of precipitation during May through September to total annual precipitation varied from 75.6 to 89.2 % for these eight northern stations. Data M for the full 12 month year was used from the remaining 10 stations located in the southern parts of China.

Calculation of the R factor at different timescales
Different timescales for RUSLE2 erosivity, R, including event, daily, month, year, average monthly, and average annual, were calculated based on the 1 min resolution data (Data M). Recall that month and year refer to individual months and years, and not averages. EI 30 (MJ mm ha −1 h −1 ) is the rainfall erosivity index for a rainfall event, where E is the total rainfall kinetic energy during an event and I 30 is the maximum contiguous 30 min intensity during an event (Wischmeier and Smith, 1978). An individual rainfall event was defined as a period of rainfall with at least six preceding and six succeeding non-precipitation hours (Wischmeier and Smith, 1978). An erosive rainfall event was defined as one with rainfall amounts greater than or equal to 12 mm, following Xie et al. (2002). We used the equation recommended by Foster 2004) for RUSLE2 to calculate the kinetic energy of the storms, which used Eq. (1) combined with where e r is the estimated unit rainfall kinetic energy (from Eq. 1) for the rth minute (MJ ha −1 mm −1 ); P r is the 1 min rainfall amount for the rth minute (mm); r = 1, 2, . . ., n represents each 1 min interval in the storm; and i r is the rainfall intensity for the rth minute (mm h −1 ). The Foster (2004) equations were chosen because they are currently used for erosion assessment for RUSLE2 in the United States and for the CSLE in China, and it appears to give results similar to the original USLE, as was discussed in the Introduction. Our evaluation included four models for events and one for daily erosivities. Event models were simply models to predict individual event erosivities, regardless of whether they occurred in 1 or more days, and regardless of whether more than one event occurred in a day. For the daily model, rainfall erosivity for each day, R day , was calculated following the method by Xie et al. (2015). When a day had only one erosive event and this event began and finished during the same day, then R day = EI 30 .
(3) Figure 1. Locations of the 18 stations with 1 min resolution rainfall data. Eleven stations marked with dots were used to calibrate 21 models. The other seven stations marked with triangles were used to validate models and conduct comparisons with previous research.
When more than one full rainfall event happened during 1 day, then where n is the number of rainfall events during the day, and E event_i and (I 30 ) event_i are the total rainfall energy and the maximum contiguous 30 min intensity, respectively, for the ith event. When only one part of a rainfall event occurred during 1 day, then where E day_d is the rainfall energy generated by the part of rainfall occurred during the dth day and (I 30 ) event is the maximum contiguous 30 min intensity for the entire event. The remaining situations were calculated by combining Eqs. (4) and (5).
Month, year, average monthly, and average annual R values were summed from the event EI 30 index by erosive storms that occurred during the corresponding period. They were calculated by using Eqs. (6)-(9).
where y is the number of years of record; (EI 30 ) y, m, j is the EI 30 value for the j th event in the mth month of the yth year; R month, y, m is the R value for the mth month of the yth year; R ave_month, m is the average R value for the mth month over the years of record; R year, y is R value in the yth year; and R ave_annual represents average annual erosivity, correspondent to the annual average R factor in USLE-type models (MJ mm ha −1 h −1 a −1 ).

Model calibration using different resolutions of rainfall data
A total of 21 models were calibrated for different timescales of R, based on varying resolutions of rainfall data (Table 2). Event amount P event and peak-intensity indices were derived based on the 1 min resolution data, including I 10 , I 30 , and I 60 , which were the maximum contiguous 10, 30, and 60 min intensities, respectively, within an event. I 10 and I 60 were used because of their close correlation with the daily (I 10 ) day and (I 60 ) day values commonly reported by the Chinese Meteorological Administration (2003). Four event-based models were developed relating measured EI 30 to estimated EI 30 (Table 2). Similar models for the other timescales were also calibrated ( Table 2). Data were organized in various ways. P day , P month , P year , P ave_month , and P annual were the daily, (individual) month, (individual) year, average monthly, and average annual rainfall amounts, respectively, for a given station. (P 60 ) month and (P 60 ) year represented maximum contiguous 60 min rainfall amount observed within a specific month or year, respectively. (P 60 ) month_max represented the maximum of (P 60 ) month values for each month of the year over the entire period of record. The average of (P 60 ) month values was(P 60 ) month . Each station had 12 values of (P 60 ) month_max and (P 60 ) month , one for each month of the year. (P 60 ) year_max was the maximum value of (P 60 ) year and (P 60 ) annual was the average of (P 60 ) year values. Each station had only one value for these two parameters. P 1440 was daily rainfall amount and its related index, including (P 1440 ) month , (P 1440 ) year , (P 1440 ) month_max , (P 1440 ) month , (P 1440 ) year_max , and (P 1440 ) annual , which were defined in an analogous way as were correspondent values for P 60 . The parameters were obtained station-by-station for calibration stations first and parameters for linear relationships were compared to determine if data from all stations could be pooled together to conduct the regressions (Snedecor and Cochran, 1989). Parameters for power-law models, including Month I, Year I, Average Monthly I, and Annual I (Table 2), were obtained by using the Levenberg-Marquardt algorithm (Seber and Wild, 2003). Note that models coded as Annual refer to annual averages.

Assessment of the models
After the 21 models in Table 2 were calibrated with the data from the 11 calibration stations, the performance for these models was assessed and compared with the performance of the previously published models listed in Table 3 using data from the seven validation stations. Symmetric mean absolute percentage error (MAPE sym ) and the Nash-Sutcliffe model efficiency coefficient (ME) were utilized to reflect the deviation of the calculated values from the observation data. MAPE sym is considered to be superior to MAPE, since it corrects the problem of MAPE's asymmetry and the possible influence by outliers (Makridakis and Hibon, 1995). MAPE sym was calculated as follows (Armstrong, 1985): where R obs is the measured rainfall erosivity for the kth period of time, such as month, year, or annual, based on 1 min resolution rainfall data. R sim is the estimated value for the same period using equations in Tables 2 or 3. ME was calculated as follows (Nash and Sutcliffe, 1970): ME compares the measured values to a perfect fit (1:1 line). Hence, ME is a combined measure of linearity, bias, and relative differences between the measured and predicted values. The maximum possible value for ME is 1. The greater the value the better the model fit. An efficiency of ME < 0 indicates the single value (the mean) for the measured data's mean is a better predictor of the data than the model. MAPE sym and ME were calculated based on all the data for the seven validation stations. Individual values for all stations were also determined.

Basic data results
Average annual rainfall ranged from 449.7 to 1728.1 mm, and average annual erosivity varied from 781.9 to 8258.5 MJ mm ha −1 h −1 yr −1 (Table 1). A total of 11 801 erosive events were used in the study. The eleven stations had 6376 erosive events, which were used to calibrate the models, and the seven validation stations had 5425 erosive events.

Validation and calibration for the new models
Parameters, MAPE sym , ME, and coefficients of determination, R 2 , for calibration models are shown in Table 4. The model Event IV, with a combination of event rainfall amount P event and I 30 , when I 30 was divided into two categories, with a threshold of 15 mm h −1 , performed slightly better in terms of the MAPE sym value than did Event II, which used the same variables but did not separate the rainfall events by intensity. The performance of Daily I with daily rainfall amount and (I 10 ) daily was similar to that for Event I with event rainfall amount and I 10 .
Using only total rainfall amount as input, the models for month, year, and average monthly scales were statistically significant, with determination coefficients R 2 greater than 0.66 (Table 4 and Fig. 2). However, their capabilities in predicting erosivity were limited based on the ME values (Table 4). Data from Tengchong and Xichang, located in the southwestern part of China, were in part responsible for these low ME values. Table 5 shows the individual values of MAPE sym and ME for the seven validation stations, with the average of each using all the stations and using only the five without Tengchong and Xichang. Results were much better without those two stations. The model Annual I, which use only average annual precipitation values, performed reasonably well, considering that the only input required was annual average precipitation (Table 4). If other information is available, other models performed better, but Annual I may be used if only average annual precipitation is available at a location.
In general, we found that the finer the temporal resolution of the rainfall input data, the better the models performed for a given erosivity timescale. Models that used some expression of maximum daily rainfall amount (Month III, Year III, Average Monthly III, Average Monthly V, Annual III, and Annual Model V) predicted the R factor better than those models with only total rainfall amount as input (Table 4), for a specific timescale. Models based on rainfall amount and maximum contiguous 60 min rainfall amounts (Month II, Year II, Average Monthly II, Average Monthly IV, Annual II, and Annual IV) generally performed better than corresponding models with rainfall amount and maximum daily rainfall amount (Month III, Year III, Average Monthly III, Average Monthly V, Annual III), except for Annual Model V, which performed well. The reason for that may be due to the fact that maximum contiguous 60 min rainfall amounts may have been more highly correlated with maximum contiguous 30 min intensity in an event as compared to just the maximum daily rainfall amount. The only annual average model that did not perform well was Annual III, which uti- 1 Parameters of models for power-law models, including α 1 , β 1 , α 2 , β 2 , α 3 , β 3 , α 4 , β 4 , α 5 , β 5 , were solved by pooling data from 11 stations together. Parameters for average annual-scale models, including λ 15 , λ 16 , λ 17 , λ 18 , were calculated by fitting data from all calibration stations and for the remainder they were the average values of parameters for the 11 calibration stations. 2 R 2 is the coefficient of determination.
lized (P 1440 ) year_max , the maximum of (P 1440 ) year values for each year over the entire period of record. Tables 3 and 4 show the models only evaluated for the erosivity temporal scale that corresponds to the input data resolution. For example, the event-based models are only evaluated on the basis of events modeled. We also evaluated the models at the aggregate scale. For example, EI 30 estimated from event-based models were summed up to month and year values, in order to evaluate if fine temporal resolution data also improve the accuracy of aggregate erosivity measures (Table 6). Two important facts emerge. First, when the models are applied at the aggregated scale their predictions get better. Secondly, the models that use fine resolution of input data predict better for the same erosivity timescale compared to models using coarser resolution input data. This has important implications for model applications.

Seasonal variations of erosivity
Taking Tonghe and Tengchong as examples, it was found that Month II generated better results than Month III, which performed better than Month I, in estimating seasonal and yearly variations (Figs. 3a, b and 4a, b). Correspondingly, seasonal variations by Average Monthly II were closer to observations as compared to those by Average Monthly III and Average Monthly I (Fig. 3c and d).
Year II and Year III produced better simulations of yearly variations compared with Year I, especially for the Tengchong station (Fig. 4c, d).
Hydrol. Earth Syst. Sci., 19, 4113-4126  Seasonal variations by monthly and average monthly models (Fig. 3) and yearly variations by month and year models (Fig. 4) were demonstrated using Tonghe and Tengchong stations. Month I and Average Monthly I captured the general seasonal pattern for the Tonghe station ( Fig. 3a and c), but the simulated peak value of monthly R was in July for the Tengchong station, which was not consistent with observation.
Month I and Year I captured the general year-to-year pattern for the Tonghe station ( Fig. 4a and c), but they overestimated yearly erosivity for the Tengchong station ( Fig. 4b and d). Month I and Year I also overestimated the yearly erosivity for the Xichang station. The reason for the overestimation for the Tengchong and Xichang stations was mainly due to two aspects: (1) the percentages of erosive rainfall amount to to-  tal rainfall at those stations were lower (71.9 and 76.9 %, respectively), suggesting that more events occurred with small amount totals that do not generate soil loss (Table 5); and (2) the ratio for event EI 30 to event rainfall amount P was lower (3.6 and 4.1, respectively), inferring that rainfall intensity and erosivity generated by per amount of rainfall were both less than that of the other stations (Table 5). This result was consistent with that of Nel et al. (2013), which demonstrated that two models using annual average rainfall and average monthly rainfall substantially overestimated annual erosivity on the west coast and the central plateau of Mauritius, which also have a large amount of non-erosive rainfall. Rainfall erosivity reflected a combined effect of rainfall amount and rainfall intensity. Therefore, it was reasonable that rainfall amount only explained part of rainfall erosivity variation at these stations.

Evaluation of models from previous research with current models
Generally speaking, the finer the resolution of input data for models, the better was the performance of the model for estimating at the same temporal erosivity scale. For example, the models with daily rainfall amount and daily maximum 60 or 10 min amount as inputs performed better than models with only daily rainfall amount as input. Similarly, results from models with maximum 60 min rainfall amount (Month II, Year II, Average Monthly IV, and Annual IV) were generally better than those with maximum daily rainfall amount (Month III, Year III, Average Monthly V, and Annual V; Fig. 5). Wang et al. (1995) used a combination of event rainfall amount P event and I 10 for event-scale models. The model using the I 10 data was divided into two categories, with a threshold of 10 mm h −1 , performed best among the four models compared (Table 3). That model had similar performance with Event IV in this study (Table 4), which also divided the data by a rainfall-intensity threshold.
There were three kinds of daily-scale models, according to the number and type of inputs required. Two models used daily rainfall amount (Zhang et al., 2002b;Xie et al., 2015), two models used daily rainfall amount and daily maximum 10 min intensity (Xie et al., 2001 and Daily I), and one model used daily rainfall amount and daily maximum 60 min inten-sity (Xie et al., 2015). The model with daily rainfall amount as input in Xie et al. (2015) performed better than that of Zhang et al. (2002b) (Table 3). Daily I, which used daily rainfall amount and daily maximum 10 min intensity as inputs in this study, performed better than the model in Xie et al. (2001). Models with an additional daily 10 or 60 min intensity index performed better than those with only a total rainfall amount (Tables 3 and 4).
There were generally four groups of models for month, year, average monthly, and annual-scale models. The first group used linear regression (Sun et al., 1990;Wu, 1994;Zhou et al., 1995) or a power-law function (Zhang and Fu, 2003; Month I, Year I, Average Monthly I, and Annual I) with only rainfall amount as input, so that the data required were relatively easy to collect. Models by Sun et al. (1990), Wu (1994) and Zhou et al. (1995), when they were used to estimate the monthly scale of R, had MAPE sym values of 86.7, 60.2 and 67.3 % and ME of −0.63, 0.57 and 0.35, respectively (Table 3). When they were used to estimate annual scale of R, there was a tendency of underestimation, especially for the stations with larger erosivity (Fig. 5a, b). Four models by Zhang and Fu (2003) overestimated the R factor, with MAPE sym varying between 34.6 and 60.8 % and ME varying between −2.11 and 0.10 (Table 3, Fig. 5), which suggested the models' abilities were limited. Two models by Zhang and Fu (2003) using the modified Fournier index generated poorer results than the model by Zhang and Fu (2003) , (b) year, (c) average monthly, and (d) average annual models using 1 min resolution data for the seven independent validation stations. Month models included models in Wu (1994), Zhou et al. (1995), and Month I, II, and III from this study. Year models included models from Sun et al. (1990), Wang et al. (1995; the one with MAPE sym of 20.3 %), Zhang and Fu (2003), and Year I, II, and III from this study. Average monthly models included models from Average Monthly I, II, and III from this study. Average annual models included models from Wang et al. (1995; the one with MAPE sym of 13.2 %), Zhang and Fu (2003; the one with MAPE sym of 34.6 %), and Annual I, II, and III from this study. using average annual rainfall as input (Table 3), which was consistent with the findings of Yu and Rosewell (1996). The power-law models in this study, including Month I, Year I, Average Monthly I, and Annual I, tended to overestimate the R factor for the stations with larger erosivity (Fig. 5).
The second group of models (Wang et al., 1995; Month II, Year II, Average Monthly IV, Annual IV) used linear regression with rainfall amount (total rainfall or total rainfall with daily rainfall no less than 10 mm) and maximum 60 min rainfall as inputs. All these seven models generated statistically significant results, with MAPE sym for R with timescale intended for the model ranging from 13.2 to 36.0 % and ME from 0.79 to 0.91 (Tables 3 and 4; Fig. 5).
The third group used linear regression with rainfall amount and maximum daily rainfall as inputs (Month III, Year III, Average Monthly V, Annual V), which generated reasonable results (Table 4) and a slightly overestimated annual R (Fig. 5). Overall they did not perform as well as did the models in the second group (Table 4).
The fourth group (Wang et al., 1995) used a combination of three indices, including rainfall amount, maximum 60 min rainfall amount, and maximum daily rainfall amount as inputs and generated good simulation results; however, there was no improvement compared with the two models by Wang et al. (1995) in the second group (Table 3).

Applications and recommendations
The results of this study provide a multitude of options for dealing with the problem of variations in available temporal resolutions of rainfall data from across the world for developing erosivity maps and databases. We present a series of 21 potential equations for use in estimating erosivity at timescales from event to average annual using input data resolution ranging from maximum 10 min rainfall intensity to average annual rainfall amount. Of the 21 equations we can recommend the use of 17. Equations Month I, Year I, and Average Monthly I, which use only total rainfall amounts for the respective timescales, all had low ME values and poor prediction capability (Table 4). Annual III, which is a linear function of average annual rainfall and the maximum daily precipitation over the recording period, performed very poorly, with a negative ME value (Table 4).
We found that using finer resolution data input produced better predictions of erosivity at a given output timescale. An exception was for the event-based models, where using I 30 gave slightly better results than using I 60 or I 10 . However, we also found that using equations with the finest data resolution possible, and aggregating or summing results for finer erosivity timescales, gave the best results (Table 6). In other words, if one were interested in average annual erosivity, but had rainfall data available for using the Daily I model, then results are better using the Daily I model and summing results over the period of data record rather than using Annual I-V models. It is also evident that predictions of erosivity using Daily I improve as the timescale increases. In other words, the predictions of average annual erosivity calculated by summing the daily values from Daily I give a higher level of fit than when using Daily I to estimate daily erosivity (Table 6).
Models in this study performed better or similar to models from previous research given the same rainfall data inputs based on these independent validation data (Tables 4  and 5). Models from previous research had higher symmetric mean absolute percentage errors, MAPEsym, and lower Nash-Sutcliffe model efficiencies, ME, with the exception of models for event, year, and average annual timescales by Wang et al. (1995), which had similar MAPEsym and ME compared to the models in this study.
Much attention has been given to monitoring the erosion process and its controlling factors at various spatiotemporal scales (Poesen et al., 2003). Characteristics of topography and soils are usually relatively constant in the timescales of interest, whereas rainfall erosivity and vegetation vary greatly. Therefore, soil erosion monitoring work is often mainly focused on the dynamics of rainfall erosivity and vegetation rather than soil and topography (Vrieling et al, 2014). Different timescales of erosivity are required in areas with different resolutions of rainfall data availability. Models provided in this study have potential to play important roles in the soil erosion monitoring framework in terms of quantifying the temporal dynamics and changes in rainfall erosivity.