A review of the Gavrilović method (Erosion Potential Method) application

A detailed review of application of the Gavrilović method (Erosion Potential Method) and its modifications, with a focus on the potential surface erosion, is presented in the paper, together with the guidelines and recommendations for future analysis and research as needed for physical planning and urban development. The Gavrilović method results are based on the source of information, expert experience, accuracy, and level of detail of the model input and output data. For further analysis, the authors propose investigation of additional sources of erosion materials, such as construction plots in expanding urban areas.


Introduction
As indicated in the proposed Directive for the Protection of Soil, and the Amending Directive from 2004 [1], a significant increase in soil degradation processes has been registered in recent centuries. There are eight main processes causing soil degradation in European countries, among which erosion is considered to be the main and the most widespread [1]. According to Morgan [2], the occurrence of erosion processes, as well as their distribution and timing, are closely related to anthropogenic factors such as local, social, economic, and political conditions. It is well known that erosion processes can be triggered by heavy rainstorm events. The impact of flash floods on the erosion and sediment transport processes is considered to be quite significant. Until today, prediction activities remain difficult due to complexity of their nature [3]. The soil erosion, flooding, and channel management are considered to be interconnected environmental problems, especially at the catchment level where intensive precipitation can cause widespread gullies, mass movements, and flooding [4]. The erosion sediment production assessment, which begins at the parcel size and then broadens to the catchment level, is considered to be the basic component of an appropriate erosion management. This information is essential for decision makers when choosing future erosion mitigation and protection measures [5], and also to stakeholders involved in spatial planning and urban development. The erosion and flash flood prevention and mitigation measures accentuate the importance of building appropriate hydraulic structures (e.g. retention structures), and implementing various protection works (e.g. afforestation, removal of sediments from riverbeds) and other measures (e.g. defining an appropriate land use, informing interest groups and local population, …) in affected areas and areas of interest [6]. Various models are currently being applied in the area of water erosion risk and sediment yield assessment. They can be classified as empirical or regression models, conceptual models, and physics-based models. In addition, they can be classified as qualitative, quantitative and semi-quantitative models [2,7]. , are all semi-quantitative models whose basic description and comparison are given by De Vente [7] and Eisazadeh et al. [8]. In this paper, the authors provide a detailed overview of implementation of the Gavrilović method for the erosion risk and sediment assessment, as well as conclusions and suggestions for future development and improvement of the method and its application.

Gavrilović method (Erosion Potential Method)
The Gavrilović model, also known as the Erosion Potential Model (EPM), was developed by Slobodan Gavrilović based on erosion field research in the Morava River catchment area in Serbia in the 1960`s [9]. The method itself is based on the Method for the (1) (2)  [12], applied today in various studies.  [7,9] GRAĐEVINAR 68 (2016) 9, 715-725 A review of the Gavrilović method (erosion potential method) application Quantitative Classification of Erosion (MQCE), formally developed in 1954. During his research, Gavrilović discovered the possibility for further developing the MQCE method, which was used for defining the intensity of erosion. Extensions of this method were directed towards quantification of erosion processes by assessing the sediment transported downstream that reaches the control profiles [10]. Today, the method encompasses erosion mapping, sediment quantity estimation, and torrent classification, and has been extensively applied since 1968 for solving erosion and torrent-related problems in the Balkan countries [11]. It is currently being applied worldwide, from Croatia, Serbia, Slovenia, Italy, the Republic of Macedonia, Bosnia and Herzegovina, Montenegro, Iran to Chile (references are given in Table 4). Several distinct erosion processes can be assessed using the Gavrilović method, such as the surface erosion, downward erosion, or lateral erosion. In this paper the emphasis is placed on the application and modifications of the Gavrilović method for the assessment of potential surface erosion and suspended sediment and bed load in river network transported downstream. The most often calculated outputs of the method (equations 1-8, Table 1) are: the total annual volume of detached soil Wa (equation 1, Table 1) the erosion coefficient Z (equation 3, Table 1) the actual sediment yield Gy (equation 7, Table 1).
According to De Vente [7], this method can be characterised as a semi-quantitative method because it is based on a combination of descriptive and quantitative procedures. However, compared to other semi-quantitative methods named in the introduction, this method is the most quantitative because it uses descriptive evaluation for three parameters only: soil erodibility, soil protection and extent of erosion in the catchment. All other parameters are quantitative catchment descriptors.

Modifications to Gavrilović method
One of the first upgrades to the method was proposed by Lazarević [13], who noted in his work the need to adjust the assigned values so as to take into account the coefficient describing the type and extent of erosion ϕ, the soil protection coefficient X a and the soil erodibility coefficient Y. These three parameters, along with the slope angle, form the erosion coefficient Z. The purpose of this modification was to transform the definition of the erosion coefficient from its original meaning as soil erodibility to the today`s version as erosion intensity. Lazarević also modified the erosion intensity classification table represented by the erosion coefficient Z [13]. Tošić and Dragićević [12] proposed a new methodology for determining the erosion coefficient Z adapted for the use in GIS environments. It is based on the empirical methodology proposed by Gavrilović, and on its extensions by Lazarević. The very essence of their work involves the use of a PDA (Personal Digital Assistant) device with an integrated GPS receiver. The use of the device was combined with an appropriate software, namely ArcPad, to merge the GPS with the GIS. The aim was to directly determine an on-site coefficient for the type and extent of erosion (ϕ), and to transform the data accordingly to the erosion parcel condition. Yet another modification was proposed by Globevnik et al. [14] who suggested values for the soil protection coefficient based on the land cover classification CORINE ( Table 2). Later on, Fanetti and Vezzoli [15] suggested a change in the categorisation of the soil protection coefficient X a based on different land use categories ( Table 2). They were the first to consider urban areas as areas of potential erosion, thus assigning them a value greater than zero. They included several stages of urbanisation as well as various vegetation types, from growing cultures to pastures and forests.

Review of Gavrilović method application
This paper summarises application of the Gavrilović method (original one and modified version) based on the analysis of more than fifty different papers from relevant scientific bases that were available to the authors. The erosion risk/intensity, and sediment production and transportation, were estimated for more than fifty different catchments worldwide ( Table 4). The most commonly calculated method output (82 % of the catchments) is the total annual volume of soil W a . The value varies from 50 m 3 /km 2 /year for Rokava, Slovenia, [14,18,19] to 12,252 m 3 /km 2 /year for Khiav Chav, Iran [35,38]. The actual sediment yield, or sediment transported downstream, given for 38 % of the catchments, ranges from 37 m 3 /km 2 /year to 2495 m 3 /km 2 /year. A small number of analysed case studies (14 % of the analysed catchments) only provide an assessment of the erosion coefficient Z, thus giving insight the into erosion severity/ intensity for certain catchments, while not providing information about the expected sediment production [40].
Depending on characteristics of a catchment area, drainage density in particular, final results for the actual sediment yield can vary from quite small to the values similar to those estimated for the total annual volume of soil, or yearly amount of sediment available for detachment. In no case should the values obtained for the actual sediment yield result in values that are larger than those calculated for the total annual volume of soil [21]. One of the reasons for this outcome lies in the use of a different formula for the sediment delivery ratio that includes the drainage density parameter. In the original form of the Gavrilović method, a constant value of 4 was used instead of the drainage density formula. Later on, the model was modified, and the drainage density was taken to be the ratio between the primary and secondary river length and the contributing/catchment area. Results such as those for the Upper Sezar River can be obtained using the constant value instead of the length/area ratio. Overall, 37 % of catchment results showing the actual sediment yield are based on a constant value of the drainage density coefficient.

Gavrilović method, GIS and remote-sensing data
According to Thieken et al. and Vogt et al. [61,62], the reliability of final GIS based results is strongly correlated with the accuracy and level of detail of input data (topographic, land-use, and soil-data sources). Newer technologies, namely areal and satellite remote-sensing data, can be used to provide a substantially greater level of detail, and therefore simplify the erosion assessment procedure in the area of interest [15]. These technologies provide an improvement by enabling defragmentation of catchments to arbitrary cell sizes. For example, Bagherzadeh et al. [16,17] subdivided a catchment into eight homogeneous terrain units based on a visual interpretation of satellite images and field observations. Additionally, Globevnik et al. [14] and Milevski et al. [22] analysed applicability of the Gavrilović method in combination with a GIS technique. Their results demonstrated the decrease in predicted values for sediment production caused by erosion processes if calculated using parameters as a spatial variant, which is in contrast with the results obtained using the traditional/automatic method/catchment-oriented soil erosion map. GIS is used in 66 % out of the total of fifty-one (51) analysed catchments. In other papers, the use of GIS is not clear or GIS is not used at all, and 42 % use a remote sensing technology for the land cover parameter determination.

Land use/cover change and erosion mitigation measures
Since their development, GIS technologies have enabled analysis of land use/cover maps in greater spatial detail, while remote sensing technologies have facilitated generation of new and varied data sources for the same parameter. Solaimani et al. [24,25] analysed the effect of applying the change in land use as an erosion mitigation and land management measure, and showed that the output of the model predicts an 89.24 % decrease in erosion sediment yield. Although the authors did not analyse the sensitivity of outputs obtained by the Gavrilović method, this paper is the first one that refers to a significant oscillation in the predicted erosion sediment quantities that depends on the change in the soil protection coefficient representing the land use component in the Gavrilović [37]. One of the objectives was a cost-benefit analysis that demonstrated the economic and social impact upon soil erosion for a time period of 80 years. Dragičević et al. [51] were the first to analyse uncertainties in the magnitude and spatial distribution of annual sediment production predictions in the Dubračina catchment, Croatia, where several alternative land cover/use inputs were applied. They used three different land cover/use data sets: a CORINE land cover map, a Spatial Plan, and a Landsat 8 scene, and demonstrated the Gavrilović method sensitivity to different land cover/use inputs. Ristić et al. [39] analysed the effect of the change of hydrological conditions by restoring the catchment upon erosion and flood processes to define effective erosion mitigation and protection measures. They compared the outputs from the Kalimanska River catchment in Serbia for two time periods: 1967, before the restoration works, and 2010, after implementation of mitigation measures. The model showed a decrease in predicted volume of detached soil, and in erosion sediment transport via the river network. In another paper, Ristić et al. [43] predicted a 44.1 % decrease in annual sediment production of eroded material if a specific combination of biotechnical, technical and administrative measures were to be implemented in the Jalešnica catchment in Serbia. During their research, they noticed that the land use is closely related to erosion processes, and that it is the key to erosion mitigation and protection. Although technical structures in the riverbed are often applied as the erosion and torrent flood mitigation measures, they are not highly effective if used as the only measure in the catchment. The same analysis was conducted for the Manastirica and Kamišna catchments in Serbia [45]. The 40-year change in erosion sediment production caused by the impact of anti-erosion measures applied at the Celije reservoir in Serbia was analysed by Milovanović et al. [50]. They concluded that the implemented anti-erosion works, which included technical, biotechnical and biological activities, lead to a 49 % decrease in erosion sediment production and transported sediment yields in 40 years. Two other research activities showing results obtained by analysing the change in sediment production from past to present time were conducted at the Dragonja catchment (time change from 1971 to 1994) [14,60] and the Botonega catchment (time change 1989 to 2000) [14]. Both studies revealed a decrease in erosion sediment production due to abandonment of agricultural lands, followed by vegetation change in these areas to bush and forest lands. These changes were simultaneously backed by implementation of erosion control measures such as weirs, dams, special vegetation protection, bank stabilisation, etc., all of which has contributed to a substantial decrease in annual erosion sediment production.

Other applications of Gavrilović method
Lakicevic and Srdjevic [48] analysed connection between socialeconomic conditions and erosion processes in small catchments in Serbia, while Tošić et al. [49] analysed anthropogenic effects (demographic changes) on erosion processes in form of changes in population over time. Both papers concluded that human emigration leading to abandonment of agricultural land leads to a decrease in the intensity of erosion processes and sediment production in that area. Barmaki  Kazimierski et al. [5] analysed the impact of climate change parameters on the sediment yield production and, based on the Gavrilović method, they developed a methodology for estimating future sediment yield production for the Upper Plata Catchment in Chile, Bolivia. These authors generated projections for sediment yield production for the period until the year 2100 based on changes in temperature and precipitation, without considering potential changes in land cover/use. The significant difference between the observed and predicted erosion sediment yield values was associated with inaccurate interpretation of the observed data and deficiencies in the Regional Climate Models. Their analysis did not indicate either a significant change in annual sediment production over time, or a relatively small contribution of temperature in comparison to precipitation to the final sediment predictions. Bemporad et al. [44] used the Gavrilović method for calculating annual and monthly sediment production. In these calculations, temperature and rainfall were the only time-varied (monthly) parameters. For estimating the erosion sediment transport, they used the hydrological water discharge model and the equation for sediment continuity and motion by Hrissanthou (not Gavrilović). The values (only annual) were verified through a one-time field observation after a flood in 1992 that filled the newly built retention dam.

Comparison of Gavrilović method with other erosion assessment methods
The results obtained using the Gavrilović method have been compared with the PSIAC, MPSIAC and RUSLE methods based on all papers available to the authors. GRAĐEVINAR 68 (2016) 9, 715-725 A review of the Gavrilović method (erosion potential method) application Tangestani [30] compared the Gavrilović and PSIAC model outputs and established that the PSIAC model is more reliable for determining the areas of very high erosion potential. A visual field overview with GPS confirmed a good estimation for areas of moderate and heavy erosion with the Gavrilović method, and a poorer accuracy for areas with slight erosion potential. Another comparison with the PSIAC method [16,17] showed the same pattern for the predicted sediment yield by both methods, with a correlation coefficient of 0.95, which confirmed applicability of both methods for semi-arid and arid catchments. Ghobadi et al. [42] compared the Gavrilović method with PSIAC and MPSIAC and concluded that the Gavrilović method is not suitable for weather conditions in Iran, and that it provides much less accurate annual sediment production assessments compared to the MPSIAC method. In their assessments, they also used a simplified formula for the sediment delivery ratio. Petraš et al. [52] compared results obtained using the RUSLE and Gavrilović methods with on-site observations and concluded that the RUSLE method is more compatible with the on-site data measurements for the Abrami test field. In comparison with some other procedures, the Gavrilović method does not explore the physics of erosion processes and is therefore advantageous for areas where minimum data are available, or where there is a lack of previous erosion research. As such, the method can provide not only the amount of sediment production and sediment transport, but also the erosion intensity as a preliminary result, and indications or areas of potential erosion threats.

Field measurement and Gavrilović method verification
Out of all analysed catchments, verification of results is presented in papers for fifteen (15) catchments only (Table  5). In these papers, different verification methods are applied, depending on available equipment and accessibility of the terrain. Measurements of sediment yield on erosion plots were conducted at the Rokava (Dragonja) River basin in Slovenia [18,19], Abrami [52] and Jukani (Butonega), Croatia [14,20]. At the Bregalnica basin, Republic of Macedonia [19,20], a very good correspondence was obtained between the results obtained using the Gavrilović method and onsite measurements. Haghizadeh et al. [21] and Tangestani [30] used a comparison of model output results with field observations and a GLASGOD (Global Assessment of Soil Degradation) map as a verification method. Bagherzadeh et al. [16,17] verified model outputs by a field survey using GPS and a visual comparison of areas characterised as areas with moderate and heavy annual sediment yields. Amini et al. [10] applied the Gavrilović method to the Ekbatan Dam drainage basin in Iran and concluded that sediment yield can be overestimated by this method because it lacks the grain size distribution structure, humus concentration, slope morphology, and runoff parameters that affect erosion processes, all of which are usually included in physical models, which is not the case for empirical models such as Gavrilović.
Kouhpeima et al. [33] analysed five different catchments in Iran and used comparison results to measure sediment deposits in the reservoir as a verification method. The same method was also used in the Prescudin catchment, Italy, [44] and it showed minimum deviation between predicted and measured sediment yield values. Nuclear probes for suspended-load measurements were used at the Esino and Musone river basin in Italy [26]. The measurements exhibited some deviations in comparison with the overall sediment yield production estimated with the Gavrilović method but, overall, a good order-of-magnitude correspondence was obtained concerning yearly sediment yield. It was concluded that further measurements are necessary because the Gavrilović model considers total sediment load, whereas the measurements conducted take into account suspended load only. Other verification methods include the use of a PDA device and on-site observations [49], and certain verification methods remain unspecified in papers [41,42] but provide poor overall ratings for the Gavrilović method, which is said to overestimate sediment yield [41].

Discussion and conclusion
A detailed review of practical application of the Gavrilović method is presented in the paper. The most commonly calculated results using the model are the total annual volume of soil and the erosion coefficient. The actual sediment yield is calculated for no more than 38 % of analysed catchments.
Although several modifications of the model have been used over the years, different variations of the model continue to be applied. These variations concern assessment of actual sediment yield involved in transportation. The analyses exhibit better results and correspondence with on-site measurements when a modified formula is used for the sediment delivery ratio. A modified sediment delivery ratio uses drainage density as a ratio between the primary and secondary river length and the catchment area. If the simplified (original) formula is used and the ratio is replaced with a constant, the values obtained using the model can exceed the predicted values for the total annual volume of soil or the overall yearly amount of detached soil. Therefore, the authors conclude that the use of the formula for drainage density is recommended for all future analyses as a means to avoid incorrect results indicating larger values for the actual sediment yield compared to those of the total annual volume of soil. However, none of the analysed papers includes an explanation as to why a given formula, original or modified, is preferred to any other. Additionally, these papers do not provide a comparison that could be used to roughly estimate the error/difference between the calculated and measured values if both formulas are used. The evolution of the Gavrilović model began with the development of GIS technologies. However, this method has not as yet fully benefited from all possibilities offered by GIS. For example, the actual sediment yield or sediment involved in transportation is calculated within the method for the entire catchment/sub-catchment and refers to the value representing sediment transport, measured at the tow of the catchment. Today, GIS technologies enable assessment of each cell within the catchment and, as such, they can provide an estimation of the material transported in each cell representing the river. This approach can to some extent simplify the process of choosing the best location for field measurements in less accessible catchments, and provide multiple options as adequate positions for field measurements. Thus, the verification of the method in terms of the assessed parameter for actual sediment yield can also be simplified and conducted on any part/length of the river, which can potentially lead to more frequent calculations of this parameter. To achieve this, the analysis must be narrowed down from the catchment and sub-catchment assessment to the cell resolution, and later gradually broadened to the catchment size. Unfortunately, this procedure will continue to depend upon resolution of the available input data. Lazarević, Globevnik, Fanetti and Vezzoli significantly improved the method using changes in the assessment of descriptive parameters within the model. It is important to note that certain catchment areas are currently affected by substantial changes in type, extent and density of vegetation cover and by expansion of urban areas. Therefore, if this is considered, the land use/cover parameter can be considered as an extremely important parameter that will affect the final estimated values, as shown in a number of previously mentioned papers [24,25,27,28,29,58,60]. The changes in land use will not only affect the changes in the soil protection coefficient but subsequently the coefficient of type and extent of erosion, whose contribution to Gavrilović method and sediment production should not be neglected. For such areas with intensive urban changes, the parameter for urban areas (Table 6) is recommended in future analyses. It is often forgotten in erosion analysis that agricultural areas and areas with low or no vegetation cover are not the only source of eroded material in a catchment. Therefore, all catchments are unique and complex in their own way, and additional sources of erosion material should be considered, such as construction areas in expanding urban regions. These areas, although short lived, have a substantial impact on the yearly amount of erosion sediment production and should be considered when planning future activities in the catchment. Another source of erosion material that is rarely considered are residential areas with small green plots used mainly for agriculture. In larger towns, such areas are not considered to be significant; however, in towns and villages where such residential areas are often represented, this can be considered a problem and an additional source of erosion material that is often forgotten and simply classified as urban/rural area. Therefore, a new categorisation for the soil protection coefficient for urban/ rural areas, including undeveloped areas designated for urban development in the near future, is suggested in Table 6. Such a categorisation would change the model output information concerning erosion intensity and the total amount of erosion material. Table 6. Proposed assessment of soil protection coefficient for urban areas X a

Descriptive evaluation Numerical evaluation X a
Dense urban areas with no or little green areas 0,05 Scattered urban/rural residential areas with green plots used mainly for agriculture 0,3 Scattered urban/rural residential areas with green plots used mainly for agriculture using agronomic, soil management or mechanical anti-erosion measures 0,1 Construction areas 0,9 GRAĐEVINAR 68 (2016) 9, 715-725 A review of the Gavrilović method (erosion potential method) application Land use/cover parameters greatly influence final results of the model, and lead to predictions of decreased erosion production if appropriate land use management is applied. Dragičević et al. [58] highlighted the problem through generation of different results by simply using a different land use/cover input source. It can therefore be concluded that reliability of final results is strongly correlated with the data source, experience of the expert in charge of map production, and the accuracy and level of detail of input data. The expert conducting erosion analysis should also consider different data sets/maps available for the same parameter, compare their differences and, based on field surveys and local population information exchange, choose the best option for future analysis, as shown in [27,28,58]. In the future, it could be interesting to analyse the interconnection between the Gavrilović model outputs from small catchments and its geology, which would lead to possible further modifications of the model. It should be noted that the verification aspect of the analysis is omitted in most of the analysed papers, which leads to a shortage of information concerning the adaptability and applicability of the Gavrilović method to different areas varying primarily in terms of geology and hydrology. The lack of these data has also hindered possible extensions of the method because these data have yet to be provided. Additionally, it is noted in several papers that a strong correlation exists between the knowledge and experience of the erosion expert and the deviation of predicted and measured sediment yield. Not one of the papers containing verification actually addresses sensitivity of the model and the uncertainty of overall results regarding the source of input data. The verification of the models should be conducted with greater frequency so as to obtain a better correspondence between on-site measured values for sediment production and those obtained with the model.