Effect of Climate Change on Soil Erosion in a Mountainous Mediterranean Catchment (Central Pindus, Greece)

The aim of this study was to assess soil erosion changes in the mountainous catchment of the Portaikos torrent (Central Greece) under climate change. To this end, precipitation and temperature data were derived from a high-resolution (25 × 25 km) RegCM3 regional climate model for the baseline period 1974–2000 and future period 2074–2100. Additionally, three GIS layers were generated regarding land cover, geology, and slopes in the study area, whereas erosion state was recognized after field observations. Subsequently, the erosion potential model (EPM) was applied to quantify the effects of precipitation and temperature changes on soil erosion. The results showed a decrease (−21.2%) in annual precipitation (mm) and increase (+3.6 ◦C) in mean annual temperature until the end of the 21st century, and the above changes are likely to lead to a small decrease (−4.9%) in soil erosion potential.


Introduction
Climatologists highlight that climate change is occurring, both in terms of air temperature and precipitation patterns.According to the Intergovernmental Panel on Climate Change's (IPCC) Fifth Assessment Report, the Mediterranean basin is expected to become warmer and dryer due to an anthropogenic increase of greenhouse gas emissions (CO 2 , CH 4 , N 2 O, and F-gases) until the end of the 21st century [1].Moreover, in Mediterranean regions, future warming will probably be larger than the global mean, accompanied by a considerable decrease of total rainfall amount [2,3] and more frequent high-intensity rainfall events [4,5].Future climate projections will highly influence catchments' responses to soil erosion [6][7][8][9].
Soil erosion by water is one of the most significant forms of land degradation, as it threatens natural ecosystems, water resources, and crop productivity.The European Commission's Soil Thematic Strategy has identified soil erosion as a key priority for the protection of soils [10] and call for quantitative assessments of soil loss rates at the European level [11].Mediterranean regions are particularly vulnerable to erosion because of the highly irregular behavior of the rainfall regime, both on spatial and temporal scales [12]; inappropriate agricultural management practices [13]; overgrazing [14,15]; and wildfires [16,17].Also, erosion rates are higher in mountainous areas than in lowlands due to the steeper relief and the higher and more intense rainfalls.
Identification of critical soil erosion-prone areas is necessary for the implementation of appropriate mitigation measures to combat erosion [18][19][20].Stakeholders should also quantify potential soil erosion risk under future climate conditions for the scheduling of infrastructure projects.
Actual measurements of soil erosion rates is a costly and time-consuming process that is usually limited to small experimental sites [21][22][23][24][25].To this end, plenty of erosion prediction models were Water 2018, 10, 1469 2 of 12 developed during the last few decades and classified on different spatial/temporal scales and various levels of complexity [26].Erosion prediction models use mathematical equations so as to express the relationships between natural factors (e.g., rainfall, land cover, vegetation type, geological subsoil, and topography) and the erosion process.The most well-known models that have been widely used are USLE [27], RUSLE [28], RMMF [29], WEPP [30], and EPM (the Gavrilovič method) [31].
The quantitative erosion prediction model of Gavrilovič was developed using field work and laboratory experiments to determine the range of values for its parameters.It was widely used in former Yugoslavia since the early 70s, and later in the Balkan region and Central Europe.It has also been recognized as the most quantitative of all the semi-quantitative models [32].Several studies in the mountainous catchments of Greek territory which compare the Gavrilovič prediction model results with actual measurements shows that it gives satisfactory results [33][34][35][36].
In recent years, climate models' projections have been used to assess soil erosion under climate change [37,38].Understanding the uncertainties of several climatological parameters is of major interest in recent climate studies [39].It should be noted that the ability of models to represent climate conditions should be examined prior to their use in impact assessment studies [40][41][42][43].
Investigation of the effects of climate change on soil erosion using high-resolution Regional Climate Models (RCMs) are limited in mountainous catchments of the Mediterranean region [44], while no such research has been carried out in Greece.
The main object of the current research is to quantify the effect of climate change on soil erosion in a mountainous catchment of Central Greece using the erosion prediction model of Gavrilovič and climate simulation of the RegCM3regional climate model.

Study Area
The study was conducted in the mountainous catchment of the Portaikos torrent (Figure 1).It is located in the Thessaly Regional Unit (Central Greece) over the mountain range of Pindus, and is a tributary of the Pinios River.It covers an area of 136.4 km 2 and the relief is rather intense.The mean elevation is 963 meters above sea level (MASL) (maximum 1862 MASL and minimum 240 MASL), whereas mean catchment slope is 52.9% and main stream slope is 7.8%.
The area is considered highly important from a hydrological point of view as it is one of the rainiest areas of Greece territory [45], whereas extensive torrential phenomena (e.g., landslides, erosion) appeared due to the vulnerable geological subsoil (e.g., limestone and flysch) [46].

Climate Simulation
Currently, general circulation models (GCMs) are being widely used to estimate future projections in order to satisfy the growing interest in the impact of climate change on a global scale.GCMs are numerical models that simulate physical processes in the atmosphere, the ocean, the cryosphere, and the surface of the soil based on the laws of conservation of momentum, mass, total energy, and quantity of water vapor.The spatial analysis of GCMs gives satisfactory results for global circulation across the planet, as well as the general characteristics of various climatic parameters on a concise scale.However, it is not possible to precisely simulate phenomena related to the effect of topography on a local and regional scale due to local conditions and particularities, such as complex topography, coastlines, lakes, and small islands [48].
To this end, a number of statistical [4,5] or dynamic [49,50] downscaling techniques have been developed to bridge the gap between the large-scale GCM information and local scales.The statistical methods use the observed relationships between large-scale circulation and the local climate, whereas dynamic techniques use physically based regional climate models (RCMs).
In the frame of the EU funding project, ENSEMBLE (http://ensembles-eu.metoffice.com), a set of multi-model RCM simulations for characterizing climate change in Europe with high spatial resolution (25 × 25 km), was produced.A recent study by Stefanidis [51] evaluated the ability of seven RCMs of the ENSEMBLE project under the A1B SRES emission scenario to represent temperature and precipitation conditions over the mountainous Central Pindus (the study area) for the baseline period 1974-2000.The results concluded that the best simulations were made by the International Centre for the Theoretical Physics Regional Climate Model (RegCM3) [52].Therefore, in this study, regional climate analysis was performed using daily temperature and precipitation-simulated data derived from the RegCM3 model for the baseline period  and future period (2074-2100).
Additionally, mean monthly precipitation and temperature simulations and observations were compared over the period of 1974-2000 (Figure 2).The observational data were derived by the nearest meteorological station in the study area (Pertouli) with available temperature and precipitation data.

Climate Simulation
Currently, general circulation models (GCMs) are being widely used to estimate future projections in order to satisfy the growing interest in the impact of climate change on a global scale.GCMs are numerical models that simulate physical processes in the atmosphere, the ocean, the cryosphere, and the surface of the soil based on the laws of conservation of momentum, mass, total energy, and quantity of water vapor.The spatial analysis of GCMs gives satisfactory results for global circulation across the planet, as well as the general characteristics of various climatic parameters on a concise scale.However, it is not possible to precisely simulate phenomena related to the effect of topography on a local and regional scale due to local conditions and particularities, such as complex topography, coastlines, lakes, and small islands [48].
To this end, a number of statistical [4,5] or dynamic [49,50] downscaling techniques have been developed to bridge the gap between the large-scale GCM information and local scales.The statistical methods use the observed relationships between large-scale circulation and the local climate, whereas dynamic techniques use physically based regional climate models (RCMs).
In the frame of the EU funding project, ENSEMBLE (http://ensembles-eu.metoffice.com), a set of multi-model RCM simulations for characterizing climate change in Europe with high spatial resolution (25 × 25 km), was produced.A recent study by Stefanidis [51] evaluated the ability of seven RCMs of the ENSEMBLE project under the A1B SRES emission scenario to represent temperature and precipitation conditions over the mountainous Central Pindus (the study area) for the baseline period 1974-2000.The results concluded that the best simulations were made by the International Centre for the Theoretical Physics Regional Climate Model (RegCM3) [52].Therefore, in this study, regional climate analysis was performed using daily temperature and precipitation-simulated data derived from the RegCM3 model for the baseline period  and future period (2074-2100).
Additionally, mean monthly precipitation and temperature simulations and observations were compared over the period of 1974-2000 (Figure 2).The observational data were derived by the nearest meteorological station in the study area (Pertouli) with available temperature and precipitation data.

Methodology
In order to evaluate the performance of the RegCM3 model, observed values of precipitation and temperature from the meteorological station were compared with simulated values of the model over the baseline period 1974-2000.The Root Mean Square Error (RMSE) and Mean Bias Error (MBE) were used as evaluation indexes.The combination of these criteria has been widely used for the evaluation of RCMs [4,53,54].The mathematical description of these indexes is given below [55,56]: ) where n is the number of observations, and xi and yi are the observed and simulated values, respectively, for i = 1, 2, …, n.The RMSE gives the weighted variations (residuals) in errors between the simulated and observed values, while the MBE measures the weighted average magnitude of the errors.MBE is the most natural and unambiguous measure of average error magnitude [57,58].RMSE, on the other hand, is one of the most widely used error measures.It is mentioned that the smaller the values of RMSE and MBE, the better the performance of the model.The Erosion Potential Model (EPM) [31] was used in order to estimate annual soil loss.The mathematical description of the model is given in the following equation: where W is the mean annual soil loss (m 3 /year), H is the mean annual precipitation (mm), π is a mathematical constant approximately equal to 3.14159, F is the catchment area (km 2 ), and T is a coefficient of temperature, given by the equation below: where t0 is the average annual temperature (°C) and Z is the coefficient of erosion, given by the next equation: where x quantifies the protective nature of the land cover, y is a coefficient expressing soil resistance to erosion, φ quantifies the observed erosion process, and J is the mean slope of the catchment (%).
The values' range of the above-mentioned coefficients (x, y, and φ) are given in Table 1.

Methodology
In order to evaluate the performance of the RegCM3 model, observed values of precipitation and temperature from the meteorological station were compared with simulated values of the model over the baseline period 1974-2000.The Root Mean Square Error (RMSE) and Mean Bias Error (MBE) were used as evaluation indexes.The combination of these criteria has been widely used for the evaluation of RCMs [4,53,54].The mathematical description of these indexes is given below [55,56]: where n is the number of observations, and x i and y i are the observed and simulated values, respectively, for i = 1, 2, . . ., n.The RMSE gives the weighted variations (residuals) in errors between the simulated and observed values, while the MBE measures the weighted average magnitude of the errors.MBE is the most natural and unambiguous measure of average error magnitude [57,58].RMSE, on the other hand, is one of the most widely used error measures.It is mentioned that the smaller the values of RMSE and MBE, the better the performance of the model.The Erosion Potential Model (EPM) [31] was used in order to estimate annual soil loss.The mathematical description of the model is given in the following equation: where W is the mean annual soil loss (m 3 /year), H is the mean annual precipitation (mm), π is a mathematical constant approximately equal to 3.14159, F is the catchment area (km 2 ), and T is a coefficient of temperature, given by the equation below: where t 0 is the average annual temperature ( • C) and Z is the coefficient of erosion, given by the next equation: Water 2018, 10, 1469 where x quantifies the protective nature of the land cover, y is a coefficient expressing soil resistance to erosion, ϕ quantifies the observed erosion process, and J is the mean slope of the catchment (%).The values' range of the above-mentioned coefficients (x, y, and ϕ) are given in Table 1.

Coefficient of type and extent of erosion ϕ
Little erosion on catchment 0.1-0.2Erosion in waterways on 20 to 50% of the catchment area 0.3-0.5 Erosion in rivers, gullies and alluvial deposits, karstic erosion 0.6-0.7 50 to 80% of catchment area affected by surface erosion and landslides 0.8-0.9Whole catchment affected by erosion 1.0 Aerial orthophotos maps, with a scale of 1:5000, were used in order to evaluate the land cover coefficient (x).The coefficient of soil erodibility (y) was determined from geological maps with a scale of 1:50,000 of the Institute of Geology and Mineral Exploration of Greece (IGME), and the coefficients of type and extent of erosion (ϕ) were evaluated from field observations and recognition of sediment sourced areas and the state of erosion.Subsequently, the mean annual precipitation (H) and the temperature coefficient (T) were derived from the data of the RegCM3 model of ENSEMBLE.Finally, the mean catchment slope (J) and catchment area (F) were determined from a digital elevation model (DEM) with a cell size of 5 × 5 m.

Results
Firstly, the ability of the RegCM3 model to simulate the climate condition (precipitation and temperature) in the study area was evaluated using criteria of the RMSE and MBE.Results from the calculation of the RMSE and MBE between the simulated and observed precipitations (mm) during the period of 1974-2000 are given in the figure below (Figure 3).The model generally underestimated monthly precipitation in most cases except for January and March, where a slight overestimation can be observed.Also, it can be seen that better simulation was achieved in spring and summer.
Additionally, the above-mentioned evaluation criteria were calculated regarding temperature ( • C) data (Figure 4).It was found that the model overestimated mean monthly temperatures for almost all months, except for September and October.Moreover, as revealed from the analysis, better temperature simulation was achieved for winter and autumn.
Firstly, the ability of the RegCM3 model to simulate the climate condition (precipitation and temperature) in the study area was evaluated using criteria of the RMSE and MBE.Results from the calculation of the RMSE and MBE between the simulated and observed precipitations (mm) during the period of 1974-2000 are given in the figure below (Figure 3).The model generally underestimated monthly precipitation in most cases except for January and March, where a slight overestimation can be observed.Also, it can be seen that better simulation was achieved in spring and summer.Additionally, the above-mentioned evaluation criteria were calculated regarding temperature (°C) data (Figure 4).It was found that the model overestimated mean monthly temperatures for almost all months, except for September and October.Moreover, as revealed from the analysis, better temperature simulation was achieved for winter and autumn.Analysis of the land cover showed that dense forest cover was approximately 35.2% of the total catchment area, and thin forests were 24.7%, agricultural crops 18%, pastures 8.5%, dense shrubs 4.9%, scarce shrubs 3%, barren land 4.6%, and settlements 1.1%.In order to determine the coefficient of land cover (x), based on Table 1, the value 0.125 were was to the dense forest, 0.20 to thin forests, 0.25 to dense shrubs, 0.60 to scarce shrubs, 0.7 to pastures, 0.8 to agricultural crops, 0.9 to barren land, and 0 to settlements.
Regarding the geology of the study area, it was found that it mainly consisted of limestone (47.2%) and flysch (44.5%), followed by neogene rocks (6.5%) and alluvial deposits (1.8%).The following values of soil erodibility (y coefficient) were assigned to each petrographic formation: 0.8 to limestone, 1.15 to flysch, 1.55 to neogene rock, and 1.0 to alluvial deposits.The coefficients of land cover (x) and soil erodibility (y) were selected based on literature [33][34][35][36].
Also, following the field observations, the value 0.6 was given to theφcoefficient.The catchment area (F coefficient) was found to be equal to 136.4 km 2 after processing the DEM of the study area through GIS techniques [9,46], and the slopes were found to be rather intense, ranging from 1.5% to 89.3%.The land cover, geology, and slope (%) maps are shown in the next figure (Figure 5).Analysis of the land cover showed that dense forest cover was approximately 35.2% of the total catchment area, and thin forests were 24.7%, agricultural crops 18%, pastures 8.5%, dense shrubs 4.9%, scarce shrubs 3%, barren land 4.6%, and settlements 1.1%.In order to determine the coefficient of land cover (x), based on Table 1, the value 0.125 were was to the dense forest, 0.20 to thin forests, 0.25 to dense shrubs, 0.60 to scarce shrubs, 0.7 to pastures, 0.8 to agricultural crops, 0.9 to barren land, and 0 to settlements.
Regarding the geology of the study area, it was found that it mainly consisted of limestone (47.2%) and flysch (44.5%), followed by neogene rocks (6.5%) and alluvial deposits (1.8%).The following values of soil erodibility (y coefficient) were assigned to each petrographic formation: 0.8 to limestone, 1.15 to flysch, 1.55 to neogene rock, and 1.0 to alluvial deposits.The coefficients of land cover (x) and soil erodibility (y) were selected based on literature [33][34][35][36].
Also, following the field observations, the value 0.6 was given to the ϕ coefficient.The catchment area (F coefficient) was found to be equal to 136.4 km 2 after processing the DEM of the study area through GIS techniques [9,46], and the slopes were found to be rather intense, ranging from 1.5% to 89.3%.The land cover, geology, and slope (%) maps are shown in the next figure (Figure 5).
Also, following the field observations, the value 0.6 was given to theφcoefficient.The catchment area (F coefficient) was found to be equal to 136.4 km 2 after processing the DEM of the study area through GIS techniques [9,46], and the slopes were found to be rather intense, ranging from 1.5% to 89.3%.The land cover, geology, and slope (%) maps are shown in the next figure (Figure 5).The climate was evaluated from the precipitation (mm) and temperature ( • C) time-series data of the RegCM3 for the periods 1974-2000 and 2074-2100.The analysis highlighted a decrease (−21.2%) in annual precipitation (mm) and a significant increase (+3.6 • C) in mean annual temperature until the end of the 21st century (Figure 6).Specifically, the annual precipitation was expected to decrease from 1071.2 mm (1974-2000) to 874.1 mm (2074-2100) and the mean annual temperature was expected to increase from 9.3 • C to 13.0 • C. The biggest decrease in precipitation level will be in the months of spring and autumn, whereas the temperature increase will be greater in summer and winter.The climate was evaluated from the precipitation (mm) and temperature (°C) time-series data of the RegCM3 model for the periods 1974-2000 and 2074-2100.The analysis highlighted a decrease (−21.2%) in annual precipitation (mm) and a significant increase (+3.6 °C) in mean annual temperature until the end of the 21st century (Figure 6).Specifically, the annual precipitation was expected to decrease from 1071.2 mm (1974-2000) to 874.1 mm (2074-2100) and the mean annual temperature was expected to increase from 9.3 °C to 13.0 °C.The biggest decrease in precipitation level will be in the months of spring and autumn, whereas the temperature increase will be greater in summer and winter.Using the above-mentioned parameter and EPM formula (Gavrilovič method), the soil loss (m 3 /year) and catchment erosion rate (m 3 /year/km 2 ) were estimated for the baseline  and future (2074-2100) period.Finally, the results were compared between them so as to quantify the effects of climate change on soil erosion.It was noticed that soil loss would decrease by7920.5m 3 /year (−4.9%) and, consequently, the erosion rate by 58.1 m 3 /year/km 2 .Detailed results can be seen in Table 2.

Discussion
According to studies which compared the results of the EPM model with actual measurements, it was found that the model overestimated soil loss by 10% [32][33][34][35][36].The erosion rate in the study area was considered rather high in comparison with other areas of the Greek territory [59].Despite the high forest cover, the vulnerable geological subsoil and the steep slopes were found to favor the Using the above-mentioned parameter and EPM formula (Gavrilovič method), the soil loss (m 3 /year) and catchment erosion rate (m 3 /year/km 2 ) were estimated for the baseline  and future (2074-2100) period.Finally, the results were compared between them so as to quantify the effects of climate change on soil erosion.It was noticed that soil loss would decrease by7920.5m 3 /year (−4.9%) and, consequently, the erosion rate by 58.1 m 3 /year/km 2 .Detailed results can be seen in Table 2.

Discussion
According to studies which compared the results of the EPM model with actual measurements, it was found that the model overestimated soil loss by 10% [32][33][34][35][36].The erosion rate in the study area was considered rather high in comparison with other areas of the Greek territory [59].Despite the high forest cover, the vulnerable geological subsoil and the steep slopes were found to favor the development of erosion phenomena.
Nowadays, there is an urgent need for reliable future climate projections for rational catchment management and infrastructure project scheduling.The most modern tool used to simulate future climate conditions is the regional climate model (RCM).However, using models can introduce uncertainty in relation to a number of significant temporal and spatial scaling issues of the input data.Generally, RCMs underestimate monthly precipitation and overestimate temperature, whereas better simulation can be achieved in spring and summer for precipitation and winter and autumn for temperature data [60,61].
Concerning future climate conditions in the Mediterranean region [60,61], and especially Greece [3,48,[62][63][64][65], the climate is expected to be warmer and dryer.Additionally, the temperature in Central Europe is projected to increase between +2 • C and +5 • C.Even though annual precipitation is projected to increase up to +10%, most RCMs project a significant decrease of precipitation in summer [66].Moreover, extreme precipitations will be more frequent [5], as well as torrential flooding and debris-flow phenomena.

Conclusions
In this study, the effect of climate change on soil erosion in a mountainous catchment of Central Greece was assessed using the erosion prediction model of Gavrilovič and climate simulation of the RegCM3regional climate model.The main conclusions drawn are the following:

•
Annual soil loss was found to be equal to 161,236.5 m 3 /year, and the erosion rate 1182.1 m 3 /year/km 2 , by applying the Gavrilovič erosion prediction model (EPM) in the mountainous catchment of the Portaikos torrent.

•
Concerning future climate conditions in the study area, based on RegCM3 data, a decrease (−21.2%) in annual precipitation (mm) and increase (+3.6 • C) in mean annual temperature is expected until the end of the 21st century.The above-mentioned model emphasized a greater decrease in monthly rainfall, especially in spring.As for the temperature, a greater increase in mean monthly temperature should occur in summer.

•
Finally, considering the reported changes in climatic condition in the research area, a slight decrease (−4.9%) in erosion rate was estimated.It is suggested that similar surveys should be done all over the Greek mountainous regions so as to identify future erosion-prone areas.Moreover, investigations of the effect of climate change on other hydrometeorological hazards (e.g., flash floods, landslides, etc.) is necessary for adaption and mitigation measures.The EURO-CORDEX project's comparison of the results of the current study with recent very high-resolution (10 × 10 km) climate change projections [67] could be the target of future research in order to highlight the importance of using the finer-resolution model over complex terrain areas, such as Greece.Simulations of finer timescale data (e.g., daily) could also be examined.

Figure 1 .
Figure 1.A map of the study area.

Figure 1 .
Figure 1.A map of the study area.

Figure 2 .
Figure 2. Comparison of (a) mean monthly precipitation (mm) and (b) mean monthly temperature ( • C) between observed (meteorological station) and simulated (RegCM3) data over the period of 1974-2000.

Figure 3 .
Figure 3. Root mean square error (RMSE) and mean bias error (MBE) values between simulated (RegCM3) and observed (meteorological station) precipitation data for the baseline period 1974-2000.

Figure 3 .
Figure 3. Root mean square error (RMSE) and mean bias error (MBE) values between simulated (RegCM3) and observed (meteorological station) precipitation data for the baseline period 1974-2000.

Figure 4 .
Figure 4. RMSE and MBE values between simulated (RegCM3) and observed (meteorological station) temperature data for the baseline period 1974-2000.

Figure 4 .
Figure 4. RMSE and MBE values between simulated (RegCM3) and observed (meteorological station) temperature data for the baseline period 1974-2000.
Water 2018, 10, x FOR PEER REVIEW 7 of 11

Figure 6 .
Figure 6.Comparison of (a) precipitation (mm) and (b) temperature (°C) regimes in the study area, as revealed by the Theoretical Physics Regional Climate Model (RegCM3).

Figure 6 .
Figure 6.Comparison of (a) precipitation (mm) and (b) temperature ( • C) regimes in the study area, as revealed by the Theoretical Physics Regional Climate Model (RegCM3).

Table 1 .
Descriptive variables used in the Erosion Potential Model (EPM).

Table 2 .
Annual soil loss and erosion rate under present and future climate conditions.

Table 2 .
Annual soil loss and erosion rate under present and future climate conditions.