Development of daily rainfall erosivity model for Kelantan state, Peninsular Malaysia

Prevention of soil loss from erosion due to the impact of rainfall and changing rainfall patterns is one of the most important global issues in soil conservation. These changes in precipitation are largely located in the tropics and hence are probably associated with convection (Tan et al. 2015). The increased precipitation causes detachments of the individual soil particles and their transportation along the slope to rivers and reservoirs. Rainfall erosivity is a predominant factor referring to the kinetic energy of a raindrop's impact and the rate of associated runoff (Wischmeier & Smith 1978). Therefore, it is essential to obtain an accurate estimate of rainfall erosivity in tropical climates for the assessment of soil erosion risk.

The Universal Soil Loss Equation (USLE) devised by Wischmeier & Smith (1978) and the Revised Universal Soil Loss Equation (RUSLE) proposed by Renard et al. (1997) for the United States are widely used models to estimate annual soil loss by both interrill and rill erosion. These include six factors such as rainfall erosivity, soil erodibility, slope steepness, slope length, cover-management, and support practice. These factors are dynamic in nature which results in heterogeneous spatial patterns of soil loss. However, RUSLE requires regional values to be developed for each factor based on local data and conditions if used outside the United States (Wischmeier 1984).

The rainfall erosivity (R) factor is considered as the most important factor in estimation of soil loss due to its high temporal variability. The R factor was derived from more than 8,000 plot years by Wischmeier (1984) and can be quantified by the product of total kinetic energy of rainfall (E) and its peak 30 minute intensity (I₃₀) which computes all individual erosive storm events. An individual rainfall event was defined as a period of rainfall with at least six preceding and six succeeding non-precipitation hours (Xie et al. 2016).

Therefore, the main disadvantage to compute RUSLE is the requirement of high spatial and temporal (maximum 30 minutes) rainfall data series for the determination of the R factor, and these data are not generally available. However, kinetic energy can also be used for the calculation of the R factor and varies for different climates. The kinetic energy can be calculated for any region by using very high temporal resolution data (1 min, 5 min) but such data is generally not available. Alternatively, kinetic energy and intensity relationship can be used which need 30 min rainfall data to calculate the kinetic energy of an event. Salles et al. (2002) listed all the kinetic energy and intensity relationships developed for different locations which can be used in the absence of very high temporal resolution data for a particular location.

Apart from event based R factor calculation, other statistical models were also developed by researchers for different climates to calculate erosivity which use commonly available data, such as daily rainfall (Richardson et al. 1983; Yu & Rosewell 1996a; Petkovšek & Mikoš 2004; Angulo-Martínez & Beguería 2009; Ali 2015; Xie et al. 2016) and monthly rainfall (Renard & Freimund 1994; Yu & Rosewell 1996b; Ferro et al. 1999; de Santos Loureiro & de Azevedo Coutinho 2001; Mikoš et al. 2006; Diodato & Bellocchi 2007; Ochoa-Cueva et al. 2015). Xie et al. (2016) mentioned three aspects of the R factor that may be useful for soil erosion estimation: (1) average annual rainfall erosivity for predicting average annual soil loss, (2) seasonal distribution curve of rainfall erosivity and (3) event or daily rainfall erosivity. They also stressed the importance of event or daily rainfall erosivity in their related past studies and developed statistical models for the successful estimation of erosion index EI₃₀ from the daily rainfall amounts. Due to the lack of long-term event-based rainfall data, these statistical models are very useful and widely used to calculate daily erosivity. However, both event and daily rainfall amounts are not similar (Bullock et al. 1990) because daily rainfall amount includes only one event, multiple events, or only part of an event (Richardson et al. 1983; Xie et al. 2016).

Angulo-Martínez & Beguería (2009) estimated the R factor for the Mediterranean climate by using five daily models, viz. the exponential model by weighted least squares (Richardson et al. 1983), the Yu and Rosewell model, the modified Yu and Rosewell model and three monthly models, viz. precipitation intensity indices, the modified Fournier index and the F index (Ferro et al. 1999) in which they found that the Yu and Rosewell model for daily and precipitation intensity indices for a monthly R factor gave the best results as compared to the other models. Xie et al. (2016) also found superiority of the Yu and Rosewell over the Richardson exponential model. All these models have some regional factors which vary for different climates and these must be determined accurately to develop an accurate rainfall erosivity model.

The tropical climate has high precipitation throughout the year. The orographic precipitation is also well-known and has been identified and studied across the world (Karnieli & Osborn 1988; Goldreich 1994; Michaud et al. 1995; Al-Ahmadi & Al-Ahmadi 2013). Therefore, the regional factor values will definitely vary in high precipitation and elevation regions. An attempt was made by Yu et al. (2001) to estimate the R factor for tropical climates by using event-based data for the Yu and Rosewell model but the disadvantages of this model is that the limited rainfall stations and less temporal data (2 years only) were used which may not be accurate for the region (Yu & Rosewell 1996a). Other studies also use event-based data to calculate soil erosion by RUSLE in tropical climates, such as Shamshad et al. (2008), Leow et al. (2011) and Kamaludin et al. (2013), while some estimate rainfall erosivity by including orographic effect such as for the Mediterranean climate (Diodato & Bellocchi 2007) and the semi-arid climate (Nearing et al. 2015).

Previous research focused on event-based, daily and monthly precipitation data to calculate erosivity for climates other than tropical. Limited attempts were made to calculate regional factors for the tropical climate models which are suitable for commonly available data such as daily precipitation. The purpose of this study is to develop a daily rainfall erosivity model for a region of tropical climate by using 30 minute rainfall intensity data which will be helpful to calculate accurate average annual rainfall erosivity using RUSLE in the absence of high resolution data. Therefore, to achieve this goal, the study was conducted with the following steps: (i) use of high spatial and temporal data (30 min) data to calculate the R factor, (ii) to calculate monthly regional factors from power law relationship for Kelantan state, Peninsular Malaysia, and finally (iii) to develop a monthly power law relationship between erosivity and daily rainfall amount for the region.

The study was conducted in Kelantan state, in the north eastern part of Peninsular Malaysia between latitudes 4°33′ and 6°14′ north, and longitudes 101°19′ and 102°39′ east with an area of approximately 15,000 km². The northern part is bounded by the South China Sea and higher elevations are on the south, south-east and south-west sides (Figure 1).

The elevation ranges from 0 to the highest point of 2,187 m named as Mt. Tahan. The area is divided into upstream (Gua Musang), midstream (Kuala Krai, Jeli, Tanah Merah and Machang) and downstream (Pasir Mass, Pasir Puteh, Kota Bharu, Bachok and Tumpat), as shown in Figure 1. The climate is tropical, humid with average temperature ranges from 20 to 30°C. The period from November to January receives maximum rainfall, while June and July are the driest months. According to the Malaysian Meteorological Department, the wind strength is generally light and variable, however, there are yearly periodic changes in the wind flow patterns and based on this, the season is divided into four, namely, the southwest monsoon, northeast monsoon and two shorter periods of inter-monsoon seasons. The state is affected by the north-eastern monsoon wind accompanied by heavy precipitation beginning early November and ending in March. The easterly or north-easterly winds of 10 to 20 knots can sometimes reach 30 knots or more as a result of the strong surges of cold air from the north. Masseran & Razali (2016) mentioned that the winds which blow across the South China Sea bring heavy rain, especially to the eastern and southern coasts of Peninsular Malaysia as well as the central part of the Titiwangsa Ranges. This rain causes flooding, especially in November and December, in the east coast states. Average annual rainfall of the area is 3,017 mm while average daily annual wind speed is 1.50 m/s. The average monthly rainfall from November to January is 418 mm with a maximum of 576 mm in December.

Thirty minute intensity data from 55 stations were obtained from the Department of Irrigation and Drainage, Malaysia. Distribution of the number of rain gauges in upstream, midstream and downstream are 16, 17 and 22, respectively. Out of the 55 stations, 25 are relatively new, which have the data recorded from 2005 onwards. Some records of the precipitation data have been replaced by the nearby station data due to non-availability of the original data of the reporting station because of system failure at that station. Thirty minute records, measured from the current day at 12:00 am up to the next day at 12:00 am, for each station from January 2005 to December 2015, were selected.

As per the rule of RUSLE an event is considered erosive if the cumulative rainfall is greater than 12.7 mm. The precipitation limit 0.1 mm, 5.0 mm and 12.7 mm were taken into account in this study to estimate the R factor because precipitation is recorded almost throughout the year in tropical climates, hence, every small event may be playing an important role in soil erosion by keeping the soil moist. For cross-validation, a 70:30 ratio of calibration (40) and validation (15) stations were randomly selected in such a way as to represent the entire study area. Validation stations were for model development while the 40 calibration stations were for model parameter development. The details of calibration and validation stations are shown in Tables 1 and 2 respectively.

Out of the 55 stations, 15 stations were randomly selected as validation for model development and the remaining for calibration to develop model parameters. The cross-validation analysis was carried out to determine the difference between calibration and validation stations results. The process for which is as follows:

1. Regional parameter values (⁠ and ⁠) are estimated for 40 calibration stations (parameter set A) and for 15 validation stations (parameter set B) for each month.

2. Daily R factor is estimated for calibration stations using parameter set B and parameter and daily R factor is estimated for validation stations using parameter set A.

3. The estimated daily R factors so obtained are compared with the observed values and daily errors were calculated. Coefficient of determination (R²) was also calculated during parameter estimation for each month.

In the proposed model, maximum intensity is not used. Hence, the error between observed and predicted values would be different in the events of less precipitation with high intensity than the events of high precipitation with low intensity. To evaluate the model efficiency, the validation statistics commonly used by Yu et al. (2001), Petkovšek & Mikoš (2004), Bennett et al. (2013) and Xie et al. (2016) for the goodness-of-fit and error statistics were selected. These are:

1. Nash and Sutcliffe efficiency coefficient (NS) (Nash & Sutcliffe 1970),

2. symmetric mean absolute percentage error (SMAPE),

3. mean absolute error (MAE) and mean bias error (MBE).

MBE is an indicator of prediction bias while MAE is a measure of average error (Angulo-Martínez & Beguería 2009). NS, SMAPE, MBE and MAE were calculated for each validation station during the model development.

The thresholds of 0.1, 5.0 and 12.7 mm daily precipitations were taken individually as the minimum values of erosive event. For each threshold, the average annual rainfall value over a period of 10 years is taken for each station. Average annual rainfall for minimum limits of 0.1, 5.0 and 12.7 mm ranged from 1,885, 1,784 and 1,498 to 4,571, 4,458 and 4,140 mm, respectively. Average annual erosivity for the same limits varied from 10,264, 8,151 and 4,958 to 54,284, 48,301 and 39,938 MJ mm ha⁻¹ h⁻¹ y⁻¹, respectively.

For erosive events calculation, most of the days were of Type I in all the precipitation limits. A total of 84,283, 39,499 and 17,306 erosive events for 0.1, 5.0 and 12.7 mm respectively were used in the study; 24,402, 11,151 and 4,765 erosive events belong to the 15 validation stations while 59,881, 28,348 and 12,541 belong to the calibration stations for 0.1, 5.0, and 12.7 mm precipitation limits respectively. For the limit of 0.1 mm, taking the average percentage over 55 stations, 78.3% erosive events were found to be of Type I, followed by Type II (5.9%) and Type III (15.7%). For the 5.0 mm limit, the average percentages are: Type I (92.6%), Type II (4.1%) and Type III (3.3%). For 12.7 mm limit, the average percentages are: Type I (95.9%), Type II (2.4%) and Type III (1.8%). The high values of Type I for the limit of 0.1 mm is looking low in terms of the values of Type I for the limits of 5 and 12.7 mm because these values are in terms of percentage. An example to understand this calculation is shown in Table 3.

All stations have similar percentage distributions in terms of event types. The majority being of Type I at all stations shows the importance of individual events (Figure 2).

A decline of 44.6 and 54.6% was observed in Type I, when limits changed from 0.1 to 5.0 and 5.0 to 12.7 mm respectively. The percentage declines seen in Type II were 67.1 and 90.2 when the limits were increased from 0.1 to 5.0 and 5.0 to 12.7 mm respectively. Type III revealed smaller differences in percentage declines of 74.9 and 76.8% when the limits were raised from 0.1 to 5.0 and 5.0 to 12.7 mm. Since the decline in Type II and Type III events is more than Type I events, it shows that the area has small duration (Type I) events even for 5.0 and 12.7 mm limits (Table 4).

This indicates the individual power precipitation and high erosive capacity events within a short time interval that would increase the chances of flash flood.

It should also be noted that smaller events (less than 5.0 mm) which comes under 0.1 mm as the minimum precipitation limit, maintains the moisture in the soil and hence decreases the infiltration capacity of soil. This increases in the surface runoff which could raise the possibility of flash floods in the area.

The month wise average of observed R factor values for all stations obtained by using Equations (6), (8) and (9) were compared to determine the difference between them. R factor values determined by using Equation (8) were slightly higher than the other two R factor values except in the month of September and November while the R factor values from Equations (6) and (9) showed similar results except in the month of November (Figure 3). The high R factor values were observed in the month of November and December due to high rainfall with high intensity. As Equation (6) belongs to a wet and humid climate similar to the study area, it can be concluded that Equation (8) has showed slight overestimation for this particular climate while Equation (9) showed almost similar results.

Daily rainfall erosivity during the study period for each station was calculated as per the described methodology for the precipitation limits of 0.1, 5.0 and 12.7 mm. The data was then separated month wise. The distributions of average monthly rainfall erosivity and precipitation with all the three limits for all the stations are shown in Figures 4–6 which shows a similar pattern except for November and December which are showing the highest rainfall erosivity at the precipitation limit of 0.1 mm. The values have declined for 5.0 mm (23.2 and 30.5% respectively) and 12.7 mm limits (51.7 and 61.7% respectively). This shows that small events play an important role in maintaining soil moisture in November and December which could result in high soil erosion during high intensity rainfall.

Minimum average precipitation was recorded in April with values of 127.6, 121.4 and 106.9 mm for 0.1, 5.0 and 12.7 mm respectively while December experienced the maximum with values of average precipitation, i.e. 579.3, 567.9 and 531.0 mm for 0.1, 5.0 and 12.7 mm respectively. November comes second in terms of high average precipitation, followed by February. February generally has low precipitation throughout the month with few events of high intensity rainfall in a short period of time. The maximum 30 min intensity values were observed in October (164 mm/h) followed by January (161 mm/h) and December (137.2 mm/h) which are very few in the study period.

The monthly α and β values for all stations were calculated using Equation (9) with precipitation limits of 0.1, 5.0 and 12.7 mm. Figures 7–9 show the box-and-whisker plots for the overall behaviour of the monthly distributions of α and β values for the precipitation limits of 0.1, 5.0 and 12.7 mm. The length of each box indicates the interquartile range with the median (middle line) and the whiskers are 25th and 75th percentiles.

Monthly distributions of β values were quite consistent for all the precipitation limits (0.1, 5.0 and 12.7 mm) while α values have large variations from January to February and November to December for precipitation limits of 5.0 and 12.7 mm and small variations for the limit of 0.1 mm for all the months. The average ranges of monthly α, β and R² values for all stations are shown in the supplementary data, Tables S1–S3 (available with the online version of this paper).

The variations in the monthly α and β values from January to February may be due to the lower occurrence of precipitation with high I₃₀ values. Both parameters show significant spatio-temporal behaviour within the study period. Many studies (Richardson et al. 1983; Bagarello & D'Asaro 1994; Yu & Rosewell 1996a; Petkovšek & Mikoš 2004; Xie et al. 2016) have minimized the influence of β by keeping β constant throughout the year or throughout the study period because of its greater influence as compared to α on the estimations and sensitivity to calibration errors (Angulo-Martínez & Beguería 2009). However, our results show that both α and β parameters vary significantly which supports the hypothesis that a model incorporating such variations could yield better results. Therefore, in this study, both the parameters will be used as variables for each month in the model development.

Because of the spatial and temporal behaviour of parameters α and β, the best monthly values for each station were obtained by applying linear regression on values from the same month for all the years and for all precipitation limits. For instance, all daily January R factor and precipitation values of 2005, 2006, 2007 and up to 2015 were separated and one best value of α and β for January were obtained by applying linear regression on this set of data for a particular station. The same procedure was applied for every month for each station so that every station obtained the best monthly values of α and β. In model calibration and validation, the so obtained best monthly parameters α and β values of 40 and 15 stations respectively for each month in the precipitation limit of 0.1, 5.0 and 12.7 mm were used in the daily rainfall erosivity estimation. In terms of under- and overestimation of the daily rainfall erosivity, both sets of parameters showed mixed results. At the precipitation limit of 0.1 mm, October showed the lowest average error (1.3% underestimation) while November experienced the highest average error of 48.3% in underestimation of the calibration (average 10% underestimation). In the validation, the lowest average error was found in September (1.1% overestimation) and the highest in February (44.3% overestimation) with an overall average error of 4.7% overestimation. The average error of 161.8 and 305.1% underestimation in calibration, while 46.1 and 40.2% underestimation in validation, were calculated for the precipitation limit of 5 and 12.7 mm respectively. Large errors are due to the low precipitation and high I₃₀ of some particular events. The closest estimation in validation was observed for February with 4.8% and January with 2.8% overestimation in the precipitation limits of 5 and 12.7 mm respectively. The results of monthly average values of parameters α and β with R² for validation stations for all precipitation limits are shown in Table 5.

The rainfall erosivity model (Equation (10)) generates reasonable approximations of average monthly R factor, especially for the precipitation limit of 0.1 mm. For validation stations, at precipitation limit 0.1 mm, the best average NS were for the months of July (0.88) followed by June (0.86) and May (0.85) while the best average SMAPE were for August (1.32) followed by May (1.33) and July and September (1.35 for both). Worst average NS and SMAPE for this limit was for November (0.28) and February (1.65) respectively. In the precipitation limits of 5.0 and 12.7 mm, the best and worst average values of NS was for August (0.57) and November (−5.06) and worst values for August (0.16) and December (−4.87). The model works well for both average MBE and MAE at precipitation limit 0.1 mm, which is low for all the months except November to February. These high values are due to the role of high intensity rainfall in observed values for precipitation limit of 0.1 while for the limit of 5.0 and 12.7 mm, the very high values are due to less number of events with higher intensity. However, model efficiency for calibration stations was not good in comparison to validations stations. Overall, it can be concluded from the results that the model estimation for validation stations was good with the precipitation limit of 0.1 mm (Tables 6 and 7 respectively). The monthly NS, MBE, MAE and SMAPE values of 15 validation stations for all precipitation limits are shown in the supplementary data (Tables S4–S7, available online).

The model efficiency in estimating seasonal variation and average annual erosivity of the other two aspects of R factor (Xie et al. 2016) were also checked for validation stations which were not good in comparison to event or daily R factor (Table 8).

The best prediction was from station number S9 in which the observed and estimated average annual R factor values were 15,874 and 15,641 MJ mm ha⁻¹ h⁻¹ y⁻¹ respectively. The worst prediction was made for station number S21 in which the observed and estimated average annual values were 19,352 and 15,692 MJ mm ha⁻¹ h⁻¹ y⁻¹ respectively. The average error for all the validation stations is found as 7.8% overestimation.

To check the results of the proposed model for this study, rainfall erosivity was estimated for all the districts of the area in which the number of stations range from 16 (Gua Musang) to 1 (Machang) (Figure 1). The previously developed model (Yu et al. 2001) was also included to compare the results of the proposed model. From the proposed model, the best estimation of average annual R factor was 16,773 MJ mm ha⁻¹ h⁻¹ y⁻¹ (observed) and 16,746 MJ mm ha⁻¹ h⁻¹ y⁻¹ (estimated) for the district Gua Musang (16 stations) while the worst prediction was for Machang (one station). The overall average error of the proposed model is 8.2% overestimation while the previously developed model was estimated with 32.8% overestimation as shown in Figure 12.

It was found that 15% overestimation from the area are due to overestimation of the average R factor in February (40%), November (30%) and December (15%) because of high precipitation in November and December, but high rainfall intensity in February each year. The overall results of the proposed model and previously developed model for the study area are shown in Table 9.

It can be concluded that the proposed model works well for all the months except February, November and December in estimating the average annual R factor.

Although the study was carried out for Kelantan state of Peninsular Malaysia since the conditions used are for a general tropical climate region, the proposed model may be applied in any other region of similar climatic conditions. The study for the area of similar climatic conditions should be carried out by using the same average monthly α and β values as given in Table 5 with daily precipitation limits greater than 0.1 mm to check the results of estimating average annual rainfall erosivity. The limitation of the model is that it depends on the precipitation values while the observed R factor depends on both the precipitation and the intensity of the rainfall. The proposed model can overestimate for values of high daily precipitation with low rainfall intensity and vice versa. Future work thus may include the improvement of this model to incorporate the effects of daily intensity of rainfall.

The R factor for the proposed model can be calculated on a daily, monthly and annual basis. To apply the proposed model to other places, only daily precipitation data is needed. The regional parameter (α and β) values of 0.1 mm precipitation limit of all months will be used in R factor calculation which are given in Table 5. Examples of daily, monthly, annual and average annual R factor calculations from the proposed model are as follows:

1. For daily R factor calculations: For instance, to calculate the R factor of 1st January, the 1st January precipitation value and the January regional parameters values are put in Equation (16).

2. For monthly R factor calculations: For instance, to calculate the R factor for the month of February of a year, daily R factor of each day of this month are calculated (on similar lines as given in (i)) and summed up to obtain the monthly R factor of February for the specified year. For average monthly R factor calculations, the average of monthly R factor values of February for a specified number of years will be taken.

3. For annual R factor calculations: The monthly calculations of all the months of a year will be calculated as per the procedure of (ii) with their respective regional values from Table 5. After monthly calculation, the sum of all the months will be the annual R factor.

4. For average annual R factor calculations: Suppose Y number of years are included to calculate the average annual R factor. The average of R factor values of January for Y number of years will be taken. Similarly, the average of February, March and so on up to December R factor will be taken. Finally, the average monthly R factor values will be summed up to calculate the average annual R factor.

A total of 84,283, 39,499 and 17,306 erosive events for 0.1, 5.0 and 12.7 mm respectively from 30 min resolution rainfall data from 55 rain gauges located in the northeastern part of Peninsular Malaysia were used to analyse the empirical relationship between rainfall erosivity and daily precipitation. The study shows the dominance of individual (Type I) events in a day in the area which are 78.3, 92.5 and 95.8% for precipitation limits of 0.1, 5.0 and 12.7 mm respectively while the remaining days with erosive rainfall are the parts of Type II and Type III erosive events.

The values of coefficients α and β of the power law relationship were calculated by minimizing the sum of minimum square error and the spatio-temporal variations of both α and β values were found. Therefore, in the model development, α and β were taken as variables instead of constant values. Validation was carried out using 15 out of 55 stations using the average values of α and β for each month. R² of the validated stations ranged from 0.93 to 0.94, 0.64 to 0.84 and 0.41 to 0.70 for precipitation limits of 0.1, 5.0 and 12.7 mm respectively for all the months. During the model development, mixed results were obtained including both the overestimations and the underestimations of the rainfall erosivity for all precipitation limits. On the basis of the number of individual events, NS, SMAPE, MBE and MAE results, the 0.1 mm precipitation limit was selected for the proposed model. The average value of α was found to be lower than the average value of α for warm and cool season (proposed by a previous study) which indicates less significance of α as compared to β in the estimation of rainfall erosivity. The average value of β was found to be higher than that proposed in a previous study which may be the cause of high rainfall erosivity estimation in the area. Furthermore, the proposed model was compared to the previously developed model for similar climatic conditions and the estimation results from the proposed model were found to be better for most of the months, except for a little overestimation in February, November and December.

We gratefully acknowledge the Malaysian Meteorological Department, Department of Irrigation and Drainage, School of Physics and School of Industrial Technology of University Sains Malaysia for providing the required research facilities, University fellowship and data for this work. We are also grateful to the reviewers of this paper for reviewing the paper and providing their valuable suggestions. We would also like to acknowledge Universiti Sains Malaysia for providing financial support through grant 203/PJJAUH/6764002 and 1001/PTEKIND/8011021 to carry out this work.