Introduction

Tropical cyclones cause devastating loss of life and property globally1,2. Taiwan is one of the most intense tropical cyclones (typhoons) prone areas in the world3 with an average number of 3–4 typhoons striking the island every year4. Typhoons’ rainfall extremes over Taiwan are accountable for more sediment discharge5,6, natural hazards like floods, debris flow, and landslides7,8,9,10. Typhoons-induced rainfall kinetic energy (KE) is the driving force for surface runoff and landslide, and can be used as a proxy to determine the landslides triggering11. The rainfall KE is the key parameter in the estimation of rainfall erosivity factor (or R-factor), a substantial factor for soil erosion risk assessment in Universal Soil Loss Equation (USLE), Revised Universal Soil Loss Equation (RUSLE), and RUSLE212,13,14,15. The ability of rainfall to cause soil erosion can be described in terms of R-factor, and its variability immensely influences the agriculture, water, and land use management. The R-factor of a given ___location is defined as the annual accumulation of EI30 index (in M J mm ha−1 h−1): the product of kinetic energy of each event (E) and its maximum 30-min rainfall intensity (I30)12,14. An apprehending of R-factor can prominently improves an accurate assessment of soil erosion.

Since the direct measurement of rainfall KE requires a precise and expensive experimental setup16,17, alternative approaches have been developed to compute the rainfall KE from rainfall intensity (I)18,19, which require empirical relations between KE and I. As the KE of raindrop is proportional to the third power of diameter and square of fall velocity, the combined information of raindrop size and fall velocity of rain allows us to derive KE-I relations18,19,20,21. The KE-I relations have been developed in the form of linear22, exponential18,23, power24, and logarithmic25, which are effective for the observational site or to the regions with similar geographical and meteorological characteristics21,23,24,25, and vary with geographical ___location, local climate and precipitation microphysics18,20,23.

Although, most of the previous studies on R-factor were devoted to seasonal or annual basis at regional, national, continental and global level26,27,28,29,30,31, limited studies were reported on the basis of storm types32, even such study has the limitation due to the adaptation of elsewhere KE-I relations in evaluating the R-factor. Here ‘elsewhere’ refers to the other country’s KE-I relations. Due to the paucity of local KE-I relations, vast majority of rainfall erosivity studies have adopted elsewhere KE-I relations15,31,32 and only limited researchers used their local KE-I relations33. Similar scenario can be seen even for Taiwan in estimating the typhoons rainfall KE/erosivity11, which may lead to overestimation or underestimation of typhoons rainfall erosivity for this island19,33,34,35. Though, there were reports on local KE-I relations for Taiwan, they are limited to the seasonal rainfall or combination of monsoon and typhoon rainfall36,37. Moreover, the raindrop size distribution (RSD) characteristics of tropical cyclones and their KE-I relations were found to be different from that of the non-tropical cyclones precipitation19,38,39,40. Henceforth, it is crucial to explore Taiwan typhoons R-factor by adopting indigenous KE-I relations.

Here, we use the long-term RSD of 65 typhoons (2002–2017) and dense network of hourly rain gauges data (1958–2017). We investigate and evaluate the KE-I relations for typhoon rainfall, and demonstrate that the power form of KE-I relation is appropriate. Furthermore, the rainfall erosivity (R-factor) and erosivity density (R-factor density) are assessed for 393 typhoon rainfall events (occurred during 1958–2017) by adopting the estimated power form of KE-I relation to all rain gauge stations distributed over Taiwan. Also, trends in typhoons’ rainfall erosivity across the Island are discussed.

Results

The rainfall statistics of 75 typhoon rainfall erosive events computed from RSD measurements of disdrometer during 2002–2017 are summarized in Table 1. The rainfall depth, number of particles, duration of the event, maximum and mean rainfall intensity of these 75 typhoon events varies from 10.43–355.75 mm, 31175–1615753, 159–6900 min, 9.97–143.17 mm h−1 and 0.24–14.22 mm h−1, respectively. The number of particles range, 31175 to 1615753, represents the number of raindrops recorded by the JWD for 75 typhoon rainfall events. Out of 75 typhoon rainfall events, we used 72 events (70600 data points) to establish the empirical relations between kinetic energy (Kinetic energy expenditure - KEtime, and kinetic energy content - KEmm) and rainfall intensity (I) by using nonlinear least square regression analysis. The rest three events (Fungwong, Lupit, and Namtheun: designated with * symbol, in bold and italic font in Table 1) are used to validate the derived KE-I (KEtime-I and KEmm-I) relations.

Table 1 Rainfall statistics (Date, duration, no. of particles, rainfall depth, and intensity) of 75 typhoon rainfall erosive events observed by disdrometer.
Full size table

Establishment of kinetic energy-rainfall intensity (KE-I) relations

Figure 1a depicts the scatter plot for kinetic energy expenditure (KEtime) and intensity (I) data points (70600) for 72 typhoon rainfall events. The KEtime ranges from 0.1–3517.9 J m−2 h−1, corresponding to rainfall intensities of 0.1–123.4 mm h−1. The mean value of KEtime is 72.65 J m−2 h−1 with standard deviation of 194.75 J m−2 h−1. The best estimates for power24,41 and linear22 forms of KEtime-I relations from the regression analysis are given Fig. 1a. The regression line derived for the data points showed relatively higher coefficient of determination for both power and linear equations. Even so, for rainfall intensity >50 mm h−1, the KEtime is underestimated by the linear equation and the power equation created good estimates with its regression line approximately passing through middle of data points, and the data points are relatively least scattered with uniform distribution over power regression line.

Figure 1
figure 1

Scatter plot of typhoon rainfall (a) kinetic energy expenditure (KEtime) and Intensity (I) fitted with linear and power laws (b) kinetic energy content (KEmm) and Intensity (I) fitted with linear, power, logarithmic and exponential equations for 73 typhoon rainfall events (70600 1-min data points).

Full size image

A total number of 70600 one-min data points of kinetic energy content (KEmm) and intensity (I) are shown in Fig. 1b. The KEmm varies from 1.05–39.11 J m−2 mm−1 for the rainfall intensity values of 0.1–123.4 mm h−1 with mean and standard deviation of 11.14 J m−2 mm−1 and 5.79 J m−2 mm−1, respectively. The maximum KEmm of current typhoons is higher than the global average maximum KEmm (28.3 ± 2.9 J m−2 mm−1)20. The KEmm data points are more scattered at lower rainfall intensities (<15 mm h−1) and are getting narrowed with the increase of rainfall intensity. By excluding the higher KEmm values at lower rainfall intensities (<15 mm h−1), a maximum KEmm of 35.71 J m−2 mm−1 is noticed. For rainfall intensity greater than 75 mm h−1, the KEmm ranges from 23.61–34.11 J m−2 mm−1 with mean and standard deviation of 26.4 J m−2 mm−1, 1.61 J m−2 mm−1, respectively. The established power, linear, logarithmic25, and exponential18,23 fitting curves for KEmm and I are depicted in Fig. 1b.

Among four models (Fig. 1b), power model showed relatively higher performance, and other two models (exponential and logarithmic) are significantly underestimated the KEmm for rainfall intensities >30 mm h−1. The exponential form of KEmm-I relation in the present study is different from that of the Chang et al.37 (KEmm = 32.19 [1-0.725e−0.029I]), which is due to their collective consideration of both seasonal and typhoon rainfall RSD in deriving KEmm-I relations. Present exponential form of KEmm-I relation is also found to be different from that of the Van Dijk et al.20 generalized equation (KEmm = 28.3 [1-0.52 exp(−0.042I)]) and evidently shows that the KE-I relations of typhoon rainfall are different from that of the non-typhoon rainfall, and confirms that, there is need to adopt local tropical cyclones KE-I relations in computing the rainfall KE/erosivity of tropical cyclones of a given region.

The goodness-of-fit statistics in terms of coefficient of determination (R2), root mean square error (RMSE), and normalized RMSE (NRMSE) for power, linear, logarithmic and exponential equations of KE-I relations for 72 events are provided in Table 2. From the Table 2, amid four regression models (power, linear, logarithmic and exponential), statistically the power law showed higher predictive capability for the two erosivity indices (KEtime and KEmm, Fig. 1). Similar characteristics were reported by the previous researchers elsewhere40.

Table 2 Statistical parameters [correlation coefficient: R2, Root mean square error (RMSE), normalized RMSE] of 72 events and the three events selected for validation.
Full size table

Validation of kinetic energy-rainfall intensity (KE-I) relations

It is important to endorse the estimated KE-I relations before using them in rainfall erosivity assessment, and so, we validated the derived empirical (KE-I) relations with three typhoon events (Event no. 35, 39, and 40 in Table 1). Among these three events, two events have higher rainfall intensities with long (Event no. 35. Fungwong) and short (Event no. 40. Namtheun) duration, and the third one (Event no. 39. Lupit) has long duration and relatively low rainfall intensities, but this is quite erosive if the whole event is considered. Figure 2a–c represents KEtime versus I for three events with linear and power relations. For Fungwong and Namtheun events (Fig. 2a,c), estimation of KEtime by power law showed a good performance and underestimated by linear law at higher rainfall intensity (>35 mm h−1). In the case of the Lupit event (Fig. 2b), even if both power and linear laws overestimated, the power law shows more predictive capability than liner law. The goodness-of-fit statistics (R2, RMSE, and NRMSE) for power and liner forms of KEtime-I relations are quantified in Table 2. Even though two equations showed relatively higher coefficient of determination (R2), the power law showed lower RMSE and NRMSE, indicating that the power law has better performance than linear in estimating the KEtime from rainfall intensity.

Figure 2
figure 2

Validation of estimated KEtime-I relations (a–c), KEmm-I (d–f) in the form of linear (blue color), power (red color), exponential (green color), and logarithmic (pink color) laws with Fungwong (Event no. 35), Lupit (Event no. 39), and Namtheun (Event no. 40) typhoon events, respectively.

Full size image

Validation of power, linear, logarithmic, and exponential forms of KEmm-I relations with three events (Event no. 35, 39, and 40) are shown in Fig. 2d–f. In all the three events (Fig. 2d–f), power, logarithmic and exponential relations showed similar prediction capability for lower rainfall intensities (<5 mm h−1). Exponential and logarithmic fits severely underestimated the KEmm for event no. 35 (for I > 30 mm h−1) and event no. 40 (for I > 40 mm h−1). The KEmm is well predicted by the power equation for long duration- high intensity rainfall event (Event no. 35. Fungwong, Fig. 2d), nonetheless, KEmm is overestimated for remaining two events (Event no. 39 and 40). The goodness-of-fit statistics (R2, RMSE, and NRMSE) for power, linear, logarithmic and exponential equations of KEmm-I relations are provided in Table 2 for the three typhoon events. Further, estimated KE-I models (KEtime-I: linear, power, KEmm-I: linear, power, logarithmic, and exponential) are cross-validated by using leave p out cross validation method (where p = 3 events), and the iterations were performed for 1000 times. We found that the power model is more reliable than the other models for both KEtime-I and KEmm-I relations.

Because of spurious self-correlations between KEmm and I, and relatively higher KEmm values at lower rainfall intensities, the KEmm-I relations lead to uncertainty in predicting the rainfall KE21 for larger I values. This characteristics is clearly persistent even for present typhoon rainfall events (Fig. 1b). Hence, we prefer to adopt KEtime-I relation rather than KEmm-I relation in computing the rainfall erosivity and was strongly suggested21.

Spatial variation of typhoon rainfall, rainfall erosivity and erosivity density

During 1958–2017, with Central Weather Bureau (CWB) typhoon warning periods, a total number of 423 typhoons were recorded by dense network of 711 rain gauge stations over Taiwan, and among 423 typhoons 393 were qualified for the erosive event criteria. Nevertheless, in the current study, in order to estimate the representative typhoons rainfall erosivity, rain gauge stations (288) with minimum record period of 20 years14, whose elevation ranges from 2–3844.8 km asl are adopted. These 288 rain gauge stations have typhoon events for a total number of 7227 years with a mean value of 25 years per station, with number of events ranging from 58 to 368. The erosivity factor (EI30) of 393 typhoons over 288 rain gauge stations are computed by applying the estimated KE-I relation.

The annual mean precipitation, rainfall erosivity, and erosivity density of 393 typhoons over Taiwan are depicted in Fig. 3. The estimated stations’ annual mean precipitation ranges from 196.76 to 1427.3 mm yr−1, with mean and standard deviation of 578.97 mm yr−1 and 231.04 mm yr−1, respectively. The Island-wide gridded annual mean precipitation values are derived by applying kriging interpolation to stations’ annual mean precipitation. The gridded annual mean typhoons-induced rainfall varies from 288.57–877.79 mm yr−1, with mean and standard deviation value of 586.90 mm yr−1 and 155.44 mm yr−1, respectively (Fig. 3a). The spatial distribution of typhoons precipitation depicts relatively higher values over eastern and south eastern part of Taiwan. Different (north, east, central, and south) regions of Taiwan are depicted with different color in Supplementary Fig. 1b. The typhoons-induced event rainfall amounts contributed to a percentage of 14.3%, 25.7%, 29.1%, and 30.9%, respectively, for north, central, east, and south regions of Taiwan. The deep central mountain ranges (CMR) of Taiwan helps in enhancing the convective activity of typhoons while they cross it, resulting in relatively higher precipitation over eastern and south eastern part of Taiwan than other areas42. The strong easterly winds of typhoons are blocked by the CMR leading to less precipitation across west coast of Taiwan43. Annual mean precipitations of Taiwan Counties are detailed in Table 3.

Figure 3
figure 3

Spatial distribution of typhoons-induced (a) mean annual precipitation (b) R-factor, and (c) R-factor density map plotted with kriging over Taiwan during 1958–2017.

Full size image
Table 3 County wise mean and standard deviation values of accumulated precipitation, annual mean precipitation, EI30, R-factor and R-factor density.
Full size table

The R-factor values of 288 rain gauges with record period of 20–60 years varied from 886.75–22653.03 MJ mm ha−1 h−1 yr−1, with mean and standard deviation value of, respectively, 4880.74 MJ mm ha−1 h−1 yr−1 and 3261.97 MJ mm ha−1 h−1 yr−1. The Island-wide gridded annual mean R-factor values are derived by applying kriging interpolation to stations’ annual mean R-factors. The gridded R-factor values differs from 1005.73–9787.41 MJ mm ha−1 h−1 yr−1, with mean and standard deviation value of 4905.83 MJ mm ha−1 h−1 yr−1 and 1882.69 MJ mm ha−1 h−1 yr−1, respectively (Fig. 3b). Regional variability of typhoons R-factor showed an increasing pattern from north to south, with relatively higher values over eastern and southern part of Taiwan (Fig. 3b). As mentioned in the above paragraph, typhoons-induced event rainfall amounts contribute to a percentage of 14.3%, 25.7%, 29.1%, and 30.9%, respectively, for north, central, east, and south regions of Taiwan. Similarly, typhoon-induced rainfall erosivity contributed to a percentage of 9.3%, 25%, 28.1% and 37.2%, respectively, for north, central, east, and south regions of Taiwan. Typhoons invading region of Taiwan is surrounded by two major typhoon paths over the northwest Pacific region. One path moves south of Taiwan in westward to the south China sea, and the other turns in the north direction toward either Japan or Korea, which pass through the east side of Taiwan44. Because of steeper height of the central mountain range (CMR) of Taiwan, which is extended from north to south of the Taiwan Island, majority of the typhoons are blocked by these CMR, resulting in higher rainfall amounts over eastern and southern part of Taiwan (Fig. 3a), subsequently, relatively higher erosivity values over eastern and southern part of Taiwan. The regions with higher R-factor over Taiwan are in order with the higher precipitations areas (Fig. 3a). In addition to that, a good correlation is observed between precipitation (P) and EI30 values at event as well as annual wise (Supplementary Fig. 4), and the EI30-P relations at event and annual level are estimated as, respectively, EI30 = 0.68 P1.44 (R2 = 0.80) and EI30 = 1.69 P1.23 (R2 = 0.86). For the mostly associated tropical cyclones rainfall over Pacific coast of Mexico, García-Oliva et al.45 perceived mean annual erosivity of 6525.2 MJ mm ha−1 h−1. Recently, Laceby et al.32 demonstrated that tropical cyclones contributed 22% of precipitation and 40% of rainfall erosivity over Fukushima region of Japan, and they illustrated that the annual mean precipitation and R-factor ranged from, respectively, 68 to 639 mm yr−1 (mean: 422 mm yr−1) and 142 to 4547 MJ mm ha−1 h−1 yr−1 (mean: 1462 MJ mm ha−1 h−1 yr−1). Storm wise variation of mean precipitation and EI30 ranged from 41–240 mm and 118–1695 MJ mm ha−1 h−1, respectively, (Table S1. of Laceby et al.32). However, for Taiwan region, event mean precipitation and EI30 of typhoons varied, respectively, from 12.5–3059.5 mm and 0.47–92144.9 MJ mm ha−1 h−1, which are greater than the values over Japan. This demonstrates that the typhoons rainfall over Taiwan has much influence with intense rainfall and higher EI30/R factor values. Further, annual mean precipitation, rainfall erosivity for each county are computed and given in Table 3. Counties’ mean accumulated precipitation varied from 5517.5–19368.39 mm yr−1 with minimum precipitation over Changhua and maximum over Taitung. Annual mean R-factor ranges from 1433.59–7441.02 MJ mm ha−1 h−1 yr−1 with minimum value over Changhua and maximum value over Pingtung. Lee and Lin46 estimated rainfall erosivity of 16560 storm events for Kaohsiung City and Pingtung County in south Taiwan. They found mean annual total rainfall between 1376–4070 mm yr−1 and annual rainfall erosivity between 15000–70000 MJ mm ha−1 h−1 yr−1. The discrepancy in annual mean rainfall and rainfall erosivity values between present study and Lee and Lin46 is might be due to their adaptation of elsewhere KE-I relations or due to different sampling (10-min) of the rainfall data.

In order to identify the regions with high risk due to erosive events, typhoons rainfall erosivity density (ratio of annual R-factor to the annual mean precipitation) is computed for each station. For all the 288 rain gauges stations, erosivity density varies from 3.11–19.17 MJ ha−1 h−1 with mean and standard deviation of 7.92 MJ ha−1 h−1, and 2.79 MJ ha−1 h−1, respectively. The Island-wide gridded R-factor density values are derived by applying kriging interpolation to stations’ R-factor density values. The R-factor density map of Taiwan (Fig. 3c) ranges from 4.57–9.68 MJ ha−1 h−1 with mean and standard deviation of 7.77 MJ ha−1 h−1 and 0.99 MJ ha−1 h−1, respectively. Areas with higher erosivity density are an indication of locations with higher rainfall intensity events with shorter duration29. Higher erosivity density values are noticed over southern part of Taiwan, which imply that southern part is most prone to the soil loss and high flood risk than other regions of Taiwan. The R-factor density values of each county of Taiwan are illustrated in Table 3. County wise R-factor density show minimum for Hsinchu (5.25 MJ ha−1 h−1) and maximum for Kaohsiung (10.52 MJ ha−1 h−1).

Trends in typhoons precipitation and erosivity

Quite diverse results were reported for rainfall erosivity trends over different parts of the world on seasonal and annual scale32,33,47,48,49,50,51,52. For instance, Webster et al.53, Emanuel1, and Mei and Xie54 noticed an increasing trend in the intensity of tropical cyclones, and in contrast, other researchers55,56 claim small or no trends. For Taiwan region, Tu et al.57 demonstrated that an abrupt shift in typhoon count series from 3.3 (1970–1999) typhoons per year to 5.7 (2000–2006). Tu and Chou58 analyzed the frequency, intensity, and duration of typhoon-induced rainfall over Taiwan by considering 21 rain gauge stations distributed over this Island for July to October months during 1970–2010. They perceived a significant increase in typhoon-induced rainfall over Taiwan which is due to the increase in number of typhoon days that affecting the Island rather than number of typhoons that pass through the typhoons invading region (18N–29.5N, 116E–126E). Recently, trend analysis for typhoon-induced rainfall at six rain gauges stations (along west coast: Taipei, Taichung, and Tainan; along east coast: Hualien, Taitung, and Hengchun) located over Taiwan was carried out by Liang et al.59. They showed that variations in the typhoon-induced rainfall trends at the selected six stations are related to the poleward shift of tropical cyclones over northwest Pacific, which is due to the weakening of the steering flow and western north Pacific subtropical high. Albeit, long-term trends in global tropical cyclones properties (number, intensity, duration, and destructive potential) have been documented1,53,54,60,61, trends in erosivity triggered by tropical cyclones rests unidentified.

Because of complex orography of Taiwan43, we performed trend analysis for Taiwan by considering different regions/counties rather than complete Island. Regions of this Island are classified into north (which include Taipei, Hsinchu, Taoyuan counties), central (Miaoli, Taichung, Changhua, Nantou, Yunlin, Chiayi), south (Tainan, Kaohsiung, Pingtung), and east (Yilan, Hualien, Taitung), and these regions are depicted with different color in Supplementary Fig. 1b. Trends in typhoons-induced annual mean precipitation and the corresponding erosivity for 15 counties of Taiwan are illustrated in Figs 4 and 5, respectively. The linear trends for all 15 counties are performed by least square regression analysis. A trend is considered as significant if ‘p’ value is less than 0.1 (90% confidence level).

Figure 4
figure 4

Time series of mean annual precipitation induced by typhoon rainfall events for 15 counties of Taiwan during 1958–2017. The rain gauge stations at Taoyuan, Hsinchu, and Maiaoli counties are available from the year 1991, and for Changhua and Yunlin from 1992. The slope of the trend line is depicted with ‘s’ in the legend. The gaps in mean annual precipitations lines denotes unavailability of data for the corresponding year.

Full size image
Figure 5
figure 5

Same as Fig. 4 except for rainfall erosivity (EI30).

Full size image

Over north Taiwan, for the annual mean precipitations of typhoons-induced rainfall, except for Taipei county, which shows statistically significant decreasing trend (Fig. 4a), the other two counties (Taoyuan and Hsinchu; Fig. 4b,c) show slightly raising trends. For eastern Taiwan, Yilan and Hualien counties show raising trend (Fig. 4d,e) and the other county (Taitung) shows decreasing trend (Fig. 4f). On the other hand, all the three counties of south Taiwan show raising trends (Fig. 4g–i), with statistically significant raising trend only at Kaohsiung. Over Central Taiwan, except for Chiayi and Miaoli, which show decreasing trend (Fig. 4j,o), the other four counties (Miaoli, Taichung, Nantou, Changhua, and Yunlin (Fig. 4k–n) show raising trends, with statistically significant only at Yunlin.

Figure 5 depicts the time series of typhoons-induced rainfall erosivity for 15 counties of Taiwan. For north Taiwan, Taipei and Taoyuan show decreasing trends (Fig. 5a,b) and Hsinchu shows raising trend (Fig. 5c). However, the trends at these three counties are statistically insignificant. Increasing trends can be noticed for the three counties of eastern Taiwan (Fig. 5d–f), nevertheless, the trend is statistically significant only at Yilan county. Raising trends are also found for the three counties of south Taiwan (Fig. 5g–i) with statistically insignificant at Pingtung and significant at Kaohsiung and Tainan. Among six counties of central Taiwan, statistically significant raising trend at Yunlin (Fig. 5n), statistically insignificant decreasing trends at Miaoli (Fig. 5j) and Chiayi (Fig. 5o), and increasing trends at Taichung (Fig. 5k), Nantou (Fig. 5l) and Changhua (Fig. 5m) can be seen.

Discussion

In computation of rainfall erosivity, most of the previous studies used either long-term rainfall data of few rain gauge stations47,52,62 or short-term data with dense rain gauges26,27,29 by using elsewhere KE-I relations rather than adopting region specific relations. Conversely, some researchers used modified Fournier index30 in estimation of erosivity, if there were no records of high temporal resolution precipitation data, and only few studies adopted the region specific KE-I relation in the reckoning of erosivity33. Nonetheless, globally, least attention was paid for the tropical cyclones-induced rainfall erosivity. Taiwan being most tropical cyclones affecting country, it is paramount to investigate the R-factor for typhoons rainfall.

In this study we used long-term raindrop size distribution information of typhoons to estimate the KE-I relations, and we establish that power law is more appropriate in estimating the Taiwan typhoons R-factor. The estimated KE-I relation is adopted in computing the typhoons R-factor with 60 years (1958–2017) of hourly rainfall data over Taiwan. The annual mean typhoons precipitation for Taiwan is noticed as 586.90 (SD 155.44) mm yr−1 with annual mean typhoons rainfall erosivity of 4905.83 (SD 1882.69) MJ mm ha−1 h−1 yr−1, which is much higher than the global mean31. Higher R-factor values are found over eastern and southern part of Taiwan, and higher R-factor density values are perceived for southern Taiwan. The typhoons invading region of Taiwan is surrounded by two major typhoons tracks in east and south direction (see Fig. 1 of Tu et al.44), and the interaction of typhoons that pass through these two major tracks interact with the CMR of the island resulting in an enormous rainfall amounts over eastern and southern region of Taiwan.

The county wise trend analysis performed for Taiwan showed increasing trends in rainfall and erosivity of typhoon events for eastern and southern Taiwan. Over northwest Pacific, during 1977–2013, Mei and Xie54 performed the trend analysis for frequency of four types of typhoon clusters (cluster 1–4). Among the four clusters, they demonstrated an increase in high intense (category 4 and 5) typhoons of cluster-1 and cluster-2, which are covered respectively, by eastern and southern Taiwan (Supplementary Fig. 8 of Mei and Xie54). An increasing trends in rainfall and erosivity over the eastern and southern region of this island is due to an increase in high intense, long-lasting typhoons and their shift to northward44,54,58.

In analyzing the trends for Taiwan, we consider all the typhoons that occurred throughout the year (Jan-Dec) during 1958–2017, and these typhoons should met an erosive event criteria as mentioned in methods sections. Also, we selected the rain gauge stations that are available for greater than or equal to 20 years of observation (≥20 years) during typhoon-induced rainfall periods of 1958–2017. It should be noted that the typhoons whose rainfall amounts are not meeting the erosive event criteria are missing in the trend analysis. Moreover, because of complex topography of Taiwan, a given typhoon need not to produce the rainfall for the complete Island. Further, for a given typhoon-induced rainfall, all the rain gauges distributed over Taiwan may not get a rainfall event that satisfy the erosive event criteria, hence, some rain gauge stations record a rainfall event and some may not. Although there were reports on trends for typhoon-induced rainfall at different regions of Taiwan59, because of the above mentioned reasons, rainfall erosive events are not the representative of total rainfall of typhoons and their characteristics need not to be similar to that of the total rainfall amounts. Our findings can contribute to better assessment of soil erosion modeling, agricultural, effective-land use and flood risk assessment for Taiwan.

Methods

Study area

Taiwan is a subtropical island in the Western North Pacific (~36197 km2) with geological composition of sedimentary and metamorphic rocks, and deep topographic CMR of average height 2 km. Taiwan is influenced by highest annual frequency of tropical cyclones which are contributing to 47.5% of the total annual rainfall63. The geographical ___location of Taiwan with its topography, 711 rain gauge stations, and the Joss-Waldvogel disdrometer (JWD)64 is shown in Supplementary Fig. 1.

Typhoons tracks data

The typhoon track data is based on the World Meteorological Organization’s Joint Typhoon Warning Center (JTWC) (http://www.usno.navy.mil/NOOC/nmfc-ph/RSS/jtwc/best_tracks/). The JTWC provides 6-hourly records of typhoon track information.

Disdrometer data selection criteria for typhoons rainfall

During 2002–2017, the raindrop size distribution measurements of disdrometer (1-min sample data) are treated as typhoon-induced rainfall when typhoon center (supplementary Fig. 2) is within a radius of 500 km from the disdrometer site65 and they are within the CWB typhoon warning periods.

Definition of erosive rainfall event

A rainfall episode with minimum duration of 30-min and rainfall depth of greater than or nearly equal to 12.5 mm and a 6-hours of continuous non-rainfall gap between two consecutive rainfall episodes is considered as erosive rainfall event13. Here we use erosive rainfall event criteria with rainfall depth of greater than or nearly equal to 10 mm. With this rainfall classification criteria, a total number of 75 rainfall events from 66 typhoons are identified from the disdrometer measurements for the period 2002–2017.

Validation of Disdrometer data

The disdrometer rainfall amounts of selected typhoon rainfall events are validated by comparing with collocated rain gauge. The scatter plot of event rainfall depths between disdrometer and collocated rain gauge is depicted in Supplementary Fig. 3. A good correlation is found between the rain gauge and disdrometer measurements.

Computation of rainfall intensity (I) and kinetic energy (KEmm and KEtime)

The rainfall intensity (I, mm h−1), kinetic energy expenditure (KEtime, J m−2 h−1), and kinetic energy content (KEmm, J m−2 mm−1) for 75 typhoon rainfall events are computed from the RSD information of JWD.

The rainfall kinetic energy (KE) is half the product of raindrop mass and the square of its velocity, and can be expressed in two forms as time specific kinetic energy (KEtime, in J m−2 h−1) and volume specific kinetic energy (KEmm, J m−2 mm−1). KEtime is kinetic energy per unit area per hour and KEmm is the kinetic energy per unit area per unit depth18,19,20,21.

The rainfall intensity and kinetic energy are computed by using below equations.

Rainfall intensity (I, mm h−1),

$${\rm{I}}=(\frac{{\rm{\pi }}}{6})(\frac{3.6}{{10}^{3}})(\frac{1}{{\rm{A}}\,{\rm{T}}})\mathop{\sum }\limits_{{\rm{i}}=1}^{20}\,{{\rm{n}}}_{{\rm{i}}}{{\rm{D}}}_{{\rm{i}}}^{3}$$
(1)

Kinetic energy expenditure (KEtime, J m−2 h−1),

$${{\rm{KE}}}_{{\rm{time}}}=(\frac{{\rm{\pi }}}{12})(\frac{1}{{10}^{6}})(\frac{3600}{{\rm{T}}})(\frac{1}{{\rm{A}}})\mathop{\sum }\limits_{{\rm{i}}=1}^{20}\,{{\rm{n}}}_{{\rm{i}}}{{\rm{D}}}_{{\rm{i}}}\,{\rm{V}}{({{\rm{D}}}_{{\rm{i}}})}^{2}$$
(2)

Kinetic energy content (KEmm, J m−2 mm−1),

$${{\rm{KE}}}_{{\rm{mm}}}=\frac{{{\rm{KE}}}_{{\rm{time}}}}{{\rm{I}}}$$
(3)

where A = 0.005 m2 is the sampling area of the sensor, T = 60 s is the sampling time, ni the number of drops of diameter Di, V(Di) is the fall velocity of drops with diameter (Di)66.

In this study, data points with rainfall intensities less than 0.1 mm h−1 are discarded67. The KEtime-I and KEmm-I empirical relations are derived for linear, exponential, logarithmic, and power laws by using non-linear regression analysis.

Rain gauge data selection criteria for typhoons rainfall

For the period of 1958–2017, hourly rain gauges data during central weather Bureau (CWB) typhoon warning periods (http://rdc28.cwb.gov.tw/TDB/ntdb/pageControl/typhoon) of Taiwan are considered as typhoons attributed rainfall.

Rain gauge stations selection criteria

In this study, we selected the rain gauge stations with minimum record of 20 years for the estimation of R-factor, with this threshold, among the dense network of 711 rain gauge stations distributed over Taiwan, a total number of 288 stations are qualified for typhoons rainfall record periods of 20 years or higher during 1958–2017. At these 288 rain gauge station, during 1958–2017, a typhoon-induced rainfall amount is considered for further analysis (estimation of R-factor and R-factor density) if it satisfy the erosive rainfall event criteria as mentioned above.

Calculation of R-factor (EI30)

The average annual typhoon rainfall erosivity (R) is computed as:

$$R=\frac{1}{{\rm{n}}}\mathop{\sum }\limits_{{\rm{j}}=1}^{{\rm{n}}}\,\mathop{\sum }\limits_{{\rm{k}}=1}^{{{\rm{m}}}_{{\rm{j}}}}\,{({{\rm{EI}}}_{30})}_{{\rm{k}}}$$
(4)

where R-factor is the average annual rainfall erosivity (MJ mm ha−1 h−1 yr−1), n is the number of years of record, mj is the number of erosive events of a given year j, and EI30 is the typhoon rainfall erosivity index of a single event k.

The single event erosivity EI30 (MJ mm ha−1 h−1) is defined as:

$${\rm{EI}}={{\rm{EI}}}_{30}={{\rm{eI}}}_{30}$$
(5)

where I30 is the maximum 30-min rainfall intensity (I, mm h−1) and is calculated by adopting the method of Yin et al.68. Where e is the rainfall kinetic energy per unit depth of rain and is obtained by using power law of KEtime-I relation:

$${{\rm{KE}}}_{{\rm{time}}}={{\rm{aI}}}^{{\rm{b}}}$$

where ‘a’ and ‘b’ are empirical coefficients

$${\rm{e}}=\frac{{{\rm{KE}}}_{{\rm{time}}}}{{\rm{I}}}=\frac{{{\rm{aI}}}^{{\rm{b}}}}{{\rm{I}}}={{\rm{aI}}}^{{\rm{b}}-1}$$
(6)

Annual R-factor values are used to drive rainfall erosivity maps by using ordinary kriging.

The annual mean precipitation value at each rain gauge station is computed by using the expression as mentioned below.

$${\rm{Annual}}\,{\rm{mean}}\,{\rm{precipitation}}\,{\rm{at}}\,{\rm{each}}\,{\rm{rain}}\,{\rm{gauge}}\,{\rm{station}}=\frac{{\sum }_{j=1}^{m}{\sum }_{i=1}^{n}\,{P}_{i,j}}{m}$$
(7)

where ‘n’ represents the number of events in a given year for a given station, and ‘m’ represents the number of years that have typhoon-induced rainfall events for a given station, and “Pi,j” is typhoon-induced rainfall accumulation of nth event in mth year.

The annual mean EI30 value at each rain gauge station is computed by using the expression as mentioned below.

$${\rm{Annual}}\,{\rm{mean}}\,{{\rm{EI}}}_{30}\,{\rm{at}}\,{\rm{each}}\,{\rm{rain}}\,{\rm{gauge}}\,{\rm{station}}=\frac{{\sum }_{j=1}^{m}\,{\sum }_{i=1}^{n}\,{X}_{i,j}}{m}$$
(8)

where ‘n’ represents the number of events in a given year for a given station, and ‘m’ represents the number of years that have typhoon-induced rainfall events for a given station, and “Xi,j” is EI30 of a typhoon-induced rainfall of nth event in mth year.

Kriging interpolation is applied to annual mean precipitation and annual EI30 values available at Island wide distributed 288 rain gauge stations (whose record period ≥ 20 years during 1958–2017) to plot the spatial distribution (Island-wide gridded data) map for annual mean precipitation and annual EI30.