Rainfall Erosivity in Peru: A New Gridded Dataset Based on GPM-IMERG and Comprehensive Assessment (2000–2020)

Abstract

In soil erosion estimation models, the variables with the greatest impact are rainfall erosivity ($RE$), which is the measurement of precipitation energy and its potential capacity to cause erosion, and erosivity density ($ED$), which relates $RE$ to precipitation. The $RE$ requires high temporal resolution records for its estimation. However, due to the limited observed information and the increasing availability of rainfall estimates based on remote sensing, recent research has shown the usefulness of using observed-corrected satellite data for $RE$ estimation. This study evaluates the performance of a new gridded dataset of $RE$ and $ED$ in Peru (PISCO_reed) by merging data from the IMERG v06 product, through a new calibration approach with hourly records of automatic weather stations, during the period of 2000–2020. By using this method, a correlation of 0.94 was found between PISCO_reed and $RE$ obtained by the observed data. An average annual $RE$ for Peru of 7840 MJ · mm · ha${}^{-1}$· h${}^{-1}$ was estimated with a general increase towards the lowland Amazon regions, and high values were found on the North Pacific Coast area of Peru. The spatial identification of the most at risk areas of erosion was evaluated through a relationship between the $ED$ and rainfall. Both erosivity datasets will allow us to expand our fundamental understanding and quantify soil erosion with greater precision.

1. Introduction

Soil erosion is one of the greatest environmental threats worldwide (Nearing et al. [1], Panagos et al. [2], Karlen et al. [3], Tripathi and Singh [4]), presenting multiple issues such as reduced crop yields, deterioration of water quality due to transport of fertilizers and pesticides, the decreased storage capacity of reservoirs due to sediment production, and losses in soils for cultivation [5,6,7]. There are various categories of erosion, such as water, wind, freezing, and mixed erosion, but the most frequent category with the highest proportion is water [8]. Water erosion of the soil damages the productive surface of the soil due to separation and transport processes, exposing the subsoil of the soil [9,10]. Therefore, the quality of the soil is affected by reducing its water retention capacity and amount of organic matter [11,12], endangering its various ecosystem services, such as CO${}_2$ fixation, agricultural productivity, and flood risk reduction [13], which are expected to increase in demand due to urban expansion and changes in consumption patterns [14]. To prevent the worsening of soil erosion, it is necessary to apply public soil conservation policies based on the monitoring of regions susceptible to $RE$ to understand and mitigate its effects, such as the reduction of agricultural productivity, food and water security, and the national economy [15,16,17,18].

The Intergovernmental Panel on Climate Change reports that there has been a continuous increase in CO${}_2$ emissions in recent decades [19]. On average, the global concentration of CO${}_2$ in the atmosphere increased by 40% from a preindustrial value in 2011. Evidence suggests that such an increase has resulted in an average increase in air temperature of 0.85 ${}^{°}$C (1880–2012), and according to the outputs of global climate models, it is predicted that by the end of the 21st century, the increase will reach 2 ${}^{°}$C with respect to the second half of the 19th century [19]. In this sense, it is expected that changes will occur in the hydrological cycle and, consequently, in the availability of water resources [19,20,21]. The change in the global water supply is expected to be ±10% with more intense storm events, depending on the region of analysis [19,22,23]. In Peru, some investigations on the impact of climate change converged toward a scenario with higher rainfall rates in the summer months (December–March), which would increase the erosive potential of storms, favoring soil loss from the available agricultural area, during these months [24]. Another consequence is the increase in the occurrence of events associated with soil loss, such as landslides [25,26]. For this reason, soil erosion in Peru should be part of the urgent national policy, aiming to identify and monitor areas more vulnerable to the loss of agricultural soils and promote actions to prevent, mitigate, or reverse its effects on desertification and soil degradation processes [27,28].

Soil erosion is caused by two physical processes: (i) separation of soil particles generated by the kinetic energy of the impact of raindrops and the (ii) transport of sediments by surface flow [1]. The level of erosion depends on the regional physiographic, soil, and precipitation characteristics [29], composed of two factors: the intensity of the precipitation and its kinetic energy at the surface of the soil [18]. One of the widely used indicators to quantitatively represent and measure the level of soil erosion, sheet and rill, is the multiannual index of RE [22,30,31,32] and erosivity density ($ED$), calculated as the ratio of $RE$ and precipitation [33,34]. Generally, the $RE$ is calculated in periods of less than 15 min, or adapted by means of statistical algorithms according to the available temporal resolution [35]. To predict soil erosion using $RE$, the empirical Revised Universal Soil Loss Equation (RUSLE) [34,36,37], which combines the influences of duration, magnitude, and intensity of storm events, can be used. Although the RUSLE method is estimated at the annual average level, it can also be calculated on shorter time scales to assess its variability [38]. In its formulation, the most dynamic and reactive factor to changes in climatic conditions is $RE$; therefore, identifying temporal variability provides a more realistic and accurate assessment of soil erosion. For example, the seasonal estimate of $RE$ is used to assess the risk of erosion in various vulnerable regions [11,39,40].

The classic $RE$ equation requires precipitation time series from 1 to 15 min; unfortunately, this information is scarce globally [41,42]. However, through empirical equations, it is possible to use hourly or 30 min data. This convenient technique is commonly used in multiple regions [40,43]. More recently, a variety of research has examined the use of observed data and a satellite precipitation product (SPP) for $RE$ estimation, with their respective limitations due to the source, data derivation model, and spatial scales [2,44,45]. Based on the above, the spatial estimation of $RE$ can be grouped into three approaches: (i) observed-based $RE$: local estimates of weather stations and subsequent geospatial interpolation [44,46,47]; (ii) satellite-based $RE$: the use of satellite-based precipitation products (SPPs) [48,49,50]; and (iii) merged-based $RE$: a mix of both observed and simulated data sources, based on the correction of the $RE$ obtained by the SPPs with respect to data from observed stations, at the national [8,18,51,52,53], regional [11,54], and global scales [40,52]. In this research, the merged-based $RE$ method is used through seasonal satellite correction factors based on automatic weather stations (AWS) on a national scale. This procedure combines the advantages of AWS (accuracy at the hourly timescale) with that of the SPP (spatial variability), widely used as a complement in the analysis of various hydrological processes [55,56].

In South America, studies have been developed with the observed-based $RE$ methodology for the estimation of the $RE$. In Brazil, Sanchez-Moreno et al. [57] used this method because they had more available information, obtaining an $RE$ range of 1672 to 22,452 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ with an increase from east to west; likewise, Mello et al. [58] identified areas in the northwest with very high $RE$ (>20,000 MJ · mm · ha${}^{-1}$ · h${}^{-1}$) and in the northeast with medium $RE$ rates (>2000 MJ · mm · ha${}^{-1}$ · h${}^{-1}$). Using merged-based $RE$, in Ecuador, Delgado et al. [59] estimated the $RE$ based on observed stations, and the Integrated Multi-satellitE Retrievals for GPM (IMERG) obtained a national average of 3173 MJ · mm · ha${}^{-1}$ · h${}^{-1}$. In Chile’s central region, Bonilla and Vidal [47] obtained an $RE$ range of 50 to 6000 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ with an increase from north to south. Moreover, Lobo and Bonilla [60] based on the hourly precipitation from AWS estimates the $RE$ at a point level with a range of 68 to 3520 MJ · mm · ha${}^{-1}$ · h${}^{-1}$. In addition, they highlighted that the use of rainfall at a higher temporal resolution results in a nonlinear decrease in the $RE$.

In Peru, there are investigations that use the three methods. Based on the observed-based $RE$ approach, local studies such as that Romero et al. [61] in the north of the Andean region estimated an $RE$ of 2950 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ at a point level, while Mejía-Marcacuzco et al. [62] on the south coast in Tacna estimated an $RE$ of 1190 MJ · mm · ha${}^{-1}$ · h${}^{-1}$. Using the satellite-based $RE$ method, some global studies determined an average $RE$ in Peru of 2246 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ [63]; on the other hand, through the Global Rainfall Erosivity Database (GloREDa) product developed by Panagos et al. [40], an $RE$ range was estimated between 148 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ in the Pacific Coast region and 14,226 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ in the lowland Amazon. Using the merged-based $RE$ technique, INRENA-Pura [64] prepared a map of soil erosion intensity at a national scale, which was published by the National Institute of Natural Resources (INRENA), using cartographic information, represented by national charts, aerial photographs, and images captured by radar and satellite. Additionally, Sabino Rojas et al. [65] developed a soil erosion atlas on an annual scale from 1981 to 2014, based on information from the PISCOp V1.0 product on a monthly scale. Aybar et al. [66] found a range from 0 to ±10,000 MJ · mm · ha${}^{-1}$ · h${}^{-1}$. In summary, the $RE$ present in Peru, including the aforementioned results, would be in the range of 0 to 15,000 MJ · mm · ha${}^{-1}$ · h${}^{-1}$, with an average of 2000 to 3000 MJ · mm · ha${}^{-1}$ · h${}^{-1}$, with lower values on the Pacific Coast, in contrast with the higher values in the high Andean and the Amazon regions. However, no references were found on the climatology and yearly and monthly evolution of the $RE$, as well as the properties of storm events, at a regional or local level, relevant information for sediment modeling and water erosion studies [67,68].

In this study, the PISCO_reed product was constructed through a seasonal calibration process based on AWS in order to (i) obtain a more accurate $RE$ product on a national scale and (ii) perform a regional assessment of erosivity, which allows us to identify the areas most at risk from the negative effects of soil loss. For this reason, the specific objectives of this research are (a) realize a cross-validation of the $RE$ database and (b) to evaluate spatiotemporal $RE$ by estimating trends and identifying danger zones. Finally, this study has the utility of demonstrating the application of precipitation data based on satellite products and observed stations to estimate the $RE$ at monthly, annual, and multiannual scales.

2. Study Area

Peru is located on the west coast of South America, between 0${}^{°}$02N–17${}^{°}$50.2S and 68${}^{°}$10.2W–81${}^{°}$90.2W, with an extension of 1285 million km${}^2$. This territory is characterized by high topographic variability, with an elevation range from sea level to 6685 meters above sea level (masl), with an average of 1489 masl. Peru exhibits high variability of various climatic factors, such as precipitation and temperature, as a result of the interaction of various influences and forcing features, such as atmospheric flow, the complex orography of the Andes, the cold Humboldt Current System, and El Niño Southern Oscillation [69,70,71,72].

In general, the average annual precipitation varies in the range of ±1 mm on the southern coast, while in the lowland Amazon, it reaches higher values of 4860 mm; the average is 1412 mm. The highest rainfall is in the month of February, and the minimum is during the month of July. In addition, in the Peruvian Andes, the climate is complex and is mainly controlled by the orography that acts as a topographic barrier to the flow of moisture, causing the formation of strong precipitation gradients on the eastern flanks of the Andes [66]. The inter-Andean valleys (≳500 mm) are mainly dominated by convective processes that channel moisture intrusions from the Amazon. At the same time, the influences of the cold and dry air masses coming from the Humboldt Current System cause the driest conditions on the Pacific Coast and on the western flanks of the Andes (≲500 mm). However, during the El Niño Southern Oscillation occurrence, the Humboldt Current System weakens, and the formation of severe convective storms can occur, especially over the North Pacific Coast [66].

For better clarity on the development and evaluation of $RE$ in Peru, the study area was divided into different regions. This segmentation was based on (i) the classification of climatic sectors [73] and (ii) the availability of AWS (Figure 1); the regions were labeled as follows: North Pacific Coast (R1), Central and South Pacific Coast (R2), North Western Andes (R3), Central and South Western Andes (R4), North Eastern Andes (R5), Central and South Eastern Andes (R6), High Amazon (R7), Northern Low Amazon (R8), and Lowland Amazon Center and South (R9). These regions can be grouped into three zones: R1 and R2 (Pacific Coast); R3, R4, R5, and R6 (Andes); and R7, R8, and R9 (Amazon).

3. Materials and Methods

3.1. Overview

Stages, functions, and gridded/observed data used in this research are shown in the flowchart (Figure 2), obtaining as results the PISCO_reed product and its evaluation. First, the $RE$ from the AWS and from IMERG was estimated on an hourly scale using the conventional RUSLE method. Then, the $RE$ from satellite series was constructed through a correction based on the AWS, with a validation of the interpolation of the multiplicative factor by seasonal periods. Next, the uncertainty of the erosivity products generated was evaluated, taking the observed data as a reference. Second, the $ED$ rates and the areas of greatest danger of erosivity were estimated. Moreover, trends of the annual $RE$ series by region were calculated.

3.2. Data

The estimation and evaluation of the $RE$ in Peru required gridded products and observed precipitation data: (i) SPP IMERG, (ii) AWS from SENAMHI, and (iii) Global $RE$ gridded products. These datasets were used in this study to develop the methodology.

3.2.1. Satellite Rainfall

The increased accessibility to high spatiotemporal resolution SPPs, such as the Global Precipitation Measurement (GPM) mission, has facilitated the development of hydrological research. GPM, jointly developed by the National Aeronautics and Space Administration (NASA) and the Japan Aerospace Exploration Agency (JAXA), was put into orbit and has been in operation since 2014 as a successor to TRMM (deployed in 1997), to unite and increase precipitation estimates from a constellation of satellites [74]. NASA’s Goddard Earth Sciences Data and Information Services Center (GES DISC) provides access to GPM-derived products through multiple applications that contribute to knowledge about the water cycle and its impact from extreme natural events [45,75]. NASA’s Earth Observing System Data and Information System has defined three levels of products for the distribution and labeling of GPM products; level 3 is called IMERG. For precipitation estimation, IMERG is based on an algorithm that uses GPM microwave observations as a reference to combine, using various interpolation methods, such as CMORPH-KF [76], with multiple satellite data, such as passive microwave (PMW), geosynchronous infrared (IR), and TRMM microwave imager estimates (TMI), among others. Therefore, it presents advantages over TRMM in terms of temporal resolution (3 to 0.5 h), spatial resolution (0.25${}^{°}$ to 0.1${}^{°}$), and area covered (60${}^{°}$S–60${}^{°}$N). Furthermore, IMERG version 6 features two enhancements. The first is the length of the period covered by GPM, using TRMM estimates (2000–2014) to obtain hourly rainfall from June 2000 to the present [74]. The second is the distinction between liquid and non-liquid precipitation, to calculate the precipitation from this first ratio. IMERG classifies its products according to latency periods: Early, Late, and Final with ±4 h, ±14 h, and ±3.5 months, respectively [77]. The Early and Late versions only have a calibration for observed climatology, while the Final version uses the observed monthly rainfall from the Global Precipitation Climatology Center (GPCC) for bias correction. The GPM-IMERG-V06-Final product (IMERGF) contains the variable PrecipitationCal, calibrated precipitation with observed stations, expressed in $mm/h$. To estimate $RE$, this research uses IMERGF data from September 2000 to August 2020. The product is available in netCDF4 format, freely downloadable from the portal https://disc.gsfc.nasagov/datasets (accessed on 1 April 2023).

3.2.2. Observed Rainfall

The observed hourly precipitation used in this study was obtained from the National Hydrology and Meteorology Service of Peru (SENAMHI). The data were extracted from hourly rainfall at 322 AWS and 9 AWS with 10 min rainfall $AWS\_10min$ located at altitudes ranging from 100 to 5000 masl, with an average of 2000 masl (Figure 1). The amount of AWS has increased in recent years, from 120 in 2014 to 320 in 2021. In addition, the temporal period of the seasons is from 1.5 to 6 years, with an average of 3.5 (Figure A1); therefore, this is the most extensive hourly dataset for the estimation of the $RE$ in Peru. The analysis of the results derived from the AWS is classified in the nine climatic regions, obtaining a total of 180, 90, and 40 AWS in the Coast, Andes, and Amazon zones, respectively. SENAMHI provides this information with basic quality control (QC), with a procedure that can be divided into two stages: (i) filtering by extreme physical, national limits (401 mm/h) and internal consistency control and (ii) visual inspection of the hourly, daily, monthly, and annual series, to validate internal consistency and homogeneity. The reliability of the AWS data has been successfully tested in various investigations that evaluate the performance and accuracy of the SPPs in Peru [56].

3.2.3. Global $RE$ Products

This study uses global $RE$ products such as GloREDa and CMORPH for comparison purposes. GloREDa is based on observed data of high temporal resolution (1 to 60 min) that mainly cover the period from 2000 to 2020, collected from various countries of different regions and climates; however, in South America, the number of stations only represents 4% with an irregular geographic distribution [40,78]. The resulting global $RE$ product has a spatial resolution of 0.01${}^{°}$ and can be accessed through the European Soil Data Centre website (https://esdac.jrc.ec.europa.eu/, accessed on 30 August 2023).

On the other hand, CMORPH is a global precipitation satellite product, based on geostationary, low-orbit, and PMW satellites [18], with temporal resolution (30 min), spatial resolution (8 km × 8 km), temporal period (1998 to present), and coverage area (60${}^{°}$S–60${}^{°}$N), developed by the National Oceanic and Atmospheric Administration (NOAA) [79,80]. This product has been applied for the identification of extreme hydrological processes, such as $RE$; Bezak et al. [52] obtained this variable at pixel level using the RUSLE method [81]. The RE from CMORPH was downloaded from https://esdac.jrc.ec.europa.eu/ (accessed on 30 January 2023).

3.3. Methodology

3.3.1. Gridded Product Construction

Estimation $RE$

The estimation of the $RE$ requires a minimum period of 20 years of information in order to reduce the uncertainties and biases generated by dry and wet years [30]. Moreover, this indicator provide valuable information for storm event assessment separately; therefore, Wischmeier and Smith [30] recommend as storm event identification requirements (Figure 3) (i) the use of a time interval minimum (TMI) of 6 h between each event, (ii) hourly rainfall greater than 0.2 mm for hourly level and 0.1 mm for 30 min level, (iii) finally, the accumulated volume by each event being greater than 0.2 mm.

The precipitation intensity (I) for the hourly and sub-hourly scales was obtained through the relationship of the accumulated precipitation and its corresponding recording time (1):

$$I=\frac{P}{T}$$

(1)

where I indicates the intensity of precipitation in $mm/h$, P is the precipitation in $mm$, and T is the recording time in hours. In each storm event, according to Wischmeier and Smith [30], the unit measure of kinetic energy ($e_r$) is estimated at each chosen time resolution interval, Equation (2):

$$e_r=0.29[1-0.72\times exp(-0.05i_r)]$$

(2)

where $i_r$ is the precipitation intensity during the time interval in $mm/h$. The sum of the unit kinetic energy, multiplied by the rainfall volume for each time interval in a storm event, results in its total kinetic energy E, expressed in the next Equation (3):

$$E=\sum _{r=1}^m{e}_r{V}_r$$

(3)

where E is expressed in MJ · ha${}^{-1}$ and $V_r$ is the rainfall in mm, during an event. Subsequently, Brown and Foster [36] define $RE$ ($E{I}_{30}$) as the result of the multiplication of E with the maximum intensity in 30 min of each storm event, as indicated in Equation (4):

$$RE=E{I}_{30}=E{I}_{30max}$$

(4)

where $E{I}_{30}$ is expressed in MJ · mm · ha${}^{-1}$ · h${}^{-1}$ and $I_{30}$ is the maximum rainfall intensity by storm event in $mm/h$; then $RE$ is $E{I}_{30}$ grouped by year and is expressed as MJ · mm · ha${}^{-1}$ · h${}^{-1}$ · yr${}^{-1}$. In the case of only having temporal resolutions of 60 min, Panagos et al. [53] and Yin et al. [11] suggest multiplying E×$I_{60max}$ with a correction coefficient ($C{C}_{60}$); this value varies from 1.15 to 3.37 (Equation (5)):

$$RE=E{I}_{30}=E\times I_{60max}\times C{C}_{60}$$

(5)

where $I_{30max}$ corresponds to the maximum intensity of 30 min identified in each storm event. $E{I}_{30}$ is the equivalent to the $RE$ of RUSLE, and $C{C}_{60}$ is the correction coefficient.

As part of the methodology, $C{C}_{60}$ was estimated by means of correlations between the $RE$ obtained from $AWS\_10min$ added to 30 min ($AWS\_30min$) and those obtained from a temporal resolution of 60 min ($AWS$); the distribution of both AWS is represented in Figure 1. Fischer et al. [82] identified an underestimation of the $RE$ of RUSLE when using rainfall series with a temporal resolution greater than 30 min.

Estimation $ED$

The $RE$ needs to be evaluated with wide availability of hourly rainfall; however, this information is scarce in Peru. To solve these deficiencies, Foster et al. [83] introduced the $ED$ function. This erosivity index is more stable and independent of rainfall data availability than $RE$, since it depends on the number of erosive events [83]; moreover, it is used to evaluate erosion patterns [11]. Although the $RE$ provides information on the erosive potential of rainfall, it does not provide information on the concentration of extreme storms during the year. On the other hand, the $ED$ better represents the $RE$ patterns and the type of precipitation during erosive events [84]. Very high $ED$ values indicate high runoff, implying an area more prone to floods and intense storms [37,45,85]. The $ED$ is the ratio of $RE$ to precipitation [83], expressed through Equation (6):

$$ED=\frac{RE}{P}$$

(6)

where the annual accumulated precipitation is measured in mm·yr${}^{-1}$ and the annual $RE$ in MJ · mm · ha${}^{-1}$ · h${}^{-1}$ · yr${}^{-1}$. The equation was applied on a monthly and annual scale. The estimation RE-IMERGF were evaluated through a pixel by pixel, while the RE-AWS was calculated by each station location; both results were classified by region.

Construction and Validation of $RE$

Various studies show a high correlation between the extreme rainfall (>99th) obtained from SPPs with rain gauge adjustment and observed stations in South America [66,86]; accordingly, their spatial distribution shows a good correlation with the $RE$ estimate [38,51,87]. Therefore, there is evidence of using hourly rainfall from IMERGF to estimate RE-IMERGF, with respective correction-based RE-AWS. Correction of simulated data based on observed stations is widely used to improve the accuracy in the generation of multiple rainfalls datasets and their derived products [79,88,89].

Several studies identified a high correlation in the spatiotemporal variability between precipitation and $RE$ [23,51,63]; therefore, the correction of RE-IMERGF was performed by rescaling precipitation with observed data, for each independent pixel, using an annual factor average by season. By extending the approach of Chen et al. [51], the correction process was as follows: (i) obtainment of the calibration factor by grouping the monthly series at the seasonal level (summer, autumn, winter, and spring) in the same period, 2015–2020; (ii) employ a linear regression between the point-gridded values from RE-AWS and RE-IMERGF in the four seasonal periods to obtain the slope of each linear model defined as the seasonal multiplicative factor (FME); (iii) spatial interpolation of the FME point values by using the inverse distance weighted interpolation (IDW) method, at the same native spatial resolution of IMERGF (0.1${}^{°}$); (iv) spatial aggregation applied to reduce the spatial resolution from 0.1${}^{°}$ to 0.25${}^{°}$, with the aim of avoiding spatial inconsistencies as a consequence of the high variability of the multiply factor; (v) and finally, PISCO_reed was obtained as a result of multiplying RE-IMERGF by the FME maps. The validation was realised at the pixel level with the observed data from RE-AWS during the common period, 2015–2020.

Metrics Validation

The performance of PISCO_reed was evaluated by regions with reference to RE-AWS, and both databases were compared by statistical correlation and magnitude difference. The statistical metrics were Pearson correlation coefficient (r) and the aggregation index (dr) [90]. The coefficient r is a measure of the relationship of strength between two variables where a result of 1 indicates a perfect relationship with a positive slope, while −1 indicates a relationship with a negative slope. The dr is similar to a correlation coefficient, except that it varies between −1 and 1, a high value (>0.5) indicating both high correlation and low absolute differences between the observed and simulated time series. In addition, due to the difference between the magnitudes of the hydrological variables in the study regions, it was necessary to have statistical metrics that indicate the relative difference between the observed and simulated data, through a normalization process. In the comparison of the samples, the relative differences are normalized by the observed sample. The normalized mean gross error (NMGE) and normalized mean bias (NMB) statistics were selected. NMGE is a measure of the mean relative deviation from the observed values and is independent of the magnitude of the hydrological variable, suitable for comparison between arid and humid regions. Meanwhile, NMB is useful for evaluating the $RE$ of different monthly rates, since the mean bias is normalized by dividing it by the observed $RE$.

3.3.2. Erosivity Evaluation

Trends

Analysis of $RE$ and $ED$ trends provides information on rates and magnitudes of change over long periods of time [30]. Thus, the trends detection were evaluated using the non-parametric Mann–Kendall (M–K) test; since it is widely used to identify monotonic trends in hydro-meteorological time series. It is more resistant to the existence of outliers and does not require that the data be normally distributed [91,92,93,94].

In addition, the magnitude of the trend was estimated using the non-parametric method of Sen’s slope ($SS$), where a positive value indicates an increase, while a negative value indicates a decrease in the trend [95,96]. In this study, the areas with statistically significant positive or negative $RE$ trends were identified at the pixel level at a significance level (p-value) of 0.1, as well as the respective magnitude of the trend at the seasonal level, expressed in 10-year (decade) changes (T), as shown in Equation (7).

$$T=SS\times 10$$

(7)

Global and National Comparative Analysis

The performance of IMERGF for $RE$ estimation and storm events properties, such as the intensity of precipitation, duration, accumulated precipitation, and number of events during 2015–2020, were evaluated using a statistical approach based on the results obtained by the AWS. In addition, a comparison was made in statistical terms by pixels between the observed data, PISCO_reed, and other global products that use the RUSLE method to estimate $RE$, such as the GloREDa product (2000–2020) developed by Panagos et al. [78] and the CMORPH $RE$ product (1998–2020) obtained by Bezak et al. [52]. The temporal availability of global products is lightly different; therefore, it was contrasted by the multiyearly average of $RE$; in the case of the correlation, a symmetric line with the origin at 0 was used as a reference.

Risk Map

In this study, we used $ED$ to evaluate erosivity patterns and their effect under different ranges of rainfall, through a risk map [2]. Based on the combination of quartiles of the variables at a multiyear scale of $ED$ and rainfall, we obtained 16 classes to characterize the susceptibility to soil erosivity, according to Das et al. [45]. In this range, areas with very high $ED$ and very high to very low mean rainfall are the most vulnerable, while areas with very low $ED$, regardless of the precipitation rate, can be considered as those less prone to soil erosion, according to Panagos et al. [2].

4. Results

4.1. Spatiotemporal Distribution of $RE$ and $ED$

In Figure 4, the properties of the storm events are compared, such as the average accumulated rainfall, maximum intensity, number storm, duration, annual rainfall, $RE$, and $ED$ during the 2015–2020 period obtained from IMERGF and AWS. The RE-AWS correction was based on Equation (5) represented in Figure A2, where $E{I}_{30}$ is equal to RE-AWS based on the twenty one annual comparisons of $RE$ identified in 9 $AWS\_30min$; linear regression forced through the origin resulted in a slope ($C{C}_{60}$) of 1.541 with a standard error of 0.01. Rainfall during the events and annual accumulated rainfall were slightly underestimated in parts of the Coast and Amazon zones. We noticed that in the Andes zone, there was a lower overestimation of the lowest daily rainfall and a high underestimation of the highest daily rainfall, a similar behavior found by Derin et al. [97]. On the contrary, the duration of the events was overestimated at the national level by IMERGF, at a ratio of up to 3 to 1 in the Andean regions. However, the greatest disagreement of magnitudes was found in the identification of the maximum intensities, with a high correlation by underestimation up to 5 times in all the regions, which has a more relevant impact on the estimation of the $RE$ and $ED$. The maximum underestimations were found in the Amazon zone for $RE$ and $ED$, reaching 6.2 and 8.3 times on average, respectively. In summary, this comparison showed good performance in correlation for storm properties since an SPP with higher temporal resolution is used, with a correlation greater than 0.74, except for the low value (0.47) found in the average duration of storms.

The precision of the annual $RE$ estimates using IMERGF was also evaluated and compared with the AWS. The correlation coefficient of the mean annual RE-AWS and RE-IMERGF was 0.84. In general, taking AWS as a reference, the values based on IMERGF indicate an underestimation of the $RE$, present in all the analysis regions, with IMERGF values in R1 and maximum values in R7 and R9.

Based on the comparison of RE-IMERGF and RE-AWS at a seasonal scale for the period of 2015–2020, the multiplicative factors for each AWS were calculated (Figure 5). Values of less than 1 indicate an underestimation of IMERGF, while values greater than 1 indicate an overestimation. Figure 5 shows the variability between 0 and 28 in the histograms with a median of 5.5. There is a smaller amount of multiplicative factors < 1 (∼4% on average) that indicate the overestimation of the $RE$, mainly in the coastal regions during the winter and spring seasons. Meanwhile, the regions with the greatest underestimation of the $RE$ are R3, R4, R5, R6, and R9, precisely where the most intense storm events are recorded. The interpolation of the multiplicative factors was added to reduce the spatial resolution to 0.4${}^{°}$, in order to reduce the factors with the greatest difference and proximity (Figure 5). Subsequently, IMERGF was corrected, through seasonal multiplication with the multiplicative factor map.

PISCO_reed from 2015–2020 was compared with the AWS data in Figure 6. The results show an improvement for dr with means in the range of 0.5 in R9 to 0.65 in R1. Additionally, the underestimation of the product is improved with means of NMB from −0.06 to −0.43. On the other hand, the NMGE medians range from 0.78 to 1.11, indicating an adequate performance of the $RE$ from PISCO_reed, which will be analyzed by regional scale in the following sections. In addition, together with the $ED$ product, both on a monthly and yearly scale, PISCO_reed was created.

4.1.1. Comparison with Global Products

When comparing the PISCO_reed product with the global products (GloREDa and CMORPH obtained from the multiyear average), taking as reference the observed data from the AWS (Figure 7 and Figure A1), a higher correlation was found in the PISCO_reed product (r = 0.94), finding slight over-estimations in some regions. On the other hand, the CMORPH product had a greater underestimation in all regions, mainly in R3, R4, R5, and R6. This underestimation was also observed in GloREDa, with the particularity of a marked overestimation of the $RE$ in R1 and R2, where the $RE$ was less than 3000 MJ · mm · ha${}^{-1}$ · h${}^{-1}$.

4.1.2. National Analysis

Figure 8 shows the RE-AWS and the map of $RE$ and $ED$ from PISCO_reed. The $RE$ and $ED$ map is presented with a spatial resolution of 0.1${}^{°}$ or ∼10 km. The mean $RE$ value is 7118 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ with a high variability that can be expressed by the standard deviation of 6231 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ or a coefficient of variation of 0.88. The median $RE$ is 7161 MJ · mm · ha${}^{-1}$ · h${}^{-1}$, the first quartile is 841 MJ · mm · ha${}^{-1}$ · h${}^{-1}$, and the third is 12,377 MJ · mm · ha${}^{-1}$ · h${}^{-1}$. In the same way, the $ED$ has an average value of 3.18 MJ · ha${}^{-1}$ · h${}^{-1}$ with a coefficient of variation of 0.62; in addition, the 25th, 50th, and 75th percentiles are 1.33, 3.48, and 4.49 MJ · ha${}^{-1}$ · h${}^{-1}$, respectively.

The $RE$ climatology map in Figure 9 shows that the month with the highest $RE$ is March (954 MJ · mm · ha${}^{-1}$ · h${}^{-1}$), followed by December, January, and February (793 to 846 MJ · mm · ha${}^{-1}$ · h${}^{-1}$). The $RE$ shows its lowest values for the period from June to September (224 to 397 MJ · mm · ha${}^{-1}$ · h${}^{-1}$), coinciding with the variability of precipitation climatologies, which have a high correlation in the range of 0.719 (March) to 0.923 (July). Likewise, the spatial distribution of the $RE$ shows a notable difference between R5 and R9, with respect to the formerly R1 to R4, with high rates of $RE$ in the first years.

4.1.3. Regional Analysis

Figure A3 shows, at the regional level, the 2000–2020 annual variability of storm properties, such as total precipitation from erosive events, average storm duration, number of storm events, $RE$, and $ED$. In 2001, 2012, and 2017, when the $RE$ reached its highest annual values, the total precipitation and the number and duration of storm events were also higher than in the other years. The mean annual duration of storm events is in the range of 3 and 10 h, the mean annual number of storms ranges from 40 to 280, the total precipitation ranges from 50 to 2500 mm, and the mean annual $RE$ ranges from 40 to 11,000 MJ · mm · ha${}^{-1}$ · h${}^{-1}$.

The properties of the average annual storm events are analyzed at the regional scale in Figure 10. The average number of storm events varies gradually from 47 in R2 to 242 in R8; that gradual scale in the regions is also shown in the average duration of the storms in the range of 2.5 to 8 h. The Amazon zone has a higher $RE$ (7677 to 13,648 MJ · mm · ha${}^{-1}$ · h${}^{-1}$) compared with the coastal and Andean zones (44 to 1488 MJ · mm · ha${}^{-1}$ · h${}^{-1}$). This contrast is mainly due to the total precipitation, which on average is up to 13 times higher between the regions of the respective groups. In the same way, the analysis of the $ED$ shows average rates of 0.5 to 2 MJ · ha${}^{-1}$ · h${}^{-1}$ in the coastal and Andean zones, while in the Amazon, the average $ED$ varies from 3.8 MJ · ha${}^{-1}$ · h${}^{-1}$ in R7 to 5.6 MJ · ha${}^{-1}$ · h${}^{-1}$ in R9.

Figure 11 shows the monthly $RE$ by region. The seasonality of the $RE$ is visible in the regions, presenting the highest values of the $RE$ during the wet season from December to April (average 510 MJ · mm · ha${}^{-1}$ · h${}^{-1}$), with maximum values in the month of March. Meanwhile, the low $RE$ values are found during the dry season from June to September (average 152 MJ · mm · ha${}^{-1}$ · h${}^{-1}$), with minimum values during the months of July or August. Finally, a descriptive statistical summary is shown in Figure 11 and Figure A4.

4.2. Spatial Variation of Risk Areas

The spatial variation of the risk areas, classified into sixteen classes in Peru, is shown in Figure 12. High-risk areas for $RE$ were mainly found in R8 and R9; only some areas were identified in R7 due to their high rates of precipitation and $ED$. These transition areas to the Amazon plain are associated with areas of high soil erosion and landslides, depending on the physiography. The rest of the regions (R3 to R7) belong to the medium precipitation classification but with high rates of $ED$ associated with a medium risk of erosion. These regions are the ones that are frequently affected by soil erosion and landslides due to their physiography with steep slopes and intense rainfall. In the Andean zone, during the period of 2004–2013, there were 38 landslides registered with fatal consequences for the local population and high damage to the supply network at the national level [98]. In addition, the Andes mountains have multiple geological faults and anthropogenic soil erosion activities, and are affected by extreme hydrological events, such as El Niño.

4.3. RE Seasonal Trend

The spatial variation of seasonal trends (decadal) of $RE$, evaluated with a confidence level of 95% (p < 0.05), including the regional trend mean, is shown in Figure 13. The spring season shows only a punctual increase in the south of R8 with positive changes greater than 500 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ by decade, while in the north of R8, the change is negative at rates of −800 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ by decade. In the rest of the country, such as R1 to R4, there are no significant changes (∼150 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ by decade) during all seasons. This same pattern is observed at R5 and R6 in the winter season; however, in the autumn and summer seasons, a slight increase is identified in these regions with trends of up to −250 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ by decade. Finally, it is highlighted that R8 and R9 are the ones with the widest range of trends (−1700 to −900 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ by decade), except during the winter season, where the variation is between −600 and −500 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ by decade.

5. Discussion

5.1. Comparison with Other Studies and Analysis of Causes

Observed subdaily station data are a fundamental source of information on rainfall frequency and amount [8]. However, the scarce availability, irregular distribution, and spatial inconsistencies, in steep regions, limit the use of the observed rainfall [99]. Therefore, SPPs are an alternative solution to precipitation estimation, since they provide spatially continuous information over large geographic areas [100]. The complementary use of this information is already used in South America for the estimation of the $RE$ with the RUSLE method on a monthly and annual scale [47,59]. The evaluation of the product obtained from $RE$ in this study is compared with the hourly rainfall data, suitable for evaluating hydrological products based on hourly information, such as RE [45], which constitute the longest, most numerous hourly base and with quality control in Peru.

The main factors of uncertainty in the generation of the $RE$ are the precision of precipitation data, the estimation model of the $RE$, and the correction factor in the conversion of time scales [8]. Furthermore, Catari et al. [101], identified the sources of the error in the estimation of the $RE$: erroneous measurements of precipitation, the efficiency of the kinetic energy equation of storms as a function of intensity, and the variation of spatial patterns.

This study is consistent with the $RE$ obtained by Barurén [102]; however, the maximum average $RE$ is higher than our product (41,106 MJ · mm · ha${}^{-1}$ · h${}^{-1}$). This may be due to the fact that the analysis periods are different, altering the final result of the multiyear product. In addition, Rosas and Gutierrez [103] agree with the quantification of erosion variability, finding an increase in erosion in the direction of R4 to R3, the latter region being the region with the highest erosion in the western Andes of Peru. In R2, the results are in the range obtained by Mejía-Marcacuzco et al. [62], with an analysis period within the years of this study of 1997–2020. In the South American region, one of the global $RE$ products is GloREDa, which estimates this variable using globally observed data and the ERA5 precipitation reanalysis product; however, the largest differences between $RE$ and GloREDa estimates were observed in South America due to underestimation of precipitation in mountainous regions, such as the Andes [52]. The observed stations used for the elaboration of this product are mainly from the Brazilian Amazon.

5.2. Limitations

The IMERGF presents strengths and weaknesses in the analysis of storm events in multiple regions around the world due to difference in rainfall intensities, stations used for correction, type of sensors, resolution temporal, etc. According to Derin et al. [104] and Manz et al. [97], IMERGF is better than other products, such as TMPA-V7, in terms of accuracy in estimating the frequency of occurrence and intensity of precipitation distribution in the Andean regions. Regarding the detection of daily rainfall volumes, Das et al. [45] reported that IMERGF is sensitive in identifying light and no rain conditions. Regarding the detection of extreme precipitation, the global average is satisfactory due to the correction of biases based on the observed stations [105]; therefore, for the geographical conditions of the Andean zone of Peru, Derin et al. [97] identified slight overestimations. These biases influence the underestimation of the $RE$ in R1, R6, and R4, since IMERGF has limitations in detecting extreme storms. Therefore, seasonality from IMERGF precision undoubtedly resulted in seasonal variation in $RE$ and $ED$ biases.

IMERGF, being a high spatial resolution dataset, reduces the uncertainties in the spatial patterns caused by the interpolation of rainfall and its derivatives in geographical areas with low AWS density [45]. However, the scarcity of AWS in R8 limits the bias correction of the $RE$, although this region is more resilient to erosion due to its flat land cover and physiography. This region presents a lower level of danger in the generation of landslide events according to the map defined by Millán-Arancibia and Lavado-Casimiro [26].

The $RE$ estimate is sensitive to the temporal resolution of the rainfall [45], which originally required 30 min for its calculation. Therefore, if we use a different temporal resolution, it is necessary to apply an appropriate regression function to obtain a result equivalent to the estimate of the $RE$. Various studies develop a correlation function based on 30 min resolutions in areas where subhourly rainfall is not available. Panagos et al. [2] found a high underestimation of the $RE$ (∼56%) when using hourly rainfall without correction, but when using a correction factor, a good correlation was obtained. The use of an erroneous function, or a low and irregular distribution of rainfall stations for its development, can increase the bias of the estimate [45]. In this investigation, information from rainfall stations with 10 min temporal resolution was used for the first time to find an adequate correction factor for the hourly $RE$. We found a factor greater than 1.5, which indicates an underestimation of $RE$ compared with the original calculation with 30 min records.

5.3. Applications

RE is a key factor in estimating soil erosion [23]. Rain and its kinetic energy are the main driver of soil water erosion processes. These processes are associated with the detachment of soil particles, the generation of runoff, and the triggering of mass movements [106,107]. On the other hand, agricultural practices without soil conservation measures together with the increase in the intensity of severe storms [23,108] would cause erosion and increase the amount of solids in suspension in rivers. This could clog reservoirs, raise riverbeds, and affect water quality, due to high turbidity levels in the rainy season.

High erosion rates have effects on the ground as well as effects external to it [11]. As a consequence, they are responsible for the contamination and low productivity of farmlands [109]. It is expected that the $RE$ results developed in this research can be used to update the maps and identify regions vulnerable to this mode of erosion. We must understand the physical processes of soil erosion; for this reason, it is important to implement a monitoring system for erosion and sediment production in experimental basins. The results of this study can be considered in the planning of public policies to reduce erosion aimed at conserving soil productivity [110] and maintaining soil ecosystem services at a tolerable level [111], especially in regions where a large increase was detected.

In Peru, data from AWS are spatially scarce. Therefore, SENAMHI expects the installation of rain radars, which will make it possible to have high-resolution spatiotemporal information. This could potentially improve [112] estimates of $RE$ [113].

6. Conclusions

In this research we generated PISCOp_reed V1.0, a new dataset in Peru for estimated and analyze by region the $RE$ using the RUSLE methodology, based on a correction of IMERGF from the hourly and subhourly AWS for the period of 2000–2020. The following is concluded: through the spatial correction of RE-IMERGF based on the RE-AWS, it was possible to reduce the biases to analyze its spatial distribution at the national and regional levels, on various time scales (climatology, monthly, and annual). At national level, the $RE$ mean was 7840 MJ · mm · ha${}^{-1}$ · h${}^{-1}$, in the range of 0 to 60,000 MJ · mm · ha${}^{-1}$ · h${}^{-1}$ (R9), with a spatial distribution similar to rainfall. The results of this study indicate that the previous analyses underestimated the $RE$ due to the underestimation of the maximum intensities by the use of daily rainfall data; however, the $RE$ from PISCO_reed was in the range of regional studies in the Amazon and Pacific Coast zones, obtained with similar methodologies.

The PISCO_reed product has the advantage of quick and simple access to information for the characterization and identification of risk areas to erosion and trends by pixel level. A high risk of $RE$ is indicated at the Amazon zone based on a combination of high rainfall and high $ED$, where the number of storms, $RE$, and $ED$ are highest of all regions. In addition, the spatial analysis of risk areas suggests that it is insufficient to predict the $RE$ only based on rainfall, especially in R1 and R4. This information provides a method and database for the implementation of soil conservation, management policies, water administration, disaster prevention, agricultural/forestry planning, and other applications for the management of hydrographic basins. Especially in the Andes and Amazon zones, there are significant changes in trends of $RE$ in the last 20 years. At the national scale, we estimated that about 1.1% of all the pixels have statistically significant positive trends and 13.5% have statistically significant negative trends at the 95% confidence level.

In the coming years, through the use of radars for the identification of observed hourly precipitation, storm events will be able to be analyzed with greater precision, improving the accuracy of the PISCO_reed product.