Response to reviews on "Computing water flow through complex landscapes, part 1: Incorporating depressions in flow routing using FlowFill"

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 5 high-resolution, high-precision DEM coverage increases the likelihood that depressions are real-world features. To address this longstanding problem of emerging significance, we developed FlowFill, an algorithm that routes a prescribed amount of runoff across the surface in order to flood depressions, but only 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 10 initial runoff depths and report the resulting filled and unfilled depressions, the drainage network structure, and the required compute time. Typical FlowFill calculations take minutes to perform and permit more realistic analyses of water flow across landscapes :: For ::: the :::::: reach:: to ::::::::::::: watershed-scale :::::::: examples :::: that :: we :::: ran, ::::::: FlowFill :::::::: compute :::: times :::::: ranged ::::: from :::::::::::: approximately : 1 :: to ::: 30 ::::::: minutes, :::: with ::::::: compute ::::: times ::: per :::: cell :: of :::::: 0.0001 :: to ::::: 0.006 : s.

issue if sinks are rather a data artefact than true sinks. But DEM preprocessing techniques actually deal with and correct for these artefacts. If FlowFill is not designed for this, and line 1 on page 18 actually suggests this, than FlowFill is not really a DEM preprocessing technique. Rather, it is a highly simplified simulation tool that models water flow across landscapes without loss through infiltration or evaporation. 25 Response: I agree in part with this comment: FlowFill is a simplified simulation tool that models water flow across landscapes. It has several potential uses, one of which is filling depressions in a landscape under a certain amount of runoff, and is the focus of this paper. I also agree that using this method, sink removal is highly dependent on network topology and that issues may arise in cases where sinks are data artefacts. This is an issue inherent to the method: by assuming that all depressions 30 in the DEM are true landscape features, we ignore cases where depressions really are data errors. Even if spurious depressions are filled by FlowFill, they have 'used up' some of the water available on the landscape and this water will not be available further downstream. This is a limitation of FlowFill, but does not negate its usefulness for this application.
First of all, spurious depressions are likely to be filled with even low runoff values since they are often a result of small uncertainties (Lindsay & Creed, 2006). Secondly, differentiating real depressions from erroneous ones is challenging. Lindsay and Creed (2006) discuss a few ways of doing so. The most reliable method is ground inspection, but this is often impractical.
Other methods, e.g. modelling, require assumptions to be made about what constitutes a real depression. These methods may also be flagging a certain depression as real when it is not, or vice versa. However, it is important to note here that there is 5 nothing stopping a user from first using a different method to remove depressions that they believe to be spurious, and then using FlowFill to see how hydrologic connectivity may change under different amounts of starting runoff. A user may also choose to use the output of FlowFill in more complex ways than we have here: for example, if they prefer carving to remove depressions, they may run FlowFill to determine which of the depressions should be removed from their study area, and then use a carving algorithm that allows them to preserve they depressions that they wish to keep (those that FlowFill did not fill). 10 Other depression-filling algorithms, including the flood-fill, carving, and combined methods we have discussed, consider all depressions to be spurious. FlowFill represents the other end-member: a method that considers all depressions to be real. A perfect method would be able to reliably differentiate between the two; perfect data would only contain real depressions. Since we have neither, we inherently assume that FlowFill successfully fills spurious depressions while recognising that this may not always be true and that a user should be vigilant for any obvious errors. 15 Comment: I wonder whether the objective of determining the amount of water in each topographic sink given a runoff volume can be obtained more easily.
Response: It certainly can -at least, more easily with respect to a user's computational time. This is something that we are 20 still working on and hope to share with the community in parts 2 and 3 of our series. The method we are working on should be much more computationally efficient, but is significantly more complex to create and explain. Despite the fact that we are working toward a better way of dealing with this particular problem, we still thought it important to publish this discussion of FlowFill for a number of reasons: 1) It is important to highlight this issue in DEM preprocessing and move forward discussion about it. We hope that others 25 will also think more about this issue and potential solutions.
2) While FlowFill is not computationally efficient, it is an intuitive method that was our first step in dealing with this particular problem.
3) There are cases (in other applications) where FlowFill may be more useful than methods of the type you have suggested.
For example, since FlowFill actually routes water across the land surface, it can be used to see how flow directions and water 30 distribution on a landscape change through time. This was not the main focus of our discussion here and is closer to what you have mentioned as FlowFill being a simplified simulation tool that models water flow. FlowFill can be used as a cellular automaton with a variety of possibilities. A user can view outputs after a given amount of iterations, or at multiple intervals after a set number of iterations each time. This could provide a simplified simulation of how flood water moves across the landscape through time. A user could also view outputs from adjacent iterations to see how dominant flow directions change 35 5 through time. These or other cell-based applications would not be possible with approaches that develop flow topology based on depressions, such as the approach proposed by the reviewer.
We will mention these alternative applications in the final version of the paper.
Comment: I think that the paper would benefit from placing more emphasis of FlowFill being a modelling tool to study 5 hydrological connectivity of complex landscapes, rather than a preprocessing technique.
Response: Thank you for your comment: our intention was to focus on one specific application of FlowFill. This ties back to our answer to the question above, where we have mentioned some other possible applications of this tool. We will mention these in the final version of the paper, although we do think that the preprocessing application of FlowFill is a significant one. 10 Smaller specific line and figure comments have been dealt with as suggested, and we thank Dr. Schwanghart for his careful attention.

Changes to the manuscript
Major changes: 15 -Ran FlowFill on a small region at a variety of resolutions. Included new a table and figures comparing these resolutions, and a discussion of the impact of resampling on the results.
-Added a discussion of usefulness of FlowFill in simulation of water flow.
Minor changes: -Added an explanation of the new option to randomise flow direction in the case of ties, as opposed to using a preferential 20 direction.
-Added a mention of the new GRASS extension, r.flowfill.
-Mentioned the removal of man-made structures from initial topography.
-Mentioned the possibility of using FlowFill results in conjunction with carving or breaching.

20
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 downhillsloping topography. By building such a flow-routing surface, we selectively remove information about the complexity of the real landscape. Real hydrologic networks include both transport of water across the land surface and temporary storage of water in depressions. 25 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). 30
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, 200 Even coarse-resolution data on a global scale can resolve real internally-drained basins (e.g. Riddick et al., 2018;Wickert, 30 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' 35 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. ? and ? ::::::::::::::::::: 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 ar-5 eas 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).

10
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 ( Figure 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 that can be contained by the depression before it overflows into an adjoining pixel is left behind 15 (?) ::::::::::::::: (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 summarised in Figure 2. FlowFill produces a flow-routing surface through an unconditionally stable method that iteratively routes water from cell to cell across the domain ( Figure 3). The required inputs are a DEM, a user-selected starting runoff value, and a user-selected 20 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.
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 20 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 25 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 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 watermoving algorithm (Figures 3 and 2) 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 5 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 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 Section 4). 10 Flooded depressions have flat surfaces which can be problematic for flow routing, so flat areas were corrected using Rich-DEM 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 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. 15 4 Implementation

Results
We used FlowFill to fill depressions both at the Sangamon River reach (Figures 5 and 6) and in the Río Toro basin (Figures 7   and 8). We applied varying amounts of runoff to demonstrate differing levels of depression filling ( Figure 9) and hydrologic connectivity ( Figure 10). Both study sites contain persistent depressions that are unlikely to be permanently filled and connected 10 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 ( Figure 5), significantly lower than the 15 m runoff required for the Río Toro site ( Figure   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 15 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 ( Figure   9). We define both drainage integration and hydrologic connectivity based on Strahler stream order ( Figure 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 20 that can, if unfilled, significantly break up the drainage network. The deep runoff required to completely fill all depressions ( Figure 9) or produce higher-order drainage networks ( Figure 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 ( Figure 10 and Table 4).
The processing time for FlowFill varies depending on the selected starting runoff depth, the number of cells in the domain, 5 and the topographic structure of the site. Runtimes for our test cases varied from 0.97 minutes to 32.42 minutes (Table 2).
Due to the slight overfilling of some depressions in the outputs from FlowFill, a correction was performed, as discussed in Section 3. In the two sample study regions discussed in this paper, the volume of the adjustment for overfilling was insignificant, 10 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.
The high number of depressions seen may also be partially due to the finer resolution of this data relative to the Río Toro study site -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.

Limitations
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;Barne Secondly, it creates flat areas in DEMs, which present their own challenges for flow routing. Thirdly, the threshold value for |∆h max | is distinct in different landscapes, and as such is a user-defined parameter. Finally, if an input topography contains 30 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 Figure  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 appealling flow networks. Fortunately, tools to do 5 so efficiently already exist. To produce our drainage networks for the stream-order calculations (Figure 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 10 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 15 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: Figure 11) in order to create a maximum water depth moved vs. iteration curve like those shown in Figure 4. From this, 20 users may pick a threshold value for |∆h max |. 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).

25
All depression-filling algorithms can produce 'lakes' as artefacts where 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 Section 3. This common problem further motivates work to remove these artificial blockages from rivers in DEMs for flow routing (?) :::::::::::::::::: (Abdullah et al., 2012).

Conclusions
Common and efficient downslope flow-routing algorithms must be run across surfaces that properly represent a true surfacewater-potential surface. As modern DEM resolution and accuracy increases, 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 5 landscape.