Temporal and Elevation Trend Detection of Rainfall Erosivity Density in Greece †

This paper presents certain characteristics of trends in rainfall erosivity density (ED), that have not been so far investigated in depth in the current literature. Raw pluviograph data were acquired from the Greek National Bank of Hydrological and Meteorological Information for 108 stations. Precipitation time series values were cleared from noise and errors, and the ratio of missing values was computed. Erosive rainfalls were identified, their return period was determined using intensity–duration–frequency (IDF) curves and erosivity values were computed. A Monte Carlo method was utilized to assess the impact of missing values ratio to the computation of annual erosivity (R) and ED values. It was found that the R values are underestimated in a linear way, while ED is more robust against the presence of missing precipitation values. Indicatively, the R values are underestimated by 49%, when only 50% of the erosive rainfall events are used, while at the same time the estimation error of ED is 20%. Using predefined quality criteria for coverage and time length, a subset of stations was selected. Their annual ED values, as well as the samples' autocorrelation and partial autocorrelation functions were computed, in order to investigate the presence of stochastic trends. Subsequently, Kendall's Tau was used in order to yield a measure of the monotonic relationship between annual ED values and time. Finally, the hypothesis that ED values are affected by elevation was tested. In conclusion: (a) It is suggested to compute ED for the assessment of erosivity in Greece instead of the direct computation of R; (b) stationarity of ED was found for the majority of the selected stations, in contrast to reported precipitation trends for the same time period; and (c) the hypothesis that ED values are not correlated to elevation could not be rejected.


Introduction
Global warming is expected to increase the intensity of rainfall in Europe [1] and; consequently, to increase the soil erosion rates [2].This potential may have a significant impact, especially on Greece, which is inflicted by the phenomenon of desertification, as a combined result of its biogeoclimatic characteristics and the overexploitation of its natural resources [3], taking into account that the most significant process responsible for soil loss in the country is related to rainfall [4].
The second revised version of the Universal Soil Loss Equation, RUSLE2 [5], introduced the erosivity density (ED), as a measure of rainfall erosivity per unit rainfall to develop erosivity values for the USA, due to the fact that ED requires shorter record lengths (as 10 years lead to acceptable results), allows more missing data than R, and is independent of the elevation.
The problem of precipitation trends in Greece has been dealt with in the literature.In general, annual precipitation presents a downward trend for the period 1955-2001 [6].Concerning the ED values in the country, Panagos et al. [7] used interpolated values of R and also interpolated monthly precipitation, both coming from different datasets, to produce maps of seasonal ED values and plotted the average values per three decades and the nine-year moving average for eight stations.However, surveys [8,9] of the above pluviograph data revealed significant proportions of missing values that affect the calculations of R.
This study aims to assess the impact of missing values ratio to the computation of R and ED values in a numerical way, as RUSLE2 uses a theoretical justification.Additionally, its intention is to test the hypothesis that ED is independent of elevation and investigate its temporal trends in Greece using the latest methodologies developed and presented with RUSLE2, taking into account the presence of missing values in precipitation records.

Materials and Methods
The data utilized in the analysis were taken from the Greek National Bank of Hydrological and Meteorological Information [10] and came from 108 meteorological stations.Due to the presence of missing values, a subset of the stations was used for the analysis using two criteria: (a) The stations must have a common time length of at least 30 years; and (b) during these years, the coverage must be at least 45%.
Initially, and after clearing the data from errors, the product of the kinetic energy of a rainfall and the maximum 30 min intensity, EI30, was computed using the pluviograph records [11,12]: where is the kinetic energy per unit of rainfall (MJ/ha/mm), the rainfall depth (mm) for the time interval of the hyetograph, which has been divided into = 1,2, … , sub-intervals, and 30 is the maximum rainfall intensity for a 30 min duration.On the grounds that the use of fixed-time intervals to measure maximum rainfall amounts can lead to an underestimation of the true value [13][14][15], the Hershfield factor equal to 1.14 was used, as Weiss proposed [13].The quantity was calculated for each using the kinetic energy equation of Brown and Foster [16], as corrected and used in RUSLE2 [5,17]: where is the rainfall intensity (mm/hr).A rainfall event was divided into two parts, if its cumulative depth for duration of 6 h at a certain location was less than 1.27 mm.A rainfall was considered erosive if it had a cumulative value greater than 12.7 mm and these were used in the calculations.All rainfalls with extreme EI30 values and a return period greater than 50 years were deleted using the intensityduration-frequency curves for each station, as they have recently been published [18].After the computation of EI30 values, the annual rainfall erosivity density (MJ/ha/h) per station was calculated: where is the number of storms during year , ( 30 ) the erosivity of storm , and the annual precipitation height.The numerator in Equation ( 3) is the annual rainfall erosivity (MJ.mm/ha/h).
A Monte Carlo procedure was used to assess the effect of missing values on the calculation.In this procedure a subset of the calculated EI30 values was extracted based on the data coverage and the water divisions for the selected stations.For 1000 iterations, a random sample per station and year was extracted to simulate different missing values ratios, and the mean absolute percentage error (MAPE) was computed using the initial and the sampled values of ED and R: where = [1, … , ] is the year, is the computed annual value using all rainfall events per station, and , is the computed value coming from the random sample.The autocorrelation coefficient function and the partial correlation coefficient function were compiled [19,20] to investigate the presence of serial correlation in the annual ED values per station.For every selected station, the hypothesis that ED does not change over time was tested using the Kendall's Tau [21] rank correlation value, and the resulting p-values per station were adjusted using the Benjamini and Hochberg method, in order to control the false discovery rate due to multiple statistical testing [22].Finally, the same method was utilized to test the hypothesis that the ED values per station were not affected by the elevation.The data importing, analysis, and presentation were done using the R language for statistical computing and graphics [23] using the packages: hydroscoper [24], hyetor [25], and ggplot2 [26].

Results and Discussion
Based on the criteria about common time length and coverage, 18 meteorological stations (Figure 1; Table 1) were selected for a common time period of 31 years, from 1966 to 1996.Using their pluviograph records, 29,333 rainfalls were extracted and 7570 of them were erosive.Utilizing intensity-duration-frequency curves, 20 rainfalls were removed as outliers, because their return period was from 50 up to 661 years.These return periods cause extreme annual R and ED values and would, disproportionately, affect the calculations.The computed mean value of R for the selected stations 930.05 MJ.mm/ha/h was close to the average value reported for Greece by Panagos et al. (807 MJ.mm/ha/h).In contrast, the computed mean value of ED, 2.41 MJ/ha/h, was two times the value reported in the same study (1.22 MJ/ha/h).The reason for this difference is that Panagos et al. did not use the same precipitation data for the calculation of R, which came from pluviograph records, and ED, in which erosivity came from pluviograph records, and precipitation had different origins: one-km-global-spatial-resolution monthly values [27].The Monte Carlo procedure results showed that ED was more robust against the presence of missing precipitation values.Using only 5% of the data, annual R values were underestimated on average by 85% when the average estimation error of ED value was 50%.R was proportionally underestimated as the missing ratio increased, while ED's estimation error followed a parabolic curve.In the presence of 50% of the data, R values were underestimated by 49%, while at the same time the estimation error of ED was 20%.
The findings regarding the samples' autocorrelation coefficient functions and the partial correlation coefficient functions did not reveal any practical meaning of the statistically significant values that were found at specific lags in the time-series of a small number of stations.On account of the previous fact, it is safe to suppose that no stochastic trends existed for the examined time-series.The Kendall's Tau rank correlation test results indicated that, for the majority of the stations, the null hypothesis that annual ED values change over time could not be rejected for a significance level α = 5% (Table 1).Thus, it is reasonable to suppose that these time-series are realizations of stationary processes.Concerning the relation between elevation and ED, the p-value = 0.053, using the same test, indicates that the null hypothesis that annual ED values is affected by the elevation could not be rejected for the same significance level α = 5%.Table 1.Location and analysis results for the stations with a common time length during 1966-1996.ID is an abbreviation for the station ID, as reported in the Greek National Bank of Hydrological and Meteorological Information, WD for the Greek Water Divisions, Lon for longitude, Lat for latitude, El for elevation, MCV for mean coverage per station, padj is the adjusted p-value from the test using the Benjamini and Hochberg method.With a star are marked the test results where the null hypothesis is rejected for a significance level α = 5%.

Figure 1 .
Figure 1.Locations of the stations.The selected stations with a common time length for the time period 1966-1996 are symbolized with red and the corresponding number.The stations not used in the trends analysis are symbolized with gray.