Spatial Assessment of the Effects of Land Cover Change on Soil Erosion in Hungary from 1990 to 2018

Abstract

As soil erosion is still a global threat to soil resources, the estimation of soil loss, particularly at a spatiotemporal setting, is still an existing challenge. The primary aim of our study is the assessment of changes in soil erosion potential in Hungary from 1990 to 2018, induced by the changes in land use and land cover based on CORINE Land Cover data. The modeling scheme included the application and cross-valuation of two internationally applied methods, the Universal Soil Loss Equation (USLE) and the Pan-European Soil Erosion Risk Assessment (PESERA) models. Results indicate that the changes in land cover resulted in a general reduction in predicted erosion rates, by up to 0.28 t/ha/year on average. Analysis has also revealed that the combined application of the two models has reduced the occurrence of extreme predictions, thus, increasing the robustness of the method. Random Forest regression analysis has revealed that the differences between the two models are mainly driven by their sensitivity to slope and land cover, followed by soil parameters. The resulting spatial predictions can be readily applied for qualitative spatial analysis. However, the question of extreme predictions still indicates that quantitative use of the output results should only be carried out with sufficient care.

1. Introduction

Soil erosion by water continues to be a significant degradation process for soil resources around the world. Due to the expected increased occurrence of extreme precipitation events, soil erosion remains a key threat to soil resources [1,2,3]. Therefore, soil erosion by water is still in the focus of a number of global and regional policies [4].

In the past decades, a significant number or model-based soil erosion susceptibility mapping studies have been carried out worldwide [5]. Batista et al. [6] provided a review of recently applied soil erosion models, evaluation methods, and limitations.

A large portion of soil erosion modeling projects focused on the use of the empirical Universal Soil Loss Equation (USLE) [7] or its revised versions MUSLE [8,9] and RUSLE [10] in variable spatial and temporal scales. A growing number of applications are utilizing remotely sensed input data [11,12,13,14,15].

The Pan-European Soil Risk Assessment (PESERA) model was primarily developed for the regional scale estimation of rill and inter-rill erosion. As a process-based model, PESERA requires a significantly larger number of input data layers and a larger amount of processing power than the USLE model. Most recently Li et al. [5] have applied the PESERA model along with the RUSLE model. National level soil erosion maps in Hungary have first been founded on field-based expert knowledge [16]. In 2000, a 1:100,000 scale USLE based national erosion risk map has been compiled [17]. In 2016, Pásztor et al. [18] developed a new soil erosion risk map for Hungary, based on the USLE and PESERA models, focusing on the extreme wet year of 2010. This map was later adopted for the National Atlas of Hungary [19]. While the primary focus of soil erosion research in Hungary is erosion by water, there have also been recent studies focusing on susceptibility to wind erosion [18,20].

Since Hungary is still lacking a comprehensive national soil erosion monitoring network [21], proper validation and/or calibration of such maps has been a challenge. A semi-quantitative evaluation study has concluded that the map from the combined USLE–PESERA approach has been in line with in situ observations of farmers [22].

Land use and land cover (LULC) play a critical role in soil erosion processes. The changes in LULC can have significant effect on susceptibility to soil erosion [23,24,25], especially in case of agricultural land use [26,27,28], and, thus, can contribute to reductions in soil organic carbon (SOC) [29]. Therefore, understanding LULC dynamics and its effects on soil erosion can provide key information to decision makers [11,30,31,32,33].

Considering the effects of land use change on soil erosion, several studies have been focused on the assessment of LULC changes. One typical approach is the simulation of LULC change by different models. For example, Shrestha et al. [34] quantified the individual and integrated impacts of climate and land use change in stream flows in the Songkhram River. They used dynamic conversion of land use and its effects (Dyna-CLUE) as a land use change model to define land use change scenarios. Zhang et al. [35] used a combination of a Markov chain model and Dyna-CLUE model to simulate future land uses. Lamichhane and Shakya [36] also utilized a CLUE-S based approach to assess the impact of LULC and climate change on watershed hydrology. Another common methodology is based on using legacy information and/or satellite data for the assessment of LULC changes [37,38,39].

As the abovementioned studies show, both the USLE and the PESERA models are still being applied for evaluating the effects of changes in climate or in land use/land cover, even though there are clear differences in their predictions [40,41]. In their recent study, Ciampalini et al. [42] have performed a sensitivity analysis with the PESERA model, focusing on the effects of climatic parameters. However, there is still a lack of information regarding the effects of other input parameters, specifically the effects of land use/land cover.

In modeling situations, where sufficient validation data is not available, a potential method is the application of random forest variable importance measures in order to assess the sensitivity of models to particular input variables [43,44]. To our knowledge, this method has not yet been applied for the comparison of USLE and PESERA models.

The aim of the work presented here was to assess the extent of land cover changes in Hungary between 1990 and 2018, and to estimate the effects these changes had on soil erosion potential, while also providing a comparison of the two models and their sensitivity to input parameters through the application of random forest variable importance. Our basic assumptions include the established nature of the applied models (i.e., they provide valuable information regarding potential loss due to water erosion [18,22]), and reliability of the input datasets. Our base hypothesis is that LULC changes from 1990 to 2018 have influenced potential erosion rates. In order to make different years comparable, we have used the climate data of 2010 as a benchmark year.

2. Materials and Methods

The current study has adopted the methodology and the data utilized for the development of the Soil Erosion Map of Hungary [18]. While a national erosion monitoring network does not exist, this map (utilizing the 2006 CORINE Land Cover data) has been semi-quantitatively evaluated and is currently accepted at the national level [22,45]. In order to assess the effect of LULC change, all variables unaffected by such changes have been kept as stationary, focusing on the wettest observed year 2010 as climatic baseline.

2.1. Study Area

Hungary is situated in central Europe (45°48′ and 48°35′ N, 16°05′, and 22°58′ E), in the Carpathian Basin. Its highest peak is at 1014 m, while its lowest point is at 78 m. The whole area of the country is in the basin of the Danube River [19].

The relief of Hungary is mostly expressed as low elevation and low vertical dissection (Figure 1). The area is dominated by lowland regions (82.4% below 200 m), with only about 0.5% of terrain above 500 m. Medium-height mountains (200–500 m) cover 2.1% of the area, while hills and foothills 15.5%. [46]

Hungary is located in the northern temperate zone, but is influenced by three climatic zones (oceanic, continental and Mediterranean). Its climate is characterized by four seasons, with a great temporal variability. Typical climatic parameters for the whole of Hungary are presented in

Table 1. Average values of basic climatic parameters of Hungary [41].

Climatic Parameter	Typical Values
Annual sunshine duration	1900–2100 h
Annual mean temperature	10–11 °C
Mean annual temperature range	22.1 °C
Average annual precipitation	580 mm
Prevailing wind direction	northwestern
Annual mean wind speed	2.5 m/s

. Summer is the warmest season with the highest seasonal precipitation, while winter is typically the coldest and driest season. However, precipitation can be especially variable in the region both temporally and spatially. The wettest observed year was 2010, also having an average of 9 days with high (more than 20 mm) precipitation [47]. Recent studies have shown that while the spatial distribution of extreme (more than 40 or 60 mm) rainfall events has not changed significantly over the past decades, their intensity and frequency is showing an increasing tendency [48,49]. Figure 2 presents the spatial distribution of annual precipitation for the present study’s benchmark year 2010.

Amongst the different natural hazards in Hungary, mass movements (primarily landslides) are concentrating only in small areas, while soil erosion by wind can amount to as much as 80–110 million m³ annually, with 10% of the total area of the country being susceptible to wind erosion [45,50,51]. Soil erosion by water affects about 24.7% of the total area of Hungary [18,21].

Hungary has a wide variety of soils. In hilly or mountainous areas, Luvisols and Cambisols are typical, while the Great Plain area has a range of Chernozems, Vertisols, Solonchaks and Solonetz soils. Alluvial plains often contain Luvisols, while certain areas with sand deposits have Arenosols dominating [52]. Figure 3 presents the spatial distribution of the USLE K factor in Hungary.

The vegetation of Hungary is mostly part of the Pannonian region, a forest steppe surrounded by the Turkey oak forest zone. The Carpathian Mountains are dominated by beech and spruce. Influenced by various biogeographic effects, the Pannonian region has a highly diverse vegetation with variable spatial distribution. Besides grasslands and dry oak forests, there is a high proportion of Sub-Mediterranean, continental and Balkan species, with Eurasian elements still dominant in most plant communities. [53]

2.2. Applied Models

2.2.1. The Universal Soil Loss Equation (USLE) Model

The empirical USLE model was the first widespread method to calculate soil loss. The model uses the following formula to calculate average annual soil loss A (t/ha^-):

$$A=R\times K\times L\times S\times C\times P$$

(1)

where R is the rainfall erosion index (MJ mm/ha/h/y), K is the soil erodibility factor (t × ha × h /ha/ MJ/mm), L is slope length factor and S is slope gradient factor, C is the cropping cover management factor, and P is the agricultural practice factor [7,21]. In case of the agricultural practice factor is not available, a default value of 1 can be used [21].

2.2.2. Pan-European Soil Risk Assessment (PESERA) Model

The process-based model partitions the precipitation into overland flow, evapotranspiration and soil moisture storage. It also includes a plant growth model to calculate biomass, leaf fall and vegetation cover. It generates daily rainfall using a monthly Gamma distribution. To compute total erosion, sediment yield Y is calculated by the expression:

$$Y= \varsigma kH\sum r^2$$

(2)

where ς is the ratio of slope base to average gradient, k is the soil erodibility (kg liter^-2 m day), H is the total slope length (m) multiplied by the average slope gradient and r is the local runoff in (mm) for each event. Detailed equations and model description are presented by Kirkby et al. [54].

The PESERA model requires a large number of input layers of which 96 relates to climate (rainfall, temperature, potential evapotranspiration), 25 to land use, crops, and planting dates (land cover type, dominant arable crops, planting dates, initial ground cover, initial surface storage, surface roughness reduction per month, root depth), 6 to soil parameters (crust storage, sensitivity to erosion, effective water storage capacity, soil water available to plants in top 300 mm, soil water available to plants, scale depth) and 1 to topography (standard deviation of elevation in a 1.5 km radius). Some climatic layers are optional and are only used in case of simulating climatic changes (not applicable for the current study). In case dominant crops are not available, a default assumption of maize can be applied to all arable land. Detailed description of model input parameters and application is presented in the PESERA Manual [55].

2.3. Data Sources and Processing

All source and input datasets (see

Table 2. The main datasets used in the analysis.

Input Data	Source	Original Grid Resolution
Land use and land cover	CORINE Land cover [57]	100 × 100 m
Climate data	CARPATCLIM database [3]	0.1° × 0.1°
	Agri4Cast MARS [58]	25 × 25 km
Topographic information	EU-DEM [59]	30 × 30 m
Soil data	DOSoReMI.hu [60]	100 × 100 m

) have been in a raster format, with their original grid size resampled and/or interpolated where necessary. The target grid size of the study was a 100 × 100 m (1 ha) grid.

Data processing and application of the USLE model (multiplication of its factors) was primarily carried out in ArcMap 10.1, while application of the PESERA model utilized ArcInfo Workstation 9.3 for pre- and post-processing. The PESERA model itself is a set of standalone executable and. aml files, including the model itself and data input/output [55]. Analysis, tables and graphs have been prepared with MS Excel and R software [56].

2.3.1. Land Use and Land Cover information

Information on LULC for the years 1990, 2000, 2006, 2012, and 2018 has been obtained from the respective CORINE Land Cover datasets at 100 × 100 m resolution.

The CORINE Land Cover (CLC) inventory includes 44 land cover classes for the abovementioned years. CLC mapping units are delineated via visual interpretation as well as by semi-automatic classification of satellite imagery. The geometric accuracy of the datasets is generally better than 100 m, while the thematic accuracy is typically above 85%, depending on countries [57].

All land use related layers for the PESERA model have been derived by the reclassification of the CORINE data based on the PESERA Manual [55]. As annual information on crops at arable land at a hectare level was not available, maize crop has been selected as a worst-case scenario for both the PESERA and the USLE (C factor) models, with an April sowing date. For the same reason, the USLE P factor has been set to 1, as detailed information was not available at a national scale.

C factor for the USLE model has been calculated based on the CLC values according to Podmaniczky et al. [61]

2.3.2. Climatic Information

The year 2010 has been selected as a baseline for comparing the effects of LULC changes. Therefore, all LULC “scenarios” have used the same climatic input variables. Where available, climate data was obtained from the CARPATCLIM database [3]. However, since a part of western Hungary is not covered by CARPATCLIM, AGRI4CAST MARS data [58] has been used to supplement it. The two datasets have been merged and re-interpolated for a 100 m grid via Ordinary Kriging in SAGA GIS. The PESERA model used monthly data (precipitation, PET, and temperature), while the USLE only required yearly precipitation data, where the extreme nature of the 2010 precipitation was represented as a 20-years return frequency. Calculation of R was based on the method proposed by Renard and Freimund [62], by the combined use of the following equations:

$$R=0.04830 P^{1.610}$$

(3)

and

$$R=587.8-1.219 P+0.004105 P^2$$

(4)

where P is the mean annual precipitation. Equation (3) has been applied where P < 850 mm, while Equation (4) has been applied for P > 850 mm.

2.3.3. Topography

Elevation-based characteristics for both models have been derived from the EU-DEM dataset [59]. The Standard deviation of elevation within 1.5 km radius (for PESERA) was calculated in ArcMap 10.1. The USLE L (slope length) and S (slope gradient) factor was calculated in SAGA GIS based on Moore’s method [63], according to the following equation:

$$LS= \left(0.4+1\right){\left(\frac{A_s}{22.13}\right)}^{0.4}{\left(\frac{\sin \beta }{0.0896}\right)}^{1.3}$$

(5)

where β is the average slope and A_s is the specific catchment area.

2.3.4. Soil Information

Soil parameters (particle size distribution, organic matter content and carbonate content) have been provided from the DOSoReMI.hu initiative [60], based on the Digital Kreybig Soil Information System (DKSIS) [64] and the Hungarian Soil Information and Monitoring System [65]. The USLE K factor has been calculated based on the method of Sharply and Williams [66]:

$$K= \left(0.2+0.3e^{\left[-0.256SAN\left(1-SIL/100\right)\right]}\right)\times {\left(\frac{SIL}{CLA+SIL}\right)}^{0.3}\times \left(1-\frac{0.25OM}{OM+e^{\left(3.72-2.95OM\right)}}\right)\times \left[1-\frac{0.7S{N}_1}{S{N}_1+e^{\left(22.9S{N}_1-5.51\right)}}\right]$$

(6)

where SAN = sand, SIL = silt, CLA = clay and OM = organic matter content of the soil while SN1 = 1-SAN/100.

The PESERA soil erodibility and surface crusting factors were calculated based on Fryrear et al. [67], while soil water available to plants, effective soil water capacity, and scale depth have been downscaled from the original pan-European dataset [55], using regression kriging and conditional generalization [18].

2.4. Harmonization and Combination of Results

The applied harmonization/combination procedure followed the method, which had been used earlier in developing and evaluating the Soil Erosion Map of Hungary [18,22,45]. Some predicted values proved to be negative. These values have resulted from the modeling process. In case of the USLE, they represent urban areas or waterbodies that were later masked out during the procedure and thus set to 0. As the PESERA model also calculates sediment accumulation, its negative values are actually meaningful. However, since sediment accumulation was not in the scope of the current study, in order to avoid any confusing results, predicted values below 0 have been set to 0. This does not affect the results, as soil erosion in deposition sites is considered zero.

The combination of outputs from the two models involved calculating a cell-by-cell mean of the two model outputs. Where only one of the models produced outputs (due to slight differences in methodology), the existing results were accepted.

2.5. Evaluation and Analysis

Statistical analysis of the results primarily included basic descriptive analysis and cross-correlation tables. Year to year comparison of erosion estimates followed the general rule of subtracting the output of the earlier scenario later one. This produced positive values where values have increased and negative ones where they decreased. Comparison of the USLE and PESERA outputs was carried out by calculating the difference between the two output layers for each modeled year.

In order to assess the influence of different input variables, Random Forest regression models have been created for each modeled year, for both models as well as the differences in the predicted results of the two model, from a randomized representative subsample of 70,236 points, using the randomForest package in R. For each model, 500 regression trees have been grown, with 3 variables tried at each node. Since the since the specific input variables were in many ways different between the two models, a set of 10 key variables have been identifies representing the potential effects of precipitation (total annual rainfall in 2010, “prec”), slope (‘slope’), land cover (CORINE Land Cover, “clc”) and soil properties (USLE K factor, “k”; PESERA crust storage, “crust”; PESERA sensitivity to erosion, “erod”; PESERA effective soil water storage capacity, “swsc_eff”; PESERA soil water available to plants in top 300 mm, “p×1”; PESERA soil water available to plants between 300 and 1000 mm, “p×2”; PESERA scale depth, “zm”).

3. Results

3.1. Land Cover Change

The identified changes in LULC are summarized in Figure 4, with original CORINE categories merged into major groups for better presentation. The figure shows an almost continuous reduction in “Non-irrigated arable land”, coupled with an increase of “Artificial surfaces” and “Transitional woodland-shrub”, while other categories seem to have only slight variations.

Figure 5 presents the changes as a percentage of the total area covered by each category. The most observable trends are the reduction of non-irrigated arable land and the increase of artificial surfaces and transitional woodland–shrub areas. These indicate that such changes are primarily driven by land use and natural succession can play a part only when arable land is (either temporarily or long term) not being cultivated. The changes are generally well-distributed in the country, and occur mostly sporadically, without any visible spatial trends at the national level. A summary of the changes from 1990 to 2018 is presented in

Table 3. Area affected by LULC change in Hungary between 1990 and 2018.

LULC 1990–2018	1990 (ha)	2018 (ha)	Changes (ha)	Changes (%)
Urban fabric	415,567	452,283	36,716	0.39
Industrial, commercial, transport	56,848	89,027	32,179	0.35
Mine, dump, and construction sites	11,869	16,489	4620	0.05
Green urban areas, sport-leisure facilities	36,373	40,379	4006	0.04
Non-irrigated arable land	4,958,900	4,708,948	–249,952	–2.69
Rice fields	14,775	8048	–6727	–0.07
Permanent crops	215,013	174,079	–40,934	–0.44
Pastures	680,005	692,421	12,416	0.13
Heterogeneous agricultural areas	483,859	452,706	–31,153	–0.34
Forests	1,684,294	1,744,546	60,252	0.65
Natural grasslands	225,809	230,786	4977	0.05
Transitional woodland-shrub	241,934	425,027	183,093	1.97
Sparsely vegetated areas	2413	2844	431	0.01
Inland wetlands	103,428	87,162	–16,266	–0.18
Inland waters	170,175	176,579	6404	0.07

, while more detailed comparison of all five years is presented in Appendix A.

3.2. Soil Erosion

As output of the erosion modeling and evaluation process, a total of 30 maps have been created (see Appendix B). Estimated soil loss maps have been prepared for the LULC distribution for the years 1990, 2000, 2006, 2012, and 2018. Results included PESERA and USLE based estimates as well as maps of the combined (mean) estimates.

Figure 6 presents the estimated soil loss based on the latest available (2018) land cover information based on the PESERA model. Output of the USLE model for the same year is presented in Figure 7, while the combined (mean) output is shown on Figure 8. Results for other years are provided in Appendix B (Figure A1 and Figure A2).

Spatial coverage of soil erosion classes showed no significant variation between the years. Most of the country area (PESERA: 83.7–83.6%; USLE: 63.6–65.1%) showed below 5 t/ha/y soil loss; while coverage of higher classes has varied between the two models for 5–10 t/ha/y (PESERA: 6.0–6.7%; USLE: 13.3–13.7%); and above 10 t/ha/y (PESERA: 10.3–10.7%; USLE: 21.6–22.7%).

In order to assess the effects of LULC changes and their spatial distribution, a number of maps have been created by subtracting predictions of the older map from the newer one. Changes in erosion rates from the year 1990 to 2018 are displayed in Figure 9, while changes for other years are presented in Appendix B.

Figure 10 presents the changes in mean soil loss over the observed time period for the different estimation (USLE, PESERA, combined) methods. The graphs clearly show that the effects of changes are more prominent in case of the USLE model, and thus the combined approach.

In order to evaluate the effect of land use on the prediction difference between the two models, Figure 11 displays the minimum and maximum values (thus also ranges) of prediction differences within CORINE categories. Some categories (211, 311, and 324) had higher variation than others. These results show that for certain classes (111—continuous urban fabric, 213—rice fields, 221—vineyards, 333—sparsely vegetated areas, 411—inland marshes, 511—water courses) USLE was more likely to provide higher estimates, while for other categories the over-/underestimation of the two methods were more balanced, with the possible exception of pastures (231), where PESERA was slightly more likely to produce higher estimates.

Table 4. Descriptive statistics on the effects of LULC changes on erosion estimates (in t/ha/y).

	1990 to 2000	2000 to 2006	2006 to 2012	2012 to 2018	1990 to 2018
Min.	–3976.6	–936.3	–600.9	–1134.7	–1134.7
Max.	1092.2	929.7	691.0	760.7	760.7
Mean.	0.005	–0.344	–0.085	–0.200	–0.286
Std. Dev.	9.787	8.948	4.768	9.244	10.148

presents the basic descriptive statistics on the effects of LULC changes. The mean shows that at a national level, with the exception of the period 1990 to 2000, the potential erosion has been continuously reduced.

Table 5. Descriptive statistics on the effects of LULC changes and different modeling methods on erosion estimates (in t/ha/y).

Year	Method	Min.1	Max.	Mean.	Std. Dev.
1990	PESERA	0.00	1835.00	4.87	20.70
	USLE	0.00	6306.58	13.86	38.18
	combined	0.00	3989.92	9.36	23.25
2000	PESERA	0.00	1755.00	5.03	20.79
	USLE	0.00	1335.61	13.71	37.64
	combined	0.00	1195.49	9.37	23.15
2006	PESERA	0.00	1835.00	5.06	22.11
	USLE	0.00	1278.95	13.01	35.80
	combined	0.00	1195.49	9.02	22.43
2012	PESERA	0.00	1835.00	5.05	22.52
	USLE	0.00	1278.95	12.85	35.51
	combined	0.00	1195.49	8.94	22.38
2018	PESERA	0.00	1834.65	4.99	21.21
	USLE	0.00	1354.66	12.51	34.69
	combined	0.00	1022.86	8.74	21.64

¹ Negative estimates (sedimentation) have been changed to 0.

compares the changes estimated by the different methods. It clearly shows that the USLE method generally provided a higher mean estimate, while from 2000 onwards, the highest estimate was predicted by the PESERA method.

It has to be noted that the maximum values presented in most of this chapter are extremely high values that are not at all typical under Hungarian conditions. In fact, in all likelihood these extreme values are outliers. However, at present we have insufficient observed information to provide an appropriate baseline from which calculation and removal of extreme outliers would be possible. Observation of the spatial distribution of estimates above 100 t/ha/y indicate that these mostly occur in the hilly regions. In case of the PESERA, they mostly show up in the Transdanubian Mountains (north of Lake Balaton), while the USLE estimates are generally evenly distributed along the mountains/hills of Hungary.

3.3. Importance of Input Variables

In order to assess the importance of the input variables, affecting both the output of the models as well as the differences in predictions between the two, 15 Random Forest (RF) regression models have been calculated. Detailed results of these models are presented in Appendix C. An example visualization of variable importance for the differences between the two models for 2018 is presented on Figure 12. The left side of the figure indicates the importance of each variable as a percentage of increase in mean squared error (MSE) upon leaving out the variable, while the right-hand side presents node purity based on the Gini index.

The ranking of variables by importance is presented in

Table 6. Variable ranking based on Random Forest regression analysis.

Ranking	1	2	3	4	5	6	7	8	9	10
PESERA	zm	clc	slope	prec	swsc_eff	p2x	erod	k	p1x	crust
USLE	slope	clc	k	p2x	prec	swsc_eff	erod	p1x	crust	zm
U-P	slope	clc	zm	p2x	swsc_eff	prec	erod	p1x	crust	k

(P—PESERA; U—USLE; %IncMSE—percentage of increase in MSE; prec—total annual rainfall in 2010; slope—slope; clc—CORINE Land Cover; k—USLE K factor; crust—PESERA crust storage; erod—PESERA sensitivity to erosion; swsc_eff—PESERA effective soil water storage capacity; p×1. PESERA soil water available to plants in top 300 mm; px2—PESERA soil water available to plants between 300 and 1000 mm; zm—PESERA scale depth).

.

Table 7. Variable importance based on Random Forest regression analysis.

		Mean of Squared Residuals:	% Variance Explained:	%IncMSE
				slope	prec	clc	erod	crust	k	p1x	p2x	zm	swsc_eff
1990	P	364.29	47.26	23.82	26.15	29.94	19.25	18.28	21.62	18.04	22.97	43.16	35.50
	U	452.24	73.80	100.29	34.57	68.14	32.36	24.32	58.62	31.79	45.57	22.30	29.46
	U-P	767.34	64.89	93.42	30.40	64.44	32.19	27.02	11.67	34.02	48.37	45.61	47.89
2000	P	331.22	51.90	27.41	39.42	53.92	27.00	16.63	32.03	26.00	24.48	58.04	32.59
	U	454.33	73.20	104.18	40.30	96.80	31.03	24.50	60.88	29.25	36.72	24.92	35.01
	U-P	695.86	66.60	89.29	47.11	86.93	37.85	21.90	44.07	39.30	45.10	56.85	39.34
2006	P	371.83	52.40	34.33	28.89	36.34	42.13	22.39	15.49	14.73	33.41	44.94	25.49
	U	422.93	71.49	84.17	36.80	95.83	35.71	28.99	58.16	21.62	42.31	26.82	34.28
	U-P	711.38	64.56	79.44	36.58	88.94	33.79	21.23	9.90	28.81	42.06	46.01	44.41
2012	P	389.83	51.52	32.03	34.14	34.74	29.31	19.97	11.86	24.50	30.79	40.60	27.75
	U	420.13	71.42	89.85	42.33	96.81	33.79	24.77	53.31	30.52	41.65	26.75	33.87
	U-P	729.05	64.13	77.07	34.81	75.40	36.78	22.71	10.17	23.07	46.08	43.10	44.76
2018	P	727.61	64.20	78.04	41.95	73.73	35.31	20.55	9.41	27.37	52.46	43.99	48.45
	U	442.73	67.40	88.86	37.17	85.72	32.21	35.01	53.53	31.85	47.75	27.77	38.59
	U-P	740.58	60.68	78.96	30.64	70.42	23.32	20.04	12.00	26.16	49.89	59.21	39.62
Mean		PESERA		39.12	34.11	45.74	30.60	19.56	18.08	22.13	32.82	46.14	33.96
		USLE		93.47	38.24	88.66	33.02	27.52	56.90	29.00	42.80	25.71	34.24
		U minus P		83.64	35.91	77.23	32.79	22.58	17.56	30.27	46.30	50.16	43.20

(P—PESERA; U—USLE; %IncMSE—percentage of increase in mean squared error (MSE); prec—total annual rainfall in 2010; slope—slope; clc—CORINE Land Cover; k—USLE K factor; crust—PESERA crust storage; erod—PESERA sensitivity to erosion; swsc_eff—PESERA effective soil water storage capacity; p×1. PESERA soil water available to plants in top 300 mm; px2—PESERA soil water available to plants between 300 and 1000 mm; zm—PESERA scale depth).

presents the importance (see above) of each input variable, as well as the mean squared residuals and the % of variance explained by the Random Forest regression. It is important to note that the importance of these variables is presented based on the RF regression only, since the models themselves do not share the exact same input variables.

4. Discussion

Analysis of the results have demonstrated that the effects of LULC changes between 1990 and 2018 are clearly present in Hungary, with a general reduction in soil erosion estimates. This indicates that the applied changes have somewhat reduced soil erosion risk in Hungary. Based on the comparison of estimated mean erosion values for the 1990 and the 2018 land cover, a reduction of 0.28 t/ha/year has been calculated. Considering the 93,030 km² total territory of Hungary, this would mean the potential conservation of approximately 2.6 million tons of soil per year, under the same climatic conditions as 2010.

However, the spatial distribution of these changes is not even and it is clear that there are small areas all around the country where erosion risk has increased due to the changes in land cover. While the applied methodology of utilizing the CORINE datasets as input is well established [29], it is best applied for analyzing long-term changes. However, the same methodology could also be utilized to evaluate short-term changes based on Sentinel-2 imagery as demonstrated by Karydas et al. [68].

Spatial pattern of changes (Figure 9) in predicted erosion rates reveal that the effect of LULC changes are present throughout the country. However, visual interpretation of the results also suggests that the effect of these changes is much more significant in hilly and mountainous regions. This effect is clearly visible in (but not limited to) the North Hungarian Mountains and the Transdanubian Mountains. Occurrence as well as intensity of such changes (as variables influencing soil erosion) is less expressed in the low-lying areas, such as the Hungarian Great Plain.

The successful application of the two models supported the validity of their application as found in similar studies recently [69]. However, longer, more demanding processing for the PESERA model indicates the growing demand for higher level processing utilities, such as parallel processing capabilities for physically based models [70].

Comparison of the two models showed that as in previous studies, the USLE model produced generally higher estimates than the PESERA. This is in line with the findings of other studies [40,41,70] where applying the RUSLE and the PESERA model provided similar results. Their findings indicate that PESERA has performed better responding to changes in vegetation rather than to slope. Evaluation of the Random Forest regression analysis has confirmed this. It is also worth noting that PESERA was also more sensitive to soil parameters than the USLE model. The PESERA model is also likely to produce lower predictions at finer spatial resolutions [41], which clearly demonstrated its limitations when applying at smaller scales. Changes in the mean maximum values of the PESERA model have been less influenced by the land cover changes, indicating that the USLE model was more sensitive to such effects. The differences in predictions do not seem to be specific for any particular LULC classes. While some land cover categories indeed produced higher variations in prediction difference, these categories can also be associated with a higher spatial coverage. Table 6 and Table 7 indicate that the differences between the two model are mainly driven by their sensitivity to slope and land cover, followed by soil parameters. The differences in predictions underline the findings of [6] regarding the importance of communicating uncertainty in erosion models. While the existence of new and updated spatially explicit erosion estimates is undeniably beneficial, end users—including decision makers and the general public—should also be more aware of the limitations of these products.

Future research should focus on the appropriate determination method for outlier predictions and their removal from the datasets. A comparison of model estimates with field-based monitoring data would also be invaluable for a proper valuation of the methodologies.

5. Conclusions

During our research, we have successfully applied a combined approach of the PESERA and USLE models for estimating the effects of recent LULC changes in Hungary. Based on our results and discussion above, we can establish the following conclusions:

• Changes in land use/land cover from 1990 to 2018 have reduced potential soil erosion by water in Hungary by up to 0.28 t/ha/year on average.
• LULC changes in Hungary have primarily affected the extent of non-irrigated arable land (decreased), artificial surfaces and transitional woodland-shrub areas (increased).
• The USLE model generally provided higher estimates, with higher sensitivity to slope and to LULC changes.
• The PESERA models has proven to be more sensitive to soil parameters.

As our results have indicated, the identification of outlier (extreme) predicted values are a critical issue and further analysis should be carried out to properly identify these. While our study has provided a combined approach for the use of the two models, the methodology could still be advanced by a potential implementation of spatially differential weighting based on the strengths/weaknesses of the applied models.