An assessment of tropical cyclones rainfall erosivity for Taiwan

Rainfall erosivity (or water erosion) has severe implications on agriculture, water, and land use management. Though, there were Rainfall erosivity studies on regional and global scale, tropical cyclones’ Rainfall erosivity is poorly assessed and have not been documented for one of the most cyclones affecting regions of the world like Taiwan. Here, using 15-years of raindrop size distributions (RSD) and 60-years of hourly rain gauges data, we estimated cyclones (also called typhoons) rainfall erosivity over Taiwan, and establish that typhoons’ mean rainfall erosivity is higher than the global mean rainfall erosivity. Moreover, regional variability of typhoons rainfall erosivity showed an increasing pattern from north to south (Taipei to Pingtung), with relatively higher values over eastern and southern parts of Taiwan. The annual mean erosivity of typhoons rainfall showed raising trends over eastern and southern Taiwan during 1958–2017. Our results provide an insight in assessing the land use and agricultural management for Taiwan.


Establishment of kinetic energy-rainfall intensity (KE-I) relations.
and linear 22 forms of KE time -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 KE time 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.
A total number of 70600 one-min data points of kinetic energy content (KE mm ) and intensity (I) are shown in Fig. 1b. The KE mm 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 KE mm of current typhoons is higher than the global average maximum KE mm (28.3 ± 2.9 J m −2 mm −1 ) 20 . The KE mm 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 KE mm values at lower rainfall intensities (<15 mm h −1 ), a maximum KE mm of 35.71 J m −2 mm −1 is noticed. For rainfall intensity greater than 75 mm h −1 , the KE mm 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, logarithmic 25 , and exponential 18,23 fitting curves for KE mm 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 KE mm for rainfall intensities >30 mm h −1 . The exponential form of KE mm -I relation in the present study is different from that of the Chang et al. 37 (KE mm = 32.19 [1-0.725e −0.029I ]), which is due to their collective consideration of both seasonal and typhoon rainfall RSD in deriving KE mm -I relations. Present exponential form of KE mm -I relation is also found to be different from that of the Van Dijk et al. 20 generalized equation (KE mm = 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 (R 2 ), 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 (KE time and KE mm , Fig. 1). Similar characteristics were reported by the previous researchers elsewhere 40 .

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 KE time versus I for three events with linear and power relations. For Fungwong and Namtheun events (Fig. 2a,c), estimation of KE time 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 (R 2 , RMSE, and NRMSE) for power and liner forms of Continued KE time -I relations are quantified in Table 2. Even though two equations showed relatively higher coefficient of determination (R 2 ), the power law showed lower RMSE and NRMSE, indicating that the power law has better performance than linear in estimating the KE time from rainfall intensity. Validation of power, linear, logarithmic, and exponential forms of KE mm -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 KE mm for event no. 35 (for I > 30 mm h −1 ) and event no. 40 (for I > 40 mm h −1 ). The KE mm is well predicted by the power equation for long duration-high intensity rainfall event (Event no. 35. Fungwong, Fig. 2d), nonetheless, KE mm is overestimated for remaining two events (Event no. 39 and 40). The goodness-of-fit statistics (R 2 , RMSE, and NRMSE) for power, linear, logarithmic and exponential equations of KE mm -I relations are provided in Table 2 for the three typhoon events. Further, estimated KE-I models (KE time -I: linear, power, KE mm -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 KE time -I and KE mm -I relations.  www.nature.com/scientificreports www.nature.com/scientificreports/ Because of spurious self-correlations between KE mm and I, and relatively higher KE mm values at lower rainfall intensities, the KE mm -I relations lead to uncertainty in predicting the rainfall KE 21 for larger I values. This characteristics is clearly persistent even for present typhoon rainfall events (Fig. 1b). Hence, we prefer to adopt KE time -I relation rather than KE mm -I relation in computing the rainfall erosivity and was strongly suggested 21 .
Spatial variation of typhoon rainfall, rainfall erosivity and erosivity density. During 1958During -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 years 14 , 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 (EI 30 ) 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 areas 42 . The strong easterly winds of typhoons are blocked by the CMR leading to less precipitation across west coast of Taiwan 43 . Annual mean precipitations of Taiwan Counties are detailed in Table 3.  (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 Taiwan 44 . 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 EI 30 values at event as well as annual wise (Supplementary Fig. 4), and the EI 30 -P relations at event and annual level are estimated as, respectively, EI 30 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 EI 30 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 EI 30 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 EI 30 /R factor values. Further, annual  www.nature.com/scientificreports www.nature.com/scientificreports/ 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 Lin 46 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 Lin 46 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 duration 29 . 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 ).  5.7 (2000-2006). Tu and Chou 58 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 documented 1,53,54,60,61 , trends in erosivity triggered by tropical cyclones rests unidentified.

Trends in typhoons precipitation and erosivity.
Because of complex orography of Taiwan 43 , 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  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 stations 47,52,62 or short-term data with dense rain gauges 26,27,29 by using elsewhere KE-I relations rather than adopting region specific relations. Conversely, some researchers used modified Fournier index 30 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 erosivity 33 . 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 mean 31 . 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 Xie 54 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 Xie 54 ). An increasing www.nature.com/scientificreports www.nature.com/scientificreports/ 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 northward 44,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 Taiwan 59 , 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 km 2 ) 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 rainfall 63 . 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.
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 site 65 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 event 13 . 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 (KE mm and KE time ). The rainfall intensity (I, mm h −1 ), kinetic energy expenditure (KE time , J m −2 h −1 ), and kinetic energy content (KE mm , 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 (KE time , in J m −2 h −1 ) and volume specific kinetic energy (KE mm , J m −2 mm −1 ). KE time is kinetic energy per unit area per hour and KE mm is the kinetic energy per unit area per unit depth [18][19][20][21] .
The rainfall intensity and kinetic energy are computed by using below equations. Rainfall intensity (I, mm h −1 ),   66 .
In this study, data points with rainfall intensities less than 0.1 mm h −1 are discarded 67  www.nature.com/scientificreports www.nature.com/scientificreports/ 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 (EI 30 ). The average annual typhoon rainfall erosivity (R) is computed as: 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, m j is the number of erosive events of a given year j, and EI 30 is the typhoon rainfall erosivity index of a single event k. The single event erosivity EI 30 (MJ mm ha −1 h −1 ) is defined as: 30 30 where I 30 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 KE time -I relation: 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. 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 "P i,j " is typhoon-induced rainfall accumulation of n th event in m th year. The annual mean EI 30 value at each rain gauge station is computed by using the expression as mentioned below. 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 "X i,j " is EI 30 of a typhoon-induced rainfall of n th event in m th year. Kriging interpolation is applied to annual mean precipitation and annual EI 30 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 EI 30 .