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

Abstract

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.

1. 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 bio- geoclimatic 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.

2. 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, EI₃₀, was computed using the pluviograph records [11,12]:

$$E{I}_{30}=\left(\sum _{r=1}^m{e}_r\cdot v_r\right)\cdot I_{30}$$

(1)

where e_r is the kinetic energy per unit of rainfall (MJ/ha/mm), v_r the rainfall depth (mm) for the time interval r of the hyetograph, which has been divided into r = 1,2, … , m sub-intervals, and I₃₀ 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 e_r was calculated for each r using the kinetic energy equation of Brown and Foster [16], as corrected and used in RUSLE2 [5,17]:

$$e_r=0.29\cdot \left(1-0.72e^{-0.82i_r}\right)$$

(2)

where i_r 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 EI₃₀ values and a return period greater than 50 years were deleted using the intensity–duration–frequency curves for each station, as they have recently been published [18]. After the computation of EI₃₀ values, the annual rainfall erosivity density ED_j (MJ/ha/h) per station was calculated:

$$E{D}_j=\frac{{\sum }_{k=1}^{m_j}{\left(E{I}_{30}\right)}_k}{P_j}$$

(3)

where m_j is the number of storms during year j, (EI₃₀)_k the erosivity of storm k, and P_j the annual precipitation height. The numerator in Equation (3) is the annual rainfall erosivity R_j (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 EI₃₀ 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:

$$MAPE=\frac{1}{n}\sum _{t=1}^n\left|\frac{Y_t-Y_{t,miss}}{Y_t}\right|$$

(4)

where t = [1, … , n] is the year, Y_t is the computed annual value using all rainfall events per station, and Y_t,miss 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].

3. Results and Discussion

Based on the criteria about common time length and coverage, 18 meteorological stations (Figure 1;

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%.

	ID	Name	WD	Lon (°)	Lat (°)	El (m)	MCV (%)	Tau	padj
1	200003	GRABIA	GR07	22.43	38.67	381	73.4	0.12	0.612
2	200011	LIDORIKI	GR04	22.20	38.53	548	69.2	−0.09	0.612
3	200015	PYRA	GR04	22.27	38.74	1137	74.8	−0.11	0.612
4	200018	AG. TRIADA	GR07	22.92	38.35	400	65.4	0.31	0.081
5	200021	DISTOMO	GR07	22.67	38.43	458	60.3	−0.02	0.919
6	200024	LEIBADIA	GR07	22.87	38.44	176	56	−0.27	0.132
7	200059	BASILIKO	GR05	20.59	40.01	747	75.8	−0.11	0.612
8	200092	ELASSONA	GR08	22.19	39.89	276	71.7	0.02	0.919
9	200135	KALYBIA	GR02	22.30	37.92	822	65.3	0.29	0.123
10	200142	NEMEA	GR02	22.66	37.83	306	63.8	−0.26	0.132
11	200144	SPATHOBOUNI	GR02	22.80	37.85	150	48.1	−0.08	0.612
12	200181	LESINIO	GR04	21.19	38.42	2	59.9	0.45	0.055
13	200190	POROS REG.	GR04	21.75	38.51	182	67.8	−0.11	0.612
14	200243	NEOCHORIO	GR03	22.48	37.67	704	63.2	0.14	0.595
15	200291	A. ARCHANES	GR13	25.16	35.24	392	51.6	0.09	0.612
16	200309	DRAMA	GR11	24.15	41.14	100	69.6	0.10	0.612
17	200311	PARANESTE	GR12	24.50	41.27	122	66.1	−0.46	0.005*
18	200346	KATERINE	GR09	22.51	40.28	30	64.2	−0.15	0.595

) 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%.

4. Conclusions

In summary, the main conclusions of our study are:

• It is suggested to compute ED for the assessment of erosivity in Greece instead of the direct computation of R, due to the large proportion of missing values in the pluviograph records.
• Stationarity of ED was found for the majority of the selected stations, in contrast to reported precipitation trends for the same time period.
• The hypothesis that ED values are not correlated to elevation could not be rejected.