Soil Erosion Assessment and Prediction in Urban Landscapes: A New G2 Model Approach

Soil erosion is a global problem that negatively affects the quality of the environment, the availability of natural resources, as well as the safety of inhabitants. Soil erosion threatens the functioning of urban areas, which was the reason for choosing the territory of the Master Plan of Belgrade (Serbia) as the research area. The calculation of soil erosion loss was analyzed using the G2 erosion model. The model belongs to a group of empirical models and is based on the synthesis of the equation from the Revised Universal Soil Loss Equation (RUSLE) and the Erosion Potential Method (EPM). The estimation of soil degradation was analyzed in two time periods (2001 and 2019), which represent the time boundaries of the management of the Master Plan of Belgrade. The novel approach used in this research is based on using the land cover inventory as a dynamic indicator of the urbanization process. Land cover was identified using remote sensing, machine learning techniques, and the random forest algorithm applied to multispectral satellite images of the Landsat mission in combination with spectral indices. Climatic parameters were analyzed on the basis of data from meteorological stations (first scenario, i.e., 2001), as well as on simulations of changes based on climate scenario RCP8.5 (representative concentration pathways) concerning the current condition of the land cover (second scenario). A comparative analysis of the two time periods identified a slight reduction in total soil loss. For the first period, the average soil loss value is 4.11 t·ha−1·y−1. The analysis of the second period revealed an average value of 3.63 t·ha−1·y−1. However, the increase in non-porous surfaces has led to a change in the focus of soil degradation. Increased average soil loss as one of the catalysts of torrential flood frequencies registered on natural and semi-natural areas were 43.29% and 16.14%, respectively. These results are a significant contribution to the study of soil erosion in urban conditions under the impact of climate change.


Introduction
Erosion is one of the most important contemporary problems of soil degradation and a significant factor endangering the environment all around the world [1]. Soil erosion and the consequent production of erosion material are affected by time-static, conditionally invariant parameters, such as relief characteristics, geological and pedological properties, and also time-dynamic (variable) parameters: climate, hydrological environmental conditions, groundwater oscillations, land use, land cover, etc. [2][3][4]. With the development of anthropogenic systems, through the processes of urbanization and inadequate use of natural resources, the process of soil erosion has been intensified in various structural and functional entities [5][6][7]. It is estimated that about two billion hectares of the total area in the world are affected by soil erosion as a result of human activities [8,9]. According to some research, about 450,000,000 ha arable land in the world was found to be unproductive created through a symbiosis of GIS with the application of remote sensing in the assessment of soil loss resulting from erosion. The application of the G2 model itself is very simple and fast, and so far it has proven to be a very flexible model with available databases and it has been used as a tool for supporting decision making in watersheds.
Until now, research on the application of the G2 model has implied larger spatial dimensions and the use of global and pan-European data (for example, GlobCover, LC-CCI, MODIS products, CORINE), which do not have an adequate level of detail for the analysis of erosion processes of urban landscapes. The resolution of these databases is in the range of 100 to 500 m. The obtained soil losses are shown in a raster base of 30 m for two time periods. The first period is the year 2001, when the MPB entered into force. The second period is the present state from 2019 on the basis of available data. Although there are satellite images with a more detailed spatial resolution (Sentinel, SPOT, PROBA, etc.), satellite images from the Landsat mission were used, since it is a rich database with half a century-long, continuous series of images enabling a unison display of land cover for 2001 and 2019.

Study Area
The MPB is part of national spatial planning legislation which shapes the processes of spatial development and represents a strategic development plan for large cities in the Republic of Serbia. The MPB region for the capital of Serbia, Belgrade, covers an area of 778.52 km 2 and extends to two distinctly different, clearly separated relief units ( Figure 1). The northern part is a distinct lowland area, while the southern part is a rugged, typical hilly-mountainous relief. The natural boundary separating the two relief units are the Danube and Sava Rivers, which determine the microclimate of the city and represent the recipients of urban watercourses. The territory of the MPB has an average altitude of 120 m above sea level (range from 70 to 500 m above sea level) with a dense hydrographic network in the hilly-mountainous area and a complex system of drainage canals in the lowland area. Numerous occasional watercourses in the southern part of the MPB are one of the factors influencing the genesis of frequent torrential floods. The geological structure of the terrain is complex and represents the dominant direction of pedogenesis and the occurrence of certain soil types, with the largest share of chernozem, cambisols and alluvial sediments [38,39].

Data
The basic factors that affect soil erosion can be divided into static and dynamic input parameters, which are obtained on the basis of different databases and from different sources.
Climate: For the needs of the first period (2001), the climatic data were obtained empirically from meteorological stations for the observation period from 1971 to 2001. The second period (climate-based) is generated by a climate model from the EURO-CORDEX database, reference, and future periods, as well as scenarios of greenhouse gas emissions [40]. The simulation of climate parameters is based on the EURO-CORDEX database, analyzed using the RCP8.5 climate scenario (scenario of a constant increase in greenhouse gas emissions) expected for the mid-21st century. The main reason for applying this model is that a period of 20 years (beginning and end of the validity of the Master Plan) is too short for registering changes in the parameters related to climate change.
Physical-geographical characteristics: The influence of relief and slope of the terrain significantly affects the occurrence of erosion processes. For the needs of this study, a digital terrain model was applied in the form of a raster. The EU-DEM was implemented, which is a pan-European database of elevation available under the Copernicus program with a spatial resolution of 25 m.
Soil: Soil resistance to erosion processes is primarily conditioned by pedological characteristics which are an almost obligatory element of erosionist models. Soil resistance to erosion is correlated with soil characteristics, such as organic matter content, textural composition, structure, and water permeability. A digitized pedological map with a 1:20,000 scale [38,39] was used for the purpose of this study, while relevant data on soil characteristics were taken on the basis of relevant literature sources [38,39], as well as on the basis of soil samples used in projects by the University of Belgrade-Faculty of Forestry [41,42]. The characteristics required in the model were developed on the basis of soil samples collected from 47 sites in study area. The soil cover in the MPB is dominated by loamy soils, including sandy clay loam, loam, silt loam, clay loam and sandy loam.
Landsat satellite data: The Landsat is the longest satellite Earth observation program providing data about the Earth's surface since the early 1970s. Different types of data are available through the program, with Landsat 1, 2 and 3 being equipped with the multispectral scanner (MSS) and Landsat 4 and 5 operating thematic mappers (TM). Landsat 7 collects data with enhanced thematic mappers (ETM+) while Landsat 8 uses operational land imager (OLI) [43,44]. Landsat 9 will be the latest satellite with further improved capabilities and is scheduled for launch in late 2021 [43]. A large number of spectral channels with high spatial resolution (mostly 30 m and up to 15 m for panchromatic channels) make data from the Landsat program useful for a wide range of applications. These include, but are not limited to, mapping land cover and land use [44][45][46][47] and vegetation monitoring [48][49][50][51][52][53]

Data
The basic factors that affect soil erosion can be divided into static and dynamic parameters, which are obtained on the basis of different databases and from diff sources.
Climate: For the needs of the first period (2001), the climatic data were obtained

Random Forest Method for Land Cover Classification
Random forest (RF) is one of the popular and vastly used ensemble approach machine learning methods, which can be used to solve both regression and classification problems [54,55]. The main idea of such an approach is that an ensemble of "weak classifiers" can be combined in order to produce one "powerful classifier". In case of an RF model, the basic ensemble unit a is a decision tree [56]. RF is based on a bagging approach, which uses uniform sampling with a replacement (also known as a bootstrap sampling). Each tree in a forest is grown from a bootstrap sample of input data, while the rest of the data are used for a performance evaluation. After a forest of decision trees is grown, the output predication can be obtained by combining the predictions of each tree in the ensemble. This is done by averaging all predictions in case of a regression or determining the most common class in case of a classification [56]. The implementation of an RF classification algorithm available in a SNAP (Sentinel application platform) was used in this study [57,58].

Soil Erosion Using the G2 Model
The G2 model belongs to a group of empirical models and enables the estimation of soil loss in tons per hectare on a monthly or annual basis. The structure of the input data that are processed using this model is very similar to the RUSLE method [29,30] and the Erosion Potential Method (EPM) [31]. The G2 model consists of two main models: one for sediment production calculation (G2los) and the other for sediment transport (G2sed). The first G2los model relies on the basic principles of sediment production and one part of the input parameters is calculated according to the RUSLE model, while the other part of the input parameters is calculated according to the Erosion Potential Method. The second model G2sed relies on the principles of sediment transport by the Erosion Potential Method, by applying the retention coefficient (Ru). An innovation in the G2 model are the two parameters obtained using remote sensing data. The very structure of the model and the application of input parameters have undergone various modifications, and is still being developed. Monthly or average annual soil losses are determined on the basis of five erosion factors of which the size depends on space and time, according to the formula: E-soil loss for month (t·ha −1 ·); R-rainfall erosivity (MJ·mm·h −1 ·ha −1 ·y −1 ); Vvegetation retention (dimensionless); S-soil erodibility (t·ha·h·ha −1 ·MJ −1 ·mm −1 ); terrain influence (dimensionless); L-landscape effect (dimensionless).

Rainfall Erosivity (Factor R)
Rainfall erosivity represents the potential of rain intensity to affect the erosion of a certain type of soil in a certain period [29]. The factor is a measure of the erosion intensity of rain, i.e., an index that defines the interdependence of precipitation and the production of erosion material. The annual value of factor R represents the function of EI 30 which is obtained on the basis of the mathematical dependence between energy (E) and a thirtyminute rain intensity (I 30 ) [30,59]. For the purposes of this study, due to the lack of data on thirty-minute values of precipitation, average annual precipitation amounts were used by applying a mathematical relation [60,61]: R-rainfall erosivity (MJ·mm·h −1 ·ha −1 ·y −1 ); b 0 -empirical coefficient (MJ·h −1 ); Pmaverage annual precipitation (mm).
The values of the empirical coefficient range from 1.1 to 1.5 [60], and for the purposes of this paper, the adopted value of the coefficient is 1.1. This value has been successfully applied in the territory of the Republic of Serbia [62]. Factor R is presented in the form of a raster with a spatial resolution of 30 m, using the IDW (inverse distance weighting) interpolation method based on the relevant rain gauge cells that have an impact on the study area. The obtained values of the R factor before the calculation of soil loss were compared with the extrapolated European data for the Balkan Peninsula at spatial resolutions of 500 m, which are available from the European Soil Data Centre (ESDAC) of the Joint Research Centre (JRC) [63]. In addition to this database, also consulted data obtained from studies by Bezak et al. [64], who used the Rainfall Erosivity Database at European Scale (REDES) and Uncertainties in Ensembles of Regional Reanalyses (UERRA) rainfall data for the period of 1961-2018 with a spatial resolution of 1000 m. Based on the obtained data, the results of the R factor are consistent values with small deviations.. Compared to the European data for the Balkan Peninsula the R factor in the MPB territory was reduced by 9%. Compared to the REDES database, the value of the R factor in the study area was reduced by 12%.

Vegetation Retention (Factor V)
Vegetation retention factor (V) has similar characteristics to factor C in the RUSLE model. Factor V, by definition, represents the degree to which the presence of vegetation and its density are expected, as well as the way in which the soil is managed to protect it from erosion [29,65]. Unlike factor C, which is obtained on the basis of empirical tables or by remote sensing, factor V in the G2 model is obtained by combining both. The input parameters are a time series of vegetation indices and constant empirical parameters of land use according to the formula: V-vegetation retention (normalized, dimensionless); LU-empirical parameter for land use (dimensionless); FCover-fraction of vegetation cover (dimensionless).
Factor V is dimensionless and has a value from 1 to +∞. If the value is V = 1, it represents an unfavorable way of land management, i.e., the occurrence of bare land, while the values V > 1 represent a favorable way of land management, from the aspect of reducing erosion processes. Empirical parameter LU was created based on the interpretation of empirical data sets from the Erosion Potential Method [31] and the coefficient of watershed management X·a. Parameter LU has a value from 1, for soil that is bare and subject to erosion processes, to 10, for soils that are protected from natural erosion processes. The values of the empirical parameter of land use were assigned to the land cover obtained by the method of remote detection according to the values developed by Karydas and Panagos [32]. Fractional vegetation cover (FCover) represents a vegetation bio-geophysical parameter with normalized values from 0 to 1 (completely bare to dense vegetation). Fractional vegetation cover was obtained by applying the spectral vegetation index NDVI according to the formula [66,67]:

Soil Erodibility (Factor S)
The ability of soil to resist erosion agents by the sum of its properties is defined by the soil erodibility factor [29]. The formula for the RUSLE model was applied to calculate soil erodibility [29,30]: where S-soil erodibility (t·ha·h·ha −1 ·MJ −1 ·mm −1 ); M-textural factor defined as percentage of silt plus very fine sand fraction content (0.002-0.1 mm) multiplied by the factor: 100-clay fraction; OM-organic matter content in percent (%); s-soil structure class (s = 1: very fine granular, s = 2: fine granular, s = 3: medium or coarse granular, s = 4: blocky, platy or massive); and p-permeability class (p = 1: very rapid, . . . , p = 6: very slow); 0.1317coefficient for converting the final value into the value of the SI unit. Class permeability is determined based on soil textural classes and tables according to Renard et al. [30].

Terrain Influence (Factor T)
Terrain influence factor T is a derivative of the digital terrain model and is equivalent to the terrain factor LS in the RUSLE model. Karydas and Panagos [65] suggest an equation developed by Moore and Burch [68] for the calculation of the LS terrain factor, i.e., factor T in the G2 model according to the formula: where T-terrain influence (dimensionless, ≥0); As-unit contributing area (or flow accumulation), defined as the surface upstream flowing into a specific unit surface (dimensionless, >0); and β-slope gradient at the unit surface (rad).

Landscape Effect (L)
Landscape effect (L) has the role of quantifying the influence of natural and seminatural linear landscape elements (corridors) based on the movement of surface runoff and the development of erosion processes. Landscape ecological analyses emphasize the importance of applying edge metrics at the level of integral landscapes as indicators of landscape processes and functions to which soil erosion also belongs [69]. On the other hand, factor L can be interpreted as part of anthropogenic activities and can be compared with the factor of conservation measures P in the RUSLE methodology. The retention coefficient of sediment in the basin (Ru) in the Erosion Potential Method is based on similar motives [31] According to this method, a part of the eroded material does not reach the lowest level of the watershed, but remains in bays and depressions due to linear landscape elements (natural or introduced) that prevent its migration. To quantify this parameter, an algorithm is used for edge detection in the landscape structure, i.e., the Sobel filter, which is applied to the spectral channel from the near-infrared (NIR) part of the spectrum according to the formula [65]: S f -Sobel filter applied to the NIR spectral channel; DN max -maximum value of digital number on the spectral channel.
The values of the Sobel filter applied to the NIR spectral channel can have values from 0 to the maximum value of the digital number (DN max ). Digital number (DN) represents a pixel value that ranges from 0 to 255 in the case of an 8-bit record. The spatial distribution of digital numbers defines the image and the image area. The values of the landscape effect range from 1 to 2 [65].

Land Cover Accuracy Assessment and Change Analysis
The land cover maps over the study area for the epochs of 2001 and 2019 were produced at a 30-m spatial resolution using Landsat satellite data and the RF method. The produced maps differentiate six land cover classes: water, artificial surfaces, forest, shrubs, grassland, and agriculture. All available Landsat channels were used in the classification model. The usefulness of the spectral indices in the land cover classification process were proved by several studies [70,71]; therefore, an additional set of spectral indices was included in the classification as well. These include Normalized Difference Vegetation Index (NDVI) [72], Enhanced Vegetation Index (EVI) [73] and from a group of water spectral indices proposed by Xu [74] the Modified Normalized Difference Water Index (MNDWI).
The training polygons for each LC class were created by manual vectorization, separately, for each epoch. The additional independent set of polygons was also vectorized for validation purposes. The validation polygons are created for each epoch so that they are equally distributed over the study area ( Figure 2). Table 1 shows the basic information for the created training and validation datasets.
ital number on the spectral channel.
The values of the Sobel filter applied to the NIR spectral channel can have values from 0 to the maximum value of the digital number (DNmax). Digital number (DN) represents a pixel value that ranges from 0 to 255 in the case of an 8-bit record. The spatial distribution of digital numbers defines the image and the image area. The values of the landscape effect range from 1 to 2 [65].

Land Cover Accuracy Assessment and Change Analysis
The land cover maps over the study area for the epochs of 2001 and 2019 were produced at a 30-m spatial resolution using Landsat satellite data and the RF method. The produced maps differentiate six land cover classes: water, artificial surfaces, forest, shrubs, grassland, and agriculture. All available Landsat channels were used in the classification model. The usefulness of the spectral indices in the land cover classification process were proved by several studies [70,71]; therefore, an additional set of spectral indices was included in the classification as well. These include Normalized Difference Vegetation Index (NDVI) [72], Enhanced Vegetation Index (EVI) [73] and from a group of water spectral indices proposed by Xu [74] the Modified Normalized Difference Water Index (MNDWI).
The training polygons for each LC class were created by manual vectorization, separately, for each epoch. The additional independent set of polygons was also vectorized for validation purposes. The validation polygons are created for each epoch so that they are equally distributed over the study area ( Figure 2). Table 1 shows the basic information for the created training and validation datasets.   RF models were trained and applied to the study area separately for the 2001 and 2019 epochs ( Figure 3). Accuracy assessment is a key component of any properly conducted remote sensing study [75][76][77][78]. Visual inspection of the produced LC maps showed meaningful distributions of the classes over the study area, which was additionally confirmed through statistical analyses. This was done by comparing outputs maps with the vectorized validation polygons, pixel-wise, and producing confusion matrixes [70,79,80]. The overall accuracy (OA), producer's accuracy (PA), user's accuracy (UA) and kappa coefficient of agreement (Kappa) metrics were calculated and analyzed as well as these are commonly used metrics for evaluation of land cover classification results [77,80,81].
The results shown in Table 2 show an accuracy assessment for the 2001 land cover classification. User's accuracy for each category separately ranges from 70.39% (grasslands) to 99.82% (water bodies). Analyzing the production accuracy for each category separately, the values range from 76.65% (grasslands) to 100% (water bodies). The overall accuracy for all categories in the analyzed image is 94.52%, while the Kappa is 0.92. An assessment of the accuracy of the classification of the image from 2019 is shown in Table 3. According to the given table and confusion matrix, the user's accuracy shows the lowest values for the grassland category (71.26%), while the highest values were recorded for water bodies (99.82%). The producer's accuracy according to each category ranges from 66.67% (artificial surfaces) to 100% (water bodies). The overall accuracy according to all categories is 96.91%, while the Kappa is 0.95. When Kappa statistics were analyzed for both periods, the classified data agreed perfectly with the representative samples.  (Figure 3). Accuracy assessment is a key component of any properly conducted remote sensing study [75][76][77][78]. Visual inspection of the produced LC maps showed meaningful distributions of the classes over the study area, which was additionally confirmed through statistical analyses. This was done by comparing outputs maps with the vectorized validation polygons, pixel-wise, and producing confusion matrixes [70,79,80]. The overall accuracy (OA), producer's accuracy (PA), user's accuracy (UA) and kappa coefficient of agreement (Kappa) metrics were calculated and analyzed as well as these are commonly used metrics for evaluation of land cover classification results [77,80,81]. The results shown in Table 2 show an accuracy assessment for the 2001 land cover classification. User's accuracy for each category separately ranges from 70.39% (grasslands) to 99.82% (water bodies). Analyzing the production accuracy for each category separately, the values range from 76.65% (grasslands) to 100% (water bodies). The overall accuracy for all categories in the analyzed image is 94.52%, while the Kappa is 0.92. An assessment of the accuracy of the classification of the image from 2019 is shown in Table  3. According to the given table and confusion matrix, the user's accuracy shows the lowest values for the grassland category (71.26%), while the highest values were recorded for water bodies (99.82%). The producer's accuracy according to each category ranges from 66.67% (artificial surfaces) to 100% (water bodies). The overall accuracy according to all   The created land cover maps were used to calculate the areas under change. In Table 4, the changes of each category of land cover for the observed period are shown. The biggest negative changes, i.e., loss of land, were recorded under agricultural land and grasslands. Positive changes in the areas were recorded under the category of urbanized (i.e., artificial and non-porous surfaces) and under the category of shrub vegetation. By analyzing the change matrix, for the observed period, agricultural areas had a loss of 18.66%, which was followed by changes in artificial areas, grasslands and shrub vegetation. Constant population growth has increased the demand for housing, which has put great pressure on agricultural land, but also on "arable and uncultivated land" in peripheral urban areas [11,15].

G2 Model Results
In order to deal with different representations and spatial resolutions of the used input data, first the conversion to 30-m spatial resolution rasters of some datasets was necessary. The 30-m spatial resolution is chosen to match the spatial resolution of the Landsat images and corresponding LC maps. Table 5 lists the overview of all preprocessing procedures performed on input data and corresponding G2 model inputs calculated using them. The spatial distribution of climatic characteristics over the study area was obtained by applying the Inverse Distance Weighted (IDW) interpolation method. The IDW method is chosen because it is fast, simple and commonly used for interpolation of climate variables [82][83][84]. The EU-DEM used to describe physical-geographical characteristics also had to be resampled to a 30-m spatial resolution. This was done using the bilinear resampling technique. The last input dataset which required preprocessing is pedological map, available in scale 1:20,000. The vectorized pedology polygons were rasterized, where the output raster is also in 30-m spatial resolution. For detailed determination of the value of factor S, the IDW method was applied to soil samples. It can be considered that polygons can be converted to such spatial resolutions because various sources of errors present in the map (error of primary data acquisition, scanning and georeferencing errors, etc.) do not exceed output spatial resolution. After the input data have been successfully preprocessed, the inputs of G2 model are calculated using the previously described methodology. The results for both time periods were defined and analyzed in the following paragraphs.
The spatial distribution of factor R is shown in a raster form with a spatial resolution of 30 m. The average value of factor R is 731. 19 MJ·mm·h  By simulating the parameters according to climate scenario RCP 8.5 for the middle of the century and the available mean annual precipitation, rainfall erosivity was calculated for the second period (climate-based scenario). The values of rainfall erosivity in the study area also increased from 5% to 10%. The spatial distribution of rainfall erosivity has a range of values from 760.084 up to 841,398 MJ•mm•h −1 •ha −1 •y −1 , with an average value of 801.43 MJ•mm•h −1 •ha −1 •y −1 (Figure 4b).
Vegetation retention (Factor V) represents the result of remote detection in satellite images of different spatial and temporal resolutions, using a combination of spectral channels. In order to obtain data that is as relevant as possible, a set of 12 satellite images was used, one image for each month, in order to obtain average annual values.  By simulating the parameters according to climate scenario RCP 8.5 for the middle of the century and the available mean annual precipitation, rainfall erosivity was calculated for the second period (climate-based scenario). The values of rainfall erosivity in the study area also increased from 5% to 10%. The spatial distribution of rainfall erosivity has a range of values from 760.084 up to 841,398 MJ·mm·h −1 ·ha −1 ·y −1 , with an average value of 801.43 MJ·mm·h −1 ·ha −1 ·y −1 (Figure 4b).
Vegetation retention (Factor V) represents the result of remote detection in satellite images of different spatial and temporal resolutions, using a combination of spectral channels. In order to obtain data that is as relevant as possible, a set of 12 satellite images was used, one image for each month, in order to obtain average annual values. images of different spatial and temporal resolutions, using a combination of spectral channels. In order to obtain data that is as relevant as possible, a set of 12 satellite images was used, one image for each month, in order to obtain average annual values. Satellite images from the vegetation and non-vegetation periods were used. For each month, one satellite image was selected which had the lowest percentage of cloudiness.  The average value of factor S is 0.041 t·ha·h·ha −1 ·MJ −1 ·mm −1 which classifies these soils in the group of highly erodible soils [29]. The values of factor S range from 0.033 up to 0.07 t·ha·h·ha −1 ·MJ −1 ·mm −1 for the investigated area in which water bodies, lakes and urbanized areas were not taken into account (Figure 6a). The obtained values of coefficient S were used for both time periods. The average value of factor S is 0.041 t•ha•h•ha −1 •MJ −1 •mm −1 which classifies these soils in the group of highly erodible soils [29]. The values of factor S range from 0.033 up to 0.07 t•ha•h•ha −1 •MJ −1 •mm −1 for the investigated area in which water bodies, lakes and urbanized areas were not taken into account (Figure 6a). The obtained values of coefficient S were used for both time periods. Terrain influence factor T was used for both time periods and is shown in Figure 6b. The range of terrain influence values ranges from 0.01 to 45.5 with an average value of 1.34. The landscape effect (L) has a range of values from 1 to 2 for both periods. The average value of factor L for 2001 for the entire study area is 1.22 (Figure 7a), while for 2019 the average value is 1.29 (Figure 7b). Terrain influence factor T was used for both time periods and is shown in Figure 6b  Terrain influence factor T was used for both time periods and is shown in Figure 6b. The range of terrain influence values ranges from 0.01 to 45.5 with an average value of 1.34. The landscape effect (L) has a range of values from 1 to 2 for both periods. The average value of factor L for 2001 for the entire study area is 1.22 (Figure 7a), while for 2019 the average value is 1.29 (Figure 7b).  The defined parameters (R, V, S, T, L) are the starting point for the calculation of soil loss according to the G2 model. Potential soil losses for both time periods are shown in the form of a raster databases with a spatial resolution of 30 m. The raster database represents the loss of soil for each pixel per unit t·ha −1 ·y −1 . Figure 8a shows the spatial distribution of soil loss for 2001. The average value in the investigated area is 4.11 t·ha −1 ·y −1 , and the value range is from 0 to 200.86 t·ha −1 ·y −1 . Soil loss assessment was done only in areas that are subject to intense erosion processes. Soil losses were estimated using the G2 model based on the current land cover from 2019, according to the adopted scenario RCP 8.5 for the period of the near future and the simulation of climate parameters. The average value for the investigated area is 3.63 t·ha −1 ·y −1 , and the values range from 0 to 193.11 t·ha −1 ·y −1 (Figure 8b). In the case of both periods, the investigated area belongs to the category of medium erosion, although there was a slight decrease in the average value of soil loss of 11.68% in 2019. According to the spatial layout shown in Figure 8a,b, and according to Table 6, all forms of the category of destructiveness are present in the research area. In addition to the analysis of the shares of different areas, a comparison was performed in relation to the land cover as one of the most dynamic parameters. According to Figure 9, the average values of loss are shown in relation to areas with a high degree of porous surfaces. Water bodies and urbanized artificial surfaces with a high degree of non-porous surfaces are excluded from the analysis of soil loss and comparison. 11.68% in 2019. According to the spatial layout shown in Figure 8a,b, and according to Table 6, all forms of the category of destructiveness are present in the research area. In addition to the analysis of the shares of different areas, a comparison was performed in relation to the land cover as one of the most dynamic parameters. According to Figure 9, the average values of loss are shown in relation to areas with a high degree of porous surfaces. Water bodies and urbanized artificial surfaces with a high degree of non-porous surfaces are excluded from the analysis of soil loss and comparison.

Discussion
The intense expansion of urban functions and the need for new space for the growth of cities leads to changes in the structure of land cover, i.e., to the transformation of agricultural, as well as natural and semi-natural elements, into non-porous areas. Such surfaces cause the cessation of natural processes and have a highly negative impact on the land surface and groundwater regime causing the intensification of rapid surface runoff, the frequency of high waters, urban floods, and the intensification of erosion processes [85]. These changes lead to an increased risk for the infrastructure (road network and bridges), residential and commercial zones, sewerage and drainage systems, water supply sources and reservoirs, thus disrupts vital urban functions [11,86,87]. On the other hand, soil erosion can presents a vector of negative impacts and lead to endangered public health through sources pollution affecting water supply [88], the loss of fertility of agricultural land by nutrient removal [89,90], as well as to participate in factors that cause climate change, given that eroded areas have increased carbon emissions into the atmosphere [91] Based on these facts, there is an increased need for the research of modalities for determining relevant research procedures, methods and techniques that will indicate the

Discussion
The intense expansion of urban functions and the need for new space for the growth of cities leads to changes in the structure of land cover, i.e., to the transformation of agricultural, as well as natural and semi-natural elements, into non-porous areas. Such surfaces cause the cessation of natural processes and have a highly negative impact on the land surface and groundwater regime causing the intensification of rapid surface runoff, the frequency of high waters, urban floods, and the intensification of erosion processes [85]. These changes lead to an increased risk for the infrastructure (road network and bridges), residential and commercial zones, sewerage and drainage systems, water supply sources and reservoirs, thus disrupts vital urban functions [11,86,87]. On the other hand, soil erosion can presents a vector of negative impacts and lead to endangered public health through sources pollution affecting water supply [88], the loss of fertility of agricultural land by nutrient removal [89,90], as well as to participate in factors that cause climate change, given that eroded areas have increased carbon emissions into the atmosphere [91] Based on these facts, there is an increased need for the research of modalities for determining relevant research procedures, methods and techniques that will indicate the spatial and temporal scale of the occurrence and development of erosion processes and the environmental degradation forms. Contemporary cities are a dynamic conglomeration of urban functions that primarily depend on the built infrastructure and facilities, but the quality and safety of life of city dwellers also depends on the stability of natural and semi-natural elements as well as on the productivity of agricultural land. The functionality and safety of life of the inhabitants is related to the quality of the environment, so urban landscapes can be interpreted as an outdoor laboratory that has the potential to contribute to determining the modalities of the degradation processes monitoring.
Soil loss estimation using the G2 model was analyzed for two time periods, i.e., for the beginning of the entry into force of the Master Plan of Belgrade (2001) and for the period that is close to the end of its validity (2019). For the first period, the annual soil loss ranged from 0 to 200.86 t·ha −1 ·y −1 , with an average value 4.11 t·ha −1 ·y −1 . The analysis of the second period revealed an annual soil loss of which values range from 0 to 193.11 t·ha −1 ·y −1 , with an average value of 3.63 t·ha −1 ·y −1 . Comparative analysis identified a slight reduction in soil losses. The total area endangered by soil erosion was reduced from 526.46 to 481.39 km 2 , which is a decrease of 8.56%. Despite this, all seven categories are present for both observed periods. In other words, the change is not so pronounced that there is a disappearance of a certain category of loss, but there is a distribution of quantities within categories, as well as a certain change in area categories ( Table 6). The share of areas in the category of low medium erosion increased from 2.43% in 2001, to 2.95% in 2019, which represents a growth of 10.95%. In addition to this category, the areas in the category of medium erosion increased from 5.38% to 8.05%, and an increase of 36.78% was recorded for that area. In contrast, the areas of high and very high category of erosion were reduced by 30.02% and 21.92%, respectively, while the categories of high-medium erosion and very low erosion were reduced by 8.04% and 9.49%, respectively. Low-intensity erosion, which is classified in the category of low erosion, has experienced the lowest degree of modification, so that it still plays a dominant role in the dynamics of soil erosion. According to Figure 8a,b, the spatial distribution is shown by the categories of soil erosion intensity. Areas with a loss of over 1 t·ha −1 ·y −1 were identified for both periods in the southern areas of the area of hilly and mountainous regions with developed intermittent and permanent streams with torrential characteristics. The highest soil loss was identified on agricultural land located on steep slopes. The results of the G2 model indicate that the highest losses were recorded in places where grasslands and shrub vegetation change into agricultural areas. The categories of weak and very weak erosion with the soil loss of below 1 t·ha −1 ·y −1 were identified in the northern part, mainly in the plains and under agricultural areas. According to the land use defined in the MPB, these soil losses occur in planned economic, agricultural and green areas.
In relation to the spatial distribution of the category of erosion processes, soil loss was also analyzed for different land covers. Land use change as part of urban strategies and the consequent change in land cover are dynamic factors that change under the influence of political, social and bio-physical elements [92]. These changes were recognized as an essential group of driving forces that determine the in the quality of the environment of urban areas, of which the range of impacts reaches different spatial and temporal dimensions [93,94]. Urbanization processes mainly develop on areas that do not have a clearly defined purpose or are part of a city's developmental strategy, such as land cover classes that belong to grasslands and partly to agricultural lands. As these classes of land cover show the largest soil losses during the analyzed period, their decrease in the average value of soil loss was recorded. A decrease of 0.63 t·ha −1 ·y −1 was recorded for agricultural areas, while there was a decrease of 1.90 t·ha −1 ·y −1 for grasslands. Considering that the current state of land cover and simulated climatic parameters according to scenario RCP 8.5 were applied for the second period in the analysis (scenario of constant increase of greenhouse gas emissions), the reduction of average soil losses on agricultural areas and grasslands is justified. On the other hand, as a consequence of this process, the pressure on natural and semi-natural classes of land cover, i.e., on the remaining fragments of forest and shrub vegetation, increased. In these areas, there was an increase in the average values of soil loss according to climate scenario RCP 8.5. The average value of soil loss in forest areas increased by 0.54 t·ha −1 ·y −1 , while in the area of shrub vegetation it increased by 0.67 t·ha −1 ·y −1 . The obtained results indicate that in areas where urban functions are dominant natural and semi-natural elements as carriers of ecological characteristics are endangered by soil erosion processes. It is necessary to change the attitude towards these areas in urban plans, considering the fact that they represent a vital component of the urban landscape.

Conclusions
Soil erosion is a global environmental problem that negatively affects ecosystem services, safety and the economy of urban landscapes. Urbanization processes, as a consequence of the need for new spaces and solving of the population growth problem, cause the acceleration of natural soil erosion processes, which are becoming a serious problem for modern cities in light of climate change. Soil erosion modeling using authoritative and detailed databases, such as the G2 model applied with remote sensing methods, enable the precise identification of soil erosion intensity categories and the use of a climate model adds a dimension of reality to their prediction.
The results obtained by applying the G2 model indicate that the processes of urbanization cause the distribution of soil losses and sediment production in such a way that the classes of natural and semi-natural elements stand out as particularly fragile categories. The increase in porous surfaces adds additional pressure to these areas which increases the intensity of water and pluvial erosion. In addition, the classes of forest and shrub elements linearly become significant hazard zones in urbanized conditions and through an increased surface runoff and a more dynamic genesis of torrential floods, they transfer their impact to other surrounding areas. This constellation of altered environmental parameters causes other problems that degrade the quality of life in urban areas, such as the increased heat island effect, the loss of agricultural land fertility, acidification and salinization of soils, reduction of carbon content and compaction, increased greenhouse effect, etc.
The identified zones responsible for the excessive production of eroded material can be recognized as hazard zones and their spatial determination plays a significant role in the environmental risk assessment of urban landscape. The location of hazard zones and prediction of their development in the context of urbanization and climate change enables the elaboration of city development scenarios, analysis of consequences as well as configuration of a system of guidelines for imitating their further impacts. Adding a spatial context to erosion processes is the starting point for regional policies related to soil resources and soil conservation, which enables recognition of their place in spatial planning documents and city strategies. Techniques that support spatial analysis and synthesis of heterogeneous data have formed a powerful nexus (geographic information systemsmodelling methods-remote sensing) to create a database of soil erosion categories that allows a realistic assessment of the ecological stability in urban areas.