Detailed Analysis of Spatial–Temporal Variability of Rainfall Erosivity and Erosivity Density in the Central and Southern Pannonian Basin

: Estimation of rainfall erosivity ( RE ) and erosivity density ( ED ) is essential for understanding the complex relationships between hydrological and soil erosion processes. The main objective of this study is to assess the spatial–temporal trends and variability of the RE and ED in the central and southern Pannonian Basin by using station observations and gridded datasets. To assess RE and ED , precipitation data for 14 meteorological stations, 225 grid points. and an erosion model consisting of daily, monthly, seasonal, and annual rainfall for the period of 1961–2014 were used. Annual RE and ED based on station data match spatially variable patterns of precipitation, with higher values in the southwest (2100 MJ · mm · ha − 1 · h − 1 ) and southeast (1650 MJ · mm · ha − 1 · h − 1 ) of the study area, but minimal values in the northern part (700 MJ · mm · ha − 1 · h − 1 ). On the other hand, gridded datasets display more detailed RE and ED spatial–temporal variability, with the values ranging from 250 to 2800 MJ · mm · ha − 1 · h − 1 . The identiﬁed trends are showing increasing values of RE (ranging between 0.20 and 21.17 MJ · mm · ha − 1 · h − 1 ) and ED (ranging between 0.01 and 0.03 MJ · ha − 1 · h − 1 ) at the annual level. This tendency is also observed for autumn RE (from 5.55 to 0.37 MJ · mm · ha − 1 · h − 1 ) and ED (from 0.05 to 0.01 MJ · ha − 1 · h − 1 ), as for spring RE (from 1.00 to 0.01 MJ · mm · ha − 1 · h − 1 ) and ED (from 0.04 to 0.01 MJ · ha − 1 · h − 1 ), due to the inﬂuence of the large-scale processes of climate variability, with North Atlantic Oscillation (NAO) being the most prominent. These increases may cause a transition to a higher erosive class in the future, thus raising concerns about this type of hydro-meteorological hazard in this part of the Pannonian Basin. The present analysis identiﬁes seasons and places of greatest erosion risk, which is the starting point for implementing suitable mitigation measures at local to regional scales


Introduction
Climate change is considered to have the most challenging worldwide environmental impact and soil erosion is expected to vary in response to changes in climate variables [1]. As highlighted by Azari et al. [2], it has several important implications for soil resources that are reflected in overall variations in the rainfall amount, intensity, spatial-temporal distribution, and increment in extreme precipitation events. The most direct effects of climate change on soil erosion are changes in rainfall erosivity due to changes in rainfall intensity [3] and the occurrence of extreme precipitation [4][5][6].
Water erosion is one of the most prominent causes of land degradation [7,8]. As highlighted by Grillakis et al. [9], the interplay between soil and rainfall has an important role regarding the ecological, hydrological, and biogeochemical cycles, since rainfall-induced soil erosion significantly affects its functions, and alters the soil structure and organic matter content. Therefore, soil erosion can be perceived as a widespread ecological and environmental crisis that poses a serious threat to the future quality of human life [10]. Xu et al. [11] highlight that this type of hydro-meteorological hazard is highly related to many natural factors, including rainfall (amount and variability), topography, soil, and vegetation properties, as well as human activities that intensify erosion processes. Water erosion has many on-and off-site effects [12], but off-site effects may have greater social, economic, and environmental impacts [13]. Soil erosion impacts deleteriously on land resources and develops flood-related risks. It also affects water quality degradation through the transportation of pesticides, fertilisers, and nutrients by the sediment [8]. The Pannonian Basin is known as a major agricultural region in Europe [14], so the impact of soil erosion and land degradation may severely affect regional agricultural production and the national economy of the states in this region. Therefore, more attention should be oriented towards soil erosion caused by rainfall in this area.
Since numerous studies around the world have shown that rainfall is the main climatological factor leading to soil erosion, with raindrop splash erosion and runoff being the main characteristics that influence this process, contemporary studies have focused on observing parameters described as rainfall erosivity-RE, e.g., [15][16][17][18][19][20][21][22][23][24][25][26]. Certain authors recognise this parameter as the potential of raindrops to trigger soil erosion, where its estimation is fundamental for the understanding of the climatic vulnerability within a given region, e.g., [11,27]. This geomorphic process can be illustrated in Figure 1.
In the previous conceptual framework, a quantitative assessment of potential soil erosion caused by rainfall is used for predicting soil erosion and for improving soil and water conservation and management. According to Bezak et al. [26], RE is an index that quantitatively describes the rainfall-driven force for soil erosion and is often regarded as one of the most important factors affecting the spatial and temporal variability of various soil erosion displacement processes. RE is expressed as the product of the total kinetic energy of a rainfall event times the 30-min maximum rainfall intensity [28]. It is calculated using high-resolution pluviographic data from rain gauges [10,11,17,[29][30][31] or by simpler models based on daily, monthly, seasonal, and annual rainfall [10,11,17,[23][24][25][26][27][28][29][30]. RE is one of the key input parameters of the universal soil loss equation (USLE) and the revised universal soil loss equation (RUSLE) [11,19,21,26,32]. Additionally, erosivity density (ED) is observed as the proportion of RE to precipitation [11].
Climate change is expected to affect future rainfall attributes due to a rise in the humidity in the atmosphere during a warmer climate [33]. The Intergovernmental Panel on Climate Change (IPCC) document "Special Report on Climate Change and Land", presents evidence of how these changes substantially affect other processes. These processes encompass desertification (water scarcity), land degradation (soil erosion, vegetation loss, wildfire, permafrost thaw), and food security (crop yield and food supply uncertainty) [23,34]. Moreover, an increase in severe storm intensity will likely enhance runoff and decrease the soil water infiltration process, thus causing even greater soil losses than at the beginning of the 21st century [19,20,[23][24][25]30]. Soil erosion already accounts for most soil loss, compared to other erosion processes in Europe, even though RE had only increased by 4% during the last decades of the 20th century [23,25]. In the previous conceptual framework, a quantitative assessment of potential soil erosion caused by rainfall is used for predicting soil erosion and for improving soil and water conservation and management. According to Bezak et al. [26], RE is an index that quantitatively describes the rainfall-driven force for soil erosion and is often regarded as one of the most important factors affecting the spatial and temporal variability of various soil erosion displacement processes. RE is expressed as the product of the total kinetic energy of a rainfall event times the 30-min maximum rainfall intensity [28]. It is calculated using high-resolution pluviographic data from rain gauges [10,11,17,[29][30][31] or by simpler models based on daily, monthly, seasonal, and annual rainfall [10,11,17,[23][24][25][26][27][28][29][30]. RE is one of the key input parameters of the universal soil loss equation (USLE) and the revised universal soil loss equation (RUSLE) [11,19,21,26,32]. Additionally, erosivity density (ED) is observed as the proportion of RE to precipitation [11].
Climate change is expected to affect future rainfall attributes due to a rise in the humidity in the atmosphere during a warmer climate [33]. The Intergovernmental Panel on Climate Change (IPCC) document "Special Report on Climate Change and Land", presents evidence of how these changes substantially affect other processes. These processes encompass desertification (water scarcity), land degradation (soil erosion, vegetation loss, wildfire, permafrost thaw), and food security (crop yield and food supply uncertainty) [23,34]. Moreover, an increase in severe storm intensity will likely Numerous respective researchers worldwide discussed the spatial-temporal variability of rainfall erosivity in study regions: from world study conducted by Borelli et al. [35] to Australia [36,37] and New Zeland [38], Africa [39][40][41], Brasil [42,43], Columbia [44], China [10,11,17,30,45], South Korea [45,46], Iran [2,47], and Uzbekistan [48]. Small-scale studies in the Pannonian Basin [8], Tuscany [49], Rhone region [50], Greece [21], Czech Republic [51], and the Netherlands [52] also suggest a positive trend in RE during the last 50 years. What all the studies have in common is that rainfall erosivity poses a significant environmental threat that needs to be thoroughly addressed.
Since the Pannonian Basin is an intensive agriculture area, the present study extends previous research in the Basin [8] to gain a better insight into soil loss in the light of precipitation variability. Despite limited data in the Western Balkans, this research strives to fill a gap in information for the central and southern parts of the Pannonian Basin in Southeast Europe (northern Serbia-Vojvodina, Hungary, and the eastern part of Croatia).
The main objectives of this study are: (1) to investigate spatial-temporal trends and variability of the RE and ED in the central and southern Pannonian Basin by using different data repositories (ECA&D and E-OBS), and (2) to emphasise the importance of the influence of the large-scale processes of climate variability. Particular attention is placed on multidecadal trends and variability in these two indices. The findings will add to European erosivity studies by filling in gaps in the erosivity map provided by Panagos et al. [19,20], The investigated area has a high diversity of physical, biophysical, and socioeconomic conditions. During the past few decades, an increasing number of extreme hydrometeorological events (heat waves, droughts, heavy precipitation, etc.) have been reported [8,[56][57][58][59][60][61], favouring increased rainfall erosivity for most of the investigated area, except in the northwestern parts. Using four pluvial indices, Lukić et al. [8] showed that precipitation extremes have been rising in terms of totals and intensities. Climate change projections suggest that there will be consequences within multiple sectors of the observed Pannonian Basin, counting as well the lethal impacts on its ecosystem services. Being that this part of Europe represents one of the regional hotspots, it is crucial to say that it is regarded as the region with the highest number of severely affected sectors in the whole of Europe [62]. Since the Pannonian Basin is located in the Northern Hemisphere mid-latitudes, the lowest air temperature values occur in January and the highest in July and August. The average annual air temperature fluctuates between 9.7 • C in Hungary to 11 • C in Vojvodina and 12.5 • C in eastern Croatia. The rainfall season ranges from mid-April to September, and cloud cover is highest from November to January. The average amount of precipitation is approximately 600 mm per year [8]. Wind speed was found to decrease in every season throughout the year, while relative humidity is increased in autumn and decreased in Sustainability 2021, 13, 13355 5 of 31 spring, summer, and winter during the period of 1960-2010 [54]. On the other hand, sunshine duration has the opposite trend [8].
The investigated area has a high diversity of physical, biophysical, and socioeconomic conditions. During the past few decades, an increasing number of extreme hydrometeorological events (heat waves, droughts, heavy precipitation, etc.) have been reported [8,[56][57][58][59][60][61], favouring increased rainfall erosivity for most of the investigated area, except in the northwestern parts. Using four pluvial indices, Lukić et al. [8] showed that precipitation extremes have been rising in terms of totals and intensities. Climate change projections suggest that there will be consequences within multiple sectors of the observed Pannonian Basin, counting as well the lethal impacts on its ecosystem services. Being that this part of Europe represents one of the regional hotspots, it is crucial to say that it is regarded as the region with the highest number of severely affected sectors in the whole of Europe [62].

Data
Daily rainfall data from meteorological stations were obtained from the European Climate Assessment and Dataset (ECA&D) to estimate RE and ED for the period of 1961-2014. In this data repository, daily precipitation values (RR) are calculated as the sum of RR of 18 UT of the current day and RR of 06 UT of the next day. The snow was not included in RR since the indices used in the study do not require this parameter. Only validated RR values were used, indicated by the flag 0 in the column labelled quality control (QC) from the ECA&D data repository [63]. Daily ECA&D RR series at each site were tested and homogenised by Wijngaard et al. [64].
As pointed out by Thornton et al. [65], gridded weather products are generally used to support ecological, agricultural, water resources management, and climate change studies, more particularly in regions with sparse meteorological stations. Being that gridded data exhibit a high resolution on spatial and temporal distributions of weather products, the European Observations (E-OBS) datasets were used in this study. The latest version v23.1e, in regular latitude-longitude coordinates and regular horizontal resolution of 0.25 • × 0.25 • were used as well. The produced database, given as part of the EU ENSEMBLES project, is comprised of daily mean, minimum, and maximum temperature values and daily precipitation, respectively. Total daily amounts of precipitation (RR) in mm (measured as the height of the equivalent liquid water in a square meter) are included within the observed datasets. The data sources for the precipitation are rain gauge data which do not possess a uniform way of defining the 24 h, after which the precipitation measurements were made. Therefore, there is no uniform time for the measurement period, which could be attributed to the daily precipitation. The Global 30 Arc-Second Elevation Dataset (GTOPO30) is used for the elevation file. It represents a global raster digital elevation model (DEM) with a horizontal grid spacing of 30 arc seconds (approximately 1 kilometre) that is developed by USGS [63,66]. The above-described dataset is comprised of 225 grids and was used to produce a high spatial resolution of RE and ED maps on multiannual and seasonal scales.
The monthly, seasonal, and annual series were derived from daily series. Standard season definitions are used, namely: winter (DJF), spring (MAM), summer (JJA), and autumn (SON). A list of the meteorological stations used in the study is given in Table 1.

Methodology
All steps, including procedures and different quantitative approaches used in this research, are presented in the following workflow ( Figure 3).
Numerous indices have been used to quantify water erosion, among which the most prominent are Fourier index (FI) and modified Fourier index (MFI). These indices are based on monthly rainfall data averages since high-quality datasets on a regional scale are very scarce [8,55]. R factor reflects the potential capability of precipitation to cause soil loss and is one of the most important factors for estimating soil erosion. In its humblest formula, the R factor stands for an average annual value and is quantified by summing the event-based energy-intensity values (EI 30 ) for a specific site divided by the number of observed years [10,11]. The relationship between FI and MFI with the USLE R has been discussed [8,55]. The main difference between the previously mentioned indices and the one used here is in the erosional rate estimation that is best represented by RE and ED [10,11,17,67]. The value of RE was calculated following Zhang et al. [68]: where RE i is the rainfall erosivity in i-th half-month (MJ·mm·ha −1 ·h −1 ). Each month was separated into two half-months: from 1st-15th and from 16th to the end of the month; dividing a year into 24 half-months. P j represents the erosive rainfall on the j-th day (mm) when actual rainfall was >12 mm. If the actual daily rainfall was <12 mm, then the P j was counted as zero. Parameter k represents the number of days in the i-th half-month [68]. The values of α and β, as empirical parameters, were determined as follows: where P d12 and P y12 present the average daily/annual erosive rainfall when daily rainfall is >12 mm. ED is the ratio of rainfall erosivity to precipitation [67]. It measures the erosivity per rainfall (in mm) and is expressed as MJ·ha −1 ·h −1 . The equation can be applied on a monthly or yearly scale (i).  The direct measurement of kinetic energy (KE) requires sophisticated and costly instruments. Methods utilized for the estimation of KE, based on power functions, are proven to be valid in practice. Due to the limited availability of hyetograph and disdrometer data, statistical models that use available data on daily precipitation were created. Three main aspects of the R factor are considered when performing soil erosion estimations. One is the average annual rainfall erosivity for the prediction of average annual soil loss. The second one is based on the seasonal distribution values of rainfall erosivity, and the third one represents the event or daily rainfall erosivity, which is crucial for the analysis of the soil loss recurrence and its assessment. Models that use events or daily rainfall as measurements for the calculation of erosivity, exhibited at all temporal scales, are significant. Moreover, those can satisfy the three usage requirements mentioned previously in the context. A power-law function has been generally used for the estimation of erosion index EI 30 gathered from the event or daily rainfall amounts. A connection between EI 30 index, obtained for the half-month period, and the sum of the estimated daily erosivity for the same period, has been tested by Xie et al. [69]. This approach was considered the average annual rainfall erosivity and its seasonal variations, while attention was not particularly given to the estimations of daily rainfall erosivity. It is crucial to say that this approach is suitable for the analysis of seasonal differences in rainfall characteristics. As highlighted by Teng et al. [70], a power function model observed by Zhang et al. [68] showed that the estimated R factor, used in this model, aligns well with the EI 30 index (R 2 = 0.72). Zhang et al. [68] also reported a strong correlation between the parameters α and β (R 2 = 0.95). Additionally, analysis of the maximum daily precipitation data indicates that a defined 12.0 mm threshold (used for a daily rainfall amount model) is similar to the threshold of 12.7 mm, as suggested by Wishmeier and Smith [71]. Therefore, this proposal can be used as a practical threshold for the separation of erosive and nonerosive storms within a given area [72]. Since (1) empirical parameters α and β are closely related to climate characteristics and topography of the study area, and (2) the presence of good agreement between R factor and EI 30 due to high-resolution data scarcity, the authors decided to apply Zhang's power-law function model for the estimation of RE and ED.
Trends in time series of RE and ED were analysed via the nonparametric Mann-Kendall (M-K) test. The M-K test is widely used because it does not require a normal distribution of data and is straightforward to calculate [17]. Two hypotheses are set here: the null hypothesis (H 0 ), which states that there is no trend in the time series; the alternative hypothesis (H a ), which states that there is a significant trend in the series, for a given level of significance. Here, the trend is regarded as statistically significant at the following confidence levels (p): 0.1, 0.05, and 0.01, respectively. If the p-value is lower than the specified level, H 0 should be rejected, accepting H a and vice versa. To prevent the effects of the serial correlation in observed time series, the Yue-Pilon method was applied [73].
To examine the relationship between the number of heavy precipitation days (R10mm) and potential climate driver North Atlantic Oscillation (NAO), Pearson correlation (r) with 0.01 and 0.05 confidence levels was utilized following the approach of Lukić et al. [8]. To examine the relationship between extreme precipitation index (R10mm) and NAO, linear correlations (the Pearson correlation coefficient) were utilised. Within this methodological approach, a correlation of +1 or −1 occurs when each of the variables is perfectly correlated to the other. On the other hand, a correlation close to zero means that there is a weak relationship between the examined variables. The values of NAO were obtained from the Climatic Research Unit of the University of East Anglia following the approach of Jones et al. [74].
The spatial distribution of RE and ED was produced using the radial basis function interpolation method with the inverse multiquadric (RBF IMQ) kernel. This technique tends to estimate the wide rainy areas and smooth rainfall fields using a small number of observed data to obtain higher accuracy. Contrary to the conventional inverse distance weighting (IDW) interpolation method or RBF with inverse quadratic (IQ) kernel, which tends to approximate the wide rainy areas and smooth rainfall field when using a small number of observations, the RBF IMQ does not flatten the field, even when the same smoothing parameter is used as in the usual RBF and IDW [75]. The same method for visualisation of the results was used in previous studies regarding rainfall erosivity (e.g., Lukić et al. 2018 [52], and Lukić et al. 2019 [8]).

Intra-Annual Distribution of RE and ED
The annual distribution of RE for the central and southern Pannonian Basin was binned and plotted by decade to better highlight temporal variability in RE (Figure 4a). The boxen plots for RE (seaborn module, boxen plot) are conducted for gridded data, due to the data volume (225 grids). As boxen plots are showing minimum, maximum, and median in the dataset boxes, and plot more quantiles, they provide a better insight into the shape of the distribution, particularly in the tails (Figure 4d-f). Overall, maximum mean long-term RE values were recorded for stations Pécs (2100 MJ·mm·ha −1 ·h −1 ) and Osijek (2093 MJ·mm·ha −1 ·h −1 ). Interestingly, the maximum value for the last decade occurred in Vršac (3739 MJ·mm·ha −1 ·h −1 ), which had a similar RE distribution to other stations until the most recent decade. Minimal average value through the decades (Figure 4a), as well as variation in RE, is recorded for Budapest meteorological station (703 MJ·mm·ha −1 ·h −1 ).
Paying attention to the median on the boxen plot ( Figure 4d), we can highlight that RE values are slightly rising since the 1980s, with their maximum detected in the last decade. This finding is consistent with station data, as well as with the previous studies [8,55,[76][77][78][79][80].
The seasonal variability of RE is pronounced. Summer has the highest mean RE values (709 MJ·mm·ha −1 ·h −1 ), followed by the means for autumn (351 MJ·mm·ha −1 ·h −1 ), spring (314 MJ·mm·ha −1 ·h −1 ), and winter (139 MJ·mm·ha −1 ·h −1 ) (Figure 4b). Szombathely is the station with the highest summer RE (975 MJ·mm·ha −1 ·h −1 ), and Budapest the lowest (260 MJ mm·ha −1 ·h −1 ). Only three meteorological stations have autumn mean RE higher than 500 MJ·mm·ha −1 ·h −1 : Osijek (548), Pécs (548), and Szombathely (566). Although similar to summer as a wet part of the year, the spring has a modest RE (314 MJ·mm·ha −1 ·h −1 ) compared with the autumn mean. The meteorological station with the lowest winter RE is Szombathely, with an average of 67 MJ·mm·ha −1 ·h −1 , making it the site with the largest range in seasonal mean RE between summer and winter. The highest winter mean RE values occurred in Osijek (217.4 MJ·mm·ha −1 ·h −1 ). Gridded data are following the station season distributions (Figure 4e).
Monthly RE station distributions also exhibit variability between the wettest (April to October) and driest (October to March) parts of the year [8]. Figure 4c shows a rapid increase in RE from April, reaching a peak in July and then progressively decreasing until October. After October, RE values are equalised as a result of the dry season. The highest RE (1783 MJ·mm·ha −1 ·h −1 ) was found for Vršac station in the eastern part of the Vojvodina region. On the boxen plots for gridded data (Figure 4e,f), this value is well above the median for July (see also Figure 4a,b in the last decade of station data for summer). This value should not be interpreted as an outlier, since it stands as an extreme event recorded by a few previous studies [8,61], as well as by the Republic Hydrometeorological Institute of Serbia [81]. Interesting to notice is the boxen plot shape from May to October. This elongation of the third quartile has more sense when referring to the decade plot (Figure 4d), where the increasing trend was observed from the 1980s to 2010s.

RE and ED Based on Station Data
The development of seasonal erosivity maps contributes to the estimation of the seasonal soil loss ratio [82]. Different seasonal spatial patterns of RE can be observed in Figure 5. During spring (March-May) and autumn (September-November), the study area is broadly divided into regions of higher RE in the southwest and lower RE in the northeast (Figure 5a,c). This gradient is less apparent during autumn, as Vojvodina (northern Serbia) has higher RE during autumn compared with spring. Both spring and autumn have high RE values (>400 MJ·mm·ha −1 ·h −1 ) in Osijek and Pécs, but the autumn maximum RE is in Szombathely (566 MJ·mm·ha −1 ·h −1 ). The lowest RE is seen in Budapest (128 MJ·mm·ha −1 ·h −1 ) in spring and Debrecen (226 MJ·mm·ha −1 ·h −1 ) in autumn. Monthly RE station distributions also exhibit variability between the wettest (April to October) and driest (October to March) parts of the year [8]. Figure 4c shows a rapid area is broadly divided into regions of higher RE in the southwest and lower RE in the northeast (Figure 5a,c). This gradient is less apparent during autumn, as Vojvodina (northern Serbia) has higher RE during autumn compared with spring. Both spring and autumn have high RE values (>400 MJ·mm·ha −1 ·h −1 ) in Osijek and Pécs, but the autumn maximum RE is in Szombathely (566 MJ·mm·ha −1 ·h −1 ). The lowest RE is seen in Budapest (128 MJ·mm·ha −1 ·h −1 ) in spring and Debrecen (226 MJ·mm·ha −1 ·h −1 ) in autumn. More than 70% of the annual RE in the central and southern Pannonian Basin occurs during the wet period, from April to October ( Figure 6). It can be seen that the RE values in the eastern and western parts of the Pannonian Basin are larger than those in the northern and central parts. The lowest average monthly RE values are observed in February (33 MJ·mm·ha −1 ·h −1 ) and January (39 MJ·mm·ha −1 ·h −1 ), with the highest in June (285 MJ·mm·ha −1 ·h −1 ) and July (246 MJ·mm·ha −1 ·h −1 ). The monthly spatial patterns of RE, in general, display a continuous increase in RE from winter to spring, followed by a sharp high transition in the summer and a gradual decrease from autumn to winter. Identification of RE hotspots for specific times and places can guide investments in soil protection and help to avoid tillage during the most at-risk times [21]. However, it should be noted that extreme erosive rainfall events are not the only drivers of soil erosion. Other factors, such as topography, soil characteristics, and vegetation cover, also affect erosion rates, and, thus, should be addressed in future studies [31]. Sustainability 2021, 13, x FOR PEER REVIEW 14 of 36 Figure 6. Spatial distribution of monthly RE from (a-l) January-December in the study area based on station data.
In the central and southern Pannonian Basin, the rainy season lasts from May to July, with high spatial variability in annual precipitation. Up to 1600 mm of annual rainfall can occur in the Ukrainian Carpathians and the Tatra Mountains, compared with 550 mm in the plains of eastern Hungary. The lowest summer precipitation is registered in northern Serbia (80 mm), while the lowest winter precipitation values are observed for northern Romania (70 mm) [83]. Increasing trends in total annual rainfall over the last 50 years have been reported for the region, except for northern Hungary [8].
Annual ED and RE follow the precipitation distribution with strong spatial variability (Figure 7a Regions with high ED are exposed to the risk of flooding and water scarcity as a result of their infrequent, but very intense and erosive, rainstorms [31]. In the central and southern Pannonian Basin, the rainy season lasts from May to July, with high spatial variability in annual precipitation. Up to 1600 mm of annual rainfall can occur in the Ukrainian Carpathians and the Tatra Mountains, compared with 550 mm in the plains of eastern Hungary. The lowest summer precipitation is registered in northern Serbia (80 mm), while the lowest winter precipitation values are observed for northern Romania (70 mm) [83]. Increasing trends in total annual rainfall over the last 50 years have been reported for the region, except for northern Hungary [8].
Annual ED and RE follow the precipitation distribution with strong spatial variability (Figure 7a-c). High ED values (>3 MJ ha −1 h −1 ) indicate a greater risk of erosive rainstorms, soil erosion, and flooding [17,31,84]. The highest annual ED values for the study area are obtained for the southwest, in Osijek and Pécs (3.3 MJ ha −1 h −1 ), east in Debrecen (2.9 MJ ha −1 h −1 ), and southeast in Banatski Karlovac (2.8 MJ ha −1 h −1 ). Where ED values are >1, a certain precipitation amount may cause relatively higher RE. The lowest ED values occur in the northern (Budapest with 1.2 MJ ha −1 h −1 ) and central parts of the study area (e.g., Kikinda station in northern Vojvodina with 2.3 MJ ha −1 h −1 ). Regions with high ED are exposed to the risk of flooding and water scarcity as a result of their infrequent, but very intense and erosive, rainstorms [31].  The annual distribution of RE for the central and southern Pannonian Basin was also obtained for the gridded E-OBS datasets. The obtained result indicates good spatial agreement with the average annual precipitation for the investigated period (Figure 7d). The values are increasing from east to west, ranging from 250 MJ·mm·ha −1 ·h −1 (for the Great Hungarian Plain) to 2800 MJ·mm·ha −1 ·h −1 in the southwestern parts of Hungary (Figure 7e). The zones with the highest average ED (Figure 7f) correspond with the obtained RE values and range between 1.00 (the Great Hungarian Plain) and 3.75 MJ·ha −1 ·h −1 (Southern Transdanubia region and western parts of the Pannonian part of Croatia).
During the summer season, the RE exhibits its largest seasonal values for western parts of Hungary and the Pannonian part of Croatia (400-1000 MJ·mm·ha −1 ·h −1 ), as well as southeastern parts of the Vojvodina Region in northern Serbia (500-700 MJ·mm·ha −1 ·h −1 ) (Figure 8b). Higher RE values (from 300-800 MJ·mm·ha −1 ·h −1 ) are also observed for the autumn season in the western parts of the study area (Figure 8c).  Figure 7. Spatial distribution of mean annual precipitation, RE, and ED for station data (a-c) and gridded data (d-f).

RE and ED Based on Gridded Data
The annual distribution of RE for the central and southern Pannonian Basin was also obtained for the gridded E-OBS datasets. The obtained result indicates good spatial agreement with the average annual precipitation for the investigated period (Figure 7d). The values are increasing from east to west, ranging from 250 MJ·mm·ha −1 ·h −1 (for the Great Hungarian Plain) to 2800 MJ·mm·ha −1 ·h −1 in the southwestern parts of Hungary (Figure 7e). The zones with the highest average ED (Figure 7f) correspond with the obtained RE values and range between 1.00 (the Great Hungarian Plain) and 3.75 MJ·ha −1 ·h −1 (Southern Transdanubia region and western parts of the Pannonian part of Croatia).
During the summer season, the RE exhibits its largest seasonal values for western parts of Hungary and the Pannonian part of Croatia (400-1000 MJ·mm·ha −1 ·h −1 ), as well as southeastern parts of the Vojvodina Region in northern Serbia (500-700 MJ·mm·ha −1 ·h −1 ) (Figure 8b). Higher RE values (from 300-800 MJ·mm·ha −1 ·h −1 ) are also observed for the autumn season in the western parts of the study area (Figure 8c).  The winter season exhibits relatively stable RE values ranging from 0 to 100 MJ·mm·ha −1 ·h −1 for the majority of the investigated area. The areas with higher RE values (varying between 100 and 200 MJ·mm·ha −1 ·h −1 ) are spotted in the southeastern and southwestern part of the study area (Figure 8d). These values are starting to divide into higher erosivity classes during the spring season (Figure 8a) following the previously observed spatial patterns extending from southeast to west (0-400 MJ·mm·ha −1 ·h −1 ).
The most erosive month is July, followed by June, August, and September. RE and ED are generally stable in the period from December to April, and then display an increase in values during the second period (Figure 9). higher erosivity classes during the spring season (Figure 8a) following the previously observed spatial patterns extending from southeast to west (0-400 MJ·mm·ha −1 ·h −1 ).
The most erosive month is July, followed by June, August, and September. RE and ED are generally stable in the period from December to April, and then display an increase in values during the second period ( Figure 9). The performance of the spatial interpolation was validated by the root mean square error (RMSE) and standard error (SE) for both station and gridded RE values. RMSE is acceptable for the months that contribute to 70% of total RE (from May to September) ( Table 2). It is expected that gridded data provide smaller errors due to much higher spatial resolution, as shown in Table 2. In combination with the SE, mean RE values can be used for estimating the intervals prediction. The performance of the spatial interpolation was validated by the root mean square error (RMSE) and standard error (SE) for both station and gridded RE values. RMSE is acceptable for the months that contribute to 70% of total RE (from May to September) ( Table 2). It is expected that gridded data provide smaller errors due to much higher spatial resolution, as shown in Table 2. In combination with the SE, mean RE values can be used for estimating the intervals prediction.

Discussion
Spatial variations in annual precipitation, RE and ED, the R factor from Borrelli et al. [35], and FI and MFI of Lukić et al. [8] are compared in Figure 10. The spatial variations in R calculated by Borrelli et al. [35] partially coincide with RE station values obtained in this study. When comparing Figure 10b,d, it becomes evident that higher rainfall erosivity values are located in the southwest of the Pannonian Basin, but with deviation in rainfall erosivity values. In Figure 10d, Osijek meteorological station has an annual RE between 950 and 1000 MJ·mm·ha −1 ·h −1 [35], and, in Figure 10b, it falls in the last range between 2000 and 2150 MJ·mm·ha −1 ·h −1 . The modelling results derived from the investigation of the changes in rainfall erosivity from Bezak et al. [25] (who used the COPERNICUS datasets and high-resolution erosive storms from the REDES database) indicate that the applied approach is showing a tendency to polish the annual erosivity values >1500 MJ·mm·ha −1 ·h −1 .
Explanations for this deviation in rainfall erosivity estimates may be due to different (1) databases and spatial-temporal resolutions in the studies, and (2) methodologies applied for R factor calculation. The R factor was computed using sub-hourly (61%) and hourly pluviograph data (39%), obtained from the global rainfall-runoff erosivity database (GloREDa) from 2001 to 2012 [35], whereas the present study used daily data from the ECA&D database from 1961 to 2014. The methodology for computing the R factor was taken from the USDA RUSLE handbook and based on the Brown & Foster [84] equation for rainfall kinetic energy. In this study, RE distribution was computed by following [68] daily precipitation distributions. Since the European continental zone is characterised by warm summers and cold winters, it can be seen as an area prone to the high variability of RE. Due to this, and different methodological approaches, the mean RE values obtained in this study tend to be higher compared to the R factor for the Pannonian biogeographical region (660 MJ·mm·ha −1 ·h −1 ) [19]. Nonetheless, spatial patterns of the calculated indices are in relatively good agreement with other studies [19,35].
The work of Xie et al. [69] shows that there is a strong correlation between the powerlaw-derived relationship of empirical parameters α and β (within the Zhang et al. [68] RE model). However, mainly β was associated with average daily (P d12 ) and annual erosive rainfall (P y12 ). These parameters can be significantly linearly related with the mean annual rainfall (as the independent variable), yet, on other occasions, display insignificant linearity. Yu [36] and Xie et al. [69] explained this in terms of the aggregation of the above-mentioned parameters, which were obtained month by month (in their respective studies), and by all months, as presented by Zhang et al. [68] and this study. If parameter β is not properly estimated, the deviation of the model output may occur. Hence, Xie et al. [69] indicated that one of the reasons for the overestimation when using the Zhang et al. [68] model may be attributed to the utilisation of a surrogate for the measured value of EI 30 . In this model, the surrogate was perceived as the product of daily rainfall P d and the maximum 10-min intensity in a day I 10d , multiplied by a conversion factor of 0.184. Although Xie et al. [69] observed that daily rainfall and maximum 10-min intensity in a day are highly correlated, the use of P d I 10d instead of EI 30 (due to data scarcity) contributes to the uncertainty of this relationship. Therefore, the results in this study emphasise that there is a great necessity for a more proper calibration for this model to be implemented in other places.
Yang and Yu [37] outlined that rainfall erosivity can vary depending on the time of year and location, especially when it comes to latitude and elevation. On the other hand, the RUSLE manual states that the R factor is mainly quantified using a defined threshold rainfall amount of 12.7 mm [85]. The observed discrepancies in the obtained R factor can be due to different rainfall thresholds (used here and in other European studies), thus leading to the increased values, as mean annual rainfall decreases because the relative contribution of small storm events is greater in dry areas. When it comes to the contemporary analysis of the spatiotemporal rainfall erosivity and its variability in Europe, Bezak et al. [31] highlighted that the regional temporal RE patterns do not display homogenous behaviour due to the relatively large differences among neighbouring meteorological stations, and also the sparse network is problematic and places a lot of dependence on the choice of interpolation technique. This leads to the conclusion that local climatological conditions determine the local temporal rainfall erosivity distribution, as seen in this work (some meteorological stations can show the similar temporal distribution of RE and ED while having quite different absolute RE and ED values). Bezak et al. [31] explained this pattern by the differences between rainfall erosivity synchrony scale (Rsync) and computed precipitation synchrony scale. Rsync is the maximum radius (in km) surrounding a meteorological station where at least half of the other meteorological stations record an erosive event.
When observing the assessment of rainfall energy, the FI and MFI results gained by Lukić et al. [8] should be discussed, as they have very good correspondence with RE spatial distribution in this study. According to respective authors, the central and southern Pannonian Basin area belongs to the erosive class of low values , with the tendencies of a progressive increase in the values of the erosion index at the annual level, hence leading to the possible transition to a higher erosive class and increase in the sensitivity to erosion by water processes in the given area (Figure 10e,f). According to these indices, the majority of the study area belongs to the low erosivity class, except in the northwest, which, in the future, could lead to a transition to a higher erosion category related to RE variation for this part of the Pannonian Basin.
On the other hand, climate gridded datasets are valuable and useful when performing detailed analysis of RE and ED at regional scales. As highlighted by Sidău et al. [86], usage of such data is important for the development of quality results that can be used as tools for climate services, policies, and strategies for climate change adaptation and mitigation. Obtained results correspond quite well with the results of Panagos et al. [19], who used data from 42 stations from Croatia and 30 from Hungary for the period 1961-2012 and 1998-2013, respectively. Respective authors argued that, on the country level, one of the highest calculated values of the R factor is observed for Austria and Croatia (mean values of >1000 MJ·mm·ha −1 ·h −1 ). In this study, it has been found that values for the Pannonian part of Croatia are increasing from east to west (ranging from 750 to 2800 MJ·mm·ha −1 ·h −1 ). Possible reasons for the difference in R factor variability between the studies could have occurred due to data availability, analysed period, data quality, sample sizes, interpolation method, and interpretation. An overall increase in the mean R factor for our study can be seen in the most western parts of the investigated area. This feature could have been influenced by the availability of long-term rainfall data with better temporal resolution from a higher number of weather stations and a denser spatial distribution than previous studies [19,20]. Moreover, respective authors [19] produced R factors of different temporal resolutions to 30-min resolution by applying a calibration factor for comparability of data from the study area. This approach may have caused a certain decrease in obtained R factors for Croatia and Hungary. Generally speaking, Mehr et al. [87] highlighted that observed trends for Osijek meteorological station (eastern part of Croatia) could pose a danger to regional agricultural production, and thus should be addressed with special attention when analysing the emergence of new climate extremes over the territory of the Pannonian Basin. Pásztor et al. [88] used the Pan-European Soil Erosion Risk Assessment (PESERA) model to perform rill and inter-rill erosion estimates at the regional level due to the lack of available in situ soil erosion datasets in Hungary. The authors identified high erosion rates in mountainous and hilly areas in the northern parts of Hungary and the southern Transdanubia region, highlighting that more than 26% of the country's total area is susceptible to moderate and high erosion rates. This was also the case for the extremely wet year of 2010, where R factor varied between 2300 and 3400 MJ·mm·ha −1 ·h −1 . Results from our study highlight that the above-mentioned areas experience rainfall erosivity ranging between 500 and 2800 MJ·mm·ha −1 ·h −1 , thus cor- responding well with the previous study. The Great Hungarian Plain generally exibits low erosion rates (250-500 MJ·mm·ha −1 ·h −1 ), with a few regions where mostly sandy soils dominate [88]. The high extent of the lowland areas is, therefore, affected only by relatively low erosion rates, while, on the other hand, significantly higher erosion rates can be seen in the hills and mountains. For northern Serbia (Vojvodina), variability of calculated RE values display certain contrasts between lowland (500-750 MJ·mm·ha −1 ·h −1 ) and mountainous/hilly areas, which exibit the highest susceptibility to rainfall erosivity during the observed period of 1961-2014 (1000-1250 MJ·mm·ha −1 ·h −1 ). These results are in good agreement with the study of Lukić et al. [61], who outlined that highest rainfall erosivity values can be seen for southeastern parts of the Vojvodina Region (Banat Region with Vršac weather station). Authors also emphasised that, as the altitude decreases, severe and very severe erosivity class values are approximately equally distributed on a multiannual basis. Since this study revealed that the highest average ED corresponds with the obtained RE values, it can be stated that RE is driven by rainfall intensity in the investigated period. Obtained results are in agreement with the findings of Panagos et al. [19]. The authors of this study state that the combination of higher values of RE and ED could lead to various environmental risks when high rainfall amounts falling on moist or saturated soils could trigger landslide events or wetland erosion, thus intensifying soil degradation processes. Within the study area, differences in the seasonal patterns of rainfall erosivity can be seen. With the summer season displaying the highest RE values in western parts of Hungary, the Pannonian part of Croatia, and southeastern parts of the Vojvodina Region in northern Serbia, it can be highlighted that these findings correspond well with the distribution of the seasons with the largest R factor values provided by Panagos et al. [21]. Moreover, further results follow the findings of the mentioned authors who distinguish a second RE zone, with the autumn being the second most erosive period, including parts of Croatia (200-900 MJ·mm·ha −1 ·h −1 ). These observations support the statement that the R factor is highest during the summer season in central Europe due to large-scale extreme precipitation events [32]. Similar patterns can be seen in northern Serbia, where Milovanović et al. [89] indicated that increases in June precipitation amounts reach up to 3 mm/decade (with small areas increasing up to 8 mm/decade, respectively). The findings from this study follow the values presented on monthly erosivity maps of Europe (for January and July) by Panagos et al. [32]. Therefore, as the magnitude of the predicted climate change is likely to have a pronounced influence on water erosion processes, investigation of regional conditions that influence RE and ED is of great importance for future studies. The results of this research clearly outline the importance of the utilisation of quality gridded datasets in research activities related to the estimation of rainfall erosivity on a regional scale. On the other hand, considering the variables of daily precipitation, the E-OBS dataset displays a certain underestimation of the precipitation magnitude in the Pannonian (Carpathian) region. Hence, the investigation of the long-term trends should be performed exclusively on weather station data due to temporal inconsistencies within gridded datasets caused by residual inhomogeneities in the used station network [89].

Long-Term Trends of Annual and Seasonal RE and ED
Long-term trends in the annual and seasonal RE for each meteorological station were examined with the M-K test method and are shown in Figure 11. It can be seen that annual RE has a generally positive trend for every station, except for Szombathely (in the western part of Hungary). However, Szeged in the south of Hungary has a positive trend (p < 0.01). Four meteorological stations with statistically significant trends are located in Vojvodina, with two stations displaying rising tendencies detected for the confidence level of 99% (p < 0.01) and two stations for the confidence level of 95% (p < 0.05). Spring RE values have increasing trends for all stations, of which 5/14 stations have statistically significant results: Szeged and Palić at the confidence level of 99% (p < 0.01), Kikinda, Novi Sad, and Szombathely at the confidence level of 95% (p < 0.05), and Banatski Karlovac at the confidence level of 90% (p < 0.1). Summer and winter season trends of RE (with the exception of Budapest and Szombathely) display an increase, but without statistical significance at the specified confidence levels. Autumn has an increasing trend for RE values in Vojvodina and eastern Croatia, with a statistical significance of 99% (p < 0.01) and 95% (p < 0.05) for three stations, and two stations with a statistical significance of 90% (p < 0.1). For the rest of Hungary, except Szeged, the autumn RE trend is decreasing. These long-term seasonal trends are in line with the results of Tošić et al. [74] for the territory of Vojvodina, as they found an increasing trend for autumn and annual precipitation. Bezak et al. [25] also reported an increasing trend for most parts of the Pannonian Basin by using the M-K test, emphasising statistical significance for central Hungary (p < 0.02), northwestern Vojvodina (from p < 0.06 to p < 0.16), the rest of Vojvodina, and eastern Croatia (p < 0.17-p < 0.28). Figure 12 shows long-term trends in the annual and seasonal ED spatial distribution. With the observation of statistical significance for positive annual trend, it can be seen that Szeged station in southern Hungary has p < 0.01, while three stations in Vojvodina have a positive trend with the statistical significance of 95% (p < 0.1) and one with 90% (p < 0.1). Other stations all have increasing trends, but without statistical significance, except Szombathely, which has a negative trend (without statistical significance). Spring and autumn seasons have similar spatial distribution for positive and negative trend distribution ( Figure 12). During the spring, a rising trend is recorded with a confidence level of 99% for two stations, 95% for three stations, and 90% confidence level for one station. A negative ED trend is recorded during spring in the western part of Hungary (Szombathely station), with a confidence level of 95%. Both summer and winter ED trend distribution are without the presence of statistical significance ( Figure 12). Panagos et al. [19] highlighted a necessity for a joint assessment of the ED and the risk areas in areas with the combination of low precipitation and high ED values (central Hungary), as well as medium-high precipitation with medium-high ED values (western Hungary and eastern Croatia). These results show that parts of the Pannonian Basin are at high risk of erosive events, as well as from water scarcity and floods [14,19].

North Atlantic Oscillation Impact
Large-scale atmospheric circulation patterns have a considerable impact on precipitation variability [90]. As shown by Dehghani et al. [91], these patterns have a considerable impact on precipitation in different regions, where various climate indices in different phases may decrease the seasonal precipitation (even up to 100%). On the other hand, seasonal precipitation may increase >100% in different seasons due to the impact of these indices.
Many recent studies [61,76,77,92] have highlighted that North Atlantic Oscillation (NAO) has a prominent impact on the rainfall amount, and, since RE and ED are based on precipitation amounts, NAO has a certain influence on them as well. As highlighted by Bice et al. [90], NAO has a strong influence on winter precipitation in the Pannonian Basin, with negative NAO phases corresponding to periods of high precipitation. This is consistent with winter RE trend values obtained in this study ( Figure 11).
Other respective authors also argued that the regional intensity and frequency of extreme precipitation increased during the period between 1976 and 2001 [8,93], while total precipitation decreased. The increasing occurrence of the precipitation extremes in this part of the Pannonian Basin characterises the climate of the region as drier in the summer and wetter during the winter season. Other studies indicate that the amount and intensity of precipitation in Serbia have significantly increased during autumn in the Vojvodina (northern Serbia) [94]. It has been shown that the climate of southern parts of the Pannonian Basin has become wetter during the last 50 years, regarding the precipitation magnitude and its frequency [8,94]. The study of Lukić et al. [8] indicated that the number of heavy precipitation days (R10mm) is increasing, with varying degrees of statistical significance at only three stations in the southeastern part of the study area (Banatski Karlovac, Kikinda, and Palić meteorological stations). This extreme precipitation index is in good agreement with the annual spatial distribution of RE and ED values (Figure 13d,e) since the used Zhang et al. [68] model includes the erosive rainfall values when the actual rainfall was >12 mm. Moreover, from Figure 13, it can be seen that this variability is strongly influenced by NAO when observing the southern stations from the investigated area (Figure 13a-c), with negative NAO phases corresponding to periods of high precipitation.
Correlation coefficients between the NAO and R10mm indices were determined in order to investigate the possible relationships between RE/ED and atmospheric variability. Results are shown in Table 3. The correlation between the analysed indices is significant (at a 99% and 95% level of significance) for six weather stations at the annual level. Weather stations such as Kikinda and Palić (northern Serbia) have a significant correlation with values varying between −0.36 and −0.37 at the 99% level of significance. On the other hand, meteorological stations Osijek (eastern Croatia), Szeged, and Debrecen (Hungary), as well as Zrenjanin (northern Serbia), display a significant correlation varying between −0.28 and −0.32 at the 95% level of significance. The correlation between the NAO and R10mm index is significant only for Pécs meteorological station (−0.29 at a 95% level of significance) for the autumn season. The winter season displays strong linearity between the observed indices for Sremska Mitrovica and Zrenjanin (values varying between −0.36 and −0.42 at the 99% level of significance) in northern Serbia. Analysed stations, such as Osijek, Szeged, Debrecen, and Kikinda, have a significant correlation varying between −0.27 and −0.34 at the 95% level of significance. Obtained results correlate well with the findings of Lukic et al. [61], who emphasised that the negative correlation coefficient could indicate the wetting effect on the RE and ED parameters. According to the findings of Lukić et al. [8], the annual total precipitation (when it exceeds the 95th percentile (R95p)) is increasing in the investigated area, except for the Budapest meteorological station. Ossó et al. [95] revealed that from 2000 to 2016, annual values of the R95p index were 5% larger in central Europe. The authors stated that, during the winter season, R95p was positive for most of Europe, thus leading to unusually intense winter precipitation. These results are in good agreement with spatio-temporal distribution of RE and ED obtained in this study, indicating that very wet days significantly contribute to the increase in extreme rainfall events and rainfall erosivity rates in general (Figure 14a,b). Generally speaking, the calculated increasing trends in RE are in goo with this pattern on both annual and seasonal scales. Positive spring and au of RE can be the result of an increase in the amount of precipitation in shor well as due to an increase in the contribution of extreme events to the tota precipitation. Moreover, ED is influenced not only by NAO but also Atlantic/West Russia (EAWR) patterns in terms of intensity of precipita frequency of dry and wet conditions [94]. Future research should be orien the investigation of large-scale atmospheric processes, as drivers of erosivity over the central and southern Pannonian Basin. Such research sh the relationship between the above-mentioned large-scale processes variability, their impact on changes in seasonal precipitation, and thu erosivity values. This approach would enable links to be made betwe outlooks of circulation indices (such as NAO) and heightened/lowered erosion across the region. Similarly, the same atmospheric, rainfall, and ero could be derived from climate model output to assess future risks of so climate change.

Conclusions
Estimation of rainfall erosivity and density is important for the unde the climate vulnerability in the region. Analysis of RE and ED for the southern Pannonian Basin was performed, complemented by seasonal trend variability analyses for the period of 1961-2014. The findings show seasons Generally speaking, the calculated increasing trends in RE are in good agreement with this pattern on both annual and seasonal scales. Positive spring and autumn trends of RE can be the result of an increase in the amount of precipitation in short periods, as well as due to an increase in the contribution of extreme events to the total amount of precipitation. Moreover, ED is influenced not only by NAO but also by the East Atlantic/West Russia (EAWR) patterns in terms of intensity of precipitation and the frequency of dry and wet conditions [94]. Future research should be oriented towards the investigation of largescale atmospheric processes, as drivers of the seasonal erosivity over the central and southern Pannonian Basin. Such research should explore the relationship between the above-mentioned large-scale processes of climate variability, their impact on changes in seasonal precipitation, and thus associated erosivity values. This approach would enable links to be made between seasonal outlooks of circulation indices (such as NAO) and heightened/lowered risk of soil erosion across the region. Similarly, the same atmospheric, rainfall, and erosion indices could be derived from climate model output to assess future risks of soil loss under climate change.

Conclusions
Estimation of rainfall erosivity and density is important for the understanding of the climate vulnerability in the region. Analysis of RE and ED for the central and southern Pannonian Basin was performed, complemented by seasonal trend and spatial variability analyses for the period of 1961-2014. The findings show seasons and hotspots with high erosion risks that are relevant to the development of prevention activities and promotion of suitable mitigation measures at local to regional scales.
Results in this study represent only a piece in the mosaic of spatio-temporal rainfall erosivity variability in Europe. Different periods, meteorological stations, and temporal resolution of data can yield differences in erosion risk indices for this part of the Pannonian Basin. The methodology used in the present study with the usage of station observations and gridded datasets display certain discrepancies when it comes to extremely high erosive events. The annual RE and ED are following the precipitation distribution patterns with a strong spatial variability: higher values are observed for southwestern and southeastern parts; minimal values dominate over north and northeastern parts. Detected trends obtained for weather stations indicate a general increase in RE and ED at the annual level, as well as during the spring and autumn seasons. This increase may cause a transition to a higher erosive class in the near future by elevating the vulnerability for this type of hydro-meteorological hazard in the central and southern Pannonian Basin. This study confirmed the importance of NAO atmospheric circulation by observing its correlation with the number of heavy precipitation days. During the winter season, as well as on an annual scale, significant correlations (for the given confidence levels) are observed, thus indicating the possible wetting effects on the RE and ED parameters.
Despite the use of high-resolution data (E-OBS), extreme erosive events can influence the mean annual RE and ED and, thus, their probability of occurrence needs to be accounted for in future regional studies. However, this study largely enhances previously performed analysis that was focused solely on the assessment of the distribution, concentration, variability of precipitation, and rainfall aggressiveness as the contributing factors in the triggering or reactivation of erosion phenomena in this part of Europe.
The results of this work contribute to growing evidence of erosion risk in southeast Europe and have scope for further refinements. This is especially referring to the territory of northern Serbia (Vojvodina Region) and other Western Balkan countries not intensively covered by previous studies. Future work should be oriented towards the usage of daily rainfall erosivity models with regionalised parameters (in particular α and β, since they are closely related to climate characteristics and the topography of the given location). A greater density of meteorological stations or gridded datasets of high spatial resolution (e.g., ERA5, 0.1 deg.), length of the observation period, and RE model based on RUSLE guidelines could ensure quality regional model parameters. However, one should bear in mind the somewhat existing limitations of the model (limited interactions between input factors, lack of deposition of sediment prediction, model upscaling issues, etc.) and future utilisation of more sensitive indices, along with physically based models for estimating soil erosion rates (e.g., RUSLE, RUSLE 2 ), could help to fill in the gaps in the rainfall erosivity map of Europe. Analysis of the REDES database could enable a more suitable comparison of the results aggregated with the usage of the above-mentioned methodologies. This approach would allow a more dynamic soil erosion assessment and identification of extreme erosive events on both local and regional scales affected by climate change impacts.