Computing water flow through complex landscapes – Part 1: Incorporating depressions in flow routing using FlowFill
- 1Department of Earth Sciences, University of Minnesota, Minneapolis, MN, USA
- 2Saint Anthony Falls Laboratory, University of Minnesota, Minneapolis, MN, USA
Correspondence: Kerry L. Callaghan (email@example.com)
Calculating flow routing across a landscape is a routine process in geomorphology, hydrology, planetary science, and soil and water conservation. Flow-routing calculations often require a preprocessing step to remove depressions from a DEM to create a “flow-routing surface” that can host a continuous, integrated drainage network. However, real landscapes contain natural depressions that trap water. These are an important part of the hydrologic system and should be represented in flow-routing surfaces. Historically, depressions (or “pits”) in DEMs have been viewed as data errors, but the rapid expansion of high-resolution, high-precision DEM coverage increases the likelihood that depressions are real-world features. To address this long-standing problem of emerging significance, we developed FlowFill, an algorithm that routes a prescribed amount of runoff across the surface in order to flood depressions if enough water is available. This mass-conserving approach typically floods smaller depressions and those in wet areas, integrating drainage across them, while permitting internal drainage and disruptions to hydrologic connectivity. We present results from two sample study areas to which we apply a range of uniform initial runoff depths and report the resulting filled and unfilled depressions, the drainage network structure, and the required compute time. For the reach- to watershed-scale examples that we ran, FlowFill compute times ranged from approximately 1 to 30 min, with compute times per cell of 0.0001 to 0.006 s.
Flow routing based on digital elevation models (DEMs) determines the paths taken by surface water (absent human interventions) and its associated sediment and/or dissolved load. Flow-routing algorithms are applied across a broad range of fields, including hydrologic and geomorphic modelling, topographic analysis, planetary science, and palaeoclimate. They are a critical component of both hydrologic (Neal et al., 2011; Ng et al., 2018; Ray et al., 2016) and geomorphic (Adams et al., 2017; Coulthard et al., 2013) models, with the former including watershed-scale processes and flood risk. In the latter case, flow routing is often recomputed over time to simulate the feedback between evolving topography and drainage patterns (Hobley et al., 2017; Tucker et al., 2011). Flow-routing calculations and drainage network construction also form the basis for topographic analysis algorithms to automatically pick channel heads (Clubb et al., 2014; Passalacqua et al., 2010; Pelletier, 2013), segment watersheds into representative hydrological units (Czuba and Foufoula-Georgiou, 2014; Ng et al., 2018; Teng et al., 2017), and link river-channel form with rates of tectonic uplift (Duvall et al., 2004; Perron and Royden, 2013; Willgoose et al., 1991) or subsidence (Paola et al., 1992; Wickert and Schildgen, 2019). These same tools have been applied to understand valley networks on Mars (Luo and Stepinski, 2009; Molloy and Stepinski, 2007), the impacts of freshwater forcing on climate during the most recent deglaciation (Ivanovic et al., 2017, 2018; Riddick et al., 2018), and links between palaeo-drainage networks and modern economic and agricultural resources (Craddock et al., 2010). Flow-routing algorithms are thus vital to our understanding of landscapes, climates, and water resources.
Multiple methods exist to distribute flow across the landscape, and these range from simple approaches to find the path of steepest descent from one DEM pixel to another (O'Callaghan and Mark, 1984) to full shallow-water equation solvers (McGuire et al., 2013). Simple topographically driven flow-routing approaches are the most popular because they are quick to compute (e.g. Braun and Willett, 2013; Gallant and Wilson, 1996; Schwanghart and Scherler, 2014), agnostic to the amount of rainfall or runoff applied, and applicable over length scales from puddles (e.g. Chu et al., 2013) and small catchments (e.g. Ng et al., 2018; Teng et al., 2017) to continents (e.g. Coe, 2000; Riddick et al., 2018; Wickert, 2016). It is this type of flow-routing algorithm that we consider for the remainder of this paper.
Because these simple algorithms route flow down the topographic gradient, local enclosed depressions in the landscape present a problem. Flow cannot be routed across them, so they disconnect the hydrologic network. Prior to computing flow routing, a DEM may be preprocessed to create a flow-routing surface in which depressions are managed in order to reliably extract stream networks (Metz et al., 2011). This step often results in the removal of all depressions from the original DEM (e.g. Jenson and Domingue, 1988; Martz and Garbrecht, 1998, 1999; O'Callaghan and Mark, 1984; Soille, 2004; Cordonnier et al., 2018). Removing depressions produces a topographic surface over which a flow-routing calculation will produce a fully connected hydrologic network.
The outright removal of all depressions, enforcing integrated drainage, means that we are biasing our landscape or hydrologic analyses towards what we are more easily able to calculate: an integrated drainage network within a continuous downhill-sloping topography. By building such a flow-routing surface, we selectively remove information about the complexity of the real landscape. Real hydrologic networks include both the transport of water across the land surface and the temporary storage of water in depressions.
We have developed a tool, FlowFill (Callaghan, 2019), that permits flow-routing calculations across landscapes that may contain real depressions. To do so, FlowFill employs mass-conserving and hydrologically consistent depression filling, allowing a user-selected depth of runoff to be spread across the landscape and flood only those depressions that would be filled by an overland flow event of the chosen magnitude. This approach eschews the assumption of fully integrated drainage and can help to improve the fit between the computed hydrologic network and field conditions (Coe, 2000).
Most existing approaches to managing depressions in flow-routing calculations assume that these are data errors and should be removed (Lindsay and Creed, 2005). In a historical context, this was a reasonable assumption: many DEMs were constructed from sparse data, and small data errors may have been enough to produce depressions, especially in flat areas (O'Callaghan and Mark, 1984). By filling these depressions, researchers bypassed this technical challenge to begin analysing the structure of drainage networks and continue to make new discoveries based on this approach (e.g. Hooshyar et al., 2017; Seybold et al., 2017). However, all of these analyses rely on an assumption of full drainage connectivity. This assumption may be broken in regions of lakes and basins, such as formerly glaciated terrains (Lai and Anders, 2018) and regions in which tectonic deformation isolates individual internally drained basins between ranges (e.g. Ballato et al., 2017; Sobel et al., 2003). Such depressions occur on a subcontinental scale as well, especially in arid regions, and our lack of an approach to self-consistently include these in our flow-routing algorithms inhibits efforts to construct and analyse large-scale drainage patterns and their changes over time (e.g. Wickert, 2016).
Several methods have been developed to remove depressions from a DEM during the creation of a flow-routing surface or to otherwise connect drainage across the landscape. A popular choice is a flood-fill algorithm, in which the elevation of all cells in depressions is raised to the level of their outlets (e.g. Jenson and Domingue, 1988; Martz and Jong, 1988). Martz and Garbrecht (1998) presented an alternative approach to the flood fill, which they called “breaching”. The breaching approach lowers select cells at depression outlets, reducing the amount by which depression cells need to be raised. Soille et al. (2003) extended this concept with their carving method. Rather than raising cells inside depressions, surrounding cells are lowered in order to eliminate all depressions in the topography. Combined methods both raise and lower cell elevations to minimize the topographic difference between the original DEM and the flow-routing surface (Lindsay and Creed, 2005; Schwanghart and Scherler, 2017; Soille, 2004). Grimaldi et al. (2007) proposed a physically based method based on steady-state topography that adjusts the elevation of cells in a DEM to that of a continuous river long profile. Metz et al. (2011) sidestep the need for topographic adjustments by instead using a least-cost-path method to determine drainage paths; this allows water to flow uphill to escape depressions and is employed in the GRASS GIS “r.watershed” algorithm (Neteler et al., 2012).
Each of the methods mentioned above ignores all depressions in the DEM, either by removing them or by allowing water to flow across them, and discounts the significant hydrologic impact of real depressions. In addition to helping to control drainage pathways (Govers et al., 2000), depressions set the volume of water that can pond on the land surface (Abedini et al., 2006; Hansen et al., 1999), thereby enhancing infiltration and both slowing and reducing surface runoff (Darboux and Huang, 2005). For example, in the prairie wetland region of North America, natural depressions hydrologically disconnect a landscape unless a runoff event is large enough to fill and overtop them (Arnold, 2010; Shaw et al., 2012). Furthermore, the growing availability of high-resolution and high-precision topographic data makes it increasingly difficult to support the assumption that these depressions are errors that should be removed (Arnold, 2010; Li et al., 2011; Lindsay and Creed, 2006). Even coarse-resolution data on a global scale can resolve real internally drained basins (e.g. Riddick et al., 2018; Wickert, 2016). Thus, the degree to which the drainage is integrated and the flows are connected is not static but rather varies as a function of runoff. All of these arguments motivate an approach that allows depressions to be filled in a way that is hydrologically realistic.
The aim of FlowFill is to generate flow-routing surfaces with an amount of drainage integration (and hence hydrologic connectivity) that is appropriate for the amount of input runoff and the shape of the land surface. Here, we use “drainage integration” to refer to the degree to which streams are connected (via lakes) instead of terminating in depressions, and greater drainage integration leads to a greater degree of surface-water hydrologic connectivity. Our goal in generating these surfaces is not a new one: Martz and Garbrecht (1998) noted the then-unrealized importance of incorporating depression storage into derived drainage patterns. Appels et al. (2011) and Chu et al. (2013) investigated the role of microtopography in connecting small surface depressions (puddles), and Shaw et al. (2013) required that ponds be filled by rainfall in their contributing areas before they be allowed to spill over their boundaries and integrate into the remainder of the catchment. In our approach, we developed a cell-by-cell runoff-routing algorithm that fills depressions while conserving runoff volume in real landscapes. This open-source algorithm, FlowFill (Callaghan, 2019), can compute flow-routing surfaces across a wide range of landscapes and is applicable at a range of length scales (see Table 1).
We present an algorithm to create more realistic flow-routing surfaces by flooding depressions with mass-conserved surface runoff. Depressions that are small or have large catchments become completely filled, allowing flow to cross them, whereas larger depressions may be only partially filled and continue to be hydrologic sinks (Fig. 1). FlowFill works by applying a user-selected runoff depth across the landscape and moving water downslope. If a parcel of water encounters a depression, as much of that parcel as can be contained by the depression before it overflows into an adjoining pixel is left behind (Callaghan, 2019). This enables users to perform flow routing across a landscape whose level of hydrologic connectivity changes through time due to storms, seasonality, or changing climate. The process of creating this flow-routing surface is summarized in Fig. 2.
FlowFill produces a flow-routing surface through an unconditionally stable method that iteratively routes water from cell to cell across the domain (Fig. 3). The required inputs are a DEM, a user-selected starting runoff value, and a user-selected threshold for convergence. This threshold specifies the maximum amount of water that may be moved between cells between two adjacent iterations for that iteration to count towards the eventual completion of the FlowFill calculation.
In FlowFill, water moves downslope, moving water from each cell in the domain once per iteration. The downstream direction for water movement is defined as the steepest downslope direction using a D8 approach (i.e. through comparison between the elevation of a target cell and the elevation of the eight neighbours with which it shares either an edge or a corner). In cases in which two or more directions tie for being the steepest downslope direction, the user selects whether a preferential or a random direction is preferred. We provide this choice since the selection of a preferential direction may systematically impact the ultimate destination of the water, whereas a random direction will solve this problem to some extent but make the result nondeterministic. When a preferential direction is selected, we arbitrarily route water preferentially northwest, then west, southwest, south, and continue anticlockwise with the least preferred direction being north. These cases should be rare since elevations between the cells would have to be identical to several decimals. This process of moving water down the slope of the (topography + water) surface is repeated until a predefined criterion is met to indicate that the solution has converged: either a maximum number of iterations have been performed, or the maximum amount of water moved per iteration, hmax, has not changed by more than the user-defined threshold for 20 000 iterations (i) (Fig. 4). This threshold defines the maximum value of that will not reset the “exit_counter” that terminates the FlowFill calculation (Fig. 2), where
Once the iterative downslope movement of water has been completed, FlowFill ensures that lake surfaces are flat. Due to the iterative algorithm in FlowFill, small spurious depressions can remain following convergence. To correct for these, we search for any pits (cells with no downslope neighbours) in the preliminary result. If a pit has one or more neighbours that contain water, it should ultimately have received water from that neighbour had FlowFill been allowed to run for longer (which is computationally expensive), so we raise its water level to the level of the water-containing neighbour. Strictly speaking, this correction means that water mass has not been conserved. However, the change affected by this correction is small relative to the total water volume. The adjustment of the lake-level surface is thresholded to a maximum value at 1∕10000 of the supplied runoff. The total volume of this adjustment ranged from 0 % to 0.009 % of the total water stored on the landscape at our study sites.
Following this lake correction, FlowFill outputs three files. The first of two binary (32-bit floating point) files contains the flow-routing surface: topography with depressions filled or partially filled in accordance with the provided input runoff depth. The second contains only the depth of water that is lying on the landscape. The third file contains runtime messages in ASCII text format.
We implement FlowFill in parallel using Message-Passing-Interface-enabled (MPI-enabled) Fortran 90. This speeds calculations by splitting the domain into multiple horizontal bands with fringes that interact via the D8 flow-routing algorithm. Source code and compilation instructions are available at https://github.com/KCallaghan/FlowFill (last access: 15 August 2019) (Callaghan, 2019). Use is simplified through a provided text file for users to enter parameters and a run file. As an additional option, users can run FlowFill through a GRASS GIS extension, r.flowfill (Wickert, 2019).
The gradual cell-to-cell water redistribution within FlowFill, along with its asymptotic approach towards equilibrium due to its moving at most half of the head difference per iteration, can cause depressions to become “overfilled” when the water-moving algorithm (Figs. 2 and 3) terminates. This could happen when water has not been able to fully equilibrate over a depression, for example, when there is only a small path for water to escape a large area. These cases are not corrected inside FlowFill but can be corrected in an additional step through comparison with the outputs of a flood-fill algorithm applied to the same initial DEM. Flood-fill algorithms fill all depressions fully to the level of their outlets, so we correct for overfilling by (1) performing a flood fill using RichDEM's complete-depression-filling command (Barnes et al., 2014a, b; Barnes, 2016) and then (2) taking the minimum of the flood fill and the FlowFill outputs to produce the final result. This correction violates the conservation of water volume, but the size of the adjustment is very small relative to the total volume of water stored on the landscape (See Sect. 4).
Flooded depressions have flat surfaces which can be problematic for flow routing, so flat areas were corrected using RichDEM to impose a gradient on these (Barnes et al., 2014a; Barnes, 2016). The result is a completed flow-routing surface that retains depressions based on the conservation of water volume. In order to test FlowFill, we ran it with variable initial runoff depths on two landscapes. We then routed surface-water flow over the computed flow-routing surfaces and evaluated the degree of drainage integration at these locations.
4.1 Example data
We generated flow-routing surfaces using FlowFill from DEMs of two study regions. The first study region includes a reach of the Sangamon River in Illinois, USA, located at 39.97∘ N, 88.72∘ W. The low-relief plains left behind in this postglacial landscape contain closed depressions that may impact hydrologic connectivity as a function of runoff (Lai and Anders, 2018). We resampled the 2.5 ft (0.76 m) resolution lidar DEM to 15 m resolution for our analysis. At several locations, bridges and other man-made structures cross the river channel at this study site. These are clearly visible on the lidar and artificially elevate the topography, creating blockages to water flow. These were manually removed using GRASS GIS before running FlowFill by digitizing the problematic bridges, converting these to null cells, and then performing a bilinear interpolation to populate these cells with more realistic values. The second study region was the Río Toro basin, located mainly in Salta Province, Argentina, around 24.5∘ S, 65.8∘ W. Its steep topography was shaped primarily by fluvial processes but has also been impacted by mountain glaciation, landslide dams, and tectonically driven isolation of much of the basin from the foreland (though it has since re-incised) (Sobel et al., 2003; Tofelde et al., 2017; Trauth and Strecker, 1999). The DEM of this region was resampled to 120 m resolution from 12 m TanDEM-X data. The two landscapes differ in their topographic setting in terms of tectonics, glaciation, drainage integration, average slope, and spatial scale.
In these examples, we prescribed uniform runoff across each DEM in order to easily compare depression filling and hydrologic connectivity (i.e. drainage integration). We tested runoff inputs of 1 mm to 15 m in the Río Toro basin and 1 mm to 20 cm in the Sangamon River basin (Table 2). For comparison, a typical storm event in the Río Toro basin drops ∼1 cm of rain (Castino et al., 2017, the Supplement), which equals the minimum runoff value discussed here. This minimum amount of runoff was not capable of filling many of the depressions in the landscape, and therefore our calculations indicate that some significant segmentation of the Río Toro catchment remains regardless of the size of the rainfall event. The median daily rainfall on a day with rain near the Sangamon River is 3.3 mm, and the maximum recorded single-day rainfall is 99 cm (USGS 05590050: data from 1 October 2005 to 19 February 2019). This large range suggests that hydrologic connectivity in the Sangamon River basin depends on storm intensity.
We used FlowFill to fill depressions both at the Sangamon River reach (Figs. 5 and 6) and in the Río Toro basin (Figs. 7 and 8). We applied varying amounts of runoff to demonstrate differing levels of depression filling (Fig. 9) and hydrologic connectivity (Fig. 10). Both study sites contain persistent depressions that are unlikely to be permanently filled and connected via surface water, as well as smaller depressions that may be filled during modest rainfall–runoff events.
We varied the input runoff depth at the two sites, with a maximum value selected based on how much runoff was required to fill all depressions in each DEM, giving a result comparable to existing flood-fill algorithms. This maximum runoff depth was 0.2 m for the Sangamon River site (Fig. 5), significantly lower than the 15 m runoff required for the Río Toro site (Fig. 7). However, at the Río Toro site, most depressions were filled by significantly shallower runoff (0.1 m or less), and only a few large depressions persisted as we dramatically increased the initial runoff depth.
At both sites, deeper runoff fills more depressions, thus increasing hydrologic connectivity across the landscape (Fig. 9). We define both drainage integration and hydrologic connectivity based on Strahler stream order (Fig. 10 and Table 4), which we calculated using the “r.stream.order” extension to GRASS GIS (Jasiewicz and Metz, 2011). The effect of variable runoff on drainage integration is more prominent in the Sangamon River landscape due to its larger proportion of depressions that can, if unfilled, significantly break up the drainage network. The deep runoff required to completely fill all depressions (Fig. 9) or produce higher-order drainage networks (Fig. 10) implies that either heavy rainfall or a long period of rainfall – which may or may not be plausible, depending on the climate and the hydraulic conductivity of the substrate – is necessary for the landscape to become fully hydrologically connected. Both landscapes contain a few larger depressions that persist once most other depressions have been filled, though based on their locations, they have somewhat less importance in setting overall hydrologic connectivity, which saturates at runoff values below the maximum required to flood all depressions (Fig. 10 and Table 4).
The processing time for FlowFill varies depending on the selected starting runoff depth, the number of cells in the domain, and the topographic structure of the site. Runtimes for our test cases varied from 0.97 to 32.42 min (Table 2). We performed each calculation using eight processors on an Intel i7-5820K CPU (3.30 GHz) on a desktop computer running Ubuntu Linux with 64 GB DDR3 RAM and a solid-state hard drive.
Due to the slight overfilling of some depressions in the outputs from FlowFill, a correction was performed, as discussed in Sect. 3. In the two sample study regions discussed in this paper, the volume of the adjustment for overfilling was insignificant, ranging from 0.003 % to 0.29 % of the total volume of water stored on the landscape. Cases in which the supplied initial runoff was deeper tended to have a slightly higher proportion of overfilling.
In addition to these two study sites, we used FlowFill on a subset of the Sangamon study site at four different resolutions in order to assess how the resolution of the input data affects the results. The results of this analysis can be seen in Figs. 12 and 13. A small subset of the Sangamon study site was selected in order to keep runtimes manageable at higher resolutions. The resolutions selected were 0.762 m (2.5 ft, the original resolution at which we obtained these data), 3, 5, and 15 m to match the resolution used for the entire Sangamon study site.
The results show that resampling data to a different resolution has an impact both on the number and morphology of depressions in the unfilled DEMs and on the results obtained from FlowFill. When resampling to coarser resolutions at this site, the number of large depressions in the study area visually appears to increase, as seen in Fig. 12. However, Table 3 shows us that the total number of depressions actually drastically decreases. This is due to the abundance of small depressions in the higher-resolution data. Instead, the total area of depressions present increases in lower-resolution data: the total area of depressions present in the unfilled DEM at 15 m resolution is almost 70 % more than the area of depressions at 0.762 m resolution. This trend is also reflected in the intermediate 3 and 5 m resolutions. While these coarser resolutions have resulted in higher depression areas, the smoothing effect of resampling has also resulted in depressions becoming shallower, and hence total depression volumes are smaller at coarser resolutions. This effect is less consistent: 15 m resolution exhibits the lowest total depression volume and 0.762 m has the highest, but the 5 m resolution DEM has a higher depression volume than the 3 m resolution DEM. Systematic changes in the shape of the landscape occur as resolution changes, but not all of these changes relate linearly to the change in resolution. All of these differences in depression number and morphology at different resolutions are an important reminder that the results of any landscape study based on remotely sensed data such as this are limited by the accuracy of the input data, and any preprocessing steps may change this accuracy. The effects of different resolutions on depression storage are discussed in more detail by Abedini et al. (2006), and Dixon and Earls (2009) discuss the effects of different resolutions on watershed delineation and streamflow prediction.
We completed 12 runs of FlowFill using DEMs with each of the four resolutions and with starting runoff depths of 0.2, 0.01, and 0.005 m (Figs. 12 and 13; Table 3). Regardless of the resolution, 0.2 m of starting runoff filled all of the depressions, while lower amounts of runoff left some depressions unfilled. Overall patterns in depressions filled appear visually similar at all resolutions, with 15 m resolution showing the greatest difference from other resolutions. At 15 m resolution with 0.01 m of runoff, several of the larger depressions on the southern edge have been filled, while these remained unfilled at finer resolutions. Drainage patterns also follow similar patterns at different resolutions, with the dominant river channels visible at all four resolutions. Channel widths are inflated at coarser resolutions due to the larger cell size (Fig. 13). Flow-routing pathways also differ between the 15 m resolution DEMs and those at finer resolution. When depressions are fully filled in the finer-resolution DEMs, a channel in the northeast corner of the DEM flows off the southern edge of the map. At 15 m resolution, however, the head of this channel is diverted and flows off the eastern edge of the map.
Using data at the highest available resolution prevents the loss or distortion of the data and ensures that analyses are not introducing any additional errors due to downsampling. Unfortunately, runtimes for FlowFill for large datasets can become prohibitively long. This is a limitation of FlowFill, and more computationally efficient methods for dealing with depressions in flow-routing surfaces are needed. We begin to address this problem in a companion paper (Barnes et al., 2019).
Runtimes for this subset of the Sangamon site (Table 3) ranged from 2 s to 13 min for 15, 5, and 3 m resolutions with all input runoff depths, but they escalated drastically to over 33 h for 0.762 m resolution with 20 cm of initial water depth. The reason for this nonlinearity lies in the fact that FlowFill moves water from cell to cell (Fig. 3). Increasing the resolution increases the total number of cells that must be calculated and requires more iterations of cell-to-cell water exchange for the water to move the same real-world distance. Runtimes times per cell are given in Fig. 11.
The flow-routing surfaces created by FlowFill account for water stored in the landscape and disconnects in the drainage network. The importance of such an approach is apparent because depressions persist in flow-routing surfaces even when the prescribed initial runoff is deep. This indicates that preprocessing topographic data with algorithms that fill all depressions is likely to result in spuriously integrated drainage networks. We have demonstrated that this effect can occur in both high- and low-relief landscapes and that, in addition to correcting spurious depressions, true lake basins and swales must be taken into account. Spurious depressions can also be filled by runoff and will therefore also be corrected by FlowFill. Some of the available runoff will be used up in doing so. Because spurious depressions caused by data errors are likely to be small, these would be filled with even low amounts of runoff (Lindsay and Creed, 2006; O'Callaghan and Mark, 1984), though it is still possible that some depressions that remain unfilled are artefacts of data errors.
FlowFill provides users with a completed flow-routing surface; however, should a user prefer carving, breaching, or combined methods for depression removal, these can still be used in conjunction with FlowFill. The result obtained from FlowFill determines which depressions should be removed during the creation of the flow-routing surface and which should remain. Depressions that FlowFill has not completely filled can be masked out, while those which were completely filled can be selectively carved or breached. This allows a user to utilize their preferred depression removal method, while still being cognisant of the importance of retaining real-world depressions.
5.1 Hydrologic connectivity
We present FlowFill results at two locations with differing landscape characteristics (Table 1). The results include flow-routing surfaces with selectively filled depressions and the associated changes in the drainage integration of the landscapes. A visual inspection of stream networks in cases in which lower runoff values are used tells us that hydrologic connectivity is lower in these cases. This is quantitatively supported by the Strahler stream-order data shown in Table 4. More higher-order streams occur when deeper runoff is used to create the flow-routing surfaces, hence filling more depressions. Channels were extracted using the FlowAccumulation functionality in RichDEM (Barnes et al., 2014a; Barnes, 2016) with a threshold of 100 units of accumulation.
The impact of using different amounts of runoff within FlowFill on the hydrologic connectivity of the landscape was more apparent in the Sangamon River basin study site, where more depressions were present. The presence of these depressions significantly reduced connectivity between stream segments. It is likely that FlowFill will be most useful in cases in which the geologic and geomorphic history of a landscape produces a surface with many depressions, such as this postglacial landscape. The high number of depressions seen may also be partially due to the finer resolution of these data relative to the Río Toro study site – the density of depressions has been shown to relate to grid spacing as an inverse power law (Lindsay and Creed, 2005) – but this does not belie the finding that significant real depressions exist and impact hydrologic connectivity.
The ability to use FlowFill with varying user-selected starting runoff values makes it ideal for comparing network connectivity in wet vs. dry seasons (Fig. 10 and Table 4) or for analysing the effects of storms of different sizes. Shallow runoff inputs to FlowFill imitate real-world conditions with low amounts of rainfall (e.g. during the dry season). During these times, hydrologic connectivity is significantly reduced, and routing flow across a completely depression-filled landscape becomes unrealistic. Deep runoff inputs simulate wet seasons or flood conditions. Due to the associated greater hydrologic connectivity, more of the region contributes water to basin outlets.
5.2 Cellular-based modelling
FlowFill is a cellular automaton that models water flow across landscapes. While we designed FlowFill to fill depressions on digital elevation models in a way that conserves water mass, rather than to reproduce a physics-based transient flow response, we propose that its mechanism of moving flow between cells may be useful for modelling applications. The amount of water moved between cells at each iteration is gradient based, meaning that FlowFill can approximate real transient flow in the landscape as a result of both topography (body forces) and water depth (pressure forces). The outputs of FlowFill diverge significantly from a true flow solution, for example the backwater equation, in that a parcel of water in FlowFill can move at most one cell per iteration, regardless of the underlying slope. Furthermore, FlowFill cannot accommodate different roughness values that would modulate flow velocity in one region versus another. With these limitations in mind, it could still provide a useful approach for simulations with approximately constant roughness and in which differences in elevation are consistently less than the flow depth. Fortunately, such examples are common in geomorphology and include reduced-complexity approaches towards simulating the dynamics of braided rivers (Murray and Paola, 1997) and river deltas (Liang et al., 2015b, a). With these uses in mind, users can view intermediate (pre-equilibrium) result outputs from FlowFill after a set number of iterations or at frequent intervals. The compute time for an individual iteration of FlowFill ranged from 10−3 to 10−6 s for the regions discussed in this article.
While we have created a way to handle the problem of real-world depressions in a DEM, FlowFill does have some limitations. Firstly, FlowFill requires significantly more compute time than flood-fill methods that fully fill DEM sinks (Barnes et al., 2014a; Barnes, 2016; Schwanghart and Scherler, 2014). Secondly, it creates flat areas in DEMs, which present their own challenges for flow routing. Thirdly, the threshold value for is distinct in different landscapes and as such is a user-defined parameter. Finally, if an input topography contains three-dimensional structures such as bridges, FlowFill can cause water to artificially dam behind them.
FlowFill can be time-consuming to run, especially for large DEMs, high starting runoff depths, or relatively flat study sites. Runtimes for each of the results shown here are shown in Table 2, with the compute time per cell shown graphically in Fig. 11. This may make it an unappealing choice in cases with large study areas or very-high-resolution data. Therefore, while FlowFill can route runoff to create flow-routing surfaces, a more computationally efficient solution to the problem outlined in this paper will permit faster analyses of a wider range of DEMs.
Like some other depression-filling algorithms, FlowFill produces flat areas where it fills depressions. Post-processing these into a gentle slope may be required in order to create reasonable or visually appealing flow networks. Fortunately, tools to do so efficiently already exist. To produce our drainage networks for the stream-order calculations (Fig. 10), we used RichDEM to impose a gradient on flat areas (Barnes et al., 2014a; Barnes, 2016) as the final step in constructing each flow-routing surface.
While a single criterion for convergence on a final flow-routing surface would be ideal, we were unable to find one and instead have left this as a user-selected parameter. We attempted to use the maximum amount of water moved between two cells in an iteration, the total amount of water moved in an iteration, the rate of change in the maximum amount of water moved between two cells in an iteration (averaged over various time windows), and the root mean square error of the linear regression that created the aforementioned slope. We also attempted each of these methods while normalizing for the initial amount of applied runoff. None of these approaches were able to collapse the response curves of flow over the landscape. However, we do observe that the maximum amount of flow between two cells in a single iteration asymptotes to a consistent value in each landscape. We therefore selected the exit criteria based on a time when the change in the maximum amount of flow between cells from one iteration to the next is below some small, user-selected threshold. It is necessary for a user to select this threshold value since the amount of noise after the plateau has been reached varies from one landscape to the next. As a result, we suggest that users who want to test multiple initial runoff depths first run FlowFill with a modest amount of runoff (to speed compute time: Fig. 11) in order to create a maximum water depth moved vs. iteration curve like those shown in Fig. 4. From this, users may pick a threshold value for . Selecting a value that is too large will cause FlowFill to exit before reaching this plateau, while a value that is too small will cause FlowFill to continue running for its maximum limit of 1 000 000 iterations, resulting in a long compute time. Based on our two study areas, suitable threshold amounts are landscape specific but appear to be agnostic to the amount of runoff selected for a given landscape. Following this approach, we chose thresholds of 0.01 mm (Sangamon) and 1 mm (Río Toro).
All depression-filling algorithms can produce “lakes” as artefacts. In these cases, the two-dimensional topography does not represent efficient three-dimensional flow paths – such as flow under bridges or through culverts (Lindsay and Dhun, 2015; Passalacqua et al., 2012). FlowFill is especially sensitive to these, as their damming effect can also create bottlenecks that significantly increase the number of iterations required to evacuate the water behind them, even when a narrow flow path exists to bypass them. Even after convergence, these areas often require the additional step to reduce overfilling discussed in Sect. 3. This common problem further motivates work to remove these artificial blockages from rivers in DEMs for flow routing (Abdullah et al., 2012).
Common and efficient downslope flow-routing algorithms must be run across surfaces that properly represent a true surface-water potential surface. As modern DEM resolution and accuracy increase, this requires that DEM depressions be appropriately filled. We have developed an algorithm, FlowFill, that fills only those depressions on a landscape which would become filled under reasonable runoff conditions. This allows for the existence of real depressions and hydrologic disconnects in the landscape.
By adding more realistic surface-water hydrology to flow routing, FlowFill's ability goes beyond that of static flood-fill algorithms and enables scientists to examine dynamic hydrologic connectivity. While FlowFill effectively solves this important problem, its long runtimes for larger datasets can make its use inconvenient. Future advances towards a more computationally efficient methodology will aid in the longer-term goal of linking real-world data with algorithms that harness the emerging power of far-reaching and high-resolution topographic data.
The latest version of FlowFill is available at https://github.com/KCallaghan/FlowFill (last access: 15 August 2019). The most recent version at the time of publication, v1.0.0 (Callaghan, 2019), is archived at https://doi.org/10.5281/zenodo.3358110. The GRASS extension r.flowfill (Wickert, 2019) is available at https://github.com/OSGeo/grass-addons/tree/master/grass7/raster/r.flowfill (last access: 15 August 2019).
The supplement related to this article is available online at: https://doi.org/10.5194/esurf-7-737-2019-supplement.
KLC developed the algorithm behind FlowFill, wrote its source code, ran the tests, and was the primary author of the paper. ADW provided guidance in the parallelization, suggested the test sites, and co-wrote the paper.
The authors declare that they have no conflict of interest.
The Deutsches Zentrum für Luft- und Raumfahrt (DLR) provided 12 m TanDEM-X DEM coverage of the Río Toro catchment via proposal DEM_GEOL1915 awarded to Taylor Schildgen, Andrew Wickert, Stefanie Tofelde, and Mitch D'Arcy. The authors acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota (http://www.msi.umn.edu, last access: 15 August 2019) for providing resources that contributed to the research results reported within this paper. Jingtao Lai and Alison Anders provided a copy of their Sangamon River DEM. Constructive and insightful reviews by Stuart Grieve and Wolfgang Schwanghart helped us to improve the paper content and organization. Kerry L. Callaghan was supported by the University of Minnesota Department of Earth Sciences HE Wright Footsteps Award and Junior F Hayden Fellowship, as well as by start-up funds awarded to Andrew D. Wickert by the University of Minnesota.
This research has been supported by the National Science Foundation (grant no. 1903606).
This paper was edited by Giulia Sofia and reviewed by Wolfgang Schwanghart and Stuart Grieve.
Abdullah, A. F., Vojinovic, Z., Price, R. K., and Aziz, N. A. A.: Improved methodology for processing raw LiDAR data to support urban flood modelling – accounting for elevated roads and bridges, J. Hydroinform., 14, 253–269, https://doi.org/10.2166/hydro.2011.009, 2012. a
Adams, J. M., Gasparini, N. M., Hobley, D. E. J., Tucker, G. E., Hutton, E. W. H., Nudurupati, S. S., and Istanbulluoglu, E.: The Landlab v1.0 OverlandFlow component: a Python tool for computing shallow-water flow across watersheds, Geosci. Model Dev., 10, 1645–1663, https://doi.org/10.5194/gmd-10-1645-2017, 2017. a
Appels, W. M., Bogaart, P. W., and van der Zee, S. E.: Influence of spatial variations of microtopography and infiltration on surface runoff and field scale hydrological connectivity, Adv. Water Resour., 34, 303–313, https://doi.org/10.1016/j.advwatres.2010.12.003, 2011. a
Arnold, N.: A new approach for dealing with depressions in digital elevation models when calculating flow accumulation values, Prog. Phys. Geog., 34, 781–809, https://doi.org/10.1177/0309133310384542, 2010. a, b
Ballato, P., Cifelli, F., Heidarzadeh, G., Ghassemi, M. R., Wickert, A. D., Hassanzadeh, J., Dupont-Nivet, G., Balling, P., Sudo, M., Zeilinger, G., Schmitt, A. K., Mattei, M., and Strecker, M. R.: Tectono-sedimentary evolution of the northern Iranian Plateau: insights from middle–late Miocene foreland-basin deposits, Basin Res., 29, 417–446, https://doi.org/10.1111/bre.12180, 2017. a
Barnes, R.: Parallel Priority-Flood depression filling for trillion cell digital elevation models on desktops or clusters, Comput. Geosci., 96, 56–68, https://doi.org/10.1016/j.cageo.2016.07.001, 2016. a, b, c, d, e, f, g, h, i
Barnes, R., Lehman, C., and Mulla, D.: An efficient assignment of drainage direction over flat surfaces in raster digital elevation models, Comput. Geosci., 62, 128–135, https://doi.org/10.1016/j.cageo.2013.01.009, 2014a. a, b, c, d, e, f, g, h, i
Barnes, R., Lehman, C., and Mulla, D.: Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models, Comput. Geosci., 62, 117–127, https://doi.org/10.1016/j.cageo.2013.04.024, 2014b. a
Barnes, R., Callaghan, K. L., and Wickert, A. D.: Computing water flow through complex landscapes – Part 2: Finding hierarchies in depressions and morphological segmentations, Earth Surface Dynamics Discussions, 1–19, https://doi.org/10.5194/esurf-2019-34, 2019. a
Braun, J. and Willett, S. D.: A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution, Geomorphology, 180–181, 170–179, https://doi.org/10.1016/j.geomorph.2012.10.008, 2013. a
Castino, F., Bookhagen, B., and Strecker, M. R.: Rainfall variability and trends of the past six decades (1950–2014) in the subtropical NW Argentine Andes, Clim. Dynam., 48, 1049–1067, https://doi.org/10.1007/s00382-016-3127-2, 2017. a
Chu, X., Yang, J., Chi, Y., and Zhang, J.: Dynamic puddle delineation and modeling of puddle-to-puddle filling-spilling-merging-splitting overland flow processes, Water Resour. Res., 49, 3825–3829, https://doi.org/10.1002/wrcr.20286, 2013. a, b
Clubb, F. J., Mudd, S. M., Milodowski, D. T., Hurst, M. D., and Slater, L. J.: Objective extraction of channel heads from high-resolution topographic data, Water Resour. Res., 50, 4283–4304, https://doi.org/10.1002/2013WR015167, 2014. a
Coe, M. T.: Modeling terrestrial hydrological systems at the continental scale: Testing the accuracy of an atmospheric GCM, J. Climate, 13, 686–704, https://doi.org/10.1175/1520-0442(2000)013<0686:MTHSAT>2.0.CO;2, 2000. a, b
Cordonnier, G., Bovy, B., and Braun, J.: A versatile, linear complexity algorithm for flow routing in topographies with depressions, Earth Surf. Dynam., 7, 549–562, https://doi.org/10.5194/esurf-7-549-2019, 2019. a
Coulthard, T. J., Neal, J. C., Bates, P. D., Ramirez, J., de Almeida, G. A., and Hancock, G. R.: Integrating the LISFLOOD-FP 2D hydrodynamic model with the CAESAR model: Implications for modelling landscape evolution, Earth Surf. Proc. Land., 38, 1897–1906, https://doi.org/10.1002/esp.3478, 2013. a
Craddock, R. A., Hutchinson, M. F., and Stein, J. A.: Topographic data reveal a buried fluvial landscape in the Simpson desert, Australia, Aust. J. Earth Sci., 57, 141–149, https://doi.org/10.1080/08120090903416278, 2010. a
Czuba, J. A. and Foufoula-Georgiou, E.: A network-based framework for identifying potential synchronizations and amplifications of sediment delivery in river basins, Water Resour. Res., 50, 3826–3851, https://doi.org/10.1002/2013WR014227, 2014. a
Duvall, A., Kirby, E., and Burbank, D.: Tectonic and lithologic controls on bedrock channel profiles and processes in coastal California, J. Geophys. Res., 109, F03002, https://doi.org/10.1029/2003jf000086, 2004. a
Grimaldi, S., Nardi, F., Benedetto, F. D., Istanbulluoglu, E., and Bras, R. L.: A physically-based method for removing pits in digital elevation models, Adv. Water Resour., 30, 2151–2158, https://doi.org/10.1016/j.advwatres.2006.11.016, 2007. a
Hansen, B., Schjonning, P., and Sibbesen, E.: Roughness indices for estimation of depression storage capacity of tilled soil surfaces, Soil Till. Res., 52, 103–111, 1999. a
Hobley, D. E. J., Adams, J. M., Nudurupati, S. S., Hutton, E. W. H., Gasparini, N. M., Istanbulluoglu, E., and Tucker, G. E.: Creative computing with Landlab: an open-source toolkit for building, coupling, and exploring two-dimensional numerical models of Earth-surface dynamics, Earth Surf. Dynam., 5, 21–46, https://doi.org/10.5194/esurf-5-21-2017, 2017. a
Ivanovic, R. F., Gregoire, L. J., Wickert, A. D., Valdes, P. J., and Burke, A.: Collapse of the North American ice saddle 14 500 years ago caused widespread cooling and reduced ocean overturning circulation, Geophys. Res. Lett., 44, 383–392, https://doi.org/10.1002/2016GL071849, 2017. a
Ivanovic, R. F., Gregoire, L. J., Burke, A., Wickert, A. D., Valdes, P. J., Ng, H. C., Robinson, L. F., McManus, J. F., Mitrovica, J. X., Lee, L., and Dentith, J. E.: Acceleration of northern ice sheet melt induces AMOC slowdown and northern cooling in simulations of the early last deglaciation, Paleoceanography and Paleoclimatology, 33, 807–824, https://doi.org/10.1029/2017pa003308, 2018. a
Lai, J. and Anders, A. M.: Modeled Postglacial Landscape Evolution at the Southern Margin of the Laurentide Ice Sheet: Hydrological Connection of Uplands Controls the Pace and Style of Fluvial Network Expansion, J. Geophys. Res.-Earth, 123, 967–984, https://doi.org/10.1029/2017JF004509, 2018. a, b
Li, S., MacMillan, R., Lobb, D. A., McConkey, B. G., Moulin, A., and Fraser, W. R.: Lidar DEM error analyses and topographic depression identification in a hummocky landscape in the prairie region of Canada, Geomorphology, 129, 263–275, https://doi.org/10.1016/j.geomorph.2011.02.020, 2011. a
Liang, M., Geleynse, N., Edmonds, D. A., and Passalacqua, P.: A reduced-complexity model for river delta formation – Part 2: Assessment of the flow routing scheme, Earth Surf. Dynam., 3, 87–104, https://doi.org/10.5194/esurf-3-87-2015, 2015a. a
Liang, M., Voller, V. R., and Paola, C.: A reduced-complexity model for river delta formation – Part 1: Modeling deltas with channel dynamics, Earth Surf. Dynam., 3, 67–86, https://doi.org/10.5194/esurf-3-67-2015, 2015b. a
Lindsay, J. B. and Creed, I. F.: Sensitivity of Digital Landscapes to Artifact Depressions in Remotely-sensed DEMs, Photogramm. Eng., 71, 1029–1036, https://doi.org/10.14358/pers.71.9.1029, 2005. a, b, c
Martz, L. W. and Garbrecht, J.: The treatment of flat areas and depressions in automated drainage analysis of raster digital elevation models, Hydrol. Process., 12, 843–855, https://doi.org/10.1002/(SICI)1099-1085(199805)12:6<843::AID-HYP658>3.0.CO;2-R, 1998. a, b, c
Martz, L. W. and Garbrecht, J.: An outlet breaching algorithm for the treatment of closed depressions in a raster DEM, Comput. Geosci., 25, 835–844, https://doi.org/10.1016/S0098-3004(99)00018-7, 1999. a
McGuire, L. A., Pelletier, J. D., Gómez, J. A., and Nearing, M. A.: Controls on the spacing and geometry of rill networks on hillslopes: Rain splash detachment, initial hillslope roughness, and the competition between fluvial and colluvial transport, J. Geophys. Res.-Earth, 118, 241–256, https://doi.org/10.1002/jgrf.20028, 2013. a
Metz, M., Mitasova, H., and Harmon, R. S.: Efficient extraction of drainage networks from massive, radar-based elevation models with least cost path search, Hydrol. Earth Syst. Sci., 15, 667–678, https://doi.org/10.5194/hess-15-667-2011, 2011. a, b, c
Murray, A. B. and Paola, C.: Properties of a cellular braided-stream model, Earth Surf. Proc. Land., 22, 1001–1025, https://doi.org/10.1002/(SICI)1096-9837(199711)22:11<1001::AID-ESP798>3.0.CO;2-O, 1997. a
Neal, J., Schumann, G., Fewtrell, T., Budimir, M., Bates, P., and Mason, D.: Evaluating a new LISFLOOD-FP formulation with data from the summer 2007 floods in Tewkesbury, UK, J. Flood Risk Manag., 4, 88–95, https://doi.org/10.1111/j.1753-318x.2011.01093.x, 2011. a
Ng, G.-H. C., Wickert, A. D., Somers, L. D., Saberi, L., Cronkite-Ratcliff, C., Niswonger, R. G., and McKenzie, J. M.: GSFLOW–GRASS v1.0.0: GIS-enabled hydrologic modeling of coupled groundwater–surface-water systems, Geosci. Model Dev., 11, 4755–4777, https://doi.org/10.5194/gmd-11-4755-2018, 2018. a, b, c
O'Callaghan, J. F. and Mark, D. M.: The extraction of drainage networks from digital elevation data, Comput. Vision Graph., 28, 323–344, https://doi.org/10.1016/S0734-189X(84)80011-0, 1984. a, b, c, d
Paola, C., Heller, P., and Angevine, C.: The large-scale dynamics of grain-size variation in alluvial basins, 1: Theory, Basin Res., 4, 73–90, 1992. a
Passalacqua, P., Do Trung, T., Foufoula-Georgiou, E., Sapiro, G., and Dietrich, W. E.: A geometric framework for channel network extraction from lidar: Nonlinear diffusion and geodesic paths, J. Geophys. Res., 115, 1–18, https://doi.org/10.1029/2009jf001254, 2010. a
Passalacqua, P., Belmont, P., and Foufoula-Georgiou, E.: Automatic geomorphic feature extraction from lidar in flat and engineered landscapes, Water Resour. Res., 48, 1–18, https://doi.org/10.1029/2011WR010958, 2012. a
Pelletier, J. D.: A robust, two-parameter method for the extraction of drainage networks from high-resolution digital elevation models (DEMs): Evaluation using synthetic and real-world DEMs, Water Resour. Res., 49, 75–89, https://doi.org/10.1029/2012WR012452, 2013. a
Ray, R., Beighley, R., and Yoon, Y.: Integrating Runoff Generation and Flow Routing in Susquehanna River Basin to Characterize Key Hydrologic Processes Contributing to Maximum Annual Flood Events, J. Hydrol. Eng., 21, https://doi.org/10.1061/(ASCE)HE.1943-5584.0001389, 2016. a
Riddick, T., Brovkin, V., Hagemann, S., and Mikolajewicz, U.: Dynamic hydrological discharge modelling for coupled climate model simulations of the last glacial cycle: the MPI-DynamicHD model version 3.0, Geosci. Model Dev., 11, 4291–4316, https://doi.org/10.5194/gmd-11-4291-2018, 2018. a, b, c
Schwanghart, W. and Scherler, D.: Short Communication: TopoToolbox 2 – MATLAB-based software for topographic analysis and modeling in Earth surface sciences, Earth Surf. Dynam., 2, 1–7, https://doi.org/10.5194/esurf-2-1-2014, 2014. a, b
Schwanghart, W. and Scherler, D.: Bumps in river profiles: uncertainty assessment and smoothing using quantile regression techniques, Earth Surf. Dynam., 5, 821–839, https://doi.org/10.5194/esurf-5-821-2017, 2017. a
Shaw, D. A., Vanderkamp, G., Conly, F. M., Pietroniro, A., and Martz, L.: The Fill-Spill Hydrology of Prairie Wetland Complexes during Drought and Deluge, Hydrol. Process., 26, 3147–3156, https://doi.org/10.1002/hyp.8390, 2012. a
Sobel, E. R., Hilley, G. E., and Strecker, M. R.: Formation of internally drained contractional basins by aridity-limited bedrock incision, J. Geophys. Res.-Earth, 108, 1–23, https://doi.org/10.1029/2002jb001883, 2003. a, b
Strahler, A. N.: Quantitative Analysis of Watershed Geomorphology, Transactions of the American Geophysical Union, Transactions, American Geophysical Union, 38, 913–920, 1957. a
Teng, F., Huang, W., Cai, Y., Zheng, C., and Zou, S.: Application of hydrological model PRMS to simulate daily rainfall runoff in Zamask-Yingluoxia subbasin of the Heihe River Basin, Water, 9, 769, https://doi.org/10.3390/w9100769, 2017. a, b
Tofelde, S., Schildgen, T. F., Savi, S., Pingel, H., Wickert, A. D., Bookhagen, B., Wittmann, H., Alonso, R. N., Cottle, J., and Strecker, M. R.: 100 kyr fluvial cut-and-fill terrace cycles since the Middle Pleistocene in the southern Central Andes, NW Argentina, Earth Planet. Sc. Lett., 473, 141–153, https://doi.org/10.1016/j.epsl.2017.06.001, 2017. a
Trauth, M. H. and Strecker, M. R.: Formation of landslide-dammed lakes during a wet period between 40,000 and 25,000 yr B.P. in northwestern Argentina, Palaeogeography, Palaeoclimatology, Palaeoecology, 153, 277–287, https://doi.org/10.1016/S0031-0182(99)00078-4, 1999. a
Tucker, G., Lancaster, S., Gasparini, N., and Bras, R.: The Channel-Hillslope Integrated Landscape Development Model (CHILD), in: Landscape Erosion and Evolution Modeling, edited by: Harmon, R. and Doe, W., chap. 12, 349–388, https://doi.org/10.1007/978-1-4615-0575-4_12, 2011. a
Wickert, A. D.: Reconstruction of North American drainage basins and river discharge since the Last Glacial Maximum, Earth Surf. Dynam., 4, 831–869, https://doi.org/10.5194/esurf-4-831-2016, 2016. a, b, c