Exploring the possible role of satellite-based rainfall data in estimating inter- and intra-annual global rainfall erosivity

Despite recent developments in modeling global soil erosion by water, to date, no substantial progress has been made towards more dynamic inter- and intra-annual assessments. In this regard, the main challenge is still represented by the limited availability of high temporal resolution rainfall data needed to estimate rainfall erosivity. As the availability of high temporal resolution rainfall data will most likely not increase in future decades since the monitoring networks have been declining since the 1980s, the suitability of alternative approaches to estimate global rainfall erosivity using satellite-based rainfall data was explored in this study. For this purpose, we used the high spatial and temporal resolution global precipitation estimates obtained with the National Oceanic and Atmospheric Administration (NOAA) Climate Data Record (CDR) Climate Prediction Center MORPHing (CMORPH) technique. Such high spatial and temporal (30 min) resolution data have not yet been used for the estimation of rainfall erosivity on a global scale. Alternatively, the erosivity density (ED) concept was also used to estimate global rainfall erosivity. The obtained global estimates of rainfall erosivity were validated against the pluviograph data included in the Global Rainfall Erosivity Database (GloREDa). Overall, results indicated that the CMORPH estimates have a marked tendency to underestimate rainfall erosivity when compared to the GloREDa estimates. The most substantial underestimations were observed in areas with the highest rainfall erosivity values. At the continental level, the best agreement between annual CMORPH and interpolated GloREDa rainfall erosivity maps was observed in Europe, while the worst agreement was detected in Africa and South America. Further analyses conducted at the monthly scale for Europe revealed seasonal misalignments, with the occurrence of underestimation of the CMORPH estimates in the summer period and overestimation in the winter period compared to GloREDa. The best agreement between the two approaches to estimate rainfall erosivity was found for fall, especially in central and eastern Europe. Conducted analysis suggested that satellite-based approaches for estimation of rainfall erosivity appear to be more suitable for low-erosivity regions, while in high-erosivity regions (> 1000–2000 MJ mm ha⁻¹ h⁻¹ yr⁻¹) and seasons (> 150–250 MJ mm ha⁻¹ h⁻¹ month⁻¹), the agreement with estimates obtained from pluviographs (GloREDa) is lower. Concerning the ED estimates, this second approach to estimate rainfall erosivity yielded better agreement with GloREDa estimates compared to CMORPH, which could be regarded as an expected result since this approach indirectly uses the GloREDa data. The application of a simple-linear function correction of the CMORPH data was applied to provide a better fit to GloREDa and correct systematic underestimation. This correction improved the performance of CMORPH, but in areas with the highest rainfall erosivity rates, the underestimation was still observed. A preliminary trend analysis of the CMORPH rainfall erosivity estimates was also performed for the 1998–2019 period to investigate possible changes in the rainfall erosivity at a global scale, which has not yet been conducted using high-frequency data such as CMORPH. According to this trend analysis, an increasing and statistically significant trend was more frequently observed than a decreasing trend.

1 Introduction

Rainfall erosivity is among the main drivers of soil erosion, which can be characterized by large spatial and temporal variability (Angulo-Martínez and Beguería, 2012; Ballabio et al., 2017; Bezak et al., 2021; Cui et al., 2020; Panagos et al., 2017; Verstraeten et al., 2006). In order to obtain robust rainfall erosivity estimates, high temporal resolution rainfall data are needed (Panagos et al., 2015; Yin et al., 2017). However, according to Panagos et al. (2017), the availability of stations with high-frequency data that can be used to estimate rainfall erosivity is on average relatively low in many parts of the world. Therefore, in areas with scarce data availability, remotely measured precipitation data can be instrumental in estimating rainfall erosivity (Ganasri and Ramesh, 2016; Li et al., 2020). Alternatively, approaches using simpler and less data-demanding methods such as erosivity density (ED) (Nearing et al., 2017; Panagos et al., 2015, 2016b) can also represent a viable option. Another condition that determines the accuracy of the rainfall erosivity estimates is the low availability of high temporal resolution rainfall data. Ideally, high-frequency (e.g., 1 min) measurements obtained using optical disdrometers are necessary to quantify the rainfall kinetic energy (Mineo et al., 2019; Nel et al., 2010; Petan et al., 2010; Sanchez-Moreno et al., 2014) of a given storm and to calculate its rainfall erosivity. However, such measuring equipment is not commonly available in regional and national measuring networks. Therefore, due to instrumental limitations, rainfall erosivity is communally estimated using hourly or sub-hourly rainfall records (generally ranging from 5 to 60 min) collected by tipping buckets or pluviographs, which do not provide information about raindrop size distribution (Panagos et al., 2016a; Petan et al., 2010; Petek et al., 2018). These kinds of data are then used together with empirically developed equations that relate rainfall kinetic power and intensity (Brown and Foster, 1987; Carollo et al., 2017; Petan et al., 2010) to obtain rainfall erosivity estimates. Alternatively, rainfall erosivity estimates can also be performed based on the rainfall volume instead of the intensity, using daily, monthly, or annual rainfall data (Renard and Freimund, 1994; Yu and Rosewell, 1996). However, it is worth mentioning that the accuracy of rainfall erosivity estimates decreases with the increase in the temporal data resolution (i.e., from 1 min to hourly, daily, monthly, or annual data). Currently, due to data scarcity, most rainfall erosivity assessments based on long-term estimates including a period of at least 10 years are limited to a few regions (Angulo-Martínez and Beguería, 2012; Nearing et al., 2015; Panagos et al., 2015, 2017), leaving large parts of the world under-researched. In this regard, a step forward is needed to enable the generation of year-by-year and sub-annual rainfall erosivity assessments for under-researched national- or larger-scale study areas.

Recent studies have already explored the possibility of estimating rainfall erosivity using satellite-based products at regional (Li et al., 2020) and national (Chen et al., 2021; Kim et al., 2020) scales, indicating their sources of uncertainties and a generally limited accuracy (Aghakouchak et al., 2012; Ghajarnia et al., 2018; Prakash, 2019; Prakash et al., 2015; Rahmawati and Lubczynski, 2018; Seo et al., 2018; Wei et al., 2018). However, to the best of the authors' knowledge, no such study has been conducted on a global scale using high temporal resolution data. A promising alternative to the often-limited rain-gauge data may be represented by satellite-based precipitation estimates, which currently have both adequate temporal and spatial resolution (Chen et al., 2021; Kim et al., 2020; Li et al., 2020; dos Santos Silva et al., 2020; Teng et al., 2017). Moreover, once further developed and fully operational, the satellite-based methods to estimate rainfall erosivity will have lower purchasing and processing costs compared to the current ones. In addition, satellite-based rainfall erosivity estimates could be especially useful in regions where rainfall erosivity estimates are currently very limited, such as some sizable sectors of Africa, Asia, and South America.

In this study, we aim to deepen the research on the use of satellite-based rainfall data in estimating rainfall erosivity by performing a first inter- and intra-annual global-scale assessment. The Global Rainfall Erosivity Database (GloREDa) data (Panagos et al., 2017) were used to evaluate both (a) the rainfall erosivity estimates obtained by satellite-based rainfall data (i.e., CMORPH) and (b) rainfall erosivity using the ED concept. Finally, a temporal trend analysis of global rainfall erosivity is presented with corrections between data based on the CMORPH and GloREDa databases.

2 Data and methods

2.1 CMORPH

The CMORPH product is a reprocessed and bias-corrected global precipitation data set covering the area between the 60^∘ S and 60^∘ N parallels, with a 30 min time step and a spatial resolution of 8 km × 8 km (Xie et al., 2017, 2021). The CMORPH data are developed by the National Oceanic and Atmospheric Administration (NOAA) and cover the period from 1998 onwards. This method generally uses the precipitation estimates derived from the low Earth orbit satellite-based passive microwave observations (Kim et al., 2020). Additionally, the geostationary satellite infrared imagery is used to account for possible coverage issues (Kim et al., 2020). Since CMORPH provides an estimate of the 30 min precipitation, each 30 min rainfall rate was assumed to be constant during this time interval (Kim et al., 2020; Xie et al., 2021). This data set has already been applied to several practical applications, such as validating the climate model simulations, identifying climate extremes, forcing numerical weather models, and characterizing the global precipitation (Xie et al., 2021). Additionally, details about the methodology can be found in the literature (Chen et al., 2020; Xie et al., 2017, 2021).

2.2 GloREDa database

GloREDa was created with the objective of developing the first ever global rainfall erosivity map using high temporal resolution data (Panagos et al., 2017) and moving towards a new generation of Revised Universal Soil Loss Equation (RUSLE)-based soil erosion assessments for the present (Borrelli et al., 2017) and future climate change and land use dynamics (Borrelli et al., 2020). GloREDa contains annual rainfall estimates for 3625 stations from 63 countries with temporal resolution ranging from 1 to 60 min (Panagos et al., 2017). The data sample lengths ranged from 5 to 52 years with a mean value of around 17 years, with most of the data covering the period from 2000 to 2010 (Panagos et al., 2017). The number of stations in different continents greatly varied from around 5 % (i.e., South America and Africa) to around 48 % (i.e., Europe). Based on the station data and applying the Gaussian process regression model, the global rainfall erosivity map was also prepared (Panagos et al., 2017). Therefore, in the scope of this study, both the station (i.e., point) estimated annual rainfall erosivity and a global rainfall erosivity map (Panagos et al., 2017) were used. The spatial resolution of the global rainfall erosivity map prepared by Panagos et al. (2017) is 30 arcsec (i.e., around 1 km at the Equator). The Rainfall Erosivity Database on the European Scale (REDES) is the predecessor of GloREDa as it was developed in 2015 (Panagos et al., 2015). As the REDES made the monthly erosivity values available (Ballabio et al., 2017), the monthly rainfall erosivity maps of Europe were also used here for the comparison of CMORPH to station-based rainfall erosivity. All data sets are available in the European Soil Data Centre (ESDAC) (Panagos et al., 2012).

2.3 Rainfall erosivity calculation

In order to calculate the annual and monthly rainfall erosivity for each grid cell that is covered by the CMORPH product, the time series with a 30 min time step were extracted from the original CMORPH data set (Xie et al., 2021). For each grid cell covered by CMORPH, a 30 min precipitation time series (mm h⁻¹) was extracted for the 1998–2019 period. The erosive events were defined according to the procedure described in the RUSLE handbook (Renard et al., 1997). Thus, two events were separated in case of less than 1.27 mm of rain within 6 h. Only erosive rainfall events with more than 12.7 mm of rain in total or 6.35 mm in 15 min were considered in the calculations (Kim et al., 2020; Renard et al., 1997). In order to calculate the specific kinetic energy e_B (MJ ha⁻¹ mm⁻¹), the Brown and Foster (1987) equation was applied since this equation was also used by Panagos et al. (2017):

$$\begin{array}{}\text{(1)}& e_{\mathrm{B}}=\mathrm{0.29}\cdot \left[\mathrm{1}-\mathrm{0.72}\cdot \exp \left(-\mathrm{0.05}\cdot I\right)\right],\end{array}$$

where I is rainfall intensity (mm h⁻¹). In order to calculate the annual rainfall erosivity R factor (MJ mm ha⁻¹ h⁻¹ yr⁻¹), the following two equations were also used (Renard et al., 1997):

$$\begin{array}{}\text{(2)}& {\displaystyle }E=e_{\mathrm{B}}\cdot I\cdot \mathrm{\Delta }t,\text{(3)}& {\displaystyle }R={\displaystyle \frac{{\sum }_{n}E\cdot I_{\mathrm{30}}}{N}},\end{array}$$

where E is the kinetic energy of the individual erosive event (MJ ha⁻¹), Δt is the time interval (h), and I₃₀ is the maximum 30 min intensity (mm h⁻¹) of erosive event n, which occurred within a time span of N years. This procedure was repeated for all grid cells covered by the CMORPH product.

2.4 ED and ERA5

The ED concept was first introduced by Kinnell (2010) and was also used in the scope of the enhanced RUSLE approach, named RUSLE2, which led to the improvements in rainfall erosivity mapping (Dabney et al., 2012; Nearing et al., 2017). The ED is defined as the ratio between annual or monthly rainfall erosivity and annual or monthly precipitation (Panagos et al., 2016b). Thus, ED is calculated as the ratio of rainfall erosivity (R) and rainfall depth (P) (Nearing et al., 2017):

$$\begin{array}{}\text{(4)}& \mathrm{ED}={\displaystyle \frac{R}{P}}.\end{array}$$

Since the introduction of the ED, it has been applied in numerous studies (Diodato et al., 2019; Kinnell, 2019; Nearing et al., 2017; Panagos et al., 2016b). The global rainfall erosivity map obtained by Panagos et al. (2017) was used in this study to obtain a global rainfall ED map. For the calculation of rainfall volume for specific years, the ERA5 reanalysis product was used.

ERA5 is one of the latest reanalysis products produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) that provide atmospheric, land-surface, and sea-state data. ERA5 includes a large number of historical observations and provides a long-term solution for ED estimation. The reanalysis data combine the model data and observations across the globe into a complete and consistent data set based on the laws of physics (ERA5, 2021a). Therefore, the ERA5 product is widely used for different purposes (Reder and Rianna, 2021; Sutanto et al., 2020; Tang et al., 2020). The monthly temporal resolution on a single level was used and a horizontal resolution of $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$. The temporal coverage used in this study was from 1979 until 2020. Comparison with CMORPH and GloREDa was made using the 1998–2019 period. Additional information can be found in the existing literature (An et al., 2020; ERA5, 2021b; Tang et al., 2020). ERA5 is updated regularly (i.e., monthly updates), which makes it the best option for the dynamic rainfall erosivity assessment at a global scale using the ED concept. In the case of the ED concept, the annual and monthly ED maps (Ballabio et al., 2017; Panagos et al., 2017) were multiplied by mean monthly or annual precipitation estimates provided by ERA5.

2.5 Data evaluation

The performance of the rainfall erosivity derived using the CMORPH product and ED concept was evaluated using the GloREDa point data set (Panagos et al., 2017). This evaluation was performed for the period 1998–2019 at global, continental, catchment, and local scales. For the latter, point data values (stations) were compared against values derived at this location from both methodologies. At the catchment scale, the HydroSHEDS catchment boundaries at the third level were used (Lehner and Grill, 2013). The idea of using third-level catchment boundaries was to evaluate whether the accuracy of the CMORPH- and ED-derived rainfall erosivity changes with scale (i.e., from global to large regional or even point scale). Moreover, a more detailed comparison was made for Europe since monthly rainfall erosivity maps (REDES) are also available (Ballabio et al., 2017) and were also used for the comparison.

In the data evaluation process, we used the following metrics: Pearson correlation coefficient, percent bias, and Gini coefficient. The Pearson correlation coefficient is a measure of linear correlation between two data sets. The percent bias is a measure of the mean tendency of the modeled data to be smaller or larger than the observed data. The Gini coefficient is a scalar metric that can be derived based on the Lorenz curve and is frequently used in economics to describe the inequality of wealth (Gini, 1914; Lorenz, 1905; Masaki et al., 2014). The Gini coefficient ranges from 0 to 1, where a value close to 1 and 0 indicates significant inequality and no inequality, respectively (Masaki et al., 2014). Thus, the idea behind using the Gini coefficient was to use an additional metric that describes the distribution of rainfall erosivity in the selected area (e.g., the distribution of rainfall erosivity grid cells at catchment or continental scale). Therefore, the Gini coefficient can be used as an indicator of the rainfall erosivity spatial patterns. Figure 1 shows an example of different Gini coefficient values for three examples. In the first one, there are similar grid values, and the Gini coefficient is close to 0. The third example shows significant inequality where the Gini coefficient is close to 1, and the second example represents more diverse grid values with a Gini coefficient of around 0.5 (Fig. 1).

Since the spatial resolutions of the input data sets (CMORPH, GloREDa, and REDES) were not the same, the GloREDa and REDES data (and the ED) were resampled to the same grid system extent and resolution that were used by CMORPH using the mean value (cell-area-weighted) method that is included in the SAGA GIS software (SAGA GIS, 2021). The same applied to the ERA5 product that was also resampled to the same grid system using B-spline interpolation (SAGA GIS, 2021). Therefore, the above-described comparison at global, continental, regional, and point scales was made using the resampled GloREDa and REDES maps (Panagos et al., 2017). A preliminary investigation was done to estimate the resampling effect on the mean global rainfall erosivity. The global mean rainfall erosivity using the GloREDa map (1 km spatial resolution) was 2190 MJ mm ha⁻¹ h⁻¹ yr⁻¹, while in the case of resampled (i.e., mean) data at 10 km resolution, this value is equal to 2260 MJ mm ha⁻¹ h⁻¹ yr⁻¹. Thus, resampling led to an around 3 % difference in the global mean value. However, further aggregation of the GloREDa data led to larger differences.

Figure 1 Gini coefficients for three examples of 10 grid cells with the same mean value (i.e., 100) as an illustration of the added value of the Gini coefficient. Example 1 has similar grid values (i.e., low Gini value), example 2 has more diverse grid values (i.e., Gini value around 0.5), and example 3 has significant inequality (i.e., Gini value close to 1).

2.6 Trends

Based on the annual CMORPH and ED global rainfall erosivity maps for specific years in the period from 1998 to 2019, the Mann–Kendall trend test was also calculated for each grid cell. The Mann–Kendall test is one of the most widely applied tests for the detection of changes in the environmental data (Burn and Hag Elnur, 2002; Rodrigues da Silva et al., 2016). A detailed description of the Mann–Kendall test can be found in the literature (Burn and Hag Elnur, 2002; Hamed, 2008; McLeod, 2011). The objective was to identify areas where the detected trend in the annual rainfall erosivity data was positive or negative with a significance level of 0.05.

3 Results

3.1 Spatial distribution of annual rainfall erosivity

The mean global annual rainfall erosivity using the CMORPH (Fig. 2a) data is 1236 MJ mm ha⁻¹ h⁻¹ yr⁻¹, with a standard deviation of 1895 MJ mm ha⁻¹ h⁻¹ yr⁻¹. The mean global annual rainfall erosivity using the ED approach (Fig. 2b) is 2480 MJ mm ha⁻¹ h⁻¹ yr⁻¹. As can be inferred from Fig. 2 and as further indicated in Table 1, CMORPH and ED approaches both agreed that the highest values of rainfall erosivity at the continental level were estimated for South America, while the smallest ones were estimated for Europe.

Concerning the inequality of CMORPH estimates, the Gini coefficient reflects a high level of rainfall erosivity inequality for Asia, followed by Africa and North America, while the smallest value was observed for Europe (Table 1). Also, with regard to the ED concept, the largest Gini coefficient was obtained for Asia, whereas Europe has the smallest value (Table 1). Both the mean global rainfall erosivity map for the 1998–2019 period based on the CMORPH product (Fig. 2a) and the one developed by using the ED concept and ERA5 will be available in ESDAC (Panagos et al., 2012).

Table 1 Mean, standard deviation, and Gini coefficient of the global rainfall erosivity maps derived using CMORPH and ED.

Figure 2 Mean global rainfall erosivity map for the 1998–2019 period based on the CMORPH product (a) and ED concept using ERA5 (b).

3.2 Temporal trends in rainfall erosivity

Table 2 shows the mean and standard deviation of monthly rainfall erosivity derived using the CMORPH product. One can notice that the highest rainfall erosivity values were obtained in July, followed by August, and the lowest were in November (Table 2).

The temporal trends for both CMORPH- and ED-derived rainfall erosivity data sets were also calculated. The annual rainfall erosivity in the period 1998–2019 ranged from 990 to 1440 MJ mm ha⁻¹ h⁻¹ yr⁻¹ using the CMORPH data (Figs. 3, S1). Using the ED concept for global rainfall erosivity assessment, the mean value ranged from 2380 to 2602 MJ mm ha⁻¹ h⁻¹ yr⁻¹ (Figs. 3, S1) for the period 1998–2019. In addition, the fluctuation of the mean annual erosivity in relation to the ED concept was smaller compared to CMORPH (Fig. 3), a condition which can be related to the fact that the adopted ED concept used a constant ED map for the entire period, while only annual precipitation (i.e., ERA5) changed from year to year.

Figure 3 Trend analysis for the mean and standard deviation (SD) for annual rainfall erosivity (R) using CMORPH and ED.

3.3 Data evaluation

3.3.1 Comparison at a global scale

For most continents, relatively large differences in the mean long-term annual rainfall erosivity between GloREDa and CMORPH were observed, while smaller differences were observed between ED and GloREDa, which could be expected due to the selected ED input data (Table 3). The most significant differences in the case of CMORPH were detected for Africa, South America, and North America (Table 3). As for the ED concept, the most considerable differences were calculated for Asia, Europe, and South America. On the other hand, much better agreement between the CMORPH and GloREDa maps was observed for Europe and partly for Asia (Table 3). In terms of the Gini coefficient, smaller bias values were obtained compared to the mean annual rainfall erosivity (Table 3). Thus, it seems that the distribution of the rainfall erosivity of CMORPH, ED concept, and GloREDa was relatively similar (i.e., smaller bias) despite the fact that the GloREDa map was based on interpolation. It should be noted that the ED provided a better fit to GloREDa compared to CMORPH at most of the continents in terms of the Gini coefficient, which means that the spatial rainfall erosivity patterns are quite similar. This can be regarded as an expected result since the ED indirectly uses the GloREDa data.

Table 2 Global monthly rainfall erosivity values using the CMORPH product. Mean and standard deviation are shown.

Table 3 A comparison between CMORPH-, ED-concept-, and GloREDa-derived global rainfall erosivity at a continental scale.

3.3.2 Comparison at regional scale

The HydroSHEDS catchment boundaries (Lehner and Grill, 2013) at the third level were used to compare the data of CMORPH to GloREDa at the regional scale. Thus, the global land surface was divided into 288 sub-catchments at the third level with a mean catchment area of around 460 000 km². Hence, this can be regarded as an extensive regional-scale investigation. The results demonstrated that the Pearson correlation between the mean annual rainfall erosivity at the sub-catchment level (sub-catchment average values were used) between CMORPH and GloREDa was 0.81 (R²=0.66 with p value <0.01). Moreover, the mean bias was around −50 % in case GloREDa data were considered to be the observed data. In terms of the Gini coefficient, the Pearson correlation coefficient was 0.56 (R²=0.31 with p value <0.01), while the mean bias was equal to 45 %. Therefore, CMORPH yielded more unequal (i.e., larger Gini coefficient) spatial erosivity patterns compared to GloREDa, which was based on interpolation. The spatial interpolations tend to smooth the extreme values (Dodson and Marks, 1997); therefore, Gini is smaller.

The comparison between the ED concept and GloREDa revealed that the Pearson correlation coefficient was equal to 0.95 (R²=0.90 with p value <0.01), and the mean bias was 7 %. Regarding the Gini coefficient, the Pearson correlation coefficient and the mean percent bias were 0.91 (R²=0.83 with p value <0.01) and 3.4 %, respectively. Therefore, GloREDa and ED maps have similar spatial erosivity patterns, which can be regarded as an expected result since the ED indirectly uses the GloREDa information. Furthermore, two examples of good (Fig. 4) and bad (Fig. 5) agreement between CMORPH, ED, and GloREDa are presented.

Figure 4 An example of relatively good agreement between the GloREDa (a), ED (b), and CMORPH (c) maps for parts of eastern Europe and Turkey.

Figure 5 An example of worse agreement between the GloREDa (a), ED (b), and CMORPH (c) maps for the parts of South America.

In Europe, we found the best agreement between CMORPH and GloREDa (Table 3) and the smallest uncertainty in GloREDa. For those reasons and due to the availability of monthly rainfall erosivity estimates (Ballabio et al., 2017), a more in-depth assessment was made for Europe. According to GloREDa, the mean annual rainfall erosivity in Europe (i.e., without Russia) was 668 MJ mm ha⁻¹ h⁻¹ yr⁻¹ with a standard deviation of 429 MJ mm ha⁻¹ h⁻¹yr⁻¹ (Fig. 6). According to the CMORPH product, the mean and standard deviation were equal to 752 and 533 MJ mm ha⁻¹ h⁻¹ yr⁻¹, respectively (Fig. 6). Additionally, the ED concept yielded a mean annual rainfall erosivity value of 804 with a standard deviation of 541 MJ mm ha⁻¹ h⁻¹ yr⁻¹ (Fig. 6). Moreover, the calculated Gini coefficient using all grid cells was 0.31 in all cases (Fig. 6). Thus, it can be seen that the CMORPH product yielded a relatively similar erosivity distribution across Europe (without Russia) compared to GloREDa, which means that all maps have a similar level of inequality (i.e., non-uniform distribution of rainfall erosivity). Slightly larger rainfall erosivity values were obtained using the ED concept. It should be noted that part of these differences can be attributed to the fact that the GloREDa data set mostly used data in the 2000–2010 period (Panagos et al., 2017). Moreover, in some areas (e.g., Italy, Balkan Peninsula, parts of eastern Europe), spatial patterns in all cases were similar, although the CMORPH product and ED concept yielded slightly variable rates (Fig. 6). On the other hand, CMORPH yielded higher annual rainfall erosivity values compared to the GloREDa map in some parts of the United Kingdom and eastern Europe (Fig. 6). Furthermore, CMORPH tends to underestimate areas with relatively high rainfall erosivity such as the Alpine region, Spain, Italy, and other parts of the Mediterranean basin (Fig. 6).

As there are available monthly erosivity data sets in the EU (Ballabio et al., 2017), we compared them to the CMORPH- and ED-derived maps (Table 4). Better agreement between CMORPH and GloREDa for fall compared to winter and summer was found (Table 4). The ED concept yielded higher rainfall erosivity values in almost all months, which also resulted in higher differences at the annual level (Table 4). This could be attributed to the underestimation of the WorldClim V1 map (Beck et al., 2020). In addition, GloREDa has lower values compared to REDES in Europe.

Moreover, Fig. 7 shows monthly erosivity values for selected months where three cases were selected (i.e., underestimation, overestimation, and almost complete agreement between CMORPH and GloREDa). In July, the CMORPH product in the Alpine region generally yielded smaller erosivity values compared to both the monthly erosivity maps prepared by Ballabio et al. (2017) and the ED map (Fig. 7). The same conclusion is reached for other regions such as parts of western Europe or the Iberian Peninsula (Fig. 7). In December, parts of eastern Europe have better agreement among the three maps (Fig. 7). On the other hand, October is the month with the best agreement among the three tested maps (Fig. 7). For October, the best agreement is found in parts of eastern and central Europe, while the worst was detected in parts of the Iberian Peninsula (Fig. 7). The CMORPH-derived rainfall erosivity, in some cases, is more equally distributed (i.e., winter), and in other cases it is more unequally distributed (i.e., summer) compared to GloREDa, while in the case of the ED concept and GloREDa, the derived Gini coefficients are relatively similar throughout the year (Table 4).

Figure 6 Comparison between the GloREDa rainfall erosivity map prepared by Panagos et al. (2017) (a), the CMORPH-derived rainfall erosivity map (b), and the ED-concept-derived map (c) for Europe.

Table 4 Comparison between monthly rainfall erosivity characteristics for Europe using Ballabio et al. (2017), the ED concept, and the CMOPRH product.

Figure 7 Comparison between monthly erosivity maps for Europe prepared by the Ballabio et al. (2017) (a, d, g), CMORPH (b, e, h), and ED concept (c, f, i) maps for July (a, b, c), October (d, e, f), and December(g, h, i).

3.3.3 Comparison at the local scale using GloREDa stations

The station data of GloREDa were also compared to grid cell values at the same location from the derived CMORPH and ED rainfall erosivity maps. The Pearson correlation coefficient between the CMORPH and GloREDa data sets was equal to 0.74 (R²=0.55 with p value <0.01), and the mean bias was equal to −32 % (Fig. 8). In general, the CMORPH product yielded smaller rainfall erosivity estimates, especially for locations where annual rainfall erosivity exceeded 5000 MJ mm ha⁻¹ h⁻¹ yr⁻¹ (Fig. 8). Additionally, CMORPH products tended to overestimate rainfall erosivity in locations near water bodies, which are also the points located above the orange line shown in Fig. 8. A comparison between the ED concept and GloREDa yielded a Pearson correlation coefficient of 0.77 (R²=0.59 with p value <0.01) and a mean bias of 10 %. Similarly, better agreement between the ED concept and GloREDa was detected at a global, continental, or large catchment scales compared to CMORPH versus GloREDa.

Figure 8 Comparison between local (i.e., grid cell values) long-term annual rainfall erosivity from CMORPH and ED and stations' values of GloREDa.

3.4 CMORPH data correction using GloREDa point data

Considering the results and comparisons presented above, the attempt to adjust the CMORPH rainfall erosivity estimates using the estimates of the GloREDa ground station database was made. A similar attempt was also made by Kim et al. (2020) and Wang et al. (2020). In the scope of this study, we developed correction factors (or functions) for each of the Köppen–Geiger climate zones (Peel et al., 2007). The corrections were made both at a global scale using all GloREDa stations and per climate zone (Fig. 9). In Table 5, we propose the best linear function, which can be applied at CMORPH-estimated values in order to be as close as possible to the measured rainfall erosivity values of GloREDa.

Table 5 Correction factors that were developed based on the GloREDa–CMORPH relationship per climate zone. R² is the coefficient of determination. Bold text shows global area.

Therefore, a generic correction linear function that can be used to derive the corrected CMORPH data (CMORPH_COR) for the whole globe can be written as follows:

$$\begin{array}{}\text{(5)}& {\mathrm{CMORPH}}_{\mathrm{COR}}=\mathrm{1.53}\cdot \mathrm{CMORPH}.\end{array}$$

From the results of comparing GloREDa to CMORPH, it is evident that CMORPH underestimates the rainfall erosivity for a factor close to 2 (1.85) in tropical areas where we estimate a high R factor (Panagos et al., 2017). In temperate areas where the R factor is close to the global mean, the CMORPH underestimation is about 1.5, while better agreement can be seen in low-erosivity areas (arid, cold) (Fig. 9). Thus, applying this simple linear transformation can yield better agreement between GloREDa and CMORPH_COR, both at a station scale (Fig. 9) as well as at a global scale (Fig. S2). The same correction was also applied to the global rainfall erosivity map derived using the CMORPH product (Fig. S2) and yielded a global mean rainfall erosivity of 2000 MJ mm ha⁻¹ h⁻¹ yr⁻¹ with a standard deviation of 3314 MJ mm ha⁻¹ h⁻¹ yr⁻¹. Even if one notices much better agreement between CMORPH_COR and GloREDa after the correction, this is relevant only for long-term mean rainfall erosivity assessments. With regard to dynamic rainfall erosivity maps (i.e., for specific years or months), different correction factors should be applied based on the relationship between CMORPH and station rainfall erosivity for specific years. It is worth mentioning that the applied correction can be regarded as a relatively simple one, which could be suitable for global-scale modeling applications. By applying a correction factor to CMORPH, we aim to provide a simple method that uses available remote-sensing data to develop dynamic erosivity values.

Figure 9 Comparison between CMORPH and GloREDa data sets at station scale and proposed correction factors for the whole data set and per climate zone (tropical, cold, temperate, arid, polar). Blue line indicates the linear trend. Red dots are the few outliers (i.e., identified based on Cook's distance) excluded from the correlation.

3.5 Temporal global erosivity trends for the period 1998–2019

The Mann–Kendall test was applied to identify areas with statistically significant (i.e., with a 0.05 significance level) changes (Fig. S3) during the period 1998–2019. According to CMORPH erosivity output, 15 % of the globe has a statistically significant change (Fig. S3). In case of the detected changes, most of the regions show a positive trend rather than a negative one according to the CMORPH product (Fig. S3). Therefore, the positive trend covers 80 % of the area where statistically significant change was detected (i.e., around 12 % of the total area) (Fig. S3), while the remaining 20 % (3 % of the total area) has a negative trend.

On the other hand, the ED concept yielded different results, as around 13 % of the total area shows a statistically significant trend, with positive and negative trends having similar shares (i.e., around 6.5 % each). Consistent trends for using both methods (i.e., CMORPH and the ED concept) were estimated in parts of North America and Asia, while opposite trends were found in parts of Africa, Asia, and Europe. A direct comparison to the study conducted by Bezak et al. (2020) that investigated rainfall erosivity trends in Europe (1961–2018) was not possible since the investigated periods did not overlap.

4 Discussion

The global mean annual rainfall erosivity derived by Panagos et al. (2017) totals about 2190 MJ mm ha⁻¹ h⁻¹ yr⁻¹ (i.e., the initial 1 km cell size map), with a standard deviation of 2974 MJ mm ha⁻¹ h⁻¹ yr⁻¹. The area covered by CMORPH is slightly smaller (between the 60^∘ S and 60^∘ N parallels) than the one covered by the GloREDa map (between the 60^∘ N and ∼75^∘ N parallels, including some parts of Scandinavia, Siberia, and Canada). However, a GloREDa global mean rainfall erosivity value of about 1.8 times higher than the one derived based on the CMORPH data (shown in Sect. 3.1) cannot be explained by the slight difference between the two study areas. It should be noted that Kim et al. (2020) also reported that the CMORPH-derived rainfall erosivity was 1.65 lower than the GloREDa estimates (Panagos et al., 2017) for the United States of America (USA), with some USA regions showing a bias smaller than −80 % (Kim et al., 2020). More specifically, Kim et al. (2020) reported a mean annual value of 1260 MJ mm ha⁻¹ h⁻¹ yr⁻¹ for the USA, while in this study a mean value of 1173 MJ mm ha⁻¹ h⁻¹ yr⁻¹ was calculated using a slightly different methodology (e.g., a different e_b−I equation was applied) and different time period. As the station density was quite low in the case of Panagos et al. (2017) study for Africa and North America and also parts of Asia, this can partly explain the larger differences between CMORH and GloREDa. On the other hand, the largest number of stations was positioned in Europe and also parts of Asia (Panagos et al., 2017), where the agreement between CMORPH and GloREDa was the best (Table 3). Thus, part of the differences between CMORPH and GloREDa can be attributed to the station density used by GloREDa and partly to the issues related to the detection of rainfall by the satellite-based products in mountainous regions (e.g., Stampoulis and Anagnostou, 2012).

In line with what was already discussed by Kim et al. (2020), the insights gained by conducting global analysis suggest that the CMORPH satellite-based rainfall erosivity estimates provide more seamless erosivity distribution without employing interpolation, uniform and good spatial coverage, and 30-time temporal resolution. However, it is also clear that this product has important disadvantages, i.e., overestimated precipitation over water bodies, that detection accuracy in hilly terrains can be problematic, that the accuracy of annual precipitation can be low, and that relative bias at event scale can be significant.

As previously noted, several studies have indicated that the difference between satellite-based products and ground-based precipitation data can be quite significant (Habib et al., 2012; Haile et al., 2015; Jiang et al., 2018). The differences in the rainfall intensity patterns can also be transformed into rainfall erosivity patterns. There were numerous studies published that investigated the accuracy of the CMORPH product in terms of precipitation. For example, Islam et al. (2020) showed that CMORPH overestimated the daily precipitation amount in Australia. A similar conclusion was made by An et al. (2020) for the Yellow River in China or by Wei et al. (2018) for mainland China. Some studies also showed a significant underestimation of CMORPH in winter seasons (Gebregiorgis and Hossain, 2015). Additionally, Palharini et al. (2020) showed that satellite-based products tended to underestimate extreme precipitation, which can have an important effect on rainfall erosivity. Underestimation of extreme rainfall events was also reported in many other studies (Jiang et al., 2019; Rahmawati and Lubczynski, 2018; Stampoulis et al., 2013; Sunilkumar et al., 2015; Wei et al., 2018b). This kind of underestimation can also lead to a negative bias of the satellite-based products compared to the station-based rainfall erosivity. Moreover, Tian et al. (2009) also showed that, in the USA, overestimation was seen for summer (i.e., overestimation of heavy precipitation with an intensity over 40 mm d⁻¹) and underestimation for winter (i.e., missing a significant amount of light precipitation with an intensity lower than 10 mm d⁻¹). Tian et al. (2009) also found out that hit bias (i.e., with respect to satellite-based data and reference data reporting precipitation coincidently) and missed precipitation were the two dominant error sources. A similar conclusion was also made by Jiang et al. (2018), who pointed out the limited detection accuracy of summer thunderstorms by the CMOPRH product in the Shanghai region. These drawbacks of CMORPH can also lead to underestimation or overestimation of the rainfall erosivity by CMORPH. Comparing the CMORPH outputs to the GloREDa-measured erosivity values for almost 3400 ground stations, we found that CMORPH underestimates erosivity in tropical areas by a factor close to 2, while there is better agreement for low-erosivity areas (cold, arid, polar). For the temperate climatic regions, CMORPH underestimates erosivity by a factor close to 1.5.

On the other hand, the only study that, to the best of the authors' knowledge, investigated this satellite-based derived rainfall erosivity (Kim et al., 2020) showed that CMORPH underestimated rainfall erosivity in the USA compared to the GloREDa map. Thus, it is clear that underestimation of the most extreme rainfall events can lead to large differences in the derived rainfall erosivity map. Such characteristics can also lead to relatively large differences in case satellite products are used for flood investigations (Dis et al., 2018). Underestimation of the precipitation amount by the CMOPRH product in southern Europe as shown in the Results section was also indicated by some other studies (Skok et al., 2016). Furthermore, Stampoulis and Anagnostou (2012) also pointed out that satellite-based precipitation product accuracy tended to be lower over the mountainous regions such as the Alpine region (or Andes). This was especially evident during the cold season (Stampoulis and Anagnostou, 2012) and was also highlighted in some other studies (Kidd et al., 2012). Also, other studies pointed a detection problem for winter precipitation and high-intensity rainfall events in some parts of Europe (Stampoulis and Anagnostou, 2012).

Different examples of good and worse agreement among presented rainfall erosivity maps can be seen around the globe (Figs. 4 and 5). Comparing the three maps (i.e., GloREDa, CMORPH, and ED) for parts of eastern Europe and Turkey (Fig. 4), relatively good agreement between all three maps was detected. The main reasons for this good agreement are (a) the relatively large number of stations with a measured R factor which contributed to the GloREDa map in countries such as Romania or Turkey (Panagos et al., 2017), (b) the relatively flat terrain without major mountainous regions in parts of eastern Europe, and (c) the relatively low-medium erosivity (<1000 MJ mm ha⁻¹ h⁻¹ yr⁻¹). By contrast, there are many regions where differences are much larger. An example is the Andes mountain region (Fig. 5), where GloREDa includes only 15 stations in the central part of Chile (Panagos et al., 2017) and gridded precipitation products such as WorldClim also underestimate precipitation (Beck et al., 2020).

The ED (based on the ERA5 and GloREDa data) rainfall erosivity estimates showed better agreement with the GloREDa point estimates, which could be regarded as an expected result due to the selected input data. The largest differences between the ED and the GloREDa estimates were observed in Asia, Europe, and South America because of the precipitation underestimation in mountainous regions such as the Andes, Himalayas, and Alps (Beck et al., 2020). The deviation of ED compared to the GloREDa map could be explained by two main reasons: (a) the difference in the spatial resolution of the GloREDa and CMORPH maps as aggregating the 1 km GloREDa map to the 0.25^∘ that is used by ERA5 yielded a global mean value of 2329 MJ mm ha⁻¹ h⁻¹ yr⁻¹ and that (b) the WorldClim V1 map that was used as input to produce the GloREDa map underestimates the precipitation and that the updated version of the WorldClim map (i.e., V2) yields around 10 % higher annual global precipitation (Beck et al., 2020). It should be noted that the ED concept indirectly uses the GloREDa data for the estimation.

ED has the following advantages compared to the CMORPH approach: (a) the ED concept can be used to prepare dynamic rainfall erosivity maps that have better agreement with GloREDa, and (b) there are no issues with the accuracy near water bodies. On the other hand, ED also has some shortcomings: (a) most of the gridded precipitation data sets underestimate precipitation over mountain regions (Beck et al., 2020) and (b) consider the erosivity–precipitation relationship to be constant, and (c) the rainfall erosivity map is needed as input.

5 Study limitations and other products

The density of stations used to produce the GloREDa map is locally low, especially in the African and South American continents. Obviously, this could have a substantial effect on the produced global rainfall erosivity map (Panagos et al., 2017) and consequently on the results presented in this study since the GloREDa map was used here as a reference. However, to the best of the authors' knowledge, GloREDa is the only global assessment using hourly and sub-hourly rainfall data and the best performing among the global assessments currently available (Panagos et al., 2017). This is due to the coarser time step of other potential global rainfall erosivity sources (Liu et al., 2020). Since the ED concept directly uses the GloREDa map, the results produced by the ED method are directly influenced by the potential shortcomings of GloREDa, and this should be taken into account when making further applications using the ED.

On the other hand, the satellite-based precipitation products have their own sources of uncertainty, as highlighted in the previous sections, and, consequently, CMORH significantly underestimates global rainfall erosivity rates compared to GloREDa. It should be noted that there are other potential products that could have been used to produce global rainfall erosivity maps and that could perhaps yield better results than CMORPH. For example, Multi-Source Weighted-Ensemble Precipitation (MSWEP) uses gauge, reanalysis, and satellite data sources, and it was shown that it outperforms some other products such as CMORPH (Beck et al., 2019a, b). Its spatial resolution is 0.1^∘, and it is available from 1979. Moreover, the Tropical Rainfall Measuring Mission (TRMM) rainfall products can also be used to derive the rainfall erosivity (Li et al., 2020). However, it should be noted that the temporal resolution of these two products is 3 h, which requires a non-standard RUSLE approach to derive the rainfall erosivity (Renard et al., 1997). Thus, alternative approaches for rainfall erosivity estimation are needed. For example, Li et al. (2020) used a modified Brown and Foster equation to calculate the specific kinetic energy and consequently the rainfall kinetic energy. However, this equation was developed based on the case study from China and can therefore be regarded as a local (not global) equation. Thus, applying this equation to the global scale could introduce additional uncertainty to the results. Furthermore, one could also apply the correction (i.e., conversion) factor that was suggested by Panagos et al. (2016b). However, a relatively high value is obtained for the 3 h duration (i.e., a value of 6.6), and the equation used to calculate the conversion factors was only developed for durations up to 1 h. Therefore, applying the correction factors developed by Panagos et al. (2016b) could lead to uncertain results. Thus, there is no globally accepted method for the calculation of the global rainfall erosivity using the 3 h data set. Moreover, these two products also have coarser spatial resolution compared to CMORPH, which also affects the detection of the most extreme rainfall events. Other potential sources (e.g., reanalysis, satellite-based, or combined) with a different temporal and spatial resolution could be additionally tested (Beck et al., 2019a).

6 Conclusions

The global rainfall erosivity was assessed using the CMORPH product and the ED concept. To the best of the authors' knowledge, high temporal (30 min) and spatial resolution satellite-based products such as CMORPH have not yet been used for the development of global rainfall erosivity maps. Past attempts to develop a global erosivity data set based on satellite-based or reanalysis products have used either monthly or daily data. The comparison of the derived maps was performed at global and multiple regional scales using annual and monthly rainfall erosivity values.

The CMORPH product leads to a marked underestimation of annual rainfall erosivity across the globe, with an average value of 1.53 times lower than the GloREDa station-based rainfall erosivity. The agreement between CMORPH and GloREDa estimates varied significantly among continents and climatic zones. While the best agreement was detected for Europe (i.e., percent bias around 10 %), on average, it has relatively low-erosivity values, and a considerably lower performance was observed for Africa and South America (i.e., percent bias around −60 %). Besides having a higher average rainfall erosivity value than Europe, these regions also suffer from a considerably lower number of measurement stations in the GloREDa database. Interpretation of the obtained results suggested that satellite-based products such as CMORPH cannot correctly capture the most extreme rainfall events that contribute to the largest proportion of the annual rainfall erosivity in some parts of the globe. Better agreement was generally detected between the ED concept and GloREDa (i.e., percent bias up to around 20 %), which can be regarded as an expected result since the ED concept indirectly uses the information from GloREDa.

A more detailed comparison was performed for Europe, where an investigation was also performed at a monthly timescale. Some spatial erosivity patterns were well detected by the CMORPH product in some regions, and monthly erosivity values in spring and fall were relatively close to the ones reported by the monthly erosivity maps prepared by Ballabio et al. (2017) (e.g., in parts of eastern and central Europe). Additionally, underestimation and overestimation were detected in summer (percent bias up to −40 %) and winter (percent bias up to 100 %) compared to GloREDa, respectively. On the other hand, the ED concept consistently slightly overestimated GloREDa but yielded better agreement with GloREDa both temporally and spatially than CMORPH (i.e., percent bias was in the range of around 30 %). As mentioned, the ED approach indirectly uses the GloREDa information, but it is to some extent independent as it uses a completely different rainfall data set (i.e., ERA5 instead of WorldClim).

We also estimated a temporal trend analysis at a global scale using high temporal and spatial resolution data of CMORPH for the period 1998–2019. A preliminary trend investigation revealed that around 15 % of the investigated area was characterized by the statistically significant change in the annual rainfall erosivity, while around 80 % of this change was positive (i.e., 12 % of the total area) according to the CMORPH product for the 1998–2019 period. According to the ED concept, 13 % of the area was characterized by a statistically significant trend. In some regions (e.g., parts of South or North America), the detected trends were consistent, while others were not consistent (e.g., parts of Africa or Asia). Thus, detected trends according to CMORPH could indicate that rainfall erosivity has been slightly increasing in 12 % of the globe during the last 2 decades. However, a more detailed investigation using longer time series is needed to confirm or reject this preliminary result.

It should be noted that in case the CMORPH product is used for the preparation of the rainfall erosivity map, it would be further used for soil erosion modeling where an uncertainty assessment should be included in such an investigation, similar to some other scientific disciplines (e.g., Kim et al., 2016; Sun et al., 2018).

Despite the mentioned shortcomings and strong underestimation of the rainfall erosivity in some parts of the globe, the satellite-based precipitation products tend to be an interesting option for the estimation of the rainfall erosivity, especially in regions with limited ground data. However, in some regions and seasons, such products require additional correction to remove bias, which is of course related to the availability of ground-based precipitation. Thus, it is clear that such ground-based high-frequency precipitation measurements are (still) essential for accurate rainfall erosivity estimates; however, one can expect that technological development in the next decades will lead to improved accuracy (Tang et al., 2020) of satellite-based products such as CMORPH. These kinds of products could be used as an input to the dynamic soil erosion models, which could be used by relevant stakeholders. At the moment, alternative approaches such as the ED concept can provide more accurate rainfall erosivity estimates, which can be computed more easily.