An effective depression filling algorithm for DEM-based 2-D surface flow modelling

The surface runoff process in fluvial/pluvial flood modelling is often simulated employing a two-dimensional (2-D) diffusive wave approximation described by grid based digital elevation models (DEMs). However, this approach may cause potential problems when using the 2-D surface flow model which exchanges flows through adjacent cells, with conventional sink removal algorithms which also allow for flow exchange along diagonal directions, due to the existence of artificial depression in DEMs. In this paper, we propose an effective method for filling artificial depressions in DEM so that the problem can be addressed. We firstly analyse two types of depressions in DEMs and demonstrate the issues caused by the current depression filling algorithms using the surface flow simulations from the MIKE SHE model built for a medium-sized basin in Southeast England. The proposed depression-filling algorithm for 2-D overland flow modelling is applied and evaluated by comparing the simulated flows at the outlet of the catchment represented by DEMs at various resolutions (50 m, 100 m and 200 m). The results suggest that the existence of depressions in DEMs can substantially influence the overland flow estimation and the new depression filling algorithm is shown to be effective in tackling this issue based upon the comparison of simulations for sink-dominated and sink-free DEMs, especially in the areas with relatively flat topography.


Introduction
The Digital Elevation Model (DEM) is a digital representation of the ground surface topography or terrain which often is represented as a Cartesian grid, a triangulated irregular network (TIN) or as contour-based flow nets (Moretti and Orlandini, 2008).The regular DEM grid is convenient for calculation and conversion because of its relatively simple structure and can, therefore, be widely applied to topographic analysis and terrain visualisation.In hydrological applications, DEM data can also be used for catchment delineation, river network definition and catchment characteristics extraction (Gallant and Wilson, 2000).Many fully and semi-distributed hydrological model structures or parameterisations are based on the DEM, such as the overland flow module in MIKE SHE (Abbott et al., 1986;DHI, 2007), TOPMODEL (Beven et al., 1995), SWAT (Arnold et al., 1993), GBDM (Yu, 1987) models and some other applications as referred in Singh (1995).Since the DEM is an approximation to terrain characteristics, the accuracy of the topography and the related hydrological applications will depend on the quality of the DEMs, especially the characteristics such as slope analysis, river network density, flow path and the topographic index.
DEMs have to be made "hydrologically correct" before being used in hydrological models.Most DEMs contain numerous topographic depressions which are defined as areas without an outlet and often referred to as sinks or pits (Zandbergen, 2006).In regular-grid DEMs, topographic depressions are recognised as an area having one or more contiguous cells that are lower than all of its surrounding Published by Copernicus Publications on behalf of the European Geosciences Union.cells.Unsurprisingly these depressions often create difficulty in determining flow directions as the flow cannot continue downstream until the depressions are filled (Jenson and Domingue, 1988).As such the so-called depression filling or removal has become an important procedure in hydrological models that utilise DEMs.
The removal of artificial depressions is, therefore, desired due to the widespread use of high resolution DEMs produced by technologies like LiDAR and IfSAR (Zhu, 2009).MacMillan et al. (2003) revealed that high resolution LiDAR-derived DEMs often contain a very large number of (mostly small) depressions because of greater surface roughness and finer resolutions.Most DEMs available today contain a substantial amount of depressions which represent the incorrect, spurious topographic errors (Martz and Garbrecht, 1998;Olivera et al., 2000;Planchon and Darboux, 2001;Soille, 2004;Lindsay and Creed, 2005;Zandbergen, 2006).Locating and removing depressions in DEMs, therefore, remain as a necessary very first step of hydrological analysis.
Currently, most depression filling algorithms are based on the 1-D single flow direction, for instance, Rho4/Rho8 (Fairfield and Leymarie, 1991), D8 (O'Callaghan and Mark, 1984), the Lea (1992) algorithm, D∞ (Tarboton, 1997).The algorithm developed by Jenson and Domingue (1988) is one of the widely used and has been implemented in many GIS and hydrological softwares, such as Arc Hydro Tools (Maidment, 2002), GRASS (Neteler and Mitasova, 2008), HEC GEO-HMS (USACE, 2002), TOPOZ (Garbrecht and Martz, 1997), etc., all of which have been proved to be quite effective in one-dimensional surface runoff models.Nevertheless, the conventional depression filling algorithms are computationally intensive and in particular for processing high resolution DEMs.In order to improve the efficiency of identifying and filling surface depressions in DEMs, Wang and Liu (2006) proposed a depression filling method which can simultaneously determine flow paths and spatial partition of watersheds with one pass of processing, assisted by a novel concept of spill elevation and the least-cost search for optimal flow paths.And the time complexity of this method is in O(N log N ), which is superior to any other depression filling algorithms (Wang and Liu, 2006).
Accurate flow prediction over complex topography is one of the main targets of fluvial/pluvial flood models which often adopt a one-dimensional (1-D) representation of channel flow linked with a 2-D description of flow over the floodplain, for example, TELEMAC-2D (Galland et al., 1991), FLO-2D (O'Brien et al., 1993), LISFLOOD-FP (Bates and De Roo, 2000), AOFD (Maksimović et al., 2009).All these 2-D surface flow modelling usually involve a diffusion-wave treatment, which is commonly based upon determining the magnitude of flow between any two adjacent cells according to the water surface elevation difference and the Manning equation (Lane, 1998;Wheater, 2002;Yu and Lane, 2006).The differences between the 1-D and 2-D dynamic runoff processes lie in the fact that the 2-D surface flow model implies drainage along four possible directions (D4), not only receiving flow from upstream in x, y directions, but also draining to the downstream in both x and y directions.Whereas the traditional 1-D overland flow path is determined by considering eight possible directions (D8), only receiving and draining the water in one direction.As a consequence, some non-depression DEM cells within D8's concept may unfortunately become depressions during the 2-D surface flow process, due to the inconsistency between D4 drainage used in 2-D surface runoff models and sink removal algorithms designed to work in combination with D8 flow direction algorithms.
For example, the elevation in x, y directions control the flow path, as illustrated in Fig. 1 in which the arrows represent the flow direction on the surface.The depression cells as those shaded in Fig. 1a have lower elevations than the surrounding cells.The water in these depressions does not start to flow until the water level reaches 381 m and 386 m, respectively, and in this study these depressions are categorised as Type A. By comparisons, the elevations of the shaded cells in Fig. 1b are lower than the cells draining along cardinal flow directions, but higher than the cells draining along diagonal flow directions, which are recognised as Type B depression.It is clear that for Type B (as in Fig. 1b), the cells with elevations 385 m, 364 m and 295 m starts to flow when the water level is increased to 406 m, 386 m and 340 m, respectively.
The Type A depression presented in Fig. 1a is quite common and can be conveniently removed by applying any existing depression filling algorithms mentioned above.However, the Type B depression formed by the DEM presented in Fig. 1b has two different scenarios.The two algorithms D8 and D4 for determining the overland flow path are depicted by Fig. 2a, b, respectively.In this instance, all the water drained to the DEM boundary and the D8 algorithm as shown in Fig. 2a produced no sinks.By comparison, there is a sink in Fig. 2b due to the fact that the algorithm (D4) does not take into account the diagonal flow paths.Therefore, the water in the grid with a 295 m elevation will not flow until the water level increases to 340 m.
Therefore, the Type B depression can still be removed by applying filling algorithms based on the D8 assumption.But, it is less likely to be identified and removed in 2-D surface flow modellings, which employ D4 algorithm.And as a consequence, it has greater impact on the surface flow using 2-D described computations, in terms of affecting the kinematic wave celerity directly and hydraulic diffusivity indirectly (Orlandini and Rosso, 1998).However, there is a dearth of comprehensive studies examining the effects of depression removal on the 2-D surface flow modelling.
One possible approach as suggested by Zhu (2009) is to extend the river network to those Type B sinks by adding some arbitrary tributaries able to drain the ponded water back to the river.This method was shown to be very efficient in terms of retrieve water from the sinks.However, adding tributary manually means it is impractical to apply this approach to the case in which many ponds scatter across the catchment.Furthermore, subjectivity and uncertainty are inevitably brought in so as to determine the paths and locations of these tributaries.Wang and Liu (2006) proposed an efficient 1-D depression filling algorithms for high resolution DEMs.The main advantage of this algorithm is that it can simultaneously determine flow paths and spatial partition of watersheds in one pass of processing by employing the innovative concepts of spill elevation and using least-cost search for optimal flow paths.It has been proved that their method outperforms the conventional filling algorithm in terms of the running time and effectiveness.However, similar to other 1-D depression filling approaches, it cannot be adapted to the DEM sinkremove process for 2-D overland flow modelling.Therefore, we propose a new method in this paper aiming to effectively fill the Type B depression for 2-D overland flow modelling by modifying and extending the fundamental computational concept of Wang and Liu (2006)'s 1-D filling algorithms.Further, a rural catchment, the Upper Medway Catchment (220 km 2 southeast of England) is taken as a case study on the impacts of the filling algorithms to overland flow modelling and whole catchment outflow modelling.The MIKE SHE model (DHI, 2007), in particular its overland flow module, is chosen as a typical 2-D surface flow model, to analyse and assess the performance of the new depression filling algorithm.

The 2-D overland flow calculation
As one of the 2-D mathematical surface models, the MIKE SHE model has been widely used across the world and is capable of adapting DEM with various resolutions.In this study, the MIKE SHE model is used to demonstrate the issues associated with artificial depression and to study the effectiveness of the new depression filling algorithms for 2-D surface modelling.It is expected that the result would be conveniently applicable to a number of other modelling systems that share the similar 2-D surface flow approach as in the MIKE SHE.The overland flow module in the MIKE SHE employs a simplified set of the Saint Venant equations: the diffusive wave approximation describing the water movement on the surface, solved by the finite difference method (DHI, 2007).
The Saint Venant equations are also called the onedimensional unsteady open channel flow shallow water equations.If the Cartesian (x, y) coordinates are applied in the horizontal plane, the flow depth on the ground surface can be denoted by h(x, y) and the flow velocity in the x-and y-directions is u(x, y) and v(x, y), respectively.Therefore, according to the conservation of mass, the net input i(x, y) added to the overland flow would be established as (DHI, 2007): and the momentum equations are: where S f and S O stand for the friction slopes and ground surface slope, respectively, in x-and y-directions, q is lateral inflow.
It is still numerically challenging to solve the twodimensional Saint Venant equations dynamically.Thus, the diffusive wave approximation is introduced in the MIKE SHE which neglects the momentum losses (last three terms of the momentum equations) due to the local and convective acceleration and lateral inflows perpendicular to the flow direction.The diffusive wave approximation in x-and y-directions are shown below, given the relationship that Z = Z g + h, where Z g is ground surface level.Hence, However, when the slope of the water surface profile is very shallow and the water movement is very slow, numerical problems will occur between the neighbouring cells and any  Hager, 2001) is employed for each friction slope and the diffusive wave approximation equations then become: The discharge per unit length along the cell boundary can be written as: (5) The Manning's number M, usually ranging from 10 m 1/3 s −1 (rough channels) to 100 m 1/3 s −1 (smooth channels), is used in the MIKE SHE to control the flow velocity.It is equivalent to the Strickler roughness coefficient, which is also defined as the inverse of the Manning's number.
The MIKE SHE overland flow module makes use of the finite difference technique to solve the Saint Venant equations.If the overland flow on the grid is considered as shown in Fig. 3, at time t having the length x, y and water depth of h(t), then the velocity terms of the finite difference form of Eq. ( 1) can be approximated as: where the subscripts W, E, N and S denote the directions of the flow quantity, e.g., x(uh) E is the flow volume across the eastern boundary.Therefore, the total water that flows into the square across its north, south, east and west boundaries is described as h could be defined by Eq. ( 7): The cross-boundary flow is schematically depicted in Fig. 4 where Z U and Z D stand for the higher and lower water levels, respectively, and correspondingly the water depths in the grid are denoted as h U and h D .The discharge between the grids can be calculated by Eq. ( 7) combined with the Eq. ( 5).Hence, where K is the Strickler coefficient and the water depth h u is the free water that flows into the next grid, which is equal to the actual water depth minus the detention storage, since the detention storage is the ponded water that is trapped in the shallow surface depression.The overland flow into the cell can be zero if the upstream depth h u is zero or there is no difference between two water levels Z U and Z D .

The depression filling algorithm
In order to tackle the depression issues that may exist in the 2-D surface modelling, a depression filling algorithm was developed based on Wang and Liu's ( 2006) methodology, consisting of two fundamental concepts: spill elevation and optimal spill path.Spill elevation is defined as the minimum elevation that the flow level in one cell needs to be raised in order to reach the catchment boundary and the path it takes was referred to as the optimal spill path.Therefore, the elevation of a grid cell does not need to be raised if it is high enough to establish a downslope flow path to an outlet at the edge of the DEM, and its spill elevation is the same as its original elevation value of the cell of DEM.Otherwise, we need to raise that cell from its original elevation to its spill elevation.
As water flows with the gravity force from a place with a high elevation to a place with a lower elevation, a downslope flow path consists of a series of connected cells to an outlet of the DEM has to be found and constructed without the depressions.In that case, the elevations of all the cells on the flow path are ensured to be non-decreasing from the outlet on the boundary to the interior cells, as the search of flow paths are always started from outlets.
As described in the previous section, the overland flow calculation process in the MIKE SHE model has multiple flow directions.Each non-depression cell not only receives flow from upstream in x, y directions, but also drains to the downstream in x, y directions.Therefore, the depression filling algorithm is comprised of four different single flow direction processes.If the flow in one particular cell can reach the catchment outlet without depression filling, then the elevation of all the cells on the path are in non-rising order, and vice versa.Hence, if there is no such pathway to connect the cells and the catchment outlet, the depression cells need to be raised sequentially to their spill elevation values to satisfy the condition that the elevation is non-decreasing.The spill elevation is used as the "cost" to search the optimal spill paths within the catchment, also known as the "least-cost search" algorithm.
This algorithm is represented in the C++ programme by employing the concept of the "priority queue" method, in which the top priority elements are listed in the front.When it is applied to the "least-cost search" algorithm, the top priority element is the least-cost path and the order of the queue follows the same rule.
Three member functions are involved in this method, PQueue.Push(), PQueue.Pop() and PQueue.Top(), which are used respectively to insert the element into the queue, delete the first element in the queue and then return to the first element of the queue.Table 1 shows the pseudo-code for this method: lines 1-5 initialise the priority queue, in which all DEM cells on the catchment boundary or river networks are regarded as potential watershed outlets.In line 2, the spill elevation values of the outflow cells are allocated by their corresponding original elevation values.In line 3, the boundary cells are inserted into the priority queue through the member function PQueue.push().In line 4, the processed outflow cells are marked.The While loop in lines 7-19 accomplishes the search and expansion of the trees of optimal paths.Line 8 obtains the node with the least cost (the lowest spill elevation) through the member function PQueue.Top() and this least-cost node is deleted from the queue through the member function PQueue.Pop() in line 9.This series of procedures is labelled in the matrix Mark in line 10 to show that the leastcost node was process.The For loop in lines 11-18 starts to search the optimal paths from the four neighbourhood cells of the least-cost node.If any one of those four cells has not yet been processed before, the spill elevation value of the cell will be denoted by the higher value between the spill elevation value of the least-cost node and its original elevation value.In line 16, all of the other non-processed cells are inserted into the priority queue through the member function PQueue.push().Then, the computation enters the second round until the priority queue is empty, so that all the DEM cells are processed in the end.
Additionally, Fig. 5a illustrates the DEM after this depression removal algorithm having been applied, using the same  DEM example showed in Fig. 2. The order of algorithm processing the DEM is indicated in Fig. 5b.It shows that the algorithm took the cells on the boundary as the priority and listed all the surrounding cells into the process tree, then the cell with value 340 was considered as the search centre due to the least cost principle and the algorithm iterated until all the possible flow paths were searched and all the cells were processed.

Case study
Regarding the performance of this algorithm introduced above, the rainfall-runoff simulation of the Upper Medway Catchment (located in the South of London, UK) was selected as a case study, assessing the suitability and effectiveness of the algorithm applied to the MIKE SHE 2-D surface hydrological model.In this study, various DEM with different resolutions were used for the overland module of the MIKE SHE.The comparisons were then made between the model simulations (both the overland flow and total flow through the outlet) using the raw DEM data, and that from the DEM processed with the filling algorithm.
In order to gain the insight of how the resolution of DEMs may influence, three DEM datasets with different resolutions were resampled from the original 10 m resolution DEM data which initially was derived from the Ordnance Survey Land-Form Profile DTM (1 : 10 000) provided by the UK National Academic Data Centre (EDINA).The resampling process uses the bilinear interpolation to produce a coarser resolution of the raster.This interpolation algorithm takes the elevation values of the four nearest input cell centres to determine the value on the output raster, which results in a smoother looking surface than that obtained using the nearest neighbour method.This resampling process is effectively a "low-pass" filter operation and to some extent introduces smoothing to the data field that helps remove high-frequency noise generated by the interpolation algorithm.
The catchment topography varies between 30 m and 240 m above mean sea level.The majority of the slope characteristic in the Upper Medway Catchment varies from 2 degrees to 8 degrees, which represents about 70 % of the whole catchment and suggests that the landscape of the Upper Medway Catchment is made up of small hills surrounding the flat, low-lying relief of the floodplain without much variation of elevation (see Fig. 6).
As the catchment area is around 220 km 2 , it is not feasible to use the 10 m resolution DEM as a topographic with respect to the computational cost and the model calculation capability.Therefore, the 10 m resolution DEM was resampled to a coarser raster at three different resolutions (50 m, 100 m and 200 m).The depression-filling algorithm then processed these up-scaled DEMs.All these DEMs (original and the ones after sink removal) were then analysed in ArcGIS before being put into the model.
Figure 7 shows the analysis of the cumulative distributions of DEM slopes as a fraction of the catchment area, which not only clearly indicates the considerable smoothing effect that has taken place as a consequence of resampling the elevation data (from 10 m to coarser resolutions), but also suggests that the depression filling algorithm using three different scales has not altered the general characteristics of the catchment in terms of the topographic properties.
Similarly, Table 2 shows the comparison of slopes of the DEM at three different resolutions before and after the depression filling algorithm was applied, in terms of the statistics of the slope in degrees.This did not reveal a clear difference between the DEMs before and after the sink removal, but indicated that the general slope of the catchment decreased during the DEM resampling as a result of the averaging process being applied.
However, the comparisons in Table 2 in terms of elevation show little difference between the DEMs before and after the sink removal, as well as in the resampling impact because the statistics of elevation can only reflect the general characteristics of terrain, rather than indicating the changes that may impact processes on each DEM cell.
The impact of applying the depression removal process on DEM elevation is evaluated by calculating the difference between the terrain elevation before and after the filling algorithm process.The results suggest that most of the elevation values in the catchment remained quite similar after the depression filling algorithm had been applied, and most of the elevation differences occurred near the existing river network, implying that the depression filling algorithm had made an effort to form synthetic river paths in order to drain the surface water to the river network.
More importantly, the statistics of depressions on the catchment shows that number of depressions in 50 m, 100 m and 200 m resolution DEMs detected by this filling algorithm were 3407, 1778 and 768, respectively, which made the corresponding area covered by the depressions to be 8.52 km 2 , 17.78 km 2 and 30.72 km 2 and the related proportion of whole catchment area to be 3.87 %, 8.08 % and 13.96 %.It implied that resampling the DEMs to coarser resolution decreased the number of pits, but the total area covered by sinks increased, which is consistent with the conclusion drawn by Grimaldi et al. (2007).Additionally, Table 3 indicates that the coarser DEM resolution produced more elevation differences and in  turn causes more differences in the volume of depression water on the surface.In order to test the impact of the depression algorithm on the surface flow simulation, firstly, the MIKE SHE model was configured in such a way that only the overland flow module was included without other input from the infiltration and groundwater interaction modules.Therefore, the total simulated flow from the model was from its overland flow output which is driven entirely by the terrain elevation and the Manning's number.All the DEMs were then compared in the context of the MIKE SHE hydrological modelling structures in terms of the simulation flow at the catchment outlet (see the location in Fig. 6).The corresponding overland flow simulations are shown in Figs.8-10.
Figure 11 shows the accumulated water depth on the surface during the model simulation.Given the same amount of precipitation, the accumulated surface water depth derived from the three original DEM inputs indicated that the proportions of surface detained water after the flow simulation were 17.68 %, 36.03 % and 55.84 %, respectively, suggesting that the resampled coarser raster resulted in more detained water on the catchment surface than the fine resolution DEM.In addition, more than half of the precipitation was retained on the surface when using the 200 m resolution DEM simulation, meaning that far less water drained to the river network when compared to the 50 m resolution results.However, after the depression-filling algorithm was applied, the processed DEM produced similar detained accumulated surface water in the model, the corresponding proportions decreasing to 14.12 %, 17.80 % and 12.62 %.This indicated that the depression filling algorithm had efficiently eliminated ponds affecting the overland flow in the MIKE SHE and the margins between the different DEM resolution simulated flows through the model were significantly reduced.
Figure 11 also shows that the algorithm-processed DEM with 100 m resolution produced similar accumulated surface water depth with the 50 m resolution raw DEM in the MIKE SHE model, implying that these two sets of DEM data can generate a similar amount of overland flow in the model.In this case, owing to the similarity of the accumulated surface water depth, all the DEMs with different resolutions after sink removal by the algorithm produced a similar amount of water on the surface.Thus, the depression filling  Secondly, the whole catchment rainfall-runoff simulation including subsaturated zone and saturated zone components in the MIKE SHE was carried out for further investigation into the efficiency of algorithm applied in this study.The water movement through the soil profile, along with the evapotranspiration is modelled by a simplified Two-Layer ET/UZ model (Yan and Smith, 1994), which is suitable for catchments that have a shallow groundwater table.The ground water flow is calculated using the linear reservoir method which can be regarded as the balance of the data availability of the geology, the complexity of the groundwater simulation and the benefit from the model simplicity.The whole model was set up using a grid size of 100 m × 100 m.The trial-and-error parameterisation was employed and the Nash-Sutcliffe efficiency (Nash and Sutcliffe, 1970) in calibration and validation was 0.93 and 0.91 (Zhu and Cluckie, 2012;Zhu, 2009).
Figure 12 to 14 shows the difference on the peak flow simulation at catchment outlet before and after the depressionfilling algorithm was applied.The simulations using all three DEM resolutions show an increase of overland flow due to the depression filling algorithm, in particular, in the case of using coarse DEM resolution in the model.Compared to the measured peak flow, Table 4 indicates that the difference between modelled peak flow and observed peak flow on all DEM resolutions in this study were reduced, after applying the depression filling algorithm.Additionally, similar to the previous overland-flow-only simulation results, the coarser DEM resolution, the more impact on the peak flow caused by the depression filling algorithm.Moreover, as shown in   4, the increase of total outflow in turn contributes to the peak flows during the simulations as also displayed in Figs. 13 and 14, for 100 m and 200 m DEM resolution, respectively.Overall, it can be seen that the whole model simulation can also benefit from the new algorithm of depression removal via its contribution to the better result in overland flow simulation.

Conclusions
In this study, two different types of depressions in DEM data were discussed in the context of distributed hydrological modelling.It shows that the Type B depression is unlikely eliminated by the traditional 1-D depression filling algorithm which may cause problems when used in 2-D surface flow models that make use of 2-D diffusive wave approximation to propagate the surface water into the river channel.A close investigation shows that the problems in fact arises from the inconsistency between D4 drainage used in 2-D surface runoff models and the sink removal algorithms designed to work in combination with D8 flow direction algorithms.
A new depression-filling method was then developed following Wang and Liu's (2006)   The effect of DEM resolution in this context was also investigated by using the 10 m DEM as a start and then gradually resampling it to 50 m, 100 m and 200 m resolution, respectively.The flow simulations using these DEMs at different resolution suggests that the DEM up-scaling process may result in more depressions and consequently a higher degree of precipitation being accumulated on the catchment surface rather than being drained into the river network.
In terms of modelling overall catchment response, it is important to note that once water from precipitation is in the channel network then it enters a much faster flow domain where the dynamic nature of flood response is generally determined.Conversely, water retained in artificial sinks across a DEM domain will be subject to a much slower catchment response.The outflow simulations of whole catchment by the MIKE SHE model show that the increased overland flow due to the application of the depression-filling procedure has considerable impact on the peak flow in particular and, therefore, improve the performance of the hydrological model.Again it implies that the depression filling algorithm is efficient and suitable for adoption in the 2-D surface flow models, not only because it produces more accurate surface flow, but because it also reduces the uncertainties introduced by the spatial resampling of the DEM data.Furthermore, it has been shown that the overland flow simulations with coarser DEM can benefit more by using the depression-filling method as the artificial depressions are more likely to appear and, therefore, have larger impacts on the accuracy of simulations.This suggests that the proposed depression-filling algorithm has a considerable potential for improving 2-D hydrological simulation over areas with insufficiently high resolution of DEM.
It is worth noting that the aim of this study was to examine the impact of depression removal on surface flows modelling using DEM and two-dimensional diffusive code in particular.And it was not meant to improve the description of surface runoff dynamics as a whole, but rather to highlight and correct for a specific problem.The hydrological model used in this study has been well calibrated and validated beforehand.Undoubtedly the performance of rainfall-runoff simulation depends on a wide variety of factors, e.g., the catchment hydrological characteristic, the model rainfall-runoff mechanism, the simulation period and initial conditions; this study only addresses the contribution of using a better designed depression removal method, whereas in future studies other important aspects such as model bias in whole catchment simulation needs to be taken into account.
Although a re-calibration is generally required when the grid resolution is changed, a multi-resolution (MR) validation test by Vázquez et al. (2002) indicates that using the effective parameters of a calibrated MIKE SHE model for different grid resolutions only differ marginally on the simulated streamflow.Our study is focused on the impact of depression removal on the whole catchment runoff simulation given the fact that the overland flow has been affected, instead of identifying the best grid resolution for this catchment and, as such, a single set of model parameters were adopted.Undoubtedly, however, a more comprehensive study including re-calibrations on different grid resolutions would be able to give further insights into the impact of resolution per se and that will be addressed in future work.
With respect to its computational efficiency, this algorithm inherits the advantages of Wang and Liu's (2006) algorithm compared with other traditional depression filling approaches.It directly computes a spill elevation value for each grid cell without prior delineation of the catchments of depressions and progressively builds the optimal flow paths and propagates spill elevation values from outlets to interior grid cells, due to the employment of the least-cost search algorithm.Consequently, it is able to produce a DEM with D. Zhu et al.: An effective depression filling algorithm fewer depressions and identify the locations and depth of the surface depressions with computational time complexity in O(N log N ), which can be easily implemented using objectoriented programming languages like C++ and Java.

Fig. 1 .
Fig. 1.Example of two types of DEM Depression.

Fig. 3 .
Fig. 3. Demonstration of flow on a square grid system in MIKE SHE model; source: DHI (2007).

1Fig. 5 .
Fig. 5. DEM with depression removal algorithm applied and the relevant processing order.

Fig. 6 .
Fig. 6.Contours and DEM for the upper medway catchment (unit: metre).The X and Y-axis stand for the easting and northing coordinates of catchment.

Fig. 8 .
Fig. 8. Surface flow simulations for the 50 m resolution DEM in MIKE SHE.

Fig. 9 .
Fig. 9. Surface flow simulations for the 100 m resolution DEM in MIKE SHE.

Fig. 10 .
Fig. 10.Surface flow simulations for the 200 m resolution DEM in MIKE SHE.

Fig. 12 .
Fig. 12. Whole catchment flow simulations using 50 m resolution DEM in MIKE SHE.

Fig. 13 .
Fig. 13.Whole catchment flow simulations using 100 m resolution DEM in MIKE SHE.

Fig. 14 .
Fig. 14.Whole catchment flow simulations using 200 m resolution DEM in MIKE SHE.

Table 1 .
Pseudo-code for algorithm implementation.

Table 4 .
Statistics of the peak flow in whole catchment simulations.