Evaluation of GPM IMERG-FR Product for Computing Rainfall Erosivity for Mainland China

Abstract

Satellite precipitation products (SPPs) have emerged as an alternative to estimate rainfall erosivity. However, prior studies showed that SPPs tend to underestimate rainfall erosivity but without reported bias-correction methods. This study evaluated the efficacy of two SPPs, namely, GPM_3IMERGHH (30-min and 0.1°) and GPM_3IMERGDF (daily and 0.1°), in estimating two erosivity indices in mainland China: the average annual rainfall erosivity (R-factor) and the 10-year event rainfall erosivity (10-yr storm EI), by comparing with that derived from gauge-observed hourly precipitation (Gauge-H). Results indicate that GPM_3IMERGDF yields higher accuracy than GPM_3IMERGHH, though both products generally underestimate these indices. The Percent Bias (PBIAS) is −55.48% for the R-factor and −56.38% for the 10-yr storm EI using GPM_3IMERGHH, which reduces to −10.86% and −32.99% with GPM_3IMERGDF. A bias-correction method was developed based on the systematic difference between SSPs and Gauge-H. A five-fold cross validation shows that with bias-correction, the accuracy of the R-factor and 10-yr storm EI for both SPPs improve considerably, and the difference between two SSPs is reduced. The PBIAS using GPM_3IMERGHH decreases to −0.06% and 0.01%, and that using GPM_3IMERGDF decreases to −0.33% and 0.14%, respectively, for the R-factor and 10-yr storm EI. The rainfall erosivity estimated with SPPs with bias-correction shows comparable accuracy to that obtained through Kriging interpolation using Gauge-H and is better than that interpolated from gauge-observed daily precipitation. Given their high temporal and spatial resolution, and timely updates, GPM_3IMERGHH and GPM_3IMERGDF are viable data products for rainfall erosivity estimation with bias correction.

1. Introduction

Soil erosion is one of the most severe environmental challenges in China [1]. Rainfall is the primary external force responsible for the detachment and transport of soil particles from soil masses, resulting in water-induced soil erosion [2]. The potential of rainfall to cause water-induced soil erosion is defined as the rainfall erosivity index in water erosion models, and the product of rainfall energy (E) and 30-min peak intensity (I₃₀) during a storm is considered as the optimal indicator to estimate rainfall erosivity during a storm. The long-term (over 22 years) average of annual total EI₃₀ is defined as the R-factor [3,4,5,6].

Given the high variability and randomness in rainfall intensity, the utilization of high-temporal-resolution precipitation data, such as breakpoint precipitation or minutely level measurements from rain gauges, should be the primary choice for accurately estimating rainfall erosivity [4,5,6,7,8]. However, these high-resolution datasets are frequently sparsely distributed and fail to represent the spatial distribution of rainfall. The deployment of gauge stations is limited by observational costs and practical feasibility, resulting in less dense networks. Furthermore, the higher the temporal resolution of precipitation data, the more challenging and less accessible they become. While gauge observations offer detailed site-specific measurements, their limited spatial distribution and inability to accurately represent spatial variations pose significant challenges [9].

Consequently, exploring alternative methods to depict the characteristics and changes of rainfall erosivity on both temporal and spatial scales becomes necessary. Temporally, models based on lower-resolution precipitation data, ranging from hourly to annual measurements, have been developed [7,8,10,11]. However, it is important to note that a decrease in temporal resolution typically results in reduced model accuracy [8]. Spatially, gridded datasets serve as an alternative to gauge observations, offering a more continuous spatial representation of precipitation. These datasets include gridded precipitation data derived from gauge-interpolation methods [12,13,14,15], and satellite or weather radar retrievals [9,16,17]. However, these datasets also have their limitations [18]. The interpolation-gridded products are not suitable to directly apply in erosivity estimation because the extreme values observed at gauges will be smoothed during the interpolation processes [18,19,20]. As Wang et al. [19] reported, the application of gridded datasets to estimate rainfall erosivity directly in China will lead to an underestimation of 15% to 40% in the R-factor and 25% to 50% in the 10-yr storm EI, respectively.

With the advancement of satellite remote sensing technology and precipitation retrieval algorithms, remote sensing precipitation products have emerged as a novel meteorological observation method, offering expansive spatial coverage. In comparison with gauge observations and interpolated products, which are often subject to delayed updates, satellite precipitation products are frequently updated, providing more timely data. These products are increasingly available at high temporal resolutions, including datasets with 3-h, hourly, and even 30-min intervals [21,22,23]. Consequently, remote sensing precipitation products have frequently been considered in hydrological and soil erosion studies in recent years [23,24,25,26,27,28,29,30,31].

Previous studies have widely explored the applicability of remote sensing precipitation in rainfall erosivity estimation, such as the level-3 satellite precipitation product released by the Tropical Rainfall Measuring Mission (TRMM) [29,30,31,32], the IMERG series satellite inversion precipitation products produced by the Global Precipitation Measurement mission (GPM) [26,27,29], the National Oceanic and Atmospheric Administration’s Climate Precipitation Center morphing technique (CMORPH) [21,33,34], and Climate Hazards Group InfraRed Precipitation with Station Data (CHIRPS) [35,36].

On one hand, many studies have highlighted the strengths of remote sensing precipitation in capturing the spatial distribution of rainfall erosivity [26,29,30,31,34,37,38]. For example, Chen et al. [29] observed a high Pearson correlation coefficient of 0.98 and 0.97, respectively, between annual rainfall erosivity derived from gauge-interpolation results and the daily products of GPM IMERG-F and TRMM at a 0.25° spatial resolution over mainland China.

On the other hand, it has been noted that most remote sensing precipitation products, including TRMM and GPM, tend to underestimate rainfall erosivity [21,26,27,30,37,39]. However, it is crucial to recognize that estimation accuracy varies with the methodologies and datasets employed. For instance, Xing et al. [30] reported that the monthly, seasonal, and annual rainfall erosivity were underestimated by TRMM precipitation products by −26.5% to −32.2%. In contrast, Chen et al. [29] reported a much lower bias in rainfall erosivity estimates using TRMM precipitation products, ranging from 1.1% to 1.2%. This discrepancy may be attributed, at least in part, to differences in the methodologies used in the two studies. Xing et al. [30] based their comparison on rainfall erosivity calculated from 738 stations, while Chen et al. [29] utilized gridded rainfall erosivity data at a 0.25° spatial resolution, which was interpolated from 2417 stations across mainland China. The process of spatial interpolation from gauge observations to large-scale grids at 0.25° resolution can lead to some spatial smoothing of rainfall erosivity values [18,19,20]. Consequently, the evaluation results reported by Chen et al. [29] may not reflect the same level of underestimation as observed in the study by Xing et al. [30].

Overall, the potential of satellite precipitation products in estimating rainfall erosivity has been widely explored. The underestimation of rainfall erosivity estimated using satellite precipitation products has been noticed. Obviously, a bias-correction approach is needed before applying these kinds of satellite precipitation products in soil erosion studies [27]. However, corresponding methods have not been developed or reported in previous studies.

Additionally, previous studies have generally overlooked the assessment of the accuracy of the 10-yr storm EI estimated by satellite precipitation products. This index, known as extreme event rainfall erosivity (EI₃₀) with a return period of 10 years, serves as a critical input in soil erosion models [5]. In the context of global warming, significant variations in the frequency and intensity of extreme weather events have been observed globally [40,41,42,43,44,45]. These changes have led to an increase in rainfall erosivity [46], potentially influencing the erosive index of the 10-yr storm EI [20]. Hence, when evaluating a satellite product for rainfall erosivity estimation, it is essential to consider not only the long-term annual average (R-factor) but also the 10-yr storm EI.

The objective of this study is to develop a sample framework for utilizing GPM satellite precipitation products, which offer high temporal and spatial resolution, comprehensive spatial coverage, and timely updates, to estimate two key rainfall erosivity indices for soil erosion models: the R-factor and 10-yr storm EI. To evaluate their accuracy, gauge-observed hourly precipitation data were collected from 2310 meteorological stations across mainland China. Then, the systematic biases between rainfall erosivity estimated from GPM products and gauge-observed hourly precipitation were addressed by developing these relationships for different water erosion regions. The coefficients from these relationships were subsequently applied as bias-correction factors.

2. Data and Method

2.1. Study Area

This study focuses on mainland China (Figure 1). Considering the variation in erosion driving factors and erosion characteristics across the entire study region, mainland China can be categorized into two major erosion zones: the eastern and southern regions experience water-induced soil erosion, leading to the classification of these areas as the water erosion region (WE); conversely, the western part, including the Tibet-Qinghai Plateau and northwestern China, is primarily influenced by wind and freeze–thaw processes, and is referred to as the non-water erosion region (NWE).

Further subdivision of the water erosion region results in five sub-erosion zones: (1) NR, the northern rocky soil region; (2) SWR, the southwestern rocky mountain region; (3) SR, the southern red soil region; (4) NWL, the northwestern loess plateau region; and (5) NEB, the northeastern black soil region (Figure 1).

2.2. Data Collection

Two satellite precipitation products, GPM_3IMERGHH (30-min resolution at 0.1°) and GPM_3IMERGDF (daily resolution at 0.1°), were obtained from the Global Precipitation Measurement mission (GPM) and downloaded from the NASA website (https://gpm.nasa.gov/data-access/downloads/gpm, accessed on 27 March 2024). These data cover the period from 2000 to 2020. Two rainfall erosivity indices, the long-term average of annual rainfall erosivity (R-factor) and extreme storm rainfall erosivity at 10-yr return level (10-yr storm EI), were calculated using these two satellite products for all grid cells within the study area.

To evaluate their accuracy, gauge-observed hourly precipitation data were collected from 2310 meteorological stations across mainland China. Grid cells from GPM precipitation that contained co-located meteorological stations were selected and extracted, and two rainfall erosivity indices of these grids were calculated and compared with corresponding values based on gauge-observed hourly precipitation located within the grid cells.

The gauge-observed hourly precipitation data were provided by the National Meteorological Information Center (NMIC) of the China Meteorological Administration (CMA), which are observed by siphon or tipping-bucket self-recording rain gauges and should not be used during the cold season to avoid damage to the instruments. Therefore, the hourly observation for most northern stations is only available through the warm season from May to September, but available throughout the year in southern China.

Hourly precipitation data have undergone basic quality control by CMA. Additionally, following the procedures in Wang et al. [45], a second quality control step was applied. Specifically, if the number of missing hours within a given year exceeded 5%, then the entire year was excluded from the analysis. Consequently, for each station, only the years within the 2000–2020 period that met the quality criteria were retained for further analysis. In other words, the data periods used for comparison varied among stations (as shown in Figure 1 and

Table 1. Station information in erosion zones within the study area.

Water Erosion Zones	Shortened Form	Total Stations	Area (104 km2)	Density (/104 km2)	Stations with Data Lengths (Year)
					2~5	6~15	16~20
Mainland China		2310	944.9	2.44	34	502	1774
Non-Water Erosion region	NWE	252	467.1	0.54	14	150	88
Water Erosion region	WE	2058	477.8	4.3	20	352	1686
Northeastern Black Soil Region	NEB	219	108.7	2.0	0	11	208
Northern Rocky Soil Region	NR	495	67.0	7.4	4	66	425
Southern Red Soil Region	SR	643	121.8	5.3	2	97	544
Northwestern Loess Plateau Region	NWL	245	60.0	4.1	14	62	169
Southwestern Rocky Soil Region	SWR	456	120.3	3.8	0	116	340

). Only the overlapping periods between the GPM products and gauge-observed hourly precipitation were considered for comparison and evaluation.

It is important to note that the remaining data may still contain some missing hours, but the percentage of missing data was less than 5%. To mitigate potential underestimations resulting from these data gaps, linear interpolation was performed on the hourly data by comparing the data with corresponding daily precipitation data. Daily data were observed by sample rain gauges and were available throughout the year for all stations used. Daily data were used for the interpolation due to their longer series, fewer missing observations, and generally better data quality when compared to the hourly data. The interpolation procedure was as follows: (1) If the daily observed precipitation amount is zero, then the missing precipitation amount for the corresponding missing period in the hourly data is directly assigned as 0; (2) If the daily observed precipitation amount is not equal to zero, then the observed precipitation values within the hourly data are summed, and the difference between the accumulated hourly precipitation amount and the daily observed precipitation amount is calculated. For example, consider a day with two missing hours, and the daily observation records 20 mm of precipitation. If the total observed precipitation in the available hourly data for that day is 19 mm, the precipitation amount for each of the two missing hours would be set at 0.5 mm.

2.3. Rainfall Erosivity Estimations

Two rainfall erosivity indices (R-factor and 10-yr storm EI) were calculated using two satellite products for all grid cells within the study area in two ways: (1) utilizing the 30-min satellite precipitation product (GPM_3IMERGHH) and the RUSLE2 EI₃₀ index, referred to as GPM-30-EI₃₀; (2) utilizing the daily satellite precipitation product (GPM_3IMERGDF) and the daily erosivity model calibrated by Xie et al. [11], referred to as GPM-Daily-DR. Meanwhile, rainfall erosivity was also estimated using gauge-observed hourly precipitation data coupled with the RUSLE2 EI₃₀ index (Gauge-H), and daily precipitation data coupled with the daily erosivity model (Gauge-DE), respectively.

For gauge-observed hourly precipitation data and the GPM satellite 30-min precipitation product (GPM_3IMERGHH), the continuous data series for each station was divided into individual storm events by the minimum dry inter-event time of 6 h [3,47]. Rainfall erosivity in this study was calculated for each storm using the rainfall erosivity equation recommended in RUSLE2 [6]. For a period m_th with a constant rainfall intensity, the rainfall kinetic energy during the period is given as follows:

$$e_m=0.29[1-0.72\mathrm{e}\mathrm{x}\mathrm{p}(-0.082i_m\left)\right]$$

(1)

where $i_m$ is the rainfall intensity at the m_th time, with a unit of mm h⁻¹, and $e_m$ is the rainfall kinetic energy during the m_th period, with a unit of ${\mathrm{M}\mathrm{J} h}^{-1}{mm}^{-1}$. The total kinetic energy of a rainfall event is

$$E={\sum }_{m=1}^n{e}_m·P_m$$

(2)

where $P_m$ is the rainfall amount during the m_th period, and n represents all periods with a constant rainfall intensity within a rainfall event. The product of the kinetic energy $\mathrm{E}$ and the maximum 30-min rainfall intensity I₃₀ is the rainfall erosivity of a rainfall event, with a unit of MJ mm ha⁻¹ h⁻¹. According to the definition of RUSLE2 [6], the sum of EI₃₀ for all rainfall events in a year is the annual rainfall erosivity, with a unit of MJ mm ha⁻¹ h⁻¹. The average of long-term annual rainfall erosivity is defined as R-factor, with a unit of MJ mm ha⁻¹ h⁻¹ a⁻¹.

The temporal resolution of GPM_3IMERGHH is 30 min; therefore, E and I₃₀ can be estimated directly. But for hourly data, I₃₀ cannot be directly obtained from the hourly observations. Yin et al. [7] reported that the E_h (rainfall energy obtained directly from hourly data), I_1h (the maximum hourly intensity), and (EI₃₀)_h (rainfall intensity index obtained directly from hourly data) have systematic underestimations compared with E, I₃₀, and EI₃₀ that derived from breakpoint rainfall data. Therefore, E_h, I_1h, and (EI₃₀)_h derived from hourly observations were adjusted by multiplying the adjustment factors of 1.105, 1.668, and 1.73 as suggested by Yin et al. [7], respectively, to approximate E, I₃₀, and EI₃₀.

For gauge-observed daily precipitation and the daily satellite precipitation product (GPM_3IMERGDF), rainfall erosivity was calculated using the daily model calibrated by Xie et al. [11] in China:

$$R_{d,n}=0.2686\left[1+0.5412\cos \left(\frac{\mathrm{\pi }}{6}n-\frac{7\mathrm{\pi }}{6}\right)\right]{P_{d,n}}^{1.7265}$$

(3)

where n is the month from 1 to 12; d is the rainy day; and $P_{d,n}$ (MJ·mm·ha⁻¹·h⁻¹) and $R_{d,n}$ (mm) represent the daily precipitation and rainfall erosivity on the $d^{\mathrm{t}\mathrm{h}}$ day of the $n^{\mathrm{t}\mathrm{h}}$ month, respectively. The sum of $R_{d,n}$ around the year is the annual rainfall erosivity, and the long-term average over the study period is the R-factor.

When it comes to the 10-yr storm EI (MJ·mm·ha⁻¹·h⁻¹), for hourly observation and GPM_3IMERGHH, EI₃₀ of each storm event can be determined. Then, the 10-yr storm EI can be calculated by the generalized extreme value distribution. For daily observations and GPM_3IMERGDF, only daily erosivity at the 10-year return level (referred to as 10-yr-daily EI) can be determined. But it has been found to have a strong linear correlation with 10-yr storm EI [48]. According to Yin et al. [48], 10-yr storm EI can be approximated by multiplying 1.17 by the 10-yr-daily EI. To simplify the expression, “R-factor in MJ mm ha⁻¹ h⁻¹ a⁻¹” and “10-year storm EI in MJ mm ha⁻¹ h⁻¹” were simply denoted as “R-factor in SI units” and “10-year storm EI in SI units”.

2.4. Bias-Correction Method

In the Water Erosion region (WE) of mainland China, rainfall erosivity estimates derived from two GPM precipitation products, namely GPM-30-EI₃₀ and GPM-Daily-DR, exhibit systematic biases when compared to estimates based on gauge-observed hourly precipitation. A power equation can be utilized to address these discrepancies:

$$\mathrm{y}=\mathrm{a}·x^b$$

(4)

where x represents the rainfall erosivity values derived from GPM-30-EI₃₀ and GPM-Daily-DR, y represents rainfall erosivity estimated from hourly observations (denoted as Gauge-H), and a and b are parameters. Based on Equation (4), two bias-correction methods were developed and compared to mitigate these discrepancies.

The first method, referred to as “adjust1”, applies a uniform correction across the entire region by establishing and utilizing a correlation between GPM-derived rainfall erosivity and Gauge-H by Equation (4). The second method, referred to as “adjust2”, creates and applies separate correlations for each sub-zone in the region of WE by Equation (4), adapting bias corrections to individual areas.

The efficacy of the two bias-correction approaches was evaluated through a five-fold cross-validation process [49]. The study area, including 2058 stations within the WE region, was segmented into five equal folds. For each iteration, bias-correction equations were derived using four folds and subsequently applied to adjust the data in the remaining fold. The accuracy of these adjustments was assessed by calculating a series of statistical indices, as detailed in Section 2.5. This procedure was repeated across all five folds, and the mean of these statistical indices was employed to quantify the overall accuracy of the bias-correction approaches.

2.5. Statistical Analysis

R-factor and the 10-yr storm EI calculated by two GPM products through two algorithms were evaluated by comparing with that derived from Gauge-H. The evaluation indicators are the Absolute Error (AE, SI units), Relative Error (RE, %), Nash Efficiency Coefficient (NSE), the Correlation Coefficient (CC), the Percent Bias (PBIAS, %), and the Kling-Gupta Efficiency Coefficient (KGE):

$${\mathrm{A}\mathrm{E}}_i={\mathrm{P}}_i-{\mathrm{O}}_i$$

(5)

$$\mathrm{R}\mathrm{E}=100\mathrm{\%}\times ({\mathrm{P}}_i-{\mathrm{O}}_i)/{\mathrm{O}}_i$$

(6)

$$NSE=1-\frac{\sum _{i=1}^n{\left(P_i-O_i\right)}^2}{\sum _{i=1}^n{\left(O_i-\overline{O}\right)}^2}$$

(7)

$$CC=\frac{\sum _{i=1}^n\left(P_i-\overline{P}\right)\left(O_i-\overline{O}\right)}{\sqrt{\sum _{i=1}^n{\left(P_i-\overline{P}\right)}^2}\sqrt{\sum _{i=1}^n{\left(O_i-\overline{O}\right)}^2}}$$

(8)

$$\mathrm{P}\mathrm{B}\mathrm{I}\mathrm{A}\mathrm{S}=100\mathrm{\%}\times (\frac{\sum _{i=1}^n{P}_i}{\sum _{i=1}^n{O}_i}-1)$$

(9)

$$KGE=1-\sqrt{{\left(CC-1\right)}^2+{\left(\frac{{\sigma }_P/\overline{P}}{{\sigma }_O/\overline{O}}-1\right)}^2+{\left(\frac{\overline{P}}{\overline{O}}-1\right)}^2}$$

(10)

in which $O_i$ and $P_i$ in Equations (5)–(9) are the observed and predicted values of a specific station; $\overline{O}$ and $\overline{P}$ in Equation (10) are the average observed and predicted values over all samples. KGE is a comprehensive index that reflects the correlation, variability, and bias between predicted and observed values, with ranges from −∞~1. The higher NSE and KGE suggest the better model estimation. The range for CC is −1 to 1, with −1 showcasing a perfect negative correlation, one indicating a perfect positive correlation, and zero denoting the absence of a linear relationship between the variables. A higher absolute value of CC indicates a stronger correlation. A negative PBIAS represents underestimation, and a positive PBIAS represents overestimation.

The rainfall erosivity estimated from two GPM precipitation products was compared with the interpolation erosivity maps generated by Kriging interpolation methods using gauge observations. The interpolation approach that calculates rainfall erosivity at the individual stations and then interpolates the data into a predefined grid is the most commonly used method [50]. The interpolation results are referred to as Gauge-H-Intp for hourly precipitation data and Gauge-DE-Intp for daily precipitation data, respectively. Additionally, a five-fold cross-validation approach was applied [49], and three indicators of NSE, PBIAS, and KGE were utilized. The satellite precipitation products and gauge-observed precipitation data were processed using MATLAB R2022a (for academic use), while the interpolation and figures were plotted using both RStudio and MATLAB.

3. Results

3.1. Evaluation of Rainfall Erosivity Derived from GPM Precipitation Products

3.1.1. General Accuracy over Mainland China

The spatial distribution of two erosivity indices, namely the R-factor and 10-yr storm EI, using gauge-observed hourly precipitation (referred to as Gauge-H), is illustrated in Figure 2a,b. The R-factor varies significantly across 2310 stations, ranging from 1.6 to 30,052.0 SI units. Similarly, the 10-yr storm EI exhibits a range from 2.3 to 20,999.9 SI units. Both indices display a spatial pattern similar to that of the annual precipitation in mainland China, reflecting the influence of the Asian Monsoon from the coastal to the inland areas.

The evaluation results for the two rainfall erosivity indices computed using two GPM products with two different temporal resolutions and erosivity models (GPM-30-EI₃₀ and GPM-Daily-DR) are represented in Figure 2, Figure 3 and Figure 4. Overall, GPM-Daily-DR exhibits better performance over GPM-30-EI₃₀ in estimating the R-factor across mainland China, especially in the eastern and central regions primarily controlled by water-induced erosion processes (Figure 2c,e). The Nash-Sutcliffe Efficiency (NSE) and Kling-Gupta Efficiency Coefficient (KGE) of GPM-30-EI₃₀ are as low as 0.28 and 0.23 for mainland China, whereas GPM-Daily-DR achieved a significantly higher NSE and KGE of 0.81 and 0.70 (Figure 4). Additionally, the Percent Bias (PBIAS) for GPM-30-EI₃₀ was −55.48%, which is nearly five times higher than that derived from GPM-Daily-DR (PBIAS = −10.86%).

Specifically, GPM-Daily-DR outperforms GPM-30-EI₃₀ in locations where the R-factor estimation from Gauge-H data is below 10,000 SI units. It can be clearly observed from Figure 3a that the scatter points for GPM-Daily-DR align closely with the 1:1 line, particularly when Gauge-H data are below 10,000 SI units. In contrast, GPM-30-EI₃₀ displays systematic biases deviating from the 1:1 line. This observation is corroborated by results presented in Figure 2. The Absolute Errors (AEs) between GPM-30-EI₃₀ and Gauge-H display consistent negative values over the five sub-regions prone to water-induced erosion in eastern China (Figure 2c). In contrast, the AEs between GPM-Daily-DR and Gauge-H exhibit both negative and positive values (Figure 2e). For stations with R-factors exceeding 10,000 SI units, both products tend to underestimate (Figure 3a). It can be inferred from Figure 2 that these stations are mainly located in the southern part of the Southern Red Soil Region (SR).

As for the erosivity index of the 10-yr storm EI, the estimation accuracy of both products is unsatisfying. Figure 2d and Figure 3d show that both products tend to underestimate the 10-yr storm EI. This is further supported by Figure 4g, which shows that the PBIAS for both products is negative. Nevertheless, in comparison, the accuracy of GPM-Daily-DR is still superior to that of GPM-30-EI₃₀ (Figure 4). The NSE and KGE of GPM-30-EI₃₀ are as low as 0.09 and 0.25, whereas the use of GPM-Daily-DR yielded higher values of 0.46 and 0.46. The PBIAS using GPM-Daily-DR for mainland China is also smaller than that using GPM-30-EI₃₀ (Figure 4g). It is noteworthy that both datasets are able to accurately represent the spatial distribution of two erosivity indices, as evidenced by the Correlation Coefficients (CC) shown in Figure 4. For both approaches and indices, the CC values are greater than 0.8 (Figure 4f).

3.1.2. Differences between Water Erosion and Non-Water Erosion Regions

The performance of two GPM products in the Water Erosion (WE) and Non-Water Erosion regions (NWE) shows great differences, and the estimation accuracy within the region of WE is basically consistent with that achieved for mainland China (Figure 3). This is mainly due to the fact that the majority of stations (2085 out of 2340) are situated in WE, where rainfall erosivity exhibits a wide-ranging variation from south to north. In contrast, there are only 252 stations situated in the region of NWE, and their R-factor and 10-yr storm EI estimates are influenced by regional climate dynamics. The scatter points predominantly cluster around the zero point, indicating that both gauge and GPM products yielded very small R-factors and 10-yr storm EI values (refer to Figure 3c,f). These clusters are primarily situated in northwestern China. A few scattered points are sparsely distributed within the range of 5000 to 30,000 SI units, primarily found in the Qinghai-Tibet Plateau. It is important to note that the NSE values in the region of NWE are higher than those in the region of WE (Figure 4a). However, higher NSE values in the NWE region do not necessarily indicate superior prediction accuracy compared to the WE region. Instead, they primarily reflect substantial erosivity differences between the Qinghai-Tibet Plateau and arid northwestern China.

3.1.3. Regional Differences among Five Water Erosion Sub-Zones

The region of WE can be further subdivided into five sub-regions based on the erosion characteristics, and it is noteworthy that the estimation accuracies vary among these different erosion regions. Firstly, GPM-Daily-DR consistently outperforms GPM-30-EI₃₀ in all regions (Figure 4). NSE values for GPM-30-EI₃₀ are predominantly negative across all five sub-regions. In contrast, NSE values for GPM-Daily-DR mostly exceed 0.5 in the majority of regions. The superiority of GPM-Daily-DR is evident in Figure 5, where the red points representing GPM-30-EI₃₀ deviate from the 1:1 line, whereas most values for NEB, NR, SR, and SWR using GPM-Daily-DR (blue points) align closely with the 1:1 line. It is worth noting that the Northwestern Loess Plateau Region (NWL) exhibits the poorest accuracy among all regions, with the lowest NSE and CC values observed for both GPM-30-EI₃₀ and GPM-Daily-DR (Figure 4).

Compared to the estimation of the R-factor, the accuracy in estimating the 10-yr storm EI in these five erosion sub-zones is less satisfactory (Figure 4e–h). NSE values for WE using GPM-30-EI₃₀ and GPM-Daily-DR are −0.13 and 0.36, respectively. Across all sub-regions, the NSE values are either below or near zero. Although GPM-Daily-DR performs slightly better than GPM-30-EI₃₀ in estimating the 10-yr storm EI (Figure 5f–j), the accuracy remains notably low. Nevertheless, the spatial distribution of rainfall erosivity is still broadly preserved, as indicated by the Correlation Coefficients (CC) (Figure 4f).

Both GPM precipitation products tend to underestimate these two rainfall erosivity indices compared to those derived from hourly gauge observations (Figure 4 and Figure 5). The degree of underestimation is higher when using GPM-30-EI₃₀ to estimate the R-factor compared to GPM-Daily-DR. Similarly, both GPM products exhibit underestimation in estimating 10-yr storm EI. Based on results achieved here, it can be concluded that GPM-30-EI₃₀ is not suitable for the direct estimation of both the R-factor and the 10-yr storm EI, and GPM-Daily-DR is not suitable for the direct estimation of the extreme erosivity index, which will lead to significant uncertainties and biases.

3.2. Evaluation of the Bias-Correction Method

Rainfall erosivity estimated by GPM-30-EI₃₀ and GPM-Daily-DR in the Water Erosion region (WE) shows systematic biases with that derived from gauge-observed hourly precipitation (Figure 3b,e). The biases can be addressed by a power relationship in Equation (4). Based on this equation, two bias-correction approaches named “adjust1” and “adjust2”, were developed and evaluated (Figure 6 and

Table 2. Comparison of R-factors and 10-yr storm EI estimated using bias-corrected GPM-30-EI₃₀ and GPM-Daily-DR versus gauge-observed hourly precipitation in Water Erosion region (WE), using two bias-correction methods (adjust1 and adjust2) through a five-fold cross validation.

Region	Nash-Sutcliffe Efficiency				Percent Bias (%)				Kling-Gupta Efficiency
	GPM-30-EI30		GPM-Daily-DR		GPM-30-EI30		GPM-Daily-DR		GPM-30-EI30		GPM-Daily-DR
	Adjust1 a	Adjust2 b	Adjust1	Adjust2	Adjust1	Adjust2	Adjust1	Adjust2	Adjust1	Adjust2	Adjust1	Adjust2
R-factor
NEB	0.73	0.74	0.56	0.68	−7.33	0.39	−12.21	0.11	0.81	0.83	0.69	0.77
NR	0.52	0.73	0.34	0.74	−10.59	0.36	−23.82	0.23	0.65	0.77	0.55	0.79
SR	0.79	0.80	0.81	0.83	2.75	−0.34	3.08	−0.72	0.78	0.87	0.78	0.9
NWL	0.15	0.36	0.29	0.47	6.98	0.34	−1.27	0.3	0.55	0.41	0.66	0.54
SWR	0.61	0.61	0.68	0.67	−4.5	0.17	2.16	0.05	0.69	0.69	0.72	0.74
WE	0.87	0.88	0.88	0.90	−1.39	−0.06	−2.52	−0.33	0.93	0.92	0.93	0.93
10-yr storm EI
NEB	0.41	0.41	0.34	0.41	4.81	0.13	11.51	0.1	0.49	0.52	0.56	0.51
NR	0.07	0.37	0.18	0.37	−6.77	0.38	−4.55	0.33	0.55	0.42	0.58	0.43
SR	0.54	0.54	0.50	0.51	−1.01	−0.23	−4.26	0.01	0.6	0.64	0.62	0.6
NWL	−0.21	0.12	−0.30	0.19	19.26	0.15	26.11	0.14	0.21	0.07	0.3	0.2
SWR	0.40	0.40	0.41	0.44	3.21	0.1	7.24	0.22	0.49	0.48	0.56	0.52
WE	0.66	0.68	0.65	0.67	−0.4	0.01	0.01	0.14	0.74	0.76	0.73	0.75

^a Adjust1 treats the entire region uniformly, establishing and applying a correlation between rainfall erosivity from GPM data and hourly gauge observations. ^b Adjust2 creates and applies separate correlations for each sub-zone, adapting bias corrections to individual areas.

). The two methods demonstrate comparable accuracy when applied to the entire region of WE (Figure 6 and Figure 7), but “adjust2” slightly outperformed “adjust1” for both GPM products, showing higher Nash-Sutcliffe Efficiency (NSE) and Kling-Gupta Efficiency Coefficient (KGE) values, and lower Percent Bias (PBIAS) (Table 2).

After bias-correction, an obvious improvement is observed (Table 2). Specifically, with the application of “adjust2”, the NSE for the R-factor derived from GPM-30-EI₃₀ improves significantly, rising from 0.16 to an impressive 0.88. Similarly, for GPM-Daily-DR, the NSE improves from 0.79 to 0.90. The PBIAS for GPM-30-EI₃₀ decreases from −55.48% to −0.06%, and for GPM-Daily-DR, it decreases from −10.86% to −0.33%.

These improvements are notable for the 10-yr storm EI as well, where the bias-correction leads to significant enhancements. Specifically, when “adjust2” is applied, the NSE for GPM-30-EI₃₀ improves from 0.09 to 0.68, while for GPM-Daily-DR, it improves from 0.46 to 0.67. After the bias-correction procedure, GPM-30-EI₃₀ and GPM-Daily-DR achieve comparable accuracy, with KGE values of 0.92 and 0.93, respectively. The PBIAS for GPM-30-EI₃₀ and GPM-Daily-DR also decrease markedly, to as low as 0.01% and 0.14%, respectively. These results observed for R-factor and 10-yr storm EI indicate that the significant underestimation previously observed with the direct utilization of the two GPM products was greatly mitigated (Figure 4).

Improvements are observed across all five sub-zones in the Water Erosion region (Table 2, Figure 7). Before applying the bias-correction methods, the Nash-Sutcliffe Efficiency (NSE) of R-factors for all five sub-zones is below 0. These values then improve to positive NSE values for all sub-zones after bias correction. Specifically, the accuracy enhancements are notable in the Northwestern Loess Plateau Region (NWL), where NSE using GPM-30-EI₃₀ and GPM-Daily-DR rose from negative values to 0.36 and 0.47, respectively, with the application of “adjust2”.

The fitted values of parameters a and b from Equation (4) are presented in

Table 3. Equation parameters used to adjust the systematic biases of rainfall erosivity derived from two GPM precipitation products.

Region	GPM-30-EI30			GPM-Daily-DR
	a	b	R2	a	b	R2
R-factor
NEB	4.6297	0.9150	0.76	1.8165	0.9161	0.69
NR	53.3326	0.5809	0.73	6.5847	0.7892	0.74
SR	1.4962	1.0409	0.80	0.0155	1.4869	0.83
NWL	35.4372	0.5538	0.36	5.1520	0.7366	0.47
SWR	4.5418	0.9155	0.62	0.0625	1.3356	0.68
WE	3.9466	0.9278	0.87	0.1064	1.2746	0.88
10-yr storm EI
NEB	14.6082	0.7442	0.43	5.8516	0.7966	0.43
NR	228.8218	0.3553	0.37	64.4091	0.5007	0.38
SR	13.8378	0.7513	0.55	3.8977	0.8824	0.51
NWL	100.2528	0.3824	0.13	22.6267	0.5629	0.20
SWR	20.0841	0.6964	0.41	5.8582	0.8109	0.45
WE	19.1163	0.7079	0.66	3.0119	0.9095	0.65

. The fitting accuracy for both datasets is generally similar, with R² values exceeding 0.62, except for the NWL region, which exhibits a lower fitting accuracy. The fitting accuracy for the 10-yr storm EI is lower compared to that of R-factor. Among all regions, the NWL region shows the poorest fitting accuracy for both GPM-30-EI₃₀ and GPM-Daily-DR.

Based on the findings presented above, it can be concluded that both GPM-30-EI₃₀ and GPM-Daily-DR hold potential for rainfall erosivity estimation. Both datasets can effectively benefit from the bias-correction procedure, resulting in comparable accuracy levels between GPM-30-EI₃₀ and GPM-Daily-DR. Nonetheless, even after bias correction, GPM-Daily-DR consistently exhibits slightly superior performance when compared to GPM-30-EI₃₀.

4. Discussion

4.1. Comparison with Rainfall Erosivity Derived from Interpolation

In previous studies, one of the most commonly used methods to generate R-factor maps in a study area, including regions lacking gauge observations, has been to initially calculate the R-factor for individual stations and subsequently apply interpolation methods to extend these values to specific grid cells covering the entire study area [10,50,51,52,53]. This method can provide more accurate results when compared to using gridded precipitation products generated through the interpolation of daily precipitation data collected at specific locations. This is primarily because the R-factor, which represents a long-term average annual value, exhibits lower spatial variability than daily precipitation data [22,54]. However, the application of the interpolation procedure also introduces some uncertainty into the rainfall erosivity estimation. The interpolation accuracy primarily depends on the density of stations utilized and the interpolation method applied. For instance, in areas with sparse station coverage, the accuracy of estimated rainfall erosivity might be impacted [55].

This study compares rainfall erosivity estimated by GPM-30-EI₃₀ and GPM-Daily-DR with results from the traditional Kriging interpolation method. Using a five-fold cross-validation, the study assesses the interpolation accuracy of erosivity based on gauge-observed hourly precipitation (Gauge-H-Intp) (

Table 4. The interpolation accuracy of rainfall erosivity derived by gauge-observed hourly precipitation (Gauge-H-Intp) and daily precipitation (Gauge-DE-Intp) through a five-fold cross validation.

Region	R-Factor						10-yr Storm EI
	Gauge-H-Intp			Gauge-DE-Intp			Gauge-H-Intp			Gauge-DE-Intp
	NSE	PBIAS (%)	KGE	NSE	PBIAS (%)	KGE	NSE	PBIAS (%)	KGE	NSE	PBIAS (%)	KGE
NEB	0.86	0.56	0.88	0.72	12.95	0.76	0.65	1.12	0.71	0.55	−8.54	0.68
NR	0.84	0.22	0.88	0.75	−7.56	0.85	0.58	0.46	0.64	0.32	−16.46	0.62
SR	0.82	−0.24	0.85	0.59	−12.90	0.57	0.61	−0.63	0.65	0.38	−18.68	0.50
NWL	0.71	1.34	0.79	−0.51	43.07	0.47	0.51	0.88	0.59	0.47	0.78	0.53
SWR	0.69	0.30	0.75	0.59	−0.70	0.55	0.50	0.33	0.56	0.32	−18.3	0.52
WE	0.90	0.03	0.92	0.81	−6.90	0.74	0.74	−0.03	0.78	0.61	−16.77	0.67

). After bias correction, the accuracy of GPM-30-EI₃₀ and GPM-Daily-DR proves comparable to that of Gauge-H-Intp. In the Water Erosion region (WE), the Kling-Gupta Efficiency (KGE) for R-factors estimated by GPM-30-EI₃₀ and GPM-Daily-DR using the adjust2 method is 0.92 and 0.93, respectively, similar to Gauge-H-Intp’s 0.92. However, in sub-zones within WE, Gauge-H-Intp demonstrates higher accuracy than the bias-corrected estimates from GPM products, especially for the Northwestern Loess Plateau Region (NWL) (Table 4).

Although high-temporal-resolution precipitation data can accurately depict rainfall processes and capture intensity variations, they present considerable challenges in data collection, especially in areas with limited gauge records, such as in Africa [23,28]. Therefore, daily precipitation data and daily erosivity models have been widely employed for generating rainfall erosivity maps when high-spatiotemporal-resolution gauge observations are unavailable [53,56,57,58,59]. However, using daily data and daily erosivity models introduces more uncertainty into the mapping process for rainfall erosivity than does using hourly data for interpolation, as the models’ accuracy tends to decline when the adopted data resolution decreases from hourly to daily [8].

To illustrate, the study also evaluated interpolation results for daily data-derived rainfall erosivity using erosivity models [11] based on daily data (Gauge-DE-Intp) through five-fold cross-validation. Table 4 reveals that the interpolation results are notably lower when compared to the estimation accuracy derived from the two GPM datasets. The KGE of Gauge-DE-Intp is 0.74, which is significantly lower than the predictive accuracy of bias-corrected GPM-30-EI₃₀ (KGE = 0.92) and GPM-Daily-DR (KGE = 0.93), as well as the interpolation results using hourly data (Gauge-H-Intp, KGE = 0.92).

Figure 8 further compares the spatial distribution of rainfall erosivity estimated with Kriging interpolation methods based on daily precipitation (Gauge-DE-Intp) and bias-corrected GPM-30-EI₃₀ and GPM-Daily-DR. The comparison reveals that the R-factors derived from GPM products exhibit a spatial distribution pattern identical to those derived from Gauge-DE-Intp. However, when considering the 10-yr storm EI, the erosivity maps generated using GPM products exhibit a smoother spatial pattern, whereas those produced by Gauge-DE-Intp display patches. This patchiness is evidently influenced by the locations of the stations used in the interpolation process, highlighting a key limitation of interpolation methods.

Both GPM precipitation products demonstrate great potential for estimating R-factor and 10-yr storm EI, with accuracy levels comparable to those attained through the interpolation using hourly gauge observations and better than those derived from daily gauge observations.

4.2. Regional Differences

The estimation accuracy of rainfall erosivity indices derived from two GPM products shows obvious regional differences (Figure 3, Figure 4 and Figure 5). Among all regions, the Northwestern Loess Plateau Region (NWL) exhibits the poorest accuracy, both before and after bias correction (Figure 6 and Figure 7). This phenomenon may be related to the lower precision of precipitation estimates within this region as provided by GPM products [22]. The NWL region is characterized by complex topography and significant elevation differences, ranging in elevation from 92 m to 5020 m [60]. The complex topography in this region may affect the accuracy of precipitation estimation, whereas the eastern regions of China, with higher precipitation intensity and relatively flatter topography, tend to have better precipitation estimation accuracy [31,45,61].

Moreover, it is noteworthy that the interpolation results derived from daily precipitation data and daily erosivity models (Gauge-DE-Intp) consistently demonstrate the lowest accuracy among the five regions. The average Nash-Sutcliffe Efficiency (NSE) value, determined through a five-fold cross-validation, is notably low at −0.51, while the NSE values for the other regions either equal or exceed 0.59 (Table 4). This discrepancy may be attributed to the limited data available for developing the daily models, with only two stations providing 1-min precipitation data used in the models’ calibration [11]. Therefore, the biases introduced by the daily erosivity model may be more pronounced in this region compared to others.

4.3. Comparison with Previous Studies

Satellite remote sensing precipitation products have emerged as alternative sources of gauge observations in the past decade and have been widely utilized in mapping rainfall erosivity [26,27,29,30,31,32,34,36]. Most studies have reported that these satellite products effectively capture the spatial distribution of rainfall erosivity [26,29,31,34], which is consistent with the results achieved in this study. However, the evaluation accuracy of rainfall erosivity magnitude derived from satellite products varies across studies, mainly due to differences in the satellite products used and the reference data employed in previous studies [26,29,31,34,36].

Particularly, in a study by Chen et al. [29], the performance of GPM_3IMERGDF in estimating erosivity was evaluated, which is one of the satellite products utilized in our study as well. Chen et al. [29] found good performance in rainfall erosivity estimation at the annual average, seasonal, and monthly scales. The Nash-Sutcliffe Efficiency (NSE) values for all eight river basins in China were above 0.89, which is higher than the direct evaluation results obtained in our study (Figure 4). It is important to note that the comparison in Chen et al. [29] between erosivity maps derived from satellite products and reference erosivity maps was conducted at a spatial resolution of 0.25°. The reference erosivity maps were interpolated using daily precipitation data and a daily erosivity model from 2471 stations to cover mainland China at a resolution of 0.25°, while the precipitation data from GPM was upscaled from 0.1° to 0.25°. This upscaling process may have resulted in the smoothing of rainfall erosivity derived from GPM precipitation data, partially explaining the discrepancies between the findings of Chen et al. [29] and the current study. Additionally, it is important to consider that the reference data used in Chen et al. [29] were interpolated maps generated from daily precipitation data and erosivity models. As discussed in Section 4.1, this interpolation procedure with daily precipitation may introduce higher uncertainty in rainfall erosivity estimation than using hourly rainfall erosivity. This may be the second factor leading to the differences between the results of Chen et al. [29] and the current study.

In previous studies, satellite precipitation products have been widely used for mapping the spatial distribution of rainfall erosivity. However, the results of this study demonstrated that the direct utilization of GPM data with a 30-min resolution (GPM-30-EI₃₀) is inappropriate due to its low accuracy, resulting in obvious underestimation of both the R-factor and 10-yr storm EI. On the other hand, the estimation accuracy of GPM-Daily-DR is considerably higher than that of GPM-30-EI₃₀, except for the Northwestern Loess Plateau Region (NWL). Despite these differences, both products show potential for improvement and can be bias-corrected (Table 2). Moreover, after applying bias correction, both GPM-30-EI₃₀ and GPM-Daily-DR can be utilized for mapping rainfall erosivity, achieving accuracy comparable to the interpolation results based on hourly gauge observations (Table 4).

5. Conclusions

Remote sensing precipitation products have emerged as a new observation method and are widely used in rainfall erosivity estimation. However, previous studies have revealed that rainfall erosivity estimated using satellite precipitation products was underestimated, but corresponding bias-correction methods have not been reported. This study evaluated the potential of two satellite precipitation products—GPM_3IMERGHH (30-min resolution) and GPM_3IMERGDF (daily resolution), both at 0.1°—to estimate two important rainfall erosivity indices in mainland China: the R-factor and the 10-year storm EI. Due to the different temporal resolutions of these two datasets, different erosivity models were applied. For GPM_3IMERGHH, erosivity was estimated using the EI₃₀ index (referred to as GPM-30-EI₃₀); for the GPM_3IMERGDF, a daily erosivity model was used (referred to as GPM-Daily-DR). The following results were achieved:

(1) Both products can preserve the spatial pattern of rainfall erosivity when compared with that derived from gauge-observed hourly precipitation (Gauge-H). The correlation coefficients (CC) are greater than 0.93 for R-factors and greater than 0.81 for 10-yr storm EI, respectively. But both erosivity indices derived from GPM products tended to underestimate the two erosivity indices, and GPM-Daily-DR yielded higher accuracy compared to GPM-30-EI₃₀. The Percent Bias (PBIAS) was −55.48% for the R-factor and −56.38% for the 10-yr storm EI when using GPM_3IMERGHH, and decreased to −10.86% and −32.99% when using GPM_3IMERGDF. This indicated that the direct utilization of GPM-30-EI₃₀ for estimating both the R-factor and the 10-yr storm EI will lead to significant uncertainties and biases.

(2) Rainfall erosivity estimated by GPM products in the Water Erosion region (WE) of mainland China showed systematic biases compared to that derived from Gauge-H, and bias-correction methods were developed. The bias-correction methods could obviously improve the estimation accuracy of the two erosivity indices, and the differences between the two datasets were reduced. The PBIAS using GPM_3IMERGHH decreased to −0.06% and 0.01%, and that using GPM_3IMERGDF decreased to −0.33% and 0.14%, respectively, for the R-factor and 10-yr storm EI.

(3) Moreover, the rainfall erosivity estimated by both GPM products after bias correction showed accuracy comparable to those achieved through Kriging interpolation using Gauge-H (R-factor: PBIAS = 0.03%; 10-yr storm EI: PBIAS = −0.03%) and better than those derived from gauge-observed daily precipitation (R-factor: PBIAS = −6.9%; 10-yr storm EI: PBIAS = −16.77%).

Based on the results achieved here, it can be concluded that both GPM_3IMERGDF and GPM_3IMERGHH hold great potential for application in rainfall erosivity estimation after a bias-correction procedure. This research offers valuable insights and recommendations regarding the application of remote sensing precipitation products in the estimation of rainfall erosivity.