Analysis of Rainfall Erosivity Trends 1980–2018 in a Complex Terrain Region (Abruzzo, Central Italy) from Rain Gauges and Gridded Datasets

: The erosive capacity of precipitation depends on its intensity, volume, and duration. The rainfall erosivity factor (R) of the Universal Soil Loss Equation (USLE) requires high frequency (subhourly) data. When these are not available, R can be estimated from simpliﬁed indices such as the Modiﬁed Fournier Index (MFI), the Precipitation Concentration Index (PCI), and the Seasonality Index (SI), which are computed from monthly precipitation. We calculated these indices for 34 stations in the complex terrain Abruzzo region (central Italy) during 1980–2018, based on both gauge (point) and grid datasets. Using 30-min rainfall data of 14 stations, we veriﬁed that MFI and PCI are reliable predictors of R (R 2 = 0.91, RMSE = 163.6 MJ mm ha − 1 h − 1 year − 1 ). For MFI, grid data do not capture the peaks in high-altitude stations and the low values in some inland areas, detected by the point dataset. Grid data show signiﬁcant MFI positive trends in 74% of the stations, while the point data display signiﬁcant positive trends in only 26% of stations and signiﬁcant negative trends in four stations in the inland areas. The grid data complex orography requires preliminary validation work.


Introduction
Soil erosion is a process of considerable complexity that depends on several factors such as climate, soil type, topography (slope and length), vegetation cover, and cultivation systems [1]. As regards the climate, the extent of the erosive phenomenon is attributable to the erosive capacity of precipitation (aggressiveness or erosivity), which in turn depends on characteristics such as intensity, volume, and duration [2,3]. Here, we estimate the trends of rainfall erosivity detected in a complex terrain region in central Italy (Abruzzo) using monthly data at 34 locations, comparing the calculation using homogenized rain-gauge timeseries (referenced as "point" in the manuscript) and timeseries interpolated from a widely used European gridded climatic dataset (E-OBS, referenced as "grid") at the raingauge locations. The main aims are detecting the presence of possible trends and verifying the suitability of the gridded dataset for a complex orography area.
One of the most commonly applied erosivity indicators is the R factor of the Universal Soil Loss Equation (USLE) [2], or its revised version, RUSLE [4]. The erosivity of a single rainstorm (R e ) is obtained by multiplying its kinetic energy by I 30 (mm/h), i.e., the maximum rainfall intensity recorded during the rainstorm in a period of 30 min. The R factor's final value is then derived by averaging the total annual erosivity computed for an adequate number of years. High-quality and high-frequency rainfall datasets are therefore simplified models to evaluate the 1961-2018 erosivity trends on a European scale based on the grid dataset Uncertainties in Ensembles of Regional Reanalyses (UERRA). This study shows an overall increase in erosivity for most Italian regions, thus providing results not much in agreement with the analyses mentioned above based on point datasets. On the other hand, in a recent discussion, Wang et al. [36] highlighted the potential inaccuracies and limits of the grid data for this type's applications.
This paper presents and discusses an analysis of the trends detected in the rainfall erosivity indices in the Abruzzo region (central Italy) over the last 40 years (1980-2018). A flowchart of the work is given in Figure S1. The objective of the work is twofold. First, to obtain a spatial assessment of the trends in the MFI, PCI, and SI series in relation to the orographic and climatic zones of the region; the second and more innovative objective is to compare the results obtained based on the point and grid datasets. This will allow us to understand the extent to which it is possible to rely on grid data to obtain fairly reliable space-time evaluations of the indices mentioned above in regions characterized by a marked spatial variability of climatic characteristics.

Rainfall Data
The study considered the daily rainfall data of 34 stations over the period 1980-2018. The stations are quite uniformly distributed over the territory (the red points in Figure 1a). The study was conducted using data from two different precipitation data sources: a high-quality homogenized dataset of daily observed data [37] and a widely used gridded dataset, namely the E-OBS product at 0.1 • degrees resolution [38], which represents the grid product with the best resolution in the region considered. In this work, we refer to the two data sources as the "point" and "grid" datasets, respectively. Although of modest extension (10,763 km 2 ), the Abruzzo region has a significant climate spatial variability, mainly due to a complex orography and adjacency to the sea [39]. Grid daily data were downloaded from the public E-OBS database [38] through the European Copernicus Climate Change Service website (doi.org/10.24381/cds.151d3ec6) and interpolated on the location of the 34 gauge stations by bilinear interpolation. The database collects time series from European national weather services, subjects them to extensive quality checks and blending procedures, and finally delivers a gridded daily dataset of several atmospheric variables onto a 0.1° × 0.1° regular grid. To compare the MFI, PCI, and SI indices with that of the rainfall erosivity factor R of the USLE, we also searched for 30-min rainfall data. These data are available for few stations and only after Gauge "point" daily data were mainly provided by the Regional Hydrographic Service of Abruzzo and only for the station of Avezzano are data derived from different sources (Regional Hydrographic Service of Campania, Regional Agrometeorological Center of Abruzzo, and Telespazio). This database was subjected to extensive quality check and homogenization procedures based on the Climatol library [40] to ensure the removal of outliers and the possibility of analyzing long-term trends. Full detail is given in [37].
Grid daily data were downloaded from the public E-OBS database [38] through the European Copernicus Climate Change Service website (doi.org/10.24381/cds.151d3ec6) and interpolated on the location of the 34 gauge stations by bilinear interpolation. The database collects time series from European national weather services, subjects them to extensive quality checks and blending procedures, and finally delivers a gridded daily dataset of several atmospheric variables onto a 0.1 • × 0.1 • regular grid. To compare the MFI, PCI, and SI indices with that of the rainfall erosivity factor R of the USLE, we also searched for 30-min rainfall data. These data are available for few stations and only after 1999. Considering the availability and completeness of the data, we selected 14 suitable stations at the Regional Hydrographic Service in the period 1999-2020. These 14 stations (as illustrated in Figure 1b) include 11 stations already considered in the daily datasets (i.e., Avezzano, Barisciano, Campotosto, Caramanico, Catignano, Chieti, L'Aquila, Pescara, Sulmona, Teramo, and Vasto) and 3 further stations (i.e., Arsita, Atri, and Bomba). Years with more than 5% missing records were discarded. The average number of years is 16.5, with minimum and maximum values of 8 and 22, respectively.
The median annual precipitation amounts of the stations were considered and obtained from both precipitation dataset is given in Figure 2. An illustration of the seasonality at each location is given in Figure S2. The two datasets depict quite similar spatial patterns with low rainfall amounts in the stations closest to the Adriatic coast. However, the grid data, derived from interpolation algorithms, present a much more reduced variability (range between 504-1085 mm) compared to point observations (range between 594-1407 mm). In particular, grid data describe a relatively regular gradient of increasing precipitation proportional to the Adriatic coast's distance. Therefore, the grid data do not capture either precipitation peaks (corresponding to Pescasseroli, Isola del Gran Sasso and Caramanico) or the low-precipitation stations in the inland valleys of the region (as illustrated in Figure 2a). According to previous analyses [39,41], these inland stations are characterized by the typical stau-foehn effects: the air masses of Adriatic and Tyrrhenian origin, pushed inwards by depression structures, discharge their humidity on the reliefs on the arrival side, and reach the intramountain areas depleted of precipitation.

Calculation of the Indices MFI and PCI, and of the USLE Erosivity Factor R
For each year, the Modified Fournier Index (MFI) [6] was computed as: (1)

Calculation of the Indices MFI and PCI, and of the USLE Erosivity Factor R
For each year, the Modified Fournier Index (MFI) [6] was computed as: where p i is the cumulated precipitation of the month i (i = 1, . . . , 12) and P is the cumulated annual precipitation. As explained in the introduction, there are several studies where MFI is adopted as an erosivity index. The index was also used in large-scale soil erosion risk assessments, such as in the Corine project [42]. One of the most commonly adopted classification criteria is reported in Table 1 [21,42]. MFI increases with both total annual precipitation and seasonality (also called "concentration"). The Precipitation Concentration Index (PCI) proposed by [7] can be computed according to the following equation: Therefore, PCI can also be computed as MFI/P. PCI ranges from 8.3% (completely uniform precipitation during the months) to 100% (precipitation concentrated in a single month). The classification of PCI given in Table 1 is taken from Lukic et al. [24] and is very similar to that originally proposed by [7].
The Seasonality Index (SI) [8] is defined as: The commonly used classification criteria of SI are given in Table 1 [43]. Both PCI and SI are often considered together with MFI in the literature to improve the interpretation of the results. The USLE erosivity factor R (MJ mm ha −1 h −1 year −1 ) was computed for the 14 stations indicated in Section 2.1 following the procedure indicated by Wischmeier and Smith [2]: where (I 30 ) k (mm h −1 ) is the maximum rainfall intensity of the storm k in a period of 30 min, E (MJ ha −1 ) is its total kinetic energy, m is the number of storms in the year j, and N is the number of years. An individual storm is preceded and followed by 6 h or more of no rain, moreover, in Equation (4), only storms with a total rainfall amount ≥ 12.7 mm or with at least 6.35 mm in 15 min are considered. The total storm kinetic energy is computed as: where n is the number of time segments in the storm, h i (mm) is the rainfall depth for each increment i, and e i (MJ ha −1 mm −1 ) is the energy per unit rainfall [44,45]: where I i (mm h −1 ) is the rainfall intensity in the increment i.

Relationships between Rainfall Erosivity Factor R and MFI, PCI, and SI
Before analyzing the spatio-temporal behavior of the MFI, PCI, and SI indices, we evaluated the relationships between the mean annual rainfall erosivity factor R (MJ mm ha −1 h −1 year −1 ) and the mean annual value of those indices for the 14 stations indicated in Figure 1b. Figure 3 shows the scatter plots of the mean R vs. the mean value of the indices (MFI, PCI, and SI, respectively); the regression lines with the corresponding coefficient of determination are also reported inset. From Figure 3a, it is evident that the mean values of MFI and R are highly correlated in the Abruzzo region. Based on the Student's t test, the regression line's slope presented in Figure 3a is significantly different from 0 (α = 0.01). The coefficient of determination (R 2 = 0.83) indicates that the regression line explains a large part of R variance. The Root Mean Square Error (RMSE) of this model is 220.2 MJ mm ha −1 h −1 year −1 . A high correlation between MFI and R was already observed in Abruzzo region [16]. However, the linear model obtained in the present analysis is more reliable due to the higher number of stations considered (14 instead of 9) and the more extended time series. Moreover, in [16], From Figure 3a, it is evident that the mean values of MFI and R are highly correlated in the Abruzzo region. Based on the Student's t test, the regression line's slope presented in Figure 3a is significantly different from 0 (α = 0.01). The coefficient of determination (R 2 = 0.83) indicates that the regression line explains a large part of R variance. The Root Mean Square Error (RMSE) of this model is 220.2 MJ mm ha −1 h −1 year −1 . A high correlation between MFI and R was already observed in Abruzzo region [16]. However, the linear model obtained in the present analysis is more reliable due to the higher number of stations considered (14 instead of 9) and the more extended time series. Moreover, in [16], the R calculation was based on the RUSLE approach [4], which tends to underestimate the energy associated with medium-low intensity rainfall events [45]. In some literature, the relationship between R and MFI was modeled by means of power models [13]. For our data, the power model has equation R = 2.42·MFI 1.41 , whose exponent is equal to that obtained by [13]. However, the coefficient of determination is only slightly better than the linear model (R 2 = 0.84), and the RMSE is higher (RMSE = 256.4 MJ mm ha −1 h −1 year −1 ).
The relationships obtained for PCI and SI (Figure 3b and 3c, respectively) indicate only a slight positive correlation between these indices and R. The scatterplots are also very similar as a consequence of the very high correlation between PCI and SI (r = 0.94). Based on the Student's t test, only the relationship obtained for PCI can be considered significant at the 0.1 significance level. For the model of Figure 3b, the RMSE is 473 MJ mm ha −1 h −1 year −1 , therefore much higher than those obtained for MFI-based models.
Two models of multiple regression of R were also tested, the first as a function of MFI and PCI (Equation (7)), the other as a function of MFI and SI (Equation (8)). These multiple regression models have the following equations, respectively: In both models, the independent variables are all significant and not significantly correlated. Based on the corresponding RMSE (163.6 and 169.6 MJ mm ha −1 h −1 year −1 , respectively), the model of Equation (7) can be considered the best performing.
This preliminary analysis demonstrates that in the region considered, a model based on both MFI and PCI, which only relies on monthly rainfall data, provides reliable estimates of the mean annual R.

Spatial Variability of MFI, PCI, and SI
The median values of MFI, based on the point and E-OBS datasets, are shown in Figure 4. The colors associated with each station correspond to the five classes of MFI (as illustrated in Table 1). The MFI spatial pattern recall that of precipitation (as illustrated in Figure 2) with lower values along the coast and in inland valleys, demonstrating that in the region, the MFI index (and therefore, erosivity) is influenced more by the cumulative annual rainfall than by its distribution throughout the year. According to the classes presented in Table 1, Figure 4 shows that the erosive risk can be defined as low or moderate in about half of the years (the median values are reported). ure 4. The colors associated with each station correspond to the five classes of MFI (as illustrated in Table 1). The MFI spatial pattern recall that of precipitation (as illustrated in Figure 2) with lower values along the coast and in inland valleys, demonstrating that in the region, the MFI index (and therefore, erosivity) is influenced more by the cumulative annual rainfall than by its distribution throughout the year. According to the classes presented in Table 1, Figure 4 shows that the erosive risk can be defined as low or moderate in about half of the years (the median values are reported). The comparison between Figure 4a,b shows only an apparent agreement between the median MFI values obtained from the two datasets. The differences can be appreciated by the confusion matrix of Table 2, which indicates that the classification agreement based on the median MFI is only 41%. It was also found that the disagreement increases with the station's altitude: for the 20 stations with elevation ≥ 500 m asl, the classification agreement is only 30% as opposed to 57% for the rest of the stations. As regards PCI (as illustrated in Figure 5), all stations' median values fall in the moderate class (as illustrated in Table 1), both considering the point and the grid dataset. The modest variability of PCI at a regional level contributes to the agreement, which remains also considering the maximum (instead of median) values of the series. The comparison between Figure 4a,b shows only an apparent agreement between the median MFI values obtained from the two datasets. The differences can be appreciated by the confusion matrix of Table 2, which indicates that the classification agreement based on the median MFI is only 41%. It was also found that the disagreement increases with the station's altitude: for the 20 stations with elevation ≥ 500 m asl, the classification agreement is only 30% as opposed to 57% for the rest of the stations. As regards PCI (as illustrated in Figure 5), all stations' median values fall in the moderate class (as illustrated in Table 1), both considering the point and the grid dataset. The modest variability of PCI at a regional level contributes to the agreement, which remains also considering the maximum (instead of median) values of the series. The seasonality index SI (as illustrated in Figure 6) has median values mainly falling in the 3rd class "Rather seasonal with a short drier season" (as illustrated in Table 1) with a very high agreement (94%) between the two datasets. The elevation doesn't affect the agreement, while it varies considerably as the severity of SI values increases; in fact, the agreement between the classifications obtained from the two datasets reduces to 82% and 35% for the 3rd quartile and maximum, respectively.

Trend Analysis of MFI, PCI, and SI
The spatial results of the trend analysis based on the Mann-Kendall test (see Appendix A for reference) are given in Figures 7-9 for MFI, PCI, and SI, respectively. For MFI, the analysis of the results based on the point dataset indicates prevailing increasing The seasonality index SI (as illustrated in Figure 6) has median values mainly falling in the 3rd class "Rather seasonal with a short drier season" (as illustrated in Table 1) with a very high agreement (94%) between the two datasets. The elevation doesn't affect the agreement, while it varies considerably as the severity of SI values increases; in fact, the agreement between the classifications obtained from the two datasets reduces to 82% and 35% for the 3rd quartile and maximum, respectively. The seasonality index SI (as illustrated in Figure 6) has median values mainly falling in the 3rd class "Rather seasonal with a short drier season" (as illustrated in Table 1) with a very high agreement (94%) between the two datasets. The elevation doesn't affect the agreement, while it varies considerably as the severity of SI values increases; in fact, the agreement between the classifications obtained from the two datasets reduces to 82% and 35% for the 3rd quartile and maximum, respectively.

Trend Analysis of MFI, PCI, and SI
The spatial results of the trend analysis based on the Mann-Kendall test (see Appendix A for reference) are given in Figures 7-9 for MFI, PCI, and SI, respectively. For MFI, the analysis of the results based on the point dataset indicates prevailing increasing

Trend Analysis of MFI, PCI, and SI
The spatial results of the trend analysis based on the Mann-Kendall test (see Appendix A for reference) are given in Figures 7-9 for MFI, PCI, and SI, respectively. For MFI, the analysis of the results based on the point dataset indicates prevailing increasing tendencies in the areas closest to the coast and prevailing decreasing tendencies in the internal area. Increasing significant trends (α = 0.1) are detected in only nine stations in the northeast, while significant negative trends characterize four stations in the internal areas.
tendencies in the areas closest to the coast and prevailing decreasing tendencies in the internal area. Increasing significant trends (α = 0.1) are detected in only nine stations in the northeast, while significant negative trends characterize four stations in the internal areas.
The analysis of the results obtained from grid data (as illustrated in Figure 7b) indicates widespread MFI increments, which are significant in most cases (25 out of 34 stations). A certain climatic differentiation between the stations closest to the coast and the inland ones is appreciable, but the differences are much less marked than those detected with the point dataset. Therefore, the two datasets provide a quite different assessment of the MFI time trend in the region. This difference can also be appreciated by the average regional values of Theil-Sen slope estimators, which are 0.14 and 0.69 for the point and grid datasets, respectively. As regards PCI (as illustrated in Figure 8), the analysis indicates a prevalent stationary condition in the region, with the presence of only four significant negative trends in some inland areas for the point dataset. Therefore, the two datasets provide for similar PCI information both in terms of classification (as illustrated in Figure 5) and trends (as illustrated in Figure 8). Considerations similar to those made for PCI apply to SI (as illustrated in Figure 9). In this case, there are only two significant negative trends on the point dataset and a general agreement between the two datasets.  Considerations similar to those made for PCI apply to SI (as illustrated in Figure 9). In this case, there are only two significant negative trends on the point dataset and a general agreement between the two datasets.

Discussion and Conclusions
The analysis reveals that, among the indices considered, MFI is the one characterized by the greatest spatial (as illustrated in Figure 4) (as and temporal illustrated in Figure 7) The analysis of the results obtained from grid data (as illustrated in Figure 7b) indicates widespread MFI increments, which are significant in most cases (25 out of 34 stations). A certain climatic differentiation between the stations closest to the coast and the inland ones is appreciable, but the differences are much less marked than those detected with the point dataset. Therefore, the two datasets provide a quite different assessment of the MFI time trend in the region. This difference can also be appreciated by the average regional values of Theil-Sen slope estimators, which are 0.14 and 0.69 for the point and grid datasets, respectively.
As regards PCI (as illustrated in Figure 8), the analysis indicates a prevalent stationary condition in the region, with the presence of only four significant negative trends in some inland areas for the point dataset. Therefore, the two datasets provide for similar PCI information both in terms of classification (as illustrated in Figure 5) and trends (as illustrated in Figure 8).
Considerations similar to those made for PCI apply to SI (as illustrated in Figure 9). In this case, there are only two significant negative trends on the point dataset and a general agreement between the two datasets.

Discussion and Conclusions
The analysis reveals that, among the indices considered, MFI is the one characterized by the greatest spatial (as illustrated in Figure 4) (as and temporal illustrated in Figure 7) variability in the region. This variability greatly affects the degree of agreement between the results obtained starting from the point and grid precipitation datasets. In fact, for MFI, somewhat discordant evaluations are obtained, both for the classification (as illustrated in Figure 4 and Table 2), and for the trends present in the 1980-2018 series (as illustrated in Figure 7). For the PCI and SI indices, considering the modest variability in the region and the near absence of significant trends, there are no substantial differences in the assessments deriving from the two datasets.
The most interesting discussion therefore concerns MFI. As already highlighted in the results, the trend analysis (as illustrated in Figure 7) reveals a marked separation in the behavior of the stations closest to the coast compared to that of the internal ones. This separation is more evident for the point dataset, which should be considered more representative of the heterogeneity of the territory. In particular, the areas of the Adriatic side are characterized by increasing trends in the MFI index, while those of the Tyrrhenian side show opposite trends. This separation may be also viewed in the grid dataset, since the few stations with insignificant trends are practically all located in inland areas.
This climatic configuration of the Abruzzo region confirms previous reports. In particular, Di Lena et al. [46], whose analysis of the main components applied to the 1951-2009 series of the Standardized Precipitation Index (SPI), revealed the presence of two main patterns, the first more correlated to the coastal areas, and the second more correlated to the inland, mountainous areas. Furthermore, for the coastal pattern, there was a marked increasing precipitation trend recent decades (i.e., an increase in rainfall). This trend is also confirmed in the present study, which goes up to 2018. In fact, the significant positive trends in MFI are always associated with significant positive trends in the corresponding annual rainfall, both for the point dataset and for the grid dataset. Moreover, a previous trend analysis based on daily data over 1980-2019 [37] showed significantly positive trends of indices of extreme precipitation in the coastal part of the region, while this trend was insignificant for inland stations.
The few significant negative trends detected for MFI in inland areas (as illustrated in Figure 7a) are instead attributable to corresponding changes in the rainfall distribution mode, as reported by the corresponding negative trends of PCI (as illustrated in Figure 8a) and partly also of SI (as illustrated in Figure 9a). For these stations, there are no significant trends in rainfall, but only a reduction in the seasonal variability of the rains which therefore tends to mitigate their annual aggressiveness (i.e., MFI). This is consistent with the seasonal precipitation trends illustrated in Curci et al. [37], where a more uniform distribution along the year was found for the recent period, with the exception of summer months which are becoming drier. Despite the greater agreement between the two datasets in terms of PCI and SI, the grid dataset does not detect the aforementioned significant negative trends.
The Equation (7) model was used to evaluate the effective practical repercussions of the rain erosivity analysis conducted with the grid data. The model is used to estimate the R factor of each station based on the average values of MFI and PCI obtained from the two datasets (R point and R grid , respectively). These R values were also used to estimate with the USLE the average annual soil loss A = R K (t ha −1 year −1 ) expected for a reference bare plot (slope 9%, length 22 m and tilled in the direction of the slope), where for the soil erodibility factor K (t·h·MJ −1 ·mm −1 ), an average value of 0.032 was assumed [47].
The comparison between the R point and R grid values for each station is shown in Figure 10, where the stations are in ascending order of annual precipitation. From the comparison, it emerges that the most marked differences occur in a limited number of stations, and these are mostly situations in which R point is decidedly higher than R grid . These stations (Campotosto, Caramanico, Isola del Gran Sasso, Pescasseroli, Pietracamela, Roccaraso, S.Eufemia Maiella, and Termine) are all located in the internal mountainous areas (as illustrated in Figure 1) and are characterized by high annual precipitation values (and therefore of MFI), which are not well captured by the grid data. For these stations, the mean difference between R point and R grid is about 840 (MJ mm ha −1 h −1 yr −1 ). Excluding these stations, the regional mean values of R point and R grid are very similar (1604 and 1595 MJ mm ha −1 h −1 yr −1 , respectively). These values are consistent with the mean Italian values indicated in Panagos et al. [3].   Figure 10 also shows the interesting case of 3 stations (Goriano Sicoli, Barisciano, and Avezzano) in which Rgrid is decidedly higher than Rpoint. Also in this case, the grid data does not capture well the pluviometric anomaly characteristic of these stations, which are subject to relatively low values of average annual precipitation (and therefore of MFI) due to a typical intramountain stau-foehn phenomenon which limits the arrival of large-scale perturbations [39].
The estimate of the annual average soil loss A obviously presents the same trend shown in Figure 10, since it was obtained for a hypothesis of constant K equal to 0.032. For the same 8 stations previously highlighted, the average difference between Apoint = Rpoint · K and Agrid = Rgrid · K is 27 t ha −1 year −1 , while for the rest of the stations, the average values are aligned and equal to 51.7 and 51.4 t ha −1 year −1 for Apoint and Agrid, respectively.
The results obtained here for the Abruzzo region therefore show that the grid data allow only general indications on the spatial-temporal behavior of the indices correlated to the erosion of rain R, with consequent errors in the estimation of the potential soil loss. In many cases, these are mountain valleys and slope locations where an accurate characterization of the factors influencing erosive phenomena is particularly important.
It should also be noted that indices requiring aggregated input data on a monthly scale were used in this study. For this scale, it is likely that the lower accuracy of the grid data is mitigated compared to that of the daily data due to the compensation of random errors. In the case of erosion estimation models that require daily rainfall [19], the difference between the values of R estimated from the grid and from point data would probably be more relevant.
On the other hand, the grid data represent a very valuable resource for large-scale analyses and were already used in various circumstances for evaluations on the erosivity of rain using simplified procedures and models. Examples can be found in [32][33][34]48]. As  Figure 10 also shows the interesting case of 3 stations (Goriano Sicoli, Barisciano, and Avezzano) in which R grid is decidedly higher than R point . Also in this case, the grid data does not capture well the pluviometric anomaly characteristic of these stations, which are subject to relatively low values of average annual precipitation (and therefore of MFI) due to a typical intramountain stau-foehn phenomenon which limits the arrival of large-scale perturbations [39].
The estimate of the annual average soil loss A obviously presents the same trend shown in Figure 10, since it was obtained for a hypothesis of constant K equal to 0.032. For the same 8 stations previously highlighted, the average difference between A point = R point ·K and A grid = R grid ·K is 27 t ha −1 year −1 , while for the rest of the stations, the average values are aligned and equal to 51.7 and 51.4 t ha −1 year −1 for A point and A grid , respectively.
The results obtained here for the Abruzzo region therefore show that the grid data allow only general indications on the spatial-temporal behavior of the indices correlated to the erosion of rain R, with consequent errors in the estimation of the potential soil loss. In many cases, these are mountain valleys and slope locations where an accurate characterization of the factors influencing erosive phenomena is particularly important.
It should also be noted that indices requiring aggregated input data on a monthly scale were used in this study. For this scale, it is likely that the lower accuracy of the grid data is mitigated compared to that of the daily data due to the compensation of random errors. In the case of erosion estimation models that require daily rainfall [19], the difference between the values of R estimated from the grid and from point data would probably be more relevant.
On the other hand, the grid data represent a very valuable resource for large-scale analyses and were already used in various circumstances for evaluations on the erosivity of rain using simplified procedures and models. Examples can be found in [32][33][34]48]. As indicated in Wang et al. [36], the mismatch between grid data and gauge data as well as interpolation errors strictly depends on the differences in spatial scales. In particular, the authors [36] indicate that the two types of data lead to different probability distributions of daily rainfall. Furthermore, the grid data tends to soften extreme values and to increase the number of rainy days.
This suggests that for erosivity risk assessment and planning of interventions in regions with complex orography, it is impossible to ignore detailed and local rainfall data. The grid data may still be applicable but require preliminary validation and possibly calibration work.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/atmos12060657/s1, Figure S1. Flowchart of the data and elaborations carried out in this study. Figure S2. Distribution of precipitation aggregated by season at each of the 34 locations used in this study. Data Availability Statement: Grid daily data were downloaded from the public E-OBS database [37] through the European Copernicus Climate Change Service website (doi.org/10.24381/cds.151d3ec6). Point data [37] are available upon request.
Acknowledgments: In this section, you can acknowledge any support given which is not covered by the author contribution or funding sections. This may include administrative and technical support, or donations in kind (e.g., materials used for experiments).

Conflicts of Interest:
The authors declare no conflict of interest.

Mann-Kendall Nonparametric Trend Test
The null hypothesis (H0) of the Mann-Kendall test [49,50] is that no trend is present in the population from which the dataset under investigation is extracted, whilst the alternative hypothesis indicates the presence of an increasing or decreasing monotonous trend. Being nonparametric, the test is independent of the variable's distribution and is less influenced by the presence of outliers in the time series.
The test statistic indicated by S is given by: where N is the number of observations. y i and y j are consecutive values of the variable under study. The sign function is defined as: Under the null hypothesis and for N ≥ 8, S follows a normal distribution with 0 mean and a variance approximately equal to: The p-Value of the test is: where Φ is the cumulative probability function of the standardized normal distribution. The Mann-Kendall test allows us to identify the existence of a monotonous trend, but does not allow its measurement; for this reason, as usually done in the literature, the nonparametric Theil-Sen slope estimator b [51,52] was used to evaluate the slope of the straight lines interpolating the data. The Zyp [53] and Kendall [54] packages of the R statistical software were respectively used to calculate the Theil-Sen nonparametric estimator and the p-value of the Mann-Kendall test.