Influence of pit removal algorithms on surface runoff simulation

The digital terrain model (DTM) with removed pits is a precondition for hydrologic analysis. Two pit removal methods, the carving method and the filling method, are investigated in the paper for three different geomorphometric areas. The input data are digital elevation models (DEM) created by interpolation from photogrammetric surveying in the spatial resolution of 5x5 meters. A recommendation for the use of a more acceptable DEM creation method without pits, involving minimum geomorphometric changes, is given based on geomorphometric and statistical analyses. The importance of selecting a correct algorithm is proven.


Introduction
The digital terrain model (DTM) is the basis for all geomorphometric and hydrological analyses.The following three data sets are officially used for DTM in the republic of Croatia Geomorphometry is a scientific discipline that uses DTM as the main basis in its research.The analysis of DTMs enables differentiation of morphometric terrain parameters and groundsurface features.The significance of DTM in morphometric analyses is considerable, and the results of geomorphometric and hydrological analyses derived from this model are dependent on the DTM properties and quality, and on the interpolation function that is used in the scope of geomorphometric and hydrological analyses.In addition, geomorphometric analyses necessitate scientific and professional cooperation between the profession in charge of creating the digital elevation model (geodesy) and other related geo-sciences (hydrology, geology, forestry, etc.).The study of the possibilities of qualification and quantification of problems related to the determination of runoff within a catchment area taken as a unit of space, based on the knowledge and special DTM development and analysis techniques, holds a significant place in many scientific and professional disciplines: hydrology and meteorology, hydraulics, geology, paedology, morphology, and also enters in the framework of biological, forestry and agronomical sciences.Pragmatic significance of the related studies can be seen in the resolution of water management problems relating to drainage area regulation, regulation and use of natural watercourses, torrent regulation and, globally, in flood protection of urban and rural areas [2].Some parameters needed for creating DTM are adjusted to the needs of the users (resolution, level of generalisation, etc.), but a common need of all subjects included in spatial analyses is the creation of DTMs representing the actual relief at a satisfactory level of quality.The DTM is a significant input component of topography in flood modelling, and so advanced microcorrections in relief increase accuracy of runoff simulation and determination of flood areas [3].The influence of DTM quality on the description of drainage areas, as the input data in algorithms for the determination and monitoring of erosion and deposition trends, is of crucial significance for obtaining an appropriate quality of results [4].The DTM quality is an important and unavoidable parameter in all analyses of surface runoff in the algorithms for monitoring erosion, deposition, landslides, flood areas, environmental impacts, and water-borne pollution.The stability of infrastructure facilities (road structures, railway embankments) is indirectly affected by these occurrences.In addition, a detailed knowledge of watercourse hydrology and flow hydraulics can alleviate harmful flood effects in bridge zones, and can be of crucial significance for the safety of people [5].
In engineering practice, it is often necessary to estimate as best as possible drainage area runoff incidences for which there are no measurement data, or where such data have not been adequately studied.A particular problem are small irregularlyshaped drainage areas where runoff can not be reliably calculated on the basis of frequently used simple parametric methods (such as the rational method or SCS method [6]).Conceptual parametric models, popularly known as CRRS models ("Conceptual Rainfall -Runoff Models"), are used in hydrologic definition of drainage area runoff under the abovementioned conditions.The CRRS models are characterized by a great number of climatological and physicogeographical parameters that describe and modify the runoff process in parametric models [7].The CRRS models contain a considerable number of parameters that can not be determined from measured values describing drainage area properties, and so the use is made of calibration [8].The CRRS modelling is made of a relatively small number of structurally arranged elements where each of the elements simulates a particular phase of a process taking place in a drainage area, i.e. a particular phase of the process in the scope of which the gross and/or effective rain is transformed into the corresponding runoff hydrograph.The key non-hydrological input element for the mentioned CRRS modelling is the DTM represented through the discretional system of drainage area elements as the carrier of information about the terrain slope, water flow direction and orientation, length of the route traversed by water flowing across the drainage area, time of flow, and reservoir elements that simulate the drainage area retention effect.Algorithms for the determination of the direction of surface runoff are used to simulate the movement of water, from higher to lower points on the terrain, exclusively by means of gravity flow, and are directly related to the quality and credibility of the DTM.Depending on whether the algorithm concentrates all water in the centre of the pixel, or distributes water across the entire pixel, we can differentiate one-dimensional algorithms such as D8, Rho8 and KRA algorithms [9][10][11] and twodimensional ones, i.e.D∞, FD8 and DEMON algorithms [12][13][14] of the CRRS surface runoff model.Regardless of whether we are dealing with one-dimensional or multi-dimensional surface runoff models, the key assumption is that each analysed DTM pixel has at least one pixel in its immediate neighbourhood that is lower than it, and toward which it directs its flow.In case of a local minimum (depression), this condition is not fulfilled and numerical problems are encountered during calculation of surface runoff.That is why local minimums (depressions) must be removed from the DTM before the surface runoff is calculated.Surface depressions (sinks, pits, etc.) are treated in DTMs as disturbances hindering hydrological modelling.When studying the use of hydrological simulations in surface runoff models, such as in determining the hydrographical network, direction of surface runoff, water accumulation across the terrain, and drainage area Influence of pit removal algorithms on surface runoff simulation calculation, it was established that the pit removal of depressions is used as a preliminary step in hydrological analyses depending on the availability of the module implemented in the GIS program in which the analyses are made, rather than depending on the geomorphometric area.The removal of depressions changes the DTM geomorphometry, which influences calculation of various geomorphometric parameters that are significant for further hydrological calculations, and so the choice of pit removal method is considered to be an important preliminary step in hydrological analyses.Although in Croatia the term DTM comprises both the original and derived DTM, the term digital elevation model (DEM) is used below for such model because of the regular structure of heights (altitudes) stored in the derived DTM.Recommendations for selection of a more appropriate depression removal algorithm, considering geomorphometric characteristics of the test area, will serve as the basis for future projects based on geomorphometric and hydrological investigations, which will open the gate of cooperation between geodesy and other geo-sciences.

Pit removal algorithms in DEM
In geomorphology, a depression is a sunken or depressed landform on the Earth's surface regardless of its origin, shape, or size.It is a part of the land surface or sea bed lying below the level of the surrounding terrain due to tectonic displacements, lowering along fault lines, epeirogenic bending or folding [15].In DTMs, depressions occur when a pixel, or a group of pixels, is completely surrounded by neighbouring pixels lying at a greater altitude.Although some of them may be natural, especially in the glacial and karst areas [16], a great majority of depressions constitute a DTM irregularity, and are therefore a significant hindrance to hydrological simulation and computation of flow parameters (infiltration lowering process, reduction and significance of water storage in depressions as related to real runoff process in the area under study, direction of flow, accumulation of flow, flow rate, concentration time, and other hydrological parameters).Although this paper focuses only on the influence of the of the removal of depressions on the surface runoff process, it is significant to note that the removal of depressions greatly influences water balance in the drainage area as a whole.The DTM data are collected by photogrammetric mapping from aerial photographs, by laser scanning, and through terrestrial measurements.Nevertheless, regardless of the way in which the data are collected or further analysis is made, it is difficult to discern from the DTM data whether the depression is a real terrain feature or it occurred during the data processing process (artificial depression).Artificial depressions are often due to errors made during sampling (inadequately classified input data) or to interpolation, generalisation, rounding interpolated values to a lesser accuracy, balancing of pixels within the area, or smoothing as a result of oversampling.Especially the present-day data collection procedures for DTM -such as airborne laser scanning or image correlation methods in digital photogrammetry -unfortunately have a low selectivity level despite exhibiting a high data collection frequency.Consequently, unrealistic morphometric forms inevitably occur during the DTM interpolation, including also the occurrence of depressions, which actually do not exist in real relief [17].The occurrence of artificial depressions is also related to the vertical and horizontal resolution of the original altitude and relief diversity data.Artificial depressions are more common in areas exhibiting low diversity of relief, which can be attributed to limited vertical accuracy of DTMs [18].On the other hand, it is sometimes difficult to define natural depressions (forested areas) as most airborne surveying techniques (RADAR, LIDAR) collect data for creating a digital model of the area, rather that for digital modelling of elevations (heights) [19].As no technique aimed at differentiating natural from artificial depressions (except on-site surveys) is capable of providing fully reliable results, and as the procedures are also quite complex, the pit removal algorithms operate in an unselective manner in all GIS programs.Most natural depressions finally overflow into the downstream discretisation pixel and, despite unselective removal of depressions, such a corrected model can in fact be accurate.The earliest method for finding depressions by means of modelling is based on the procedure that is carried out in two steps.The first step involves definition of the lowest point by the row-by-row model scanning.In the second step, delimitation pixels (i.e.pixels where depression ends) and overflow point (pixel from which evacuation is possible) are identified for each depression point.The procedure is iterative, and the number of iterations depends on the number and complexity of depressions.A simpler and more efficient method involves the procedure called "from lower to higher" and "from edges Sanja Šamanović, Damir Medak, Duška Kunštek inwards".It is based on the data structure priority known as the binary heap [20].Unlike row-by-row scanning, this method requires only one search to recognise depressions and overflow points.The CRRS models, in which hydrographical network is defined based of the surface flow simulation using the GIS technology, require that all depressions are first removed from DTM and the usual practice is to identify and remove them from DEM already at the first procedural step of hydrologic analysis.The presence of depressions in DEMs prevents simulation of the flow of water on the surface of terrain, and along the high-rank and lower-rank watercourses, which consequently results in unrelated water flow patterns in the drainage area with defragmented hydrographical network, and with false sub-drainage area and drainage area space determinants (watersheds).To ensure homogeneity and full connection within the hydrographical network, the water from each DEM pixel must be re-directed toward the lines of concentrated runoff and, at that, the presence of depressions in the DEM creates interruption in the runoff line continuity, and modelling of other spatial hydrological processes (time concentration modelling, computation of height and duration of flooding, modelling retention and retarding time in drainage area, erosion, etc.).Two main methods are used to remove depressions (Figure 1): carving or breaching method and filling method.Although Soil [21] and Lindsay and Creed [22] concluded when comparing corrected DEMs that the method combining two basic methods influences the original DEM the least, most GIS programs still use the simpler filling method because of the robustness of the algorithm.In addition, isolated but deep depressions are problematic for the carving method [19] as they change the terrain configuration considerably.The modern technologies of collecting data such as the LIDAR method also encourage development of new algorithms based on optimised depression-removal tools [23].When comparing two methods, first implemented in ArcView (filling) and second implemented in TOPAZ (carving) on a 30 meter DEM in the Valley region in Virginia, Srivastava established that, for the tested area, the carving method changes less pixels in the flat-land area, while the canal network changes near the modified pixels [24].The study of depressions on the model, and comparison of two pit removal methods using DEM in various resolutions (50, 100 and 200 meters), led to the conclusion that greater resolutions result in the occurrence of greater depressions, and in greater water storage in the drainage area [25].Standard method for removing depressions using the DEM (filling method) has been criticized by the academic community because of greater influence on the modelling of the surface flow network system [26], and so great efforts are currently made do develop new algorithms.The filling is the most widely used method primarily because of its simplicity.Early algorithms for the removal of depressions used the filling methods and were created for smaller areas, low resolution, and smoothened DEMs.The basic deficiency of this method is a significant change of the original digital model [9,27].The filling is the procedure by which the elevation of the depression is increased until reaching the point from which runoff is possible [28], i.e. the filling increases until the overflowing point, and the results are pixels whose value is always the same or greater than the DEM.The procedure is repeated until all depressions have been filled or, more precisely, until runoff lines have been associated with all pixels.Most depressing-filling algorithms are based on one-dimensional algorithms for determining the runoff direction [29], where each pixel assumes one out of eight values, depending on which of the eight neighbouring pixels lies in the direction of the steepest slope and, at that, the pixel value is assumed to be the centre of each pixel.One of the first algorithms was developed by Jenson and Domingue.It creates flat areas that enable further runoff.The procedure is complex and long-lasting, and it involves filling of individual depression pixels, by increasing the elevation of each pixel to the lowest neighbouring elevation at the edge of the depression [27].Each pixel in the corrected model will be a part of the path of a monotonously falling function that ends at the edge of the model.The path is composed of neighbouring pixels in the horizontal, vertical or diagonal directions of the grid, the value of which is continuously decreasing.The procedure developed by Jenson and Trautwein is similar but slower because of a greater number of iterations [30].Each pixel situated in the depression is elevated to the lowest value at the edge of the depression, and so each pixel of the output image has at least one path leading across lower pixels to the edge of the depression.The path consists of pixels in the neighbourhood (all eight directions) that must meet the requirement involving continuous decrease of the value.If the input area is deducted from the output area, each pixel of the obtained value is equal to the depth of the depression in the unit of the input area.The removal of depressions using the filling method developed by Mark is similar to this method but, in an intermediate step, it identifies drainage areas for pixels that do not have neighbouring pixels at lower elevations.The module identifies pixels with elevations lower than the lowest elevation of the drainage area boundary, and uses these pixels as initial plane in the drainage area delimitation procedure.In the drainage area generation procedure, the drainage areas increase one by one, and flat areas are attributed to the first neighbouring drainage area [31].The depression removal procedure based on the method developed by Planchon and Darboux can be divided into two phases [32].The first phase consists of flooding the entire DEM with water, while in the second phase the excess water is evacuated from the model.The flooding is operated by associating the highest water level to all DEM pixels, after which the surplus water is iteratively evacuated from each pixel.To find the downward path for each pixel and to ensure full evacuation of surplus water, the entire DEM is analysed in all eight directions.At the final stage, the water is drained from the depression to the level of the highest overflow point, thus creating the line of flow all the way to the DEM boundary, which results in the model devoid of depressions.The first depression removal algorithm founded on the carving method was created by Rieger.It was based on pixel lowering in the vicinity of the depression.Subsequent algorithms involved depression removal with carving method based on the carving canal length, which results in development of the network that creates the runoff path from the depression to the first neighbouring lower depression [33], creating the network of lines of concentrated runoff.The method is in fact most frequently used for solving depressions Influence of pit removal algorithms on surface runoff simulation that are not to far away from one another [34].Starting from the depression, the algorithm does not fill the depression but rather carves the terrain so as to enable evacuation from the depression, i.e. the carving lowers the elevation of pixels situated along the evacuation path [35] thus evacuating water to the first next lower pixel or to the edge of the model.Depression removal according to Conrad is based on the pixel value lowering so as to find an exit from the depression, and it consists of two steps [36].In the first step, the algorithm recognises the route of the concentrated runoff line and, at that, the area without depression is marked with zero, and then the values ranging from 0 and 8 are attributed to the depression, depending on the route by which evacuation from the depression will be operated.The second step involves depression removal using the concentrated runoff calculation line.Wang and Liu have proposed a new efficient depression filling method that can simultaneously determine runoff direction and spatial distribution of drainage areas, all this in a single step [37].They have introduced a new concept of evacuation from depression and a progressive method involving an optimum evacuation route, which is based on the structure of high priority data, and on the analysis of an optimum evacuation route.The jump limiting method developed by Martz and Garbrecht limits the canal carving length to no more than two pixels, while all other depressions are solved by filling [34].The approach developed by Lindsay and Creed involves either filing or evacuation depending on which method causes the smallest change to the DEM [22].An optimum depression removal approach was developed by Soille who combined the filling and evacuation methods so as to reduce differences in elevation between the original DEM and the modified one [21].The approach is based on filling up to a certain level, after which the evacuation is used.The level has been set in such a way to reduce the number of differences in height between original data and modified ones, or to reduce the number of modified pixels.Most GIS programs contain in their modules the Planchon and Darboux method published in 2001 or the Wang and Liu method from 2006 [26].GRASS and IDRISI presented the filling method based on Jenson and Domingue (1988), while the Arc GIS uses the method developed by Tarboton et al [38] and Hutchinson [39].The carving method is used by GIS Whitebox GAT, TopoToolbox [25].SAGA enables selection between a great number of depression removal modules: Planchon/Darboux, Wang/Liu, XXL Wang/Liu or Conrad.It is precisely because of this possibility of comparing two basic pit removal methods within a single GIS software -which is not possible in case of other software programs due to implementation of one method only -the open code software SAGA was selected for defining an optimum method.The analyses were conducted using the carving method developed by Conrad and the filling method developed by Planchon and Darboux.The method developed by Planchon and Darboux starts from the initial area Z and defines the final area Wf in two steps [32].In the first step, in addition to boundary model operations, the method initialises the area W with infinite elevations.In the second phase, the area W elevations are iteratively reduced and, finally, the area Wf is obtained, keeping the two properties: -Wf ≥Z for each pixel c of the area Wf there is a path leading to model boundary and there is a slope ε leading from one pixel to another (downward path).
When a pixel is at the boundary of the model, i.e. when W(c) = Z(c), the W(c) reaches its minimum, and so W(c) = Wf(c).A consequence of this is that the pixel elevation c no longer changes (Equation 1).
Inversely (Formula 2), when Z(c) is lower than W(n), W(c) can be described as: The next step involves a new iteration that brings W to the final result (Wf).
Olaf Conrad developed a method that starts from the digital elevation model, makes use of the runoff path, and creates a modified digital model.A one-dimensional runoff algorithm D8 is used to define the runoff path.The runoff direction is determined starting from the depression point via the minimum elevation pixel that is lower than the depression point, taking into account eight neighbouring pixels [40].
If such neighbouring pixel does not exist, the direction is considered undefined.Matrix values are defined as the runoff direction.
Although this principle enables definition of multiple runoff directions for evacuation from the depression, only one direction is ultimately allowed.In order to compensate for the route increase in diagonal direction, the elevation of the neighbouring pixel is defined with the factor 1/√2.The occurrence of neighbouring pixels of equal value is thus reduced.In case of equal-value pixels, the first pixel is selected in a clockwise direction, starting from the north.The depression removal procedure starts by defining the depression pixel.The overflow point (Equation 3) is defined as the point at the depression boundary that has the lowest height difference with regard to the depression.Let us assume that: N(i,j) = {(i-1,j-1), (i-1,j), (i-1,j+1), (i,j-1), (i,j+1), (i+1,j-1), (i+1,j), (i+1,j+1)} (3) with eight neighbouring pixels with the overflow point (i,j), and are the set of points at the boundary of the depression (p,q).The overflow point (i,j) in the overflow area A(p,q) is defined as follows in Equation ( 5): where ELDIF calculates the difference in altitude of individual points.
The area increases through examination of boundary points, and the procedure stops when all depression points are determined and the overflow point is established.Once the overflow point is established, runoff directions are modified so as to become oriented toward this point.The runoff direction from the overflow point is defined as the direction of the neighbouring lowest-height depression toward the overflow point.
From the hydrological standpoint, the main difference between the two methods lies in the fact that the carving algorithm calculates the flow trajectory and sticks to this trajectory regardless of the number of pixels in the input sample of the terrain, while the depression filling algorithm increases the level of pixels (changes their height) in the area, while their number as the input signal remains the same.This enables us to conclude that -besides one-dimensional signals -the carving algorithm is more selective and hence more appropriate for defining the trajectory of flow along the terrain.

Research zone and data
Input data for testing two depression-removal algorithms include the DEM resolution of 5x5 meters at three geomorphometrically distinct areas in the Republic of Croatia (mountain massif of Biokovo, Papuk mountain, and flat areas in Eastern Croatia within the Nijemci district) (Figures 2 and 3).The data for this analysis, i.e. aerial photogrammetric surveys on the scale of 1/20000 (cyclical aerial photogrammetric surveys of "Dubrovnik and its surroundings", 2005, "Western Slavonia", 2001, and "Eastern Slavonia", 2002), external orientation parameters, and camera calibration certificates, were obtained from the State Geodetic Administration.All models were analysed using the software program MATCH-T, and the same parameters were used for each model.As input data (aerial photogrammetric surveys) relating to different cyclical surveys were recorded using different cameras, the images of Papuk and Nijemci are black-and-white, while the images of Biokovo are in colour.All were scanned in the same resolution (21 µm) using the same scanner (Intergraph, PhotoScan TD).The data (points on the ground surface) were collected semi-automatically at regular intervals using the stereoscopic method.These data were extended by adding watershed lines for the entire drainage area and sub-drainage areas, using manual methods.Altitude point grids (digital elevation models) that were used in subsequent analyses were obtained by interpolation.

Figure 2. Map of Republic of Croatia with research locations
Interpolation is a procedure in the scope of which the value of a function is determined based on two known values [41].The surface determined by regular or irregular distribution of points is defined as follows: Let us assume that some points (x i ,y i ,z i ) for i=1,2,…,n, are distributed over an area.The task is to find the function z=f(x,y) which assumes given values at given points, while logical values are assumed within a given area for all other points (x,y) [42].
The data are gathered by combining three interpolation methods: linear method, least-squares method, and finiteelement method.The basic idea behind the linear interpolation method is to use several first degree polynomials instead of a single higher degree polynomial [43].Influence of pit removal algorithms on surface runoff simulation At each subinterval [x k-1 ,x k ], the polynomial p k ( 6) is singularly defined, and is written relatively considering the initial point of the interval The interpolation polynomial (7) can be written as which is followed by coefficients (8) for each k k = 1,…n An appropriate polynomial of the order of m<(N-1) is obtained by interpolation via the least squares method for N points.
The polynomial obtained by the least squares method is identical to the form of y p (x)=pa, where p is the Nxm matrix.
The interpolation based on the finite-element method (FEM) is used for approximation by parts [44].

Description of geomorphometrically different areas
Biokovo is a coastal mountain massif with a slightly elevated coastal part, which is continued with a steep rocky zone with an undulating plateau at the mountain top.From the top, this plateau gently descends towards the hinterland.The presentday "pockmarked face" of Biokovo is formed of sinkholes, karren, panholes, caves, and pits.It is a highly complex underground hydrographical network [45] abounding in natural depressions where water is contained due to conditions created by both erosion and geological substratum.Papuk is the longest mountain in Slavonia, characterised by a highly varied relief.On the one side, one encounters relatively flat forms intersected by small watercourses while the other side features steep inaccessible mountain slopes.The drainage area is generally characterised by torrential an unstable watercourses that find their origin in steep foothills, while main watercourses meander in their valleys.The area of the Nijemci district is characterized by relief presenting negligible differences in elevation (maximum height differences amount to 32.5 m and this exclusively in the northern part of the district where slight elevations are encountered).
Although the Nijemci district is situated in a flatland area, a surface and underground watershed (formation of two drainage areas) has formed along the plateau in the northern part due to terrain morphology.However, as this watershed is not impermeable, both drainage areas form a single hydrological system.The area is densely intertwined with drainage canals for the basic and detailed drainage, in which most of the surface flows actually end.The data collected for all three selected areas were analysed using harmonized coordinate differences as shown in Table 1, and coverage data as shown in Table 2.
Considering the number of rows and columns, the total number of pixels in the tested areas amounts to 614400.

Results and discussion
Elevation point grids were converted for all three areas into the ASCII format, and stored into the open-code GIS software called SAGA, where regular network points with associated heights were once again obtained in the first step (5x5 metres).Depressions were removed by the carving method (DEEPEN DRAINAGE ROUTES) using the drainage route created in advance (SINK DRAINAGE ROUTE DETECTION) and the filling method (FILL SINK).2) In order to enable better differentiation between the data without depressions and the data with depressions, appropriate shaded representations (azimuth = 315°, declination = 45°, augmentation factor = 4) showing situation before (Figure 4) and after removal of depressions (Figure 5) -contour line representation, were made in the second step.Then typical areas with overlapping contour lines were separately presented (Figures 6 and 7).Representation of Sanja Šamanović, Damir Medak, Duška Kunštek topography in the narrow zone of canals leading from one depression to another, preserving at that proper elevations of all contour lines.The situation is quite different in case of the filling method that fills the contour lines from the elevation of 1220 m to the elevation of 1275 m and changes them into a flat area.
At the initial form of the DEM for Biokovo, as obtained by interpolation, the DEM data were corrected using the filling method and carving method, which represented the basis for hydrological simulation of flow along the surface of the terrain.This simulation resulted in hydrographical network, i.e. in the trajectories of flow along the surface of the terrain, along water tracks, and along smaller bodies of water.At first glance, the comparison of these two models of flow along the surface of the terrain, along water tracks, and along smaller bodies of water, does not show any significant difference in flow.However, taking into account the fact that Biokovo mountain is steep and highly diverse with regard to its relief, with a multitude of depressions (in this case, such depression occupy up to 1.13% of the area), the morphometry of the terrain changes considerably once the depressions are filled, and so the flow of water does not follow the lowest points on the terrain, but rather fills these points and crosses them linearly until reaching the following change in terrain level (Figure 8, red arrows).It can be seen that participation of actual depressions is greater in the area under study, and that the hydrographical image changes considerably.The shaded representation of the Biokovo area shows a highly diverse and steep terrain with a multitude of depressions.Figure 5 shows the situation after removal of depressions, with canals for evacuation of water from depressions (left side) and the depressions filled to the point of overflowing (right side).Already at first glance, the terrain morphometry greatly differs from that seen in case of the filling method (Figure 5).It can therefore be assumed that better results will be obtained by carving method in the Biokovo area.
In overlapping representation showing situation before and after pit removal, we can see at the parts selected after removal of depressions using two methods (Figue 6) that the carving method changes the Influence of pit removal algorithms on surface runoff simulation The carving method respects the occurrence of depressions, ditches and similar formations to a much greater extent than the filling method, and so it can be seen that the flow of water follows depressions and always chooses the area with a lower altitude, as is normally the case on the terrain (Figure 8, green arrows).The shaded representation of the Papuk area (Figure 9) reveals a slightly inclined terrain with very few depressions, with the terrain gently rising from south toward north.Shaded areas along which water is evacuated are shown in the model on which depressions have not been removed.The water runoff in the area characterized by continuous sloping terrain is directed from higher areas in the north to lower areas in the south, along the cuttings that are present in the original relief topography, preventing the occurrence of depressions.Figure 10, in which relief is presented by shading, shows very slight changes in relief topography when the carving method is used (left side).On the other hand, changes are considerable after removal of depressions using the filling method (right side).This representation is a good example of how one can be tricked by visualisation.
Although at first glance the change in morphometry is considerable, it can be seen that considerable filling (changes in topography) occurs only at southern parts of the area under study.The entire area under study was represented with contour lines so as to better interpret deficiencies of visualisation in the Papuk area.If we consider Figure 11 that shows relief after removal of depressions via the filling method, it can be observed that a part Sanja Šamanović, Damir Medak, Duška Kunštek of contour lines "leaves" the model, i.e. is not confined within the studied area.As such contour lines can not properly be interpreted by the program within the shading module, considerable changes in geomorphometry occur in the model with removed depressions, which do not exist on the real model with removed depressions.This theory is further supported by the fact that southern parts of the area contain contour lines confined within the model, in which no change in geomorphometry can be observed.Unlike the Biokovo area, if we leave out of consideration irregular interpretation of contour lines, these two presentations do not show with certainly which of the methods provides better results.
Overlapping presentations with contour lines before and after removal of depressions on the selected parts of the Papuk area (Figure 12) provide a better insight into the differences in geomorphometry in the studied area with confined contour lines.As can be assumed, the geomorphometry of the area in which contour lines are confined within the model changes very little, regardless of whether the carving method or the filling method is used.As there are practically no depressions in this area (in this example, they take up only 0.05% of the area), the filling method changes the geomorphometry in smaller areas by filling depressions, while the carving method changes geomorphometry by creating narrow torrential flows leading to depressions, but elevations of contour lines remained unchanged regardless of which of the two methods was used.It is difficult to define which of the two tested methods provides better results by means of this visualisation in the Papuk area where the number of depressions is small.This statement is also confirmed by results of hydrological simulations of flow along the surface of the terrain, along ditches, and along similar formations (Figure 13 and Figure 14).2) In fact, in both cases of DEM modification the overlapping of hydrographical network does not reveal significant differences in simulation, while flow trajectories are quite correspondent, i.e. they overlap in almost all cases.The shaded representation of the Nijemci area (Figure 15) shows a flatland terrain with a considerable number of drainage canals, and with some depressions.Some of the canals also act as depressions and so they have been taken into consideration in the calculation, i.e. in the module for removal of depressions.
It can be seen in Figure 16 showing the terrain after removal of depressions according to the two methods that the presentation of relief by shading according to the method of pit removal by carving (left Influence of pit removal algorithms on surface runoff simulation side) does show the change in topography.However, this change is much smaller than in the case of the filling method (right side) that fills great depressions and a part of drainage canals which evidently exist on the terrain and should certainly not be filled in the context of defining a complete and reliable hydrographical network.Although it is often difficult to define which depressions constitute natural terrain features , and which ones are in fact computing errors caused by inaccurate calculation of DEM, this differentiation is simpler in this sort of terrain.As to smaller depressions within drainage canals, we can not be sure whether they are natural terrain features or simply the DEM error.Nevertheless, drainage canals recognised by the program as depressions do exist on the terrain and as such should not be taken into account in the calculation for removal of depressions and so, after analysis of Figure 12, it is our belief that the carving method will provide better results in this sort of terrain, compared to the pit removal by filling method.
Overlapping views with contour lines showing selected parts of the Nijemci area, before and after removal of depressions, provide a better insight into differences in topography in the area under study, which is full of drainage canals, natural depressions, and depressions created by inadequate digital-model calculations.
The carving method changes the topography in a smaller narrow area of newly created canals, while the filling method changes the topography considerably by creating a completely new terrain configuration (contour line 84 m in elevation).By comparing hydrographical networks created for both DEM modification models, with simulation of flow in its original form, it can be seen that much caution should be observed in this particular example (Figure 18).In fact, none of the flow simulations according to the modified DEM respects the existing drainage canals and hence does not make the surrounding streams to flow into these canals.According to the filling method, depressions and drainage canals are filled based on the same methodology, by creating parallel ditches that actually do not exist (Figure 19, red arrows).On the other hand, in case of the carving method, about 10% of flow is forced to use drainage canals, while the remaining flow is operated according to the principle of flow into the first lower terrain level (Figure 19, green arrows).Therefore, in such a specific flatland area intersected with surface drainage canals, we would recommend that none of the studied DEM modification methods be used, so as not to change the flow image considerably, as that would result in inadequate interpretation.It should also be Sanja Šamanović, Damir Medak, Duška Kunštek noted that advantage could be given to carving method in case no drainage canals are situated in such type of flatland area.In such a scenario the carving method is preferred as the filling method assigns considerable areas to depressions, which would change the geomorphometry of the terrain considerably.The depression filling procedure increases the elevation (height) of cells in which depressions have been identified, and so greater height determinants in the observed zones result in greater positive deviations from the real altitude.Raising of such zones around depressions reduces the value of inclination for these very zones, which finally results in negative deviation of inclination (slope).
A particular negative influence of this effect can be seen in drainage areas situated in flatland areas and in lower lying flattened areas such as agricultural surfaces.In large flatland areas, such as those occurring in the Danube River drainage basin, DEMs with rougher surface resolution do not provide a good-quality hydrographical network in the domain of both flows: along the surface of the terrain and along higher-rank and lower-rank watercourses.In fact, what the width of surface flows (drainage canals included) in flatland areas is smaller than the pixel size in DEM, then it is impossible to reliably define the hydrographical network by any of the mentioned DEM algorithms.
The DEM data for all models, with removed depressions, were downloaded into the GIS software IDRISI.This software was used to create views showing the pixel difference before and after pit removal by individual methods.In addition, histograms showing the differences in metres on the abscissa and the number of pixels on the ordinate, were also created (Figures 20,21,and 22).Tables showing numerical differences between the two models are also presented.Calculation of the total number of changed pixels was made for each area under study (Table 3).In this calculation, value 1 stands for the number of changed pixels, and value 0 for the number of unchanged pixels, out of the total number of 614400 pixels, Influence of pit removal algorithms on surface runoff simulation following implementation of each pit removal method.To facilitate comparison, percentage rates were also calculated for both values.
In each model, the difference in altitude in metres was also calculated for each changed pixel, and the information obtained was used to calculate the total volume change (Table 4).Removal of depressions by the carving method changes the terrain in such a way that its height (altitude) is lowered, and so the values presented are negative.On the other had, removal of depressions by the filling method changes the terrain by increasing the height (altitude) of the terrain, and so the values presented are positive.Tabular values of histograms for Biokovo area, for carving method, show that up to one metre of difference in geomorphometry was registered on 2368 pixels out of the total of 6295 pixels (defined accuracy of digital model is 0.9 m).It can also be observed that greater errors are practically inexistent.In case of the filling method, it can be noted that the greatest number of pixels is burdened with smaller differences.Thus, the difference is up to 0.5 m in 5616 pixels, while the difference ranges from 0.5 to 1 m in 4577 pixels.
The reduction in the number of pixes with greater differences occurs gradually until the pixels with the greatest differences are reached.This distribution of heigh differences between the two newlycreated models shows that the model in which depressions were removed using the carving method changes the topography to a lesser extent (most pixels with differences are withing the allowable tolerance) when compared to the model in which depressions were removed using the filling method.Tabular values from histograms for the Papuk area show that slightly better results were obtained by the filling method.In this method, most pixels are situated within the difference range of up to 1 m (only 77 pixels exhibit the difference greater than one metre), while 2738 pixels out of the total of 4476 pixels for which difference in altitude was calculated are situated within the difference range of 10 cm only.Figure 22 and tabular values from histograms for the Nijemci area show that better results were obtained using the carving method (48408 pixels are within the difference range of 0.1 m), while in case of models where depressions were removed using the filling method the number of pixels for which altitudes were changed falls quite uniformly, but the values are greater.

Table 4. Total volume change after removal of depressions by carving method and filling method in m 3
Table 3 shows differences after removal of depressions by the carving method and filling method, as related to the model without removal of depressions.In case of Biokovo area, the filling method changes the digital model on many more pixels compared to the carving method (1.02 % as related to 11.60 % of changed pixels in the model), and the change in volume (Table 4) is much lower in case of the carving method compared to the filling method, which shows that the geomorphometry was changed much less by the model created after removal of depression by the carving method.Table 3 shows that differences between the two methods are not great for the Papuk area (0.77 % as related to 0.73 % of changed pixels in the model).As already indicated, the Papuk area does not have many depressions, and the terrain topography is not greatly  Sanja Šamanović, Damir Medak, Duška Kunštek change after removal of depressions by carving method or by filling method.Considering the number of changed pixels and comparing them with Figure 21, it can be confirmed that incorrect results are obtained by visualisation of the Papuk area (Figure 10, right side) on the model with depressions removed according to the filling method.Unlike the Biokovo area where the proportion of changed pixels was considerably lower for the carving method, the differences for the Papuk area are too small for defining a method that could be recommended.The calculation of change in volume (Table 4) shows that the volume change generated by the filling method is much lower compared to the volume change generated by the carving method, which shows that the removal of depressions using the filling method is more favourable for the Papuk area, i.e. its effect on geomorphometry is lower.
In the Nijemci area, the differences between models without removal of depressions and models with removal of depressions (Table 3) can be observed in a great number of pixels both for the carving method (11.57%) and for the filling method (36.63 %).This is due to spatial resolution of models and to small differences in altitude, which is why most typical features (drainage canals, dykes, etc.) are neglected in flatland areas.The most notable differences between the two depression removal methods, i.e. influences on morphometry of the relief, can be seen in the calculation of the change in volume (Table 4), where geomorphometry of digital model for the area of Nijemci is less changed when the removal of depressions was operated using the carving method.Based on the above considerations (calculation of the total change in pixels and volume), recommendations are given in Table 5 for the use of appropriate methods for the three areas under study.

Conclusion
Concerning the data from three selected test areas (Biokovo, Papuk and Nijemci) with the resolution of 5 x 5 metres, the dimensions of the areas under study were harmonized to 4800 x 3200 metres or 960 x 640 pixels, which enabled easier comparison of results.Although some conclusions were derived from geomorphometric analyses (model created by shading, representation via contour lines, overlapping), statistical analyses (calculation of change in the number of pixels and volume) define with certainty which pit removal method changes to a lesser extent the geomorphometry of a particular area.The area of Biokovo is covered with a great number of natural and artificial depressions and, in this setting, the carving method changes the topography in the narrow zone of canals leading from one depression to another keeping at that the altitudes of all contour lines.On the other hand, the filling method connects depressions and creates larger filled flat areas.In the Biokovo area, the carving method changes altitude on a much smaller number of pixels of the digital model, and the change in volume of the area of Biokovo is smaller for the model in which depressions have been removed by the carving method.These results, together with geomorphometric analysis results, show that the geomorphometry is less changed using the carving method for removal of depression in areas characterized by a considerable number of depressions and great differences in altitude -such as the Biokovo area.This conclusion has also been confirmed by simulation of flow along the surface of the terrain, where an advantage is given to the carving method as it takes into account the occurrence of depressions, ditches and similar formations on the ground, and this to a much greater extent compared to the filling method, where consequently the flow of water follows depressions and always chooses the area with a lower altitude, as is normally the case on the terrain.For the Papuk area, one has to be cautious when making conclusions based on geomorphometric analyses because "unconfined" contour lines at model boundaries prevent accurate calculation of the model by the filling method, thus creating an unrealistic and considerably altered model.However, even if the areas of inaccurate calculations are disregarded, geomorphometric analysis results can not unambiguously define an acceptable pit removal method.It can be seen from the view with contour lines that the filling method changed geomorphometry in smaller areas, while the carving method generated narrow drainage passages leading to depressions, but the altitudes of contour lines remained unchanged in both methods.Hydrological simulations of flow do not reveal significant differences between the two models and, from this aspect, neither of the two methods can be considered to be preferable.Unlike the Biokovo area where differences in the change of the number of pixels were significant, the hilly area of Papuk shows small differences in the change of altitudes, and it is only after volume change calculation that it can be concluded that the filling method provides better depression-removal results for the Papuk area.For the Nijemci area, depression removal by the carving method reveals the change in topography, but the change is much less significant than in the case of the filling method which fills large and nested depressions but also a part of the canals that certainly should not be filled as it creates a completely new configuration of the terrain.Although both methods exhibit a great number of pixels with change in altitude after removal of depressions, the differences in the use of these two methods are nevertheless considerable.Calculation of the total modified volume, together with geomorphometric analysis results, enable recommendations favouring the use of the carving method for the Nijemci area, which is intersected by a great number of canals and features small differences in altitude.However, hydrological simulation of flow does not favour any of these methodologies, because of the specificity of the area in which priority in forcing the flow must be given to drainage canals, which should be preferred to any other criteria of flow along the terrain.The total number of changed pixels is considerable for both methods and, despite a relative advantage of one method over the other, one should exercise caution when considering methods for simulating

Nijemci Carving
Influence of pit removal algorithms on surface runoff simulation flow in flat parts without using the data about the real position of the river flow, because river networks created by such algorithms may deviate from real on-site situation considerably.That is why it is necessary to integrate a priori knowledge about the real position of river flows so that hydrological simulations can be tested and corrected.
Removal of depressions from the digital model of altitudes is a precondition for geomorphometric analyses, and so the choice of an optimum pit removal method, dependant on geomorphometry of relief, is a key factor for creating a minimally modified digital model of heights (altitudes).As many GIS software programs, with hydrological simulations and flow calculations, plan the DEM modification by initially using the "SINK" command, the methodology to be used in DEM modification should be defined depending on the topography of the terrain, presence of actual depressions, and some specific features in the area.It is recommended to always use the carving methodology in areas where the terrain is steep, dynamic and full of depressions, ditches, etc.In more moderate-relief terrains with flatter slopes and with very few depressions, any method can in fact be used as the results are very close and similar, while in flatland areas with very small inclinations and minimum changes of level, an advantage can be given to the carving method if the area is not intertwined with an artificial surface drainage network, because in the filling method great areas are attributed to depressions, which greatly changes the terrain geomorphometry and the pattern of runoff from the drainage area.
The use of any of the proposed algorithms influences the change of actual directions of flow along the terrain.The change of DEM, and hence the change of simulation of the hydrographical network, depends on the number of specific features in the area under study (size of drainage area that is being analysed, its topography -steep terrain, flat-sloped terrain, scale in which measurements are made, and a number of other micro-location determinants that locally change the direction of flow along the surface and along the higher-rank and lower-rank watercourses.
In any case, the author recommends that an extra caution be exercised when using GIS tools in hydrological simulations, especially the ones for which the depression-filling algorithm is not given as a choice, but is automatically incorporated in the module.In other words, in such situations, it is necessary to compare results from as many as possible GIS modules for DEM modification, with historical data from the terrain, as this approach will enable us to avoid creation of false trajectories of flow.An alternative to this approach is to use new hybrid depression-filling algorithms that have in fact been derived from the knowledge about limitations of existing algorithms.

Figure 1 .
Figure 1.Example of pit removal methods

Figure 6 .
Figure 6.Selected area at Biokovo DEM represented by contour lines

Figure 8 .
Figure 8. Selected part of Biokovo DEM showing depressions and hydrographical network

Figure 12 .
Figure 12.Selected area of Papuk DEM presented with contour lines

Figure 14 .Figure 15 .
Figure 14.Selected part of Papuk DEM with depressions and hydrographical network

Figure 17 .
Figure 17.Selected part of Nijemci DEM, view with contour lines

Figure 19 .
Figure 19.Selected part of Nijemci DEM with depressions and hydrographical network

Figure 20 .Figure 21 .
Figure 20.Histograms showing differences between models with depressions removed by carving method (left side) and filling method (right side), as related to the model in which depressions were not removed (Biokovo area); differences in metres are shown on abscissa, and number of pixels is shown on ordinate

Figure 22 .
Figure 22.Histograms showing differences between models with depressions removed by carving method (left side) and filling method (right side), as related to the model in which depressions were not removed (Nijemci area); differences in metres are shown on abscissa, and number of pixels is shown on ordinate