1D and 2D Flow Routing on a Terrain

For a point $x\in\mathbb{M}$, a depression $\beta_{x}$ of $h_{<\ell}$ is said to be delimited by the point $x$ if $x$ lies on the boundary of $\beta$, which implies that $h(x)=\ell$. A depression $\beta_{1}$ is maximal if every depression $\beta_{2}\supset\beta_{1}$ contains strictly more sinks than $\beta_{1}$. A maximal depression that contains exactly one sink is called an elementary depression. Note that each maximal depression is delimited by a saddle, and a saddle that delimits more than one maximal depression is called a negative saddle. For a maximal depression $\beta$, let $\text{Sd}(\beta)$ denote the saddle delimiting $\beta$, and let $\text{Sk}(\beta)$ denote the set of sinks in $\beta$. The volume of a depression $\beta$ of $h_{<\ell}$ is

$\mathsf{T}$ can be computed in $O(n\log n)$ time (Carr et al., 2003), and it can be preprocessed in $O(n)$ additional time so that for a point $x\in\mathbb{R}^{2}$, $\text{Vol}(\beta_{x})$, the volume of the depression delimited by $x$ can be computed in $O(\log n)$ time (Carr et al., 2003). In the I/O model, $\mathsf{T}$ can be constructed using $O(\operatorname{Sort}(n))$ I/Os and as shown in (Arge and Revsbæk, 2009), the volume $\text{Vol}(\beta_{x})$ and the smallest maximal depression containing $x$ can be computed for all vertices $x$ using $O(\operatorname{Sort}(n))$ I/Os.

We now describe the flooding model, which is the same as in (Lowe and Agarwal, 2019), and define flow-rate functions.

The value of $\lambda(u,v,t)$ is, in general, based on relative heights of the downslope neighbors of $u$. If $u$ is not a negative saddle vertex, then $\lambda(u,v,t)$ remains the same for all $t$, so we will often drop $t$ and write $\lambda(u,v)$ to denote $\lambda(u,v,t)$. If $u$ is a negative saddle, then $\lambda(u,v)$ changes when one of its depressions fills up as no water flows for $u$ to that saddle; see below for further discussion.

Letting $\tau_{v}$ be the time at which a vertex $v$ becomes flooded, $\phi_{v}(t)$ and $\phi_{(u,v)}(t)$ for any $v\in\mathbb{M}$ are undefined for $t\geq\tau_{v}$. That is to say, when a vertex is flooded, the flow-rate function is undefined.

Let $v$ be a negative saddle delimiting depressions $\alpha$ and $\beta$. Until one of $\alpha$ or $\beta$ is filled, $\phi_{v}$ is defined using equation 3. Without loss of generality assume that depression $\alpha$ fills first, say at time $\tau_{\alpha}$, and water starts spilling from $\alpha$ to $\beta$ through $v$. The spill rate $S_{\alpha}$ specifies the rate at which water spills from $\alpha$ to $\beta$. It is tempting to simply add $S_{\alpha}$ to equation 3, but it double counts the amount of water that was flowing from $v$ to depression $\alpha$. For $t<\tau_{\alpha}$, $\phi_{v}$ is defined as in (3). For $t\geq\tau_{\alpha}$, $\phi_{v}$ is defined as follows,

Finally for $t\geq\tau_{\alpha}$ and for any $w\in\beta$,

We first give an overview of the algorithm. We sweep through the vertices of $\mathsf{T}$ in descending height order, maintaining the following values for each height $\ell$:

1. (1)

   the set of depressions $\alpha_{i}$ in the sublevel set $h_{<\ell}$;

2. (2)

   for each $\alpha_{i}$, maintain the fill rate $F_{\alpha_{i}}(t)$, and

3. (3)

   for each $\alpha_{i}$ and the set of edges $(u,v)\in\mathbb{M}$ with $h(v)<\ell\leq h(u)$ and $v\in\alpha_{i}$, denoted by $\mathbold{E}(\alpha_{i})$, maintain the flow-rate function $\phi_{(u,v)}(t)$.

Note that $F_{\alpha_{i}}$ depends only on the sinks contained in $\alpha_{i}$, so we compute fill-rate functions of only maximal depressions. At each non-saddle vertex $v$ we compute $\phi_{v}$ and $\phi_{(v,w)}$, for each edge $(v,w)\in\mathbb{M}$, in a straightforward manner using (2) and (3). If $v$ is a negative saddle delimiting two depressions $\alpha$ and $\beta$, we must account for any rain spilling from one depression to the other. We do so by first partitioning the edges into two sets, $\mathbold{E}(\alpha)$ and $\mathbold{E}(\beta)$. Then given the flow-rates on the edges in each set, we compute the fill-rates $F_{\alpha}$ and $F_{\beta}$. Given these fill-rates, we then determine which depression (if any) fills first. If at least one depression becomes full, assume that $\beta$ fills first without loss of generality. Then we can compute the spill-rate function $S_{\beta}(t)$. Finally, we compute $\phi_{(v,w)}$ for each edge $(v,w)\in\mathbb{M}$.

We will now describe the algorithm in detail. We begin by augmenting $\mathsf{T}$ with additional information which will be used when computing the fill-rates.

Given a rain distribution $R$, before performing the downward sweep, we first compute $R(\beta,\cdot)$ for each maximal depression $\beta$ as follows: first assign $R(v,\cdot)$ for each vertex with nonzero rainfall to the smallest maximal depression containing $v$, and then perform an upward sweep on the merge-tree $\mathsf{T}$, maintaining the sum of rainfall functions at each saddle vertex.

As we sweep top-down, for each depression $\alpha_{i}$ we store the flow-rates of the edges in $\mathbold{E}(\alpha_{i})$ in a data structure that supports the following three operations: (i) add $\phi_{(u,v)}(t)$ to the data structure; (ii) remove $\phi_{(u,v)}(t)$ from the data structure; (iii) and given a saddle vertex delimiting two depressions $\beta,\gamma\in\alpha_{i}$, partition the flow-rates into two sets $\mathbold{E}(\beta)$ and $\mathbold{E}(\gamma)$ and return the sum of flow-rates in each set.

Noting there are $O(\chi+k)$ values of $t$ at which the flow-rate functions can change (corresponding to the times at which $R$ changes and spill-times of tributaries of the depression in question) we build the data structure as follows. For each time $t_{i}$ at which the functions can change, maintain the values $\Delta_{i}$ for the flow-rate of each edge $(u,v)$ sorted according to the label $\ell(v)$. Additionally maintain the prefix sum of the values for each $t_{i}$. There will be $O(|\phi_{(u,v)}|)$ values to consider at each edge, so we can insert and remove flow-rates of edges in $O(|\phi_{(u,v)}|\log n)$-time.

To partition the edges into $\mathbold{E}(\beta)$ and $\mathbold{E}(\gamma)$, recall each vertex is labeled according to an in-order traversal of $\mathsf{T}$. Given two depressions $\beta$ and $\gamma$ delimited by a saddle vertex $v$ with label $\ell(v)$, all vertices contained in $\beta$ will have a label less than $\ell(v)$, and all vertices contained in $\gamma$ will have a label greater than $\ell(v)$. Thus, to partition the edges into the depressions they cross into we simply take the prefix of all edges with label less than $\ell(v)$ to be one set, and the remaining suffix to be the other set. Then given the partition, we can use the prefix sums to determine the sum of flow-rate functions into $\beta$ and $\gamma$ in $O(\chi+k)$ time.

With this data structure at hand, we now describe how to process each vertex $v$ as we encounter it in our sweep.

The first sum has already been computed in the upward sweep of the merge tree. However, to compute $\phi_{(v,w)}(t)$ for all values of $t$ we would need to know which depression fills first and when. So instead, we define the pseudo-flow-rate function, $\phi^{\prime}_{(v,w)}$, in the same manner that flow-rate for non-saddle vertices is computed as in (3). Using this, we compute a modified pseudo-fill-rate function

Before $\alpha$ or $\beta$ become full (i.e. $t<\min(\tau_{\alpha},\tau_{\beta})$) we have that $\phi^{\prime}_{(v,w)}(t)=\phi_{(v,w)}(t)$, which in turn implies that $F^{\prime}_{\alpha}(t)=F_{\alpha}(t)$ (resp. $F^{\prime}_{\beta}(t)=F_{\beta}(t)$) before $\alpha$ or $\beta$ becomes full. Given this, compute $\phi^{\prime}_{(v,w)}(t)$ for all edges from the saddle $v$ and add these flow-rates to the data structure. We then use the data structure to partition the edges and compute the pseudo-fill-rates of the two depressions, from which we can determine which (if any) depression fills first.

If neither $\alpha$ nor $\beta$ becomes full, we are done. If not, assume without loss of generality that $\alpha$ becomes full first and spills into $\beta$. Then given $F_{\alpha}(t)$ we compute $S_{\alpha}(t)$ and in turn compute $\phi_{v}(t)$ for $\tau_{\alpha}\leq t$ using (4) along with $F_{\beta}(t)$ using (6).

Computing $R(\beta,\cdot)$ for all maximal depressions can be done in $O(|R|+\rho k)$-time where $\rho$ is the number of sinks in the terrain and $k$ is the number of times the rain distribution changes. To compute the flow-rate functions of edges from all non-negative saddle vertices, along with the pseudo-fill-rate functions from negative saddles, we spend $O(|\phi|\log n)$-time computing the flow-rate functions and adding them to the data structure. At a negative saddle delimiting two depressions $\alpha,\beta$, if we partition the edges by walking from the front and back of the sorted list of edges we spend $O(\min(|\alpha|,|\beta|)$-time, where $|\alpha|$ is the number of vertices contained in $\alpha$. A simple recurrence shows that the total time spent partitioning edges at all negative saddles is $O(n\log n)$. Finally computing the fill and spill rates at a negative saddle takes time proportional to their complexity, which is $O(k+\chi)$. Noting that $|\phi|=O(n(\chi+k))$, we get the following.

In particular, if $k=O(n)$, we have that the terrain-flow query can be answered in $O(n^{2}\log n)$ time.

We modify the above algorithm so that its running time becomes $O(\rho k+|\phi|\log n)$ where $|\phi|$ is the total complexity of all non-zero flow-rate functions. If water flows across a small portion of $\Sigma$ then $|\phi|\ll n$ and the new algorithm will be faster. The key idea is we do not have to process vertices with zero flow-rate.

The algorithm begins in the same manner as the previous algorithm, computing $R(\beta,\cdot)$ for each maximal depression $\beta$.

However rather than sweeping over all vertices of $\mathsf{T}$, we instead maintain a priority queue of vertices with nonzero flow-rate with their heights as their priority, i.e., the larger the height, the higher the priority. Every vertex $v$ with nonzero flow-rate will either have rain falling directly on it, or have a path in $\mathcal{M}$ from a source vertex (i.e. either a vertex with rain falling directly on it, or a saddle from which water is spilling.) We initialize the priority queue with all vertices on which rain falls directly along with all negative saddles.

While the priority queue is not empty, we pop the highest priority vertex $v$ from the queue and compute the flow-rate function $\phi_{v}(t)$, as we have described in the previous algorithm. For each vertex $w$ adjacent to $v$ in $\mathcal{M}$ with $\phi_{(v,w)}>0$, add $w$ to the priority queue. For each edge, in addition to each value $\Delta_{i}$ representing the change in the flow-rate function at time $t_{i}$, we also store $V_{i}$, the integral of the flow-rate function from $t=0$ to $t=t_{i}$. These are also stored in sorted order, and we maintain the prefix sum of these values for each edge.

At negative saddles, to determine which depression fills first, when we partition the edges we also compute the sum of volumes at the last time step $t_{j}$ for the depressions $\alpha$ and $\beta$ delimited by the saddle vertex. If both have filled by this time, consider the sum at $t_{j-1}$ and so on until we find the latest time at which only one depression is filled. This will be the depression which fills first, and ensures the number of time steps considered will be $O(|S_{\beta}|)$. Then we compute the fill and spill-rates as we did in the prior algorithm.

Since we compute the flow and spill-rates as before, it suffices to prove that we consider each negative saddle with positive spill-rate and each vertex with positive flow-rate. We begin with all negative saddles in the priority queue, so the former holds trivially. For the latter, given a vertex $v$ with positive flow-rate, either rain falls directly on it, in which case we added $v$ when initializing the priority queue, or water flows to it from some higher vertex $u$. If $u$ is added to the priority queue, when it is processed $v$ will in turn be added and later processed. Since all non-source vertices with positive flow-rate have a path from source vertex, and all source vertices begin in the priority queue, we have that we will reach and process $v$.

The time spent processing non-saddle vertices $v$ and edges $(v,w)$ will be $O(|\phi_{v}(t)|)$ and $O(|\phi_{(v,w)}(t)|)$ respectively. So the total time processing these will be $O(|\phi|)$. At negative saddle vertices when we partition the edges, if we perform the search by walking from the front and back of the list we will spend a total of $O(|\phi|\log|\phi|)$ time partitioning edges at all negative saddles. Finally, at each negative saddle $v$ which is a source we spend $O(|S_{\beta}|)$-time computing which depression fills first. Since $|S_{\beta}|\leq|\phi_{v}|$, the total time processing saddle vertices is $O(\rho k+|\phi_{v}(t)|\log n)$.

We now describe our first I/O-efficient algorithm. In order to extend the internal terrain-flow query algorithm, we introduce the following notation: let $\mathbold{\hat{E}}(\alpha)\subseteq\mathbold{E}(\alpha)$ be the set of edges such that, for each edge $(u,v)\in\mathbold{\hat{E}}(\alpha)$, $\alpha$ is the smallest maximal depression containing $v$.

As in the previous section, we proceed by performing an upward sweep on the terrain followed by a downward sweep.

We now describe how each step of the sweep is performed. Let $v$ be the current vertex. First we remove the function $R(v,t)$ from $Q$. Furthermore, for each edge $(u,v)\in\Sigma$, we remove the flow-rate $\phi_{(u,v)}(t)$ from $Q$. Since the algorithms visits vertices in descending order of height, the functions must have been inserted into $Q$ at an earlier point and will be at the front of $Q$. Since each edge $(u,v)$ no longer crosses the sweep line, we update the depression flow-rate sums by subtracting $\phi_{(u,v)}(t)$ from the flow-rate sum of the smallest maximal depression containing $v$.

Let $\alpha_{i}$ be the smallest maximal depression containing $v$. If $v$ is a non-negative-saddle vertex, we compute $\phi_{v}$ and determine if or when $v$ becomes flooded, using $F_{\alpha_{i}}(t)$ and $\text{Vol}(\beta_{v})$ as described in Section 3.

If $v$ is a negative saddle delimiting two depressions $\alpha$ and $\beta$, we wish to compute whether $\alpha$ or $\beta$ becomes full and, if so, which fills first. In order to compute the flood times of the two depressions, we recall that the fill-rate of a depression $\alpha$ is equal to the sum of flow-rates across all the edges $(u,w)\in\mathbold{E}(\alpha)$ plus the amount of rain falling in $\alpha$. When we are processing the vertex $v$ we have that all edges crossing into $\alpha$ and $\beta$ respectively are also crossing the sweep line. We therefore compute the fill-rate of $\alpha$ ($\beta$ resp.) by summing the depression flow-rate sums for all maximal depressions $\alpha_{i}\subseteq\alpha$ ($\beta_{i}\subseteq\beta$ resp.):

The flood times and the flow-rate $\phi_{v}(t)$ is then computed from $F_{\alpha}(t)$ and $F_{\beta}(t)$ as described in Section 3.

Finally, we can propagate the flow-rates on the outgoing edges by pushing $\phi_{(v,w)}(t)=w_{(v,w)}\cdot\phi_{v}(t)$ to $Q$ for all vertices $w$ where $w_{(v,w)}>0$. Furthermore, we update the depression flow-rate sums with the flow-rates of the outgoing edges. Note that each flow-rate is added to only one flow-rate sum.

The flow-rate $\phi_{(v,w)}(t)$ of each edge $(v,w)$ is forwarded only once using the priority queue. Thus, we spend a total of $O(\operatorname{Sort}(|\phi|))$ I/Os and $O(|\phi|\log|\phi|)$ internal computation time processing all vertices in the sweep excluding computation of the fill-rates at negative saddles. In order to efficiently compute the fill-rates in internal memory, we store the depression flow-rate sums in a sorted list ordered by the in-order traversal index of the depressions. We then partition the sums as in (Lowe and Agarwal, 2019) and compute the fill-rates using $O(|\phi|\log|\phi|)$ internal computation time in total. Since the initial sorting of vertices and edges can be performed in $O(\operatorname{Sort}(n))$ I/Os and $O(n\log n)$ internal computation time, the algorithm uses a total of $O(\operatorname{Sort}(|\phi|))$ I/Os and $O(|\phi|\log|\phi|)$ internal computation time.

We now extend the algorithm to relax the assumption on the size of the internal memory from $\rho(\chi+k)=O(m)$ to $\rho=O(m)$, at the cost of a greater number of I/Os.

We use the same framework as described previously, however, we avoid storing the depression flow-rate sums and fill-rates in memory for each depression in the sublevel set $h_{<\ell}$; We instead the priority queue $Q$ to forward fill-rates as well as the edge flow-rates used to compute fill-rates at negative saddle vertices.

We now bound the number of edges forwarded using a similar recurrence as (Lowe and Agarwal, 2019); Let $|\alpha|$ denote the total number of edges with at least one vertex contained in the depression $\alpha$ and note that $|\alpha|\geq|\mathbold{E}(\alpha)|$. Letting $T(\beta_{v})$ be the total number of flow-rates summed to compute the fill-rates for all saddles contained in $\beta_{v}$, we have that

Noting that the set of edges crossing into the two maximal depressions $\alpha$ and $\beta$ are disjoint, it follows that $|a|+|\beta|\leq|\beta_{v}|$. Using this, the recurrence solves to $T(\beta_{v})=O(|\beta_{v}|\log|\beta_{v}|)$, and we thus need to forward only a total of $O(n\log n)$ flow-rate functions. Since the complexity of each flow-rate function is bounded by $O(\chi+k)$, we spend a total of $O(\operatorname{Sort}((\chi+k)n\log n))$ I/Os and $O((\chi+k)n\log^{2}(kn))$ internal computation forwarding edges and computing fill-rates.