Rainfall erosivity estimation using gridded daily precipitation datasets

Rainfall erosivity is one of the most important factors incorporated into the empirical soil erosion models USLE 10 (Universal Soil Loss Equation) and RUSLE (Revised Universal Soil Loss Equation). Gridded precipitation datasets have been widely used in the estimation of rainfall erosivity, whereas biases due to scale differences between gridded data and gauge data have been ignored. Based on daily precipitation observations from over 2000 stations in China, as well as four widely used gauge-based gridded daily precipitation datasets, CPC, GPCC, CN05.1 and NMIC, this study compared the probability density functions (PDFs) of the gridded and gauge datasets using the skill score method, quantified the bias of rainfall erosivity 15 (including the R-factor and 1-in-10-year event rainfall erosivity) estimated using the gridded daily precipitation datasets from that estimated using the gauge daily precipitation dataset based on the area reduction factor (ARF) method, and established correction factors for rainfall erosivity maps generated from gridded datasets. The results showed that the gridded daily data reduced the frequency of no-rain days and the intensity of heavy precipitation. In the eastern part of China, the grid-estimated R-factor values were underestimated by 15–40 % compared with the gauge-estimated values, and the grid-estimated 1-in-1020 year event rainfall erosivity values were underestimated by 25–50 %, whereas in the western part of China, noticeable random errors were introduced. The lower probability and intensity of the daily precipitation larger than the 90th percentile in the gridded datasets were mainly responsible for the underestimation. CN05.1 was the most-recommended among the four datasets, as it had the lowest mean relative error (MRE), and the accuracy was higher for the eastern part of China than for the western part of China. The MREs were 16.1 % and 25.1 % for the R-factor after applying correction factors of 1.708 and 1.010, 25 respectively, for the eastern and western part of China. The 1-in-10-year event erosivity had larger correction factors and MREs than did the R-factor, with the MREs being 22.1 % and 27.2 % after applying correction factors of 1.959 and 1.880, respectively, for the eastern and western part of China. This study pointed out that in the applications of gridded precipitation datasets, the empirical models established based on gauge precipitation data should not be used directly for the gridded data, or a bias correction process needed to be considered for the model outputs. 30 https://doi.org/10.5194/hess-2020-633 Preprint. Discussion started: 10 December 2020 c © Author(s) 2020. CC BY 4.0 License.

GCMs, the gridded precipitation outputs need to be spatially downscaled to points, from which the future rainfall erosivity values can be calculated and then interpolated to generate an R-factor map Shiono et al., 2013).
The major objectives of this study are to (1) quantify the biases of precipitation metrics and rainfall erosivity values estimated 100 with four gauge-based gridded daily precipitation datasets (CPC, GPCC, CN05.1 and NMIC) based on gauge daily precipitation observations and (2) develop correction factors of the grid-estimated rainfall erosivity to minimize the bias. First, we compared four gridded datasets with gauge daily precipitation observations for 2243 meteorological stations in China in terms of the differences in the PDFs of the daily precipitation values, precipitation metrics and rainfall erosivity values estimated from a power function daily model. Second, a set of correction factors was established by relating the rainfall 105 erosivity values from the four gridded daily datasets to those from the state-of-the-art rainfall erosivity map for mainland China generated based on hourly data from 2396 stations (Yue et al., 2020b). These correction factors can be useful when applying gridded daily precipitation datasets in the estimation of rainfall erosivity and soil erosion in China.

Gauge-observed daily precipitation data
Daily precipitation data from 2481 meteorological stations, archived by the China Meteorological Administration (CMA), were used in this study. The longest overlapping period between the gauge data and four sets of gridded data was from 1 January 1982 to 31 December 2014. If the precipitation data were missing for more than six days in a month during the overlapping period, the data from that station were excluded from this study (CMA, 2003). Finally, the gauge observations of 115 2243 stations from 1982 to 2014 were selected.

Gridded daily precipitation data
In this study, four gridded daily datasets were used: two worldwide gridded daily precipitation datasets, the CPC (Climate Prediction Center) Unified Gauge-based Analysis of Global Daily Precipitation (hereinafter referred to as CPC; Xie et al., 2007) and the GPCC (Global Precipitation Climatology Centre) First Guess Daily product (hereinafter referred to as GPCC; 120 Schamm et al., 2014); and two nationwide datasets, CN05.1, developed by China Meteorological Administration (Wu et al., 2013), and the China Gridded Daily Precipitation Product, developed by the National Meteorological Information Center (hereinafter referred to as NMIC; Shen et al., 2010). Details about the four datasets are shown in Table 1.  Xie et al., 2007Schamm et al., 2014Wu et al., 2013Shen et al., 2010 https://doi.org/10.5194/hess-2020-633 Preprint. The skill score provided a quantitative measure of the similarity of the PDFs. The similarity of the PDFs of the daily precipitation amounts between gauge and gridded data was measured with the skill score ( ) following the methods outlined in Perkins et al. (2007): where n is the number of bins used to calculate the PDF for a given station (the bin width was set to 0.5 mm, and n was determined by bin width and maximum precipitation); is the fraction of values in a given bin from the gridded data; and is the fraction of values in the same bin from the gauge-observed data. and were calculated when the daily precipitation amount was ≥ 0.1 mm. This metric, being the sum of the minimum value of two frequency distributions, measured the common area between two PDFs. The ranged from 0 to 1; the closer its value was to 1, the more similar the two PDFs were.

Definition of the precipitation metrics, rainfall erosivity, and areal reduction factors
Four metrics were used to measure the performances of the different datasets for both the average precipitation conditions and the frequency or intensity of extreme events (  (Alexander et al., 2006). 140 To estimate rainfall erosivity, we used a model of erosivity as a power function of the daily precipitation amount with a coefficient that varies seasonally as a sinusoidal function of the month (Xie et al, 2016): where m is the month, from 1 to 12; , is the rainfall erosivity value for the day in the mth month; and , is the daily precipitation amount for the day in the mth month and was no less than 10 mm, which was the threshold of erosive daily 145 precipitation. The daily rainfall erosivity values were accumulated for each month, and the R-factor was defined as the sum of the average monthly rainfall erosivity values (MJ· mm· hm -2 ·h -1 ·a -1 ) over the study period.
In addition, the 1-in-10-year event EI30, which was the event rainfall erosivity value with different return periods of 10 years (MJ· mm· hm -2 ·h -1 ), was also considered. Though from daily precipitation records, only daily rainfall erosivity could be generated, Yin et al. (2019) found that there was a good linear relationship between event and daily rainfall erosivity values 150 for corresponding return periods (the 1-in-10-year event EI30 was approximately 1.17 times as much as the 1-in-10-year daily rainfall erosivity value). The generalized extreme value distribution (Smith, 2001) was first used to fit the annual series of the maximum daily rainfall erosivity values, and the 1-in-10-year daily rainfall erosivity value was then generated based on the calibrated extreme value distribution. Finally, the 1-in-10-year daily rainfall erosivity value was multiplied by the conversion factor of 1.17 to obtain the 1-in-10-year event EI30. 155 The precipitation metrics and rainfall erosivity values were first calculated for individual stations from gauge data and then interpolated into grids consistent with the four gridded datasets in terms of the spatial resolutions (Table 1), using ordinary kriging to obtain and . The metrics and rainfall erosivity values were calculated directly from the gridded data to obtain and . Used to measure the difference between the two methods, Area Reduction Factor (ARF) was defined in this study as the ratio of over ( over ) (Fowler et al., 2005;Chen and 160 Knutson, 2008): = , for precipitation metrics; = , for rainfall erosivity, In all the spatial interpolation processes, leave-one-out cross validation was used to assess the accuracy of the interpolation method. PBIAS (percent bias), NSE (Nash-Sutcliffe coefficient of efficiency) and RMSE (root-mean-square error) were calculated to examine the differences between the actual values calculated from gauge daily precipitation data and the values 165 predicted by the kriging interpolation: https://doi.org/10.5194/hess-2020-633 Preprint. Discussion started: 10 December 2020 c Author(s) 2020. CC BY 4.0 License.
where is the rainfall erosivity value of the ith site calculated from the gauge-observed data and is the rainfall erosivity 170 value of the ith site predicted by the kriging method from the cross-validation process. The closer NSE was to 1, and the closer PBIAS and RMSE were to 0, the better the interpolation model was. A negative PBIAS value indicated that the interpolation model overestimated the values at the sites, and conversely, a positive PBIAS value indicated that the model underestimated the values at the sites.
The ARFs in the eastern part of China, where the climate is relatively humid and soil erosion is mainly caused by water, may 175 be different from those in the western part of China, where the climate is drier and soil erosion is mainly caused by wind and freeze-thaw processes. In addition, the different densities of meteorological stations could also influence the results. Therefore, the study area was divided into two parts (Chen and Zhu, 1989), and the evaluation was carried out in the two parts separately.

PRCPTOT
Mean annual total precipitation from wet days WD Wet days: mean annual total days when precipitation ≥ 1 mm

R95pTOT
Mean annual total precipitation from days when precipitation > 95th percentile on wet days in the study period

R×1 day
Mean annual maximum 1-day precipitation amount

Correction factors for rainfall erosivity from gridded daily data 180
Because of the lack of long-term precipitation data with high temporal resolution that can be used to calculate EI30, the main method for obtaining nationwide rainfall erosivity datasets is based on daily precipitation data at present (Xie et al., 2016). To improve the accuracy of the rainfall erosivity map, Yue et al. (2020b) collected hourly precipitation data from more than 2000 stations and generated a state-of-the-art rainfall erosivity map for China. In this study, we considered the rainfall erosivity map for the period of 1982-2014 generated from Yue's method as a reference to evaluate the accuracy of the rainfall erosivity 185 values estimated from gridded daily precipitation data. MAE (mean absolute error) and MRE (mean relative error) were calculated for each grid box where a rain gauge was located: https://doi.org/10.5194/hess-2020-633 Preprint. Discussion started: 10 December 2020 c Author(s) 2020. CC BY 4.0 License.
where is the rainfall erosivity extracted from the reference map. The map contained two aspects: the average annual 190 rainfall erosivity (the R-factor) and the 1-in-10-year event EI30. The spatial resolution of the rainfall erosivity map was converted to match those of the four gridded datasets through bilinear interpolation resampling.
To facilitate users in obtaining rainfall erosivity values as close to the accurate value as possible based on gridded daily data in China, we established correction factors applicable to the four gridded datasets using the reference map. Linear regression analysis was conducted based on the characteristics of the errors. It was assumed that the grid boxes that contained 195 meteorological sites that were used in the generation process of the gridded data had higher accuracy in the gridded daily precipitation datasets than grid boxes that did not contain meteorological sites; therefore, these grid boxes were selected for use in developing the correction factors. For the eastern part of China, 647 grid boxes with meteorological sites used in the generation process of the four gridded precipitation products were selected for the correction factors, and the remaining 1348 grid boxes with meteorological sites inside them were used for the evaluation of the correction factors. Similarly, for the 200 western part of China, 167 grid boxes were selected for use developing the correction factors, and 81 grid boxes were used for the evaluation. Then, the linear relationships were established as: where a is the correction factor. The rainfall erosivity value of each grid box estimated from the gridded data was corrected by applying a. Then, the MAE and MRE from Eq. (7) and Eq. (8) were calculated to assess the improved rainfall erosivity values 205 once the corrections had been made.

Comparison of PDFs between gridded and gauge data based on skill score
On the national scale, the skill score values of CN05.1 were the smallest among the four gridded datasets (   Guangzhou and Yumen (Fig. 1), which represented conditions in the northern and southern regions in the eastern part of China as well as the western part of China, respectively. In Beijing and Guangzhou, the four gridded datasets reduced no-rain days (0 mm) but increased drizzle precipitation (0-1 mm) and light precipitation (1-10 mm) compared to those of the original gauge observations. For erosive precipitation, the four gridded datasets demonstrated slight increases in the frequency of occurrence of precipitation in the 10-30 mm range and decreases in the frequency of occurrence of precipitation larger than 40 mm in Beijing and precipitation larger than 50 mm in Guangzhou. For Yumen, erosive precipitation events rarely occurred 230 (PDF < 0.3 %), and no-rain days accounted for more than 90 % of the gauge observations; CPC, GPCC and CN05.1 decreased the frequency of occurrence of no-rain days and erosive days, whereas NMIC obtained an almost perfect estimation of no-rain and light precipitation days and increased the frequency of occurrence of erosive days. For these three places, the differences between CN05.1 and the gauge observations were considerable, especially for the no-rain and low-intensity precipitation events in the northern and western regions, whereas NMIC showed similar PDFs to the gauge observations. Figure 3 also 235 shows a comparison of the extreme daily precipitation amounts (90th, 95th, 99th percentiles and the maximum). As the percentile and extreme degree increased, the reduction in the daily precipitation amounts of CPC, GPCC and CN05.1 increased, whereas these percentiles based on NMIC were close to or even exceeded those of the gauge observations. https://doi.org/10.5194/hess-2020-633 Preprint. Discussion started: 10 December 2020 c Author(s) 2020. CC BY 4.0 License.

Comparison of the precipitation metrics and rainfall erosivity values based on ARF
ARFs indicating the differences in the precipitation metrics obtained through the two approaches ( and ), as well as the rainfall erosivity ( and ), are shown in Fig. 4. The cross-validation results showed that 245 the accuracy of the spatial interpolation in the ( ) generation process was quite high (Table 3), which indicated that the differences between and ( and ) were mainly produced during the upscaling of the gauge precipitation measurements to grids; that is, the discrepancy occurred because of the gridded precipitation products. In the eastern part of China, the gridded data generally conserved the average annual precipitation amounts well, but wet days were overestimated by 10-30 %, which led to lower daily precipitation intensities than those seen in the gauge data. For erosive 250 precipitation (≥ 10 mm day -1 , not shown in Fig. 4), the total erosive precipitation amounts of CPC, GPCC, and CN05.1 were approximately 10-25 % lower than the amounts seen in the gauge data, whereas the erosive days were captured well, leading to a reduction in erosive precipitation intensity. For extreme precipitation metrics, compared with the R95pTOT, which comprehensively reflected the intensity and frequency of daily precipitation of the high percentile, the reduction in the average annual maximum daily precipitation was more obvious. The medians of the ARFs for the R95pTOT of CPC, GPCC, CN05.1 255 and NMIC were 0.85, 0.93, 0.97 and 0.97, respectively, and those for R×1 day were 0.69, 0.78, 0.76 and 0.85, respectively. In terms of rainfall erosivity, the medians of the ARFs for the R-factors of CPC, GPCC, CN05.1 and NMIC were 0.63, 0.77, 0.75 and 0.85, respectively, and for the 1-in-10-year event rainfall erosivity (1-in-10-year EI30), they were 0.49, 0.62, 0.55 and 0.71, respectively. The results revealed that rainfall erosivity was more affected by extreme precipitation intensity, and the 1-in-10year EI30 amplified the differences between gridded and gauge data in extreme precipitation conditions (Table 4). 260 The ARFs in the western part of China generally had similar patterns to those in the eastern part of China but with greater spatial variability (Fig. 4). The medians of the ARFs for the R95pTOT of CPC, GPCC, CN05.1 and NMIC were 0.82, 0.88, 0.93 and 1.08, respectively, and those for R×1 day were 0.78, 0.81, 0.66 and 1.01, respectively. The characteristics of extreme precipitation resulted in a larger discrepancy between and . The medians of the ARFs for the R-factors of CPC, GPCC, CN05.1 and NMIC were 0.57, 0.67, 0.41 and 1.06, respectively, and those for the 1-in-10-year EI30 were 0.54, 0.49, 0.31 and 265 0.70, respectively.

Correction factors for rainfall erosivity values estimated using gridded daily data
Based on Yue's rainfall erosivity map, we conducted regression analysis for the eastern and western part of China separately to establish correction factors for the R-factors and 1-in-10-year event EI30 values using stations for calibration. The linear equations and coefficients of determination, R 2 , are shown in Fig. 5 and Fig. 6. 280 For the R-factor, in the eastern part of China, the results showed that the difference between the R-factor estimated from the gridded daily data and that extracted from Yue's map was larger than the difference between the R-factors estimated from the gridded daily data and that estimated from the gauge daily data; the daily erosivity model (Eq. (2)) generally underestimated the R-factor by 10-20 % in the eastern part of China.  Taken together, the results showed that the correction factors of CN05.1 had the best performance on both the Rfactor and 1-in-10-year event EI30. Maps of rainfall erosivity in China were generated based on CN05.1 by applying the correction factors ( Fig. 5 and Fig. 6). To solve the discontinuity near the boundary area between the eastern and western part of China caused by the difference in correction factors, a buffer was used for each region following the method outlined in Yue et al. (2020a). https://doi.org/10.5194/hess-2020-633 Preprint. Discussion started: 10 December 2020 c Author(s) 2020. CC BY 4.0 License. Fig. 5, but for the 1-in-10-year event EI30.

Discussions
The gridded daily precipitation datasets have different PDFs of daily precipitation than does the gauge dataset. The main difference is that the gridded datasets reduced the frequency of both no-rain days and heavy precipitation compared to those 330 of the original station observations. In the daily rainfall erosivity model, rainfall erosivity increased exponentially with the increase in erosive daily precipitation amount, which meant that a reduction in heavy and extreme precipitation would result in a more serious underestimation of rainfall erosivity from the gridded daily precipitation datasets, especially for the event rainfall erosivity of high intensity and long return periods. The difference mainly comes from the spatial scale between the grid and the point, the station density used for generating the gridded dataset, the interpolation method and the follow-up 335 bias correction method.
Although it is very difficult to define the ground truth for the PDFs of mean precipitation over the target grid boxes without a very dense gauge network or reliable remote sensing data, it can be inferred that the spatial scale discrepancy between the grid and point is objective and is influenced by precipitation characteristics and the spatial resolutions of the grids. The skill score, which qualifies the similarity of PDFs between the grid and point scales, demonstrated larger values in the southern 340 regions in the eastern part of China. In western China and northern regions in eastern China, precipitation tends to be concentrated in summer rainstorms, which are local and short-term, whereas in the southern regions in eastern China, there tends to be continuous precipitation events covering a wide area and lasting a long time, such as the Meiyu, which results in a smaller difference between the grid and point scales. It can be inferred that the differences between the grid and point scales would decrease with decreasing grid size. One extreme case is that if the grid size is small enough to a point, a 345 difference would not exist.
The station density used to generate the gridded dataset, the interpolation method and the follow-up bias correction method influenced the accuracy of the descriptions on the discrepancy between the spatial scales. Taking an extreme example, if the station density is high enough and the areal precipitation in a grid is observed by all gauges in the grid, the true value of the areal precipitation would be obtained, and the true difference of PDFs due to the spatial scale could be quantized. However, 350 this is not practical over a wide region, and the interpolation error was generated and increased with the decrease in station density in the generation process of the gridded dataset. Among the four datasets in this study, CN05.1 and NMIC were generated based on more than 2000 stations over mainland China, which used nearly three times the number of stations compared with the other two datasets, CPC and GPCC. The increase in station density was believed to be of critical importance in the interpolation of precipitation, especially for daily scales with high spatial variability. This finding can be 355 verified from the better performance of the gridded datasets in the eastern region, which has a higher station density (> 4 stations per 10,000 square kilometers) than does the western region (< 1 station per 10,000 square kilometers).
The largest difference between CN05.1 and NMIC was the different interpolation method, and NMIC applied a bias correction process such as quantile mapping, which resulted in a much higher similarity of PDFs compared with other datasets. Since the ground truth for the PDFs of the grid was hard to know, the scale difference of the PDFs between the grid and point scale was difficult to quantify, it was doubtful whether the bias correction method adapted by NMIC was necessary and suitable. In addition, CN05.1 has the highest spatial resolution of 0.25° × 0.25°. To determine whether the spatial resolution affected the accuracy of the rainfall erosivity estimations, two spatial resolutions of 0.5° × 0.5° and 1° × 1° were derived from the CN05.1 dataset based on the bilinear interpolation method and were used to calculate rainfall erosivity. It was found that the linear relationship between the results from the gridded data and gauge data tended to be better with the 365 increasing spatial resolution; thus, the bias of rainfall erosivity estimated from CN05.1 resulting from the difference in PDFs of daily precipitation between the grid and the point scale could be better corrected. In summary, CN05.1 was the mostrecommended dataset among the four gridded datasets for the estimation of rainfall erosivity in China.
Although rainfall erosivity maps based on hourly data are currently available (Yue et al., 2020b), they are not easy to update in a timely manner since the collection of hourly data is harder than that of daily data. Gauge-based gridded daily 370 precipitation datasets, such as CN05.1 in China, are easily available and can be used conveniently for the generation of rainfall erosivity maps (Zhu and Yu, 2015). Spatial scale discrepancies exist not only in gauge-based gridded precipitation data but also in satellite, merged and climate model-simulated precipitation products. This makes the direct use of gridded precipitation datasets dangerous in applications for which the empirical estimation method of a variable (rainfall erosivity, in this study) is developed based on gauge-observed precipitation. If absolute values or absolute changes of the variable are 375 calculated from the gridded dataset, it is suggested that the empirical model is rebuilt at the grid scale (Biasutti and Seager, 2015) or the systematic bias of the results generated from the model developed at the gauge scale are corrected. In this study, the latter method was adopted.

Conclusions
Based on the daily precipitation observations of over 2000 stations in China as well as four widely used gridded daily 380 precipitation datasets, CPC, GPCC, CN05.1 and NMIC, this study compared the PDFs of the daily precipitation amounts, precipitation metrics and rainfall erosivity factors between the gridded daily datasets and the gauge daily dataset and established correction factors for the four gridded daily precipitation datasets using a high-precision rainfall erosivity map for China. The main conclusions are as follows: (1) The PDFs of gridded daily data were different from that of the gauge data, mainly reflecting reductions in no-rain days and 385 heavy precipitation days and increases in light precipitation days. NMIC had the most similar PDF with the gauge data (for 75.2 % of the stations, > 0.9), whereas CN05.1 had the most different PDF (for 2.5 % of the stations, > 0.9, and 42.2 % of the stations > 0.8).
(2) In the eastern part of China, the medians of the ARFs for the R-factors of CPC, GPCC, CN05.1 and NMIC were 0.63, 0.77, 0.75 and 0.85, respectively, and those for the 1-in-10-year EI30 were 0.49, 0.62, 0.55 and 0.71, respectively. In the western part 390 of China, the medians of the ARFs for the R-factors of CPC, GPCC, CN05.1 and NMIC were 0.57, 0.67, 0.41 and 1.06, respectively, and those for the 1-in-10-year EI30 were 0.54, 0.49, 0.31 and 0.70, respectively. The results indicated that rainfall https://doi.org/10.5194/hess-2020-633 Preprint. Discussion started: 10 December 2020 c Author(s) 2020. CC BY 4.0 License. erosivity values estimated using gridded daily data were significantly lower than those estimated using gauge daily data in most cases, which was mainly caused by reductions in extreme precipitation intensity in the gridded data, especially for the 1in-10-year EI30. 395 (3) In the eastern part of China, correction factors for the R-factor of CPC, GPCC, CN05.1 and NMIC were 1.958, 1.570, 1.708, and 1.445, respectively, and those for the 1-in-10-year event EI30 were 2.133, 1.702, 1.959, and 1.477, respectively.
With the correction, for the R-factor, the MREs based on separate validation datasets were reduced by 28.4 %, 20.3 %, 23.2 % and 13.3 %, respectively, and for the 1-in-10-year event EI30, the MREs were reduced by 34.0 %, 20.1 %, 26.0 % and 17.1 %, respectively. There was no obvious improvement after applying the correction factors for the western part of China except the 400 R-factor of NMIC and the 1-in-10-year event EI30 of CN05.1. CN05.1 was the most-recommended dataset, with the lowest MREs of 16.1 % for the R-factor and 22.1 % for the 1-in-10-year event EI30 in the eastern part of China and MREs of 25.1 % for the R-factor and 27.2 % for the 1-in-10-year event EI30 in the western part of China after the application of correction factors.