1D and 2D Flow Routing on a Terrain
11D and 2D Flow Routing on a Terrain
Aaron Lowe
Duke University
Svend C. Svendsen
Aarhus University
Pankaj K. Agarwal
Duke University
Lars Arge
Aarhus University
ABSTRACT
An important problem in terrain analysis is modeling how waterflows across a terrain creating floods by forming channels andfilling depressions. In this paper we study a number of flow-query related problems: Given a terrain Σ , represented as a triangulated xy -monotone surface with n vertices, a rain distribution R whichmay vary over time, determine how much water is flowing over agiven edge as a function of time. We develop internal-memory aswell as I/O-efficient algorithms for flow queries.This paper contains four main results:(i) We present an internal-memory algorithm that preprocesses Σ into a linear-size data structure in O ( n log n ) time that for a (possiblytime varying) rain distribution R can return the flow-rate functionsof all edges of Σ in O ( ρk + | ϕ | log n ) time, where ρ is the number ofsinks in Σ , k is the number of times the rain distribution changes,and | ϕ | is the total complexity of the flow-rate functions that havenon-zero values; | ϕ | = Θ ( n ( χ + k )) in the worst case, where χ is theheight of the merge tree of Σ and k is the number of times the raindistribution changes, but | ϕ | is much smaller in practice.(ii) We also present an I/O-efficient algorithm for preprocessing Σ into a linear-size data structure using O ( Sort ( n )) I/Os and O ( n log n ) internal computation time, so that for a rain distribution R , it cancompute the flow-rate function of all edges using O ( Sort (| ϕ |)) I/Osand O (| ϕ | log | ϕ |) internal computation time.(iii) Σ can be preprocessed in O ( n log n ) time, into a linear-sizedata structure so that for a given rain distribution R , the flow-ratefunction of an edge ( q , r ) ∈ Σ under the single-flow direction (SFD)model can be computed in O (| R | + | b q | k log n ) time, where | R | isthe number of vertices in R on which nonzero rain falls, and | b q | isthe number of tributaries of q in which rain directly falls in.(iv) We present an algorithm for computing the two-dimensionalchannel along which water flows using Manning’s equation; awidely used empirical equation that relates the flow-rate of waterin an open channel to the geometry of the channel along with theheight of water in the channel. Given the flow-rates along a path ofedges, the algorithm computes the height of water and boundaryalong the channel in O (| C | log | C |) time, where | C | is the numberof wetted faces at least partially covered by water in the channel. Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for components of this work owned by others than ACMmust be honored. Abstracting with credit is permitted. To copy otherwise, or republish,to post on servers or to redistribute to lists, requires prior specific permission and/or afee. Request permissions from [email protected].
SIGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA © 2020 Association for Computing Machinery.ACM ISBN 978-1-4503-8019-5/20/11...$15.00https://doi.org/10.1145/3397536.3422269
ACM Reference Format:
Aaron Lowe, Svend C. Svendsen, Pankaj K. Agarwal, and Lars Arge. 2020.1D and 2D Flow Routing on a Terrain. In
ACM, New York, NY, USA, 12 pages. https://doi.org/10.1145/3397536.3422269
An important problem in terrain analysis is modeling how waterflows across a terrain and creates floods by forming channels andfilling up depressions. The rate at which a depression fills up duringa rainfall depends not only on the shape of the depression and thesize of its watershed (i.e., the area of the terrain that contributeswater to the depression) but also on other depressions that arefilling up. Water falling on the watershed of a filled depressionflows to a neighboring depression effectively making the watershedof the latter larger and filling it up faster. Modeling how depressionsfill and how water spills into other depressions during a flash floodevent is therefore an important computational problem.Besides determining which areas of a terrain become flooded andwhen they become flooded, determining the 2D channels (rivers)along which water flows across the terrain is also important. Theflow queries we ask can also be used to solve related flood-riskqueries, and the algorithms developed will provide a simpler andfaster algorithm for a previously studied flood-risk queries. Weassume we are given a terrain Σ , represented as a triangulated xy -monotone surface with n vertices. As in earlier papers, [4, 13, 15]we assume that water flows along the edges of Σ . Two models ofwater flow along edges have been proposed: (i) a simple and morewidely used model called the single flow-direction (SFD) model inwhich water from a vertex flows along one of its downward edges,and (ii) a more accurate but more complex model called the multi-flow-direction (MFD) model in which water at a vertex splits andflows along all of its downward edges. There is also some earlierwork by Liu and Snoeyink [12] which allows water to flow in theinterior of edges. We consider both of these models.We study the following three problems in this paper: • Terrain-flow query : given a rain distribution (possibly vary-ing with time), compute as a function of time the flow rate(of water) for all edges of Σ . • Edge-flow query : given a rain distribution and a query edge e of Σ , compute the flow rate of e under the single flow-direction (SFD) model. •
2D flow network : Given a 1D flow network, represented asa set of edges along with their flow values, compute 2Dchannels along which water flows.Finally, as high-resolution terrain data sets are becoming easilyavailable, their size easily exceeds the size of main memory of a a r X i v : . [ c s . C G ] S e p IGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA Aaron Lowe, Svend C. Svendsen, Pankaj K. Agarwal, and Lars Arge standard computer, so movement of data between main memoryand external memory (such as disk) becomes the bottleneck incomputations.We use the
I/O-model with one disk by Aggarwal and Vitter [1],in which, the computer is equipped with a two-level memory hier-archy consisting of an internal memory , which is capable of holding m data items, and an external memory of unlimited size. All com-putation happens on data in internal memory. Data is transferredbetween internal and external memory in blocks of b consecu-tive data items. Such a transfer is referred to as an I/O-operation or an
I/O . The cost of an algorithm is the number of I/Os it per-forms. The number of I/Os required to read n items from disk isScan ( n ) = O ( n / b ) . The number of I/Os required to sort n items isSort ( n ) = Θ (cid:0) ( n / b ) log m / b ( n / b ) (cid:1) [1]. For all realistic values of n , m ,and b we have Scan ( n ) < Sort ( n ) ≪ n . Related work.
Due to its importance, the problem of modeling howwater flows on a terrain has been studied extensively, and manyapproaches have been taken to address this problem. One approachfocuses on accurately simulating fluid dynamics using non-linearpartial differential equations such as the Navier-Stokes equations.These equations have no closed form solutions and are usuallysolved using numerical methods. They often account for additionalfactors, such as the effects of different terrain types and drainagenetworks. While these models tend to be more accurate, they arecomputationally expensive and do not scale to large terrains. Batesand De Roo [6] developed one such model for simulating flooding ondigital elevation models (DEMs) using two flow models for differentregions of the terrain: the first handles flow within rivers and thesecond models flow of water as it spreads over floodplains. Whilethere has been some research into refining the representation ofchannels, such as Wood et al . [17], often the channel geometry isassumed to be a simple model (e.g. rectangular or trapezoidal.)Water-flow modeling on a terrain also has been studied in theGIS community. These approaches use simpler models, which arecomputationally efficient and suitable for large datasets. Althoughsome early work, e.g. [12] allowed water to flow in the interior offaces, recent work assumes that water flows along the edges of theterrain and the SFD and MFD models described above are used tomodel the water flow locally at a vertex, see e.g. [4, 13, 15] Arge etal. [5] described an I/O-efficient algorithm for the flow-accumulation problem in the SFD model, water falls uniformly on the terrainwhich asks how much water flows over each point in a terrainassuming the terrain has only one sink at infinity. Their algorithmperforms a total of O ( Sort ( n )) I/Os, where Sort ( n ) is the number ofI/Os required to sort n items. The flow accumulation model onlyprovides rough solution to flow modeling, since it assumes thateither the terrain does not have any local minima or that they havebeen filled in advance.In order to accurately model flow it is necessary to compute timesat which depressions fill and simulate how water spills from onedepression into others. Arge et al. [4] proposed the first I/O-efficientalgorithm that computes the fill times of all maximal depressionsin O ( Sort ( ρ ) log ( ρ / m )) I/Os, where ρ is the number of depressionsin the terrain and m is the size of the internal memory. If ρ = O ( m ) ,the algorithm can be simplified and requires only O ( Sort ( n )) I/Os. Arge et al. [2] described an I/O-efficient algorithm that computeswhich points on the terrain become flooded if a total volume ψ of rain falls according to a distribution R . It performs O ( Sort ( n ) + Scan ( χ · ρ )) I/Os, where χ is the height of the merge tree of theterrain. In the worst case χ ρ = Ω ( n ) , but it can be bounded to O ( Sort ( n )) under certain practical assumptions.Lowe et al. [13] presented efficient algorithms for several floodqueries on a terrain under the multiflow direction (MFD) model.They presented a O ( n log n ) -time algorithm to answer terrain floodqueries. They also showed that a terrain Σ can be preprocessed in O ( n log n + nρ ) -time into a data structure that can answer point-flood queries: Given a rain distribution R , a volume of rain ψ , anda point q ∈ Σ , determine whether q will be flooded. The querytime is O (| R | b q + b q ) time, where b q is the number of maximaldepressions that contain the query point q ; b q = Ω ( n ) in the worstcase, but in practice it is much smaller. Finally, they presented a O ( nb q + b ωq ) -time algorithm to determine when a query point q gets flooded, where ω is the exponent of fast matrix multiplication.To our knowledge, no I/O-efficient algorithms are known for theseflooding queries under the MFD model. Our results.
There are four main results in the paper.(i) We present a O ( n log n ) time algorithm for preprocessing Σ into a linear-size data structure for answering terrain-flow queries:for a rain distribution R , it can compute the flow rate of all edgesin O ( ρk + | ϕ | log | ϕ |) time, where | ϕ | is total complexity of nonzeroflow-rate functions, ρ is the number of sinks in Σ , and k is thenumber of times the rain distribution changes. In the worst case | ϕ | = Θ ( n ( χ + k )) , where χ , as above, is the height of the merge treeof Σ , but | ϕ | is much smaller in practice. An immediate corollary ofour result is that a flood-time query (i.e. given R and a point q ∈ Σ ,when will q be flooded) can be answered in the same time, whichis a significant improvement over the result in [13].(ii) We present two I/O-efficient algorithms for the terrain flowalgorithm. We first preprocess the terrain using O ( Sort ( n )) I/Osand O ( n log n ) internal computation time. The first algorithm that ρ ( χ + k ) = O ( m ) , where χ , k and ρ are as above, and considers aterrain-flow query in O ( Sort (| ϕ |)) I/Os and O (| ϕ | log | ϕ |) internalcomputation time. The second algorithm assumes ρ = O ( m ) andanswers a terrain-flow query in O ( Sort (( χ + k ) n log n )) I/Os and O (( χ + k ) n log ( kn )) internal computation time. We additionallynote that since the terrain-flow query is a more general problem,these algorithms also yield I/O-efficient algorithms for the terrain-flood and flood-time queries under the MFD model studied in [13].(iii) The terrain-flow query algorithm naturally can also be usedto perform edge-flow queries. Under SFD flow, we can build alinear sized data structure in O ( n log n ) time which given an edge e = ( q , r ) and rain distribution R computes ϕ e in O (| R | + b q k log n ) time, where | R | is the complexity of the rain distribution, b q is thenumber of tributaries of q in which rain is falling, and k is thenumber of times the rain distribution changes.(iv) We present an algorithm that given a 1D flow network,represented as a path of edges along with their flow values, iden-tifies the depth and width of water along the 2D flow network in O (| C | log | C |) time, where | C | is the number of “wetted” faces atleast partially covered by the water in the channel. We do so by D and 2D Flow Routing on a Terrain SIGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA utilizing Manning’s equation [14], a widely used empirical formularelating flow-rate of water in an open channel to the geometryof the channel. The algorithm utilizes the key idea that while theprofile of the channel changes continuously as we sweep alongthe 1D flow network, the combinatorial structure only changes atdiscrete events. By tracking when these events occur, the algorithmefficiently computes the depth of the river along the flow network.We note that previous work computing the 2D channel assumesthe cross section of each channel has a simple geometry (e.g. rect-angular or trapezoidal) [6]. In contrast, we do not make any suchassumption.We note that our terrain-flow query algorithm builds on someideas in [13] namely, it also performed a sweep in the (− z ) -direction.However, our new algorithm is more general, simpler, faster, andseveral new ideas are needed. First note that the terrain-flood andflood-query problems studied in [13] are special cases of the terrain-flow query problem. In the process of computing the flow rate of alledges, it computes the flood-time of all depressions of Σ . In contrast,the algorithms in [13] either compute which depression got floodedor the flood time of only one depression. Furthermore, computationof the latter required a complicated algorithm that relied on matrixmultiplication and ran in O ( n ω ) time in the worst case, while ouralgorithm takes O ( n log n ) time in the worst case. Additionally,our algorithm can handle rain distributions that vary over time,and our algorithm is output-sensitive. Finally, we also present anI/O-efficient algorithm for answering terrain-flow queries. In this section we give a number of preliminary definitions anddescribe the flooding model, most of the text here follows closely [13,15].
Terrains.
Let M be a triangulation of R , and let V be the set ofvertices of M ; set n = | V | . We assume that V contains a vertex v ∞ at infinity, and that each edge { u , v ∞ } is a ray emanating from u ;the triangles in M incident to v ∞ are unbounded. Let h : M → R be a height function. We assume that the restriction of h to eachtriangle of M is a linear map, that h approaches + ∞ at v ∞ , and thatthe heights of all vertices are distinct. Given M and h , the graph of h , called a terrain and denoted by Σ = ( M , h ) , is an xy -monotonetriangulated surface whose triangulation is induced by M . Critical vertices.
There is a natural cyclic order on the neighborvertices of a vertex v of M , and each such vertex is either an upslope or downslope neighbor. If v has no downslope (resp. upslope) neigh-bor, then v is a minimum (resp. maximum ). We also refer to a mini-mum as a sink . If v has four neighbors w , w , w , w in clockwiseorder such that max ( h ( w ) , h ( w )) < h ( v ) < min ( h ( w ) , h ( w )) then v is a saddle vertex. Level sets, contours, depressions.
Given ℓ ∈ R , the ℓ -sublevelset of h is the set h <ℓ = { x ∈ R | h ( x ) < ℓ } , and the ℓ -level set of h is the set h = ℓ = { x ∈ R | h ( x ) = ℓ } . Each connected component of h <ℓ is called a depression , and each connected component of h = ℓ iscalled a contour . Note that a depression is not necessarily simplyconnected, as a maximum can cause a hole to appear; 𝑣 𝑣 𝑣 𝑢 𝑢 𝑢 𝑢 𝛼 𝛽 𝑢 𝑢 𝑢 𝑢 𝑣 𝑣 𝑣 𝛼 𝛽 𝛼 𝛽 𝛼 𝛽 𝛼 𝛽 𝛽 𝛼 (a) (b) Figure 1:
An example terrain with saddle vertices v - v . Eachsaddle v i delimits two maximal depressions α i and β i . (a) Terrainseen from above. Sinks are marked with a square and saddles aremarked with a cross. (b) Terrain seen from the side,For a point x ∈ M , a depression β x of h <ℓ is said to be delimitedby the point x if x lies on the boundary of β , which implies that h ( x ) = ℓ . A depression β is maximal if every depression β ⊃ β contains strictly more sinks than β . A maximal depression thatcontains exactly one sink is called an elementary depression . Notethat each maximal depression is delimited by a saddle, and a saddlethat delimits more than one maximal depression is called a negativesaddle . For a maximal depression β , let Sd ( β ) denote the saddledelimiting β , and let Sk ( β ) denote the set of sinks in β . The volume of a depression β of h <ℓ isVol ( β ) = ∫ β ( ℓ − h ( v )) dv . (1) Merge tree.
The maximal depressions of a terrain form a hierarchythat is easily represented using a rooted tree, called the merge tree [7, 11] and denoted by T . Suppose we sweep a horizontal plane from −∞ to ∞ . As we vary ℓ , the depressions in h <ℓ vary continuously,but their structure changes only at sinks and negative saddles. Ifwe increase ℓ , then a new depression appears at a sink, and twodepressions merge at a negative saddle. The merge tree is a treethat tracks these changes. Its leaves are the sinks of the terrain,and its internal nodes are the negative saddles. The edges of T , themerge tree, are in one-to-one correspondence with the maximaldepressions of Σ , that is, we associate each edge ( u , v ) with themaximal depression delimited by u and containing v . See Figure 1for an example. We assume that T has an edge from the root of T extending to + ∞ , corresponding to the depression that extends to ∞ . For simplicity, we assume that T is binary, that is, each negativesaddle delimits exactly two depressions. Non-simple saddles can beunfolded into a number of simple saddles [9]. T can be computed in O ( n log n ) time [7], and it can be prepro-cessed in O ( n ) additional time so that for a point x ∈ R , Vol ( β x ) ,the volume of the depression delimited by x can be computed in O ( log n ) time [7]. In the I/O model, T can be constructed using O ( Sort ( n )) I/Os and as shown in [3], the volume Vol ( β x ) and thesmallest maximal depression containing x can be computed for allvertices x using O ( Sort ( n )) I/Os.
We now describe the flooding model, which is the same as in [13],and define flow-rate functions.
IGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA Aaron Lowe, Svend C. Svendsen, Pankaj K. Agarwal, and Lars Arge
Flow graph and flow functions.
We transform M into a directedacyclic graph M , referred to as the flow graph , by directing eachedge { u , v } of M from u to v if h ( u ) > h ( v ) , and from v to u oth-erwise, i.e., each edge is directed in the downward direction. Foreach (directed) edge ( u , v ) , we define the local flow λ ( u , v , t ) to bethe portion of water arriving at u that flows along the edge ( u , v ) to v at time t . By definition for any u ∈ V , (cid:205) ( u , v )∈ M λ ( u , v , t ) = λ ( u , v , t ) is, in general, based on relative heights ofthe downslope neighbors of u . If u is not a negative saddle vertex,then λ ( u , v , t ) remains the same for all t , so we will often drop t and write λ ( u , v ) to denote λ ( u , v , t ) . If u is a negative saddle, then λ ( u , v ) changes when one of its depressions fills up as no waterflows for u to that saddle; see below for further discussion. Rain distribution.
Let R ( v , t ) : V × R → R ≥ denote a rain dis-tribution , that is, for each vertex v ∈ V , R ( v , t ) indicates the rateat which the rain is falling on v at time t . We assume that fora each v , R ( v , ·) is piecewise-constant function of time, with thefunction changing at discrete time values { t = , t , . . . , t k } , andfor all v and t ≥ t k , R ( v , t ) =
0. For a depression β , we define R ( β , t ) = (cid:205) v ∈ β R ( v , t ) . For i ≤ k , let | R i | denote the number ofvertices for which R ( v , t i ) ≥
0, and let | R | = (cid:205) ki = | R i | . In practice | R i | ≪ n . Fill and spill rates.
For a maximal depression β , we define the fill rate F β : R ≥ → R ≥ as the rate at which water is arriving inthe depression β as a function of time. That is, the rate at whichrain is falling directly in β plus the rate at which other depressionsare spilling water into β . The fill rate F β is a monotonically nonde-creasing piecewise-constant function. Similarly, we define the spillrate S β : R ≥ → [ , ] as the rate (as a function of time) at whichwater spills from β through the saddle that delimits β . Flow rate.
Next we define flow rates ϕ e and ϕ v for edges e andvertices v of M , which is the amount of water flowing through e and v , respectively, at time t . For an edge ( u , v ) ∈ M , ϕ ( u , v ) ( t ) isthe fraction of water from u that passes along ( u , v ) as a functionof time. That is, ϕ ( u , v ) ( t ) = λ ( u , v , t ) ϕ v ( t ) . (2)The flow-rate ϕ v of a non-saddle vertex v is the sum of the flow-rates along incoming edges to v plus the rain on the vertex v . Thatis, ϕ v ( t ) = R ( v , t ) + (cid:213) ( u , v )∈ M ϕ ( u , v ) ( t ) . (3)Letting τ v be the time at which a vertex v becomes flooded, ϕ v ( t ) and ϕ ( u , v ) ( t ) for any v ∈ M are undefined for t ≥ τ v . That is to say,when a vertex is flooded, the flow-rate function is undefined.Let v be a negative saddle delimiting depressions α and β . Untilone of α or β is filled, ϕ v is defined using equation 3. Without lossof generality assume that depression α fills first, say at time τ α ,and water starts spilling from α to β through v . The spill rate S α specifies the rate at which water spills from α to β . It is temptingto simply add S α to equation 3, but it double counts the amount ofwater that was flowing from v to depression α . For t < τ α , ϕ v isdefined as in (3). For t ≥ τ α , ϕ v is defined as follows, ϕ v ( t ) = (cid:18) R ( v , t ) + (cid:213) ( u , v )∈ M ϕ ( u , v ) ( t ) (cid:19) (cid:213) w ∈ β λ ( v , w , ) + S α ( t ) . (4) Finally for t ≥ τ α and for any w ∈ β , λ ( v , w , t ) = λ ( v , w , ) (cid:205) z ∈ β λ ( v , z , ) . (5) In this section we describe an internal-memory algorithm that,given a terrain Σ and rain distribution R , computes the flow-ratefunctions ϕ ( u , v ) ( t ) for all edges ( u , v ) ∈ M .Each function is piecewise constant and represented as a se-quence of pairs ( t i , ∆ i ) , where ∆ i represents the change in thefunction value at time t i . We denote | ϕ ( u , v ) | to be the combinato-rial complexity of ϕ u , v ( t ) , i.e., the number of pieces defining ϕ ( u , v ) .Set | ϕ | = (cid:205) ( u , v )∈ E | ϕ u , v | . If the rain distribution changes at O ( ) times, then | ϕ | = O ( nχ ) , where χ is the height of the merge tree.The preprocessing step builds, in O ( n log n ) time, the merge tree T of Σ , and labels each node in T according to its in-order traversalas described in [13]. Furthermore, it computes Vol ( β v ) for eachvertex v ∈ Σ and augments each edge ( u , v ) ∈ Σ with the index ofthe smallest maximal depression containing v .We first present a simple algorithm and then describe how toadapt it to make it output-sensitive. We first give an overview of the algorithm. We sweep through thevertices of T in descending height order, maintaining the followingvalues for each height ℓ :(1) the set of depressions α i in the sublevel set h <ℓ ;(2) for each α i , maintain the fill rate F α i ( t ) , and(3) for each α i and the set of edges ( u , v ) ∈ M with h ( v ) < ℓ ≤ h ( u ) and v ∈ α i , denoted by 𝐸 ( α i ) , maintain the flow-ratefunction ϕ ( u , v ) ( t ) .Note that F α i depends only on the sinks contained in α i , sowe compute fill-rate functions of only maximal depressions. Ateach non-saddle vertex v we compute ϕ v and ϕ ( v , w ) , for each edge ( v , w ) ∈ M , in a straightforward manner using (2) and (3). If v isa negative saddle delimiting two depressions α and β , we mustaccount for any rain spilling from one depression to the other. Wedo so by first partitioning the edges into two sets, 𝐸 ( α ) and 𝐸 ( β ) .Then given the flow-rates on the edges in each set, we compute thefill-rates F α and F β . Given these fill-rates, we then determine whichdepression (if any) fills first. If at least one depression becomes full,assume that β fills first without loss of generality. Then we cancompute the spill-rate function S β ( t ) . Finally, we compute ϕ ( v , w ) for each edge ( v , w ) ∈ M .We will now describe the algorithm in detail. We begin by aug-menting T with additional information which will be used whencomputing the fill-rates.Given a rain distribution R , before performing the downwardsweep, we first compute R ( β , ·) for each maximal depression β asfollows: first assign R ( v , ·) for each vertex with nonzero rainfall tothe smallest maximal depression containing v , and then perform anupward sweep on the merge-tree T , maintaining the sum of rainfallfunctions at each saddle vertex.As we sweep top-down, for each depression α i we store the flow-rates of the edges in 𝐸 ( α i ) in a data structure that supports thefollowing three operations: (i) add ϕ ( u , v ) ( t ) to the data structure; D and 2D Flow Routing on a Terrain SIGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA (ii) remove ϕ ( u , v ) ( t ) from the data structure; (iii) and given a saddlevertex delimiting two depressions β , γ ∈ α i , partition the flow-ratesinto two sets 𝐸 ( β ) and 𝐸 ( γ ) and return the sum of flow-rates ineach set.Noting there are O ( χ + k ) values of t at which the flow-rate func-tions can change (corresponding to the times at which R changesand spill-times of tributaries of the depression in question) we buildthe data structure as follows. For each time t i at which the functionscan change, maintain the values ∆ i for the flow-rate of each edge ( u , v ) sorted according to the label ℓ ( v ) . Additionally maintain theprefix sum of the values for each t i . There will be O (| ϕ ( u , v ) |) valuesto consider at each edge, so we can insert and remove flow-rates ofedges in O (| ϕ ( u , v ) | log n ) -time.To partition the edges into 𝐸 ( β ) and 𝐸 ( γ ) , recall each vertexis labeled according to an in-order traversal of T . Given two de-pressions β and γ delimited by a saddle vertex v with label ℓ ( v ) ,all vertices contained in β will have a label less than ℓ ( v ) , and allvertices contained in γ will have a label greater than ℓ ( v ) . Thus, topartition the edges into the depressions they cross into we simplytake the prefix of all edges with label less than ℓ ( v ) to be one set,and the remaining suffix to be the other set. Then given the parti-tion, we can use the prefix sums to determine the sum of flow-ratefunctions into β and γ in O ( χ + k ) time.With this data structure at hand, we now describe how to processeach vertex v as we encounter it in our sweep. Non-negative-saddle vertex. If v is a non-negative-saddle vertexcomputing ϕ v is easy using (3). Note that ϕ ( u , v ) for all incomingedges to v has been computed. Additionally, using F α i ( t ) , for themaximal depression α i containing v , along with Vol ( β v ) , we candetermine if or when v becomes flooded. Negative saddles. If v is a negative saddle delimiting two depres-sions α and β , we first determine whether either of α or β becomesfull and, if so, which fills first. The fill rate of a depression α is thesum of the rain falling directly on vertices in α plus the flow-ratesacross all edges crossing into α , that is F α ( t ) = (cid:213) u ∈ α R ( u , t ) + (cid:213) ( v , w )∈ 𝐸 ( α ) ϕ ( v , w ) ( t ) . (6)The first sum has already been computed in the upward sweepof the merge tree. However, to compute ϕ ( v , w ) ( t ) for all values of t we would need to know which depression fills first and when. Soinstead, we define the pseudo-flow-rate function , ϕ ′( v , w ) , in the samemanner that flow-rate for non-saddle vertices is computed as in (3).Using this, we compute a modified pseudo-fill-rate function F ′ α ( t ) = (cid:213) u ∈ α R ( u , t ) + (cid:213) ( v , w )∈ 𝐸 ( α ) ϕ ′( v , w ) ( t ) . (7)Before α or β become full (i.e. t < min ( τ α , τ β ) ) we have that ϕ ′( v , w ) ( t ) = ϕ ( v , w ) ( t ) , which in turn implies that F ′ α ( t ) = F α ( t ) (resp. F ′ β ( t ) = F β ( t ) ) before α or β becomes full. Given this, compute ϕ ′( v , w ) ( t ) for all edges from the saddle v and add these flow-ratesto the data structure. We then use the data structure to partitionthe edges and compute the pseudo-fill-rates of the two depressions,from which we can determine which (if any) depression fills first.If neither α nor β becomes full, we are done. If not, assumewithout loss of generality that α becomes full first and spills into β . Then given F α ( t ) we compute S α ( t ) and in turn compute ϕ v ( t ) for τ α ≤ t using (4) along with F β ( t ) using (6).Computing R ( β , ·) for all maximal depressions can be done in O (| R | + ρk ) -time where ρ is the number of sinks in the terrain and k is the number of times the rain distribution changes. To compute theflow-rate functions of edges from all non-negative saddle vertices,along with the pseudo-fill-rate functions from negative saddles,we spend O (| ϕ | log n ) -time computing the flow-rate functions andadding them to the data structure. At a negative saddle delimitingtwo depressions α , β , if we partition the edges by walking from thefront and back of the sorted list of edges we spend O ( min (| α | , | β |) -time, where | α | is the number of vertices contained in α . A simplerecurrence shows that the total time spent partitioning edges atall negative saddles is O ( n log n ) . Finally computing the fill andspill rates at a negative saddle takes time proportional to theircomplexity, which is O ( k + χ ) . Noting that | ϕ | = O ( n ( χ + k )) , weget the following.Theorem 3.1. Given a triangulation M of R with n vertices anda height function h : M → R that is linear on each face of M , a datastructure of size O ( n ) can be constructed in O ( n log n ) time so thatfor a (time varying) rain distribution R , a terrain-flow query can beanswered in O ( ρk + n ( χ + k ) log n ) time, where ρ is the number ofsinks in M , k is the number of times at which the rain distributionchanges. In particular, if k = O ( n ) , we have that the terrain-flow querycan be answered in O ( n log n ) time. We modify the above algorithm so that its running time becomes O ( ρk + | ϕ | log n ) where | ϕ | is the total complexity of all non-zeroflow-rate functions. If water flows across a small portion of Σ then | ϕ | ≪ n and the new algorithm will be faster. The key idea is wedo not have to process vertices with zero flow-rate.The algorithm begins in the same manner as the previous algo-rithm, computing R ( β , ·) for each maximal depression β .However rather than sweeping over all vertices of T , we insteadmaintain a priority queue of vertices with nonzero flow-rate withtheir heights as their priority, i.e., the larger the height, the higherthe priority. Every vertex v with nonzero flow-rate will either haverain falling directly on it, or have a path in M from a source vertex(i.e. either a vertex with rain falling directly on it, or a saddle fromwhich water is spilling.) We initialize the priority queue with allvertices on which rain falls directly along with all negative saddles.While the priority queue is not empty, we pop the highest priorityvertex v from the queue and compute the flow-rate function ϕ v ( t ) ,as we have described in the previous algorithm. For each vertex w adjacent to v in M with ϕ ( v , w ) >
0, add w to the priority queue.For each edge, in addition to each value ∆ i representing the changein the flow-rate function at time t i , we also store V i , the integral ofthe flow-rate function from t = t = t i . These are also storedin sorted order, and we maintain the prefix sum of these values foreach edge.At negative saddles, to determine which depression fills first,when we partition the edges we also compute the sum of volumesat the last time step t j for the depressions α and β delimited bythe saddle vertex. If both have filled by this time, consider the sum IGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA Aaron Lowe, Svend C. Svendsen, Pankaj K. Agarwal, and Lars Arge at t j − and so on until we find the latest time at which only onedepression is filled. This will be the depression which fills first, andensures the number of time steps considered will be O (| S β |) . Thenwe compute the fill and spill-rates as we did in the prior algorithm.Since we compute the flow and spill-rates as before, it sufficesto prove that we consider each negative saddle with positive spill-rate and each vertex with positive flow-rate. We begin with allnegative saddles in the priority queue, so the former holds trivially.For the latter, given a vertex v with positive flow-rate, either rainfalls directly on it, in which case we added v when initializing thepriority 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 turnbe added and later processed. Since all non-source vertices withpositive flow-rate have a path from source vertex, and all sourcevertices begin in the priority queue, we have that we will reach andprocess v .The time spent processing non-saddle vertices v and edges ( v , w ) will be O (| ϕ v ( t )|) and O (| ϕ ( v , w ) ( t )|) respectively. So the total timeprocessing these will be O (| ϕ |) . At negative saddle vertices when wepartition the edges, if we perform the search by walking from thefront and back of the list we will spend a total of O (| ϕ | log | ϕ |) timepartitioning edges at all negative saddles. Finally, at each negativesaddle v which is a source we spend O (| S β |) -time computing whichdepression fills first. Since | S β | ≤ | ϕ v | , the total time processingsaddle vertices is O ( ρk + | ϕ v ( t )| log n ) .Theorem 3.2. Given a triangulation M of R with n vertices anda height function h : M → R that is linear on each face of M , a datastructure of size O ( n ) can be constructed in O ( n log n ) time so thatfor a (time varying) rain distribution R , a terrain-flow query can beanswered in O ( ρk + | ϕ | log n ) time, where ρ is the number of sinks in M , k is the number of times at which the rain distribution changes,and | ϕ | is the total complexity of all non-zero flow-rate functions. In this section we describe two I/O-efficient algorithms that given aterrain Σ and a rain distribution R determine the flow-rate ϕ ( u , v ) ( t ) for all edges ( u , v ) ∈ Σ . The algorithms use the same frameworkas described in Section 3. However, we cannot explicitly maintainthe list of edges crossed by the sweep line in memory due to thebounded internal memory size. We instead store the merge tree inmemory and rely on the time-forward processing technique, wherethe flow-rates computed at vertices are forwarded to downslopeneighbors using an I/O-efficient priority queue [8]. The I/O-efficientpriority queue supports n insertions and deletions in O ( Sort ( n )) I/Os and O ( n log n ) internal computation time [16].In the preprocessing step of both algorithms, we compute themerge tree T of Σ and label each node in T according to its in-ordertraversal as described in [13]. Furthermore, we compute Vol ( β v ) for each vertex v ∈ Σ and augment each edge ( u , v ) ∈ Σ with theindex of the smallest maximal depression containing v . This canbe computed in O ( Sort ( n )) I/Os using the algorithm described byArge et al. [3]. Since both algorithms assume ρ = O ( m ) , we assumethe T is stored in internal memory. We now describe our first I/O-efficient algorithm. In order to ex-tend the internal terrain-flow query algorithm, we introduce thefollowing notation: let ˆ 𝐸 ( α ) ⊆ 𝐸 ( α ) be the set of edges such that,for each edge ( u , v ) ∈ ˆ 𝐸 ( α ) , α is the smallest maximal depressioncontaining v .As in the previous section, we proceed by performing an upwardsweep on the terrain followed by a downward sweep. Upward sweep.
For the upward sweep of our algorithm, we com-pute R ( α , ·) for each maximal depression α . Given a rain distribution R , we start by computing the sum of rain falling directly in eachmaximal depression α . We do this by assigning R ( v , ·) for eachvertex with nonzero rainfall to the smallest maximal depressioncontaining v . This step can be implemented I/O-efficiently using asort and a scan of the terrain and R ; for each maximal depression α store R ( α , ·) in memory. We then perform the upward sweepdescribed in Section 3 in memory, maintaining the sum of rainfallfunctions at each saddle vertex. This upward sweep can be triviallyimplemented in O ( Sort (| R | + n )) I/Os and O (| R | log | R | + n log n + ρk ) internal computation time. Downward sweep.
We sweep through vertices of Σ in descendingheight order, maintaining the following values in memory for eachheight ℓ :(1) for each depression α i in the sublevel set h <ℓ , maintain thefill-rate F α i ( t ) , and(2) for each α i maintain the depression flow-rate sum (cid:205) ( u , v )∈ ˆ 𝐸 ( α ) ϕ ( u , v ) ( t ) We cannot explicitly store the flow-rates for the set of edges crossingthe sweep line in memory. Instead, we store each flow-rate ϕ u , v ( t ) in an I/O-efficient priority queue Q keyed on the height of vertex v . Furthermore, we initialize Q by inserting the functions R ( v , t ) for each vertex v with R ( v ) >
0. Whenever we process a vertex v ,for all ( v , u ) ∈ Σ we propagate the flow-rate ϕ ( v , u ) forward to u using Q such that the function can be removed from Q whenever u is processed. The values maintained for each depression a i insublevel set h <ℓ can be maintained in memory since we assume ρ ( χ + k ) = O ( m ) .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 ) ∈ Σ , we remove the flow-rate ϕ ( u , v ) ( t ) from Q . Since the algorithms visits vertices in descendingorder of height, the functions must have been inserted into Q at anearlier point and will be at the front of Q . Since each edge ( u , v ) nolonger crosses the sweep line, we update the depression flow-ratesums by subtracting ϕ ( u , v ) ( t ) from the flow-rate sum of the smallestmaximal depression containing v .Let α i be the smallest maximal depression containing v . If v isa non-negative-saddle vertex, we compute ϕ v and determine if orwhen v becomes flooded, using F α i ( t ) and Vol ( β v ) as described inSection 3.If v is a negative saddle delimiting two depressions α and β , wewish to compute whether α or β becomes full and, if so, which fillsfirst. In order to compute the flood times of the two depressions,we recall that the fill-rate of a depression α is equal to the sum offlow-rates across all the edges ( u , w ) ∈ 𝐸 ( α ) plus the amount of D and 2D Flow Routing on a Terrain SIGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA rain falling in α . When we are processing the vertex v we havethat all edges crossing into α and β respectively are also crossingthe sweep line. We therefore compute the fill-rate of α ( β resp.) bysumming the depression flow-rate sums for all maximal depressions α i ⊆ α ( β i ⊆ β resp.): F α ( t ) = R ( α , t ) + (cid:213) α i ⊆ α (cid:213) ( u , v )∈ ˆ 𝐸 ( α i ) ϕ ( u , v ) ( t ) . (8)The flood times and the flow-rate ϕ v ( t ) is then computed from F α ( t ) and F β ( t ) as described in Section 3.Finally, we can propagate the flow-rates on the outgoing edgesby pushing ϕ ( v , w ) ( t ) = w ( v , w ) · ϕ v ( t ) to Q for all vertices w where w ( v , w ) >
0. Furthermore, we update the depression flow-rate sumswith the flow-rates of the outgoing edges. Note that each flow-rateis added to only one flow-rate sum.The flow-rate ϕ ( v , w ) ( t ) of each edge ( v , w ) is forwarded onlyonce using the priority queue. Thus, we spend a total of O ( Sort (| ϕ |)) I/Os and O (| ϕ | log | ϕ |) internal computation time processing allvertices in the sweep excluding computation of the fill-rates atnegative saddles. In order to efficiently compute the fill-rates ininternal memory, we store the depression flow-rate sums in a sortedlist ordered by the in-order traversal index of the depressions. Wethen partition the sums as in [13] and compute the fill-rates using O (| ϕ | log | ϕ |) internal computation time in total. Since the initialsorting of vertices and edges can be performed in O ( Sort ( n )) I/Osand O ( n log n ) internal computation time, the algorithm uses a totalof O ( Sort (| ϕ |)) I/Os and O (| ϕ | log | ϕ |) internal computation time.Theorem 4.1. Given a triangulation of M with n vertices, a heightfunction h : M → R which is linear on each face of M and a raindistribution R , a terrain-flow query can be answered in O ( Sort (| ϕ |)) I/Os and O (| ϕ | log | ϕ |) internal computation time assuming ρ ( χ + k ) = O ( m ) , where | ϕ | is the total complexity of all flow-rate functions whichwe return, χ is the height of the merge tree, k is the number of timesat which the rain distribution changes, ρ is the number of sinks in M ,and m is the size of internal memory. We now extend the algorithm to relax the assumption on the sizeof the internal memory from ρ ( χ + k ) = O ( m ) to ρ = O ( m ) , at thecost 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 inmemory for each depression in the sublevel set h <ℓ ; We instead thepriority queue Q to forward fill-rates as well as the edge flow-ratesused to compute fill-rates at negative saddle vertices. Forwarding fill-rates.
Let v be non-negative-saddle vertex andlet α i be the smallest maximal depression containing v . Let u be thevertex visited after v in the downward sweep, where the smallestmaximal depression containing u is also α i . When performing thesweep, we forward F α i ( t ) from v to u using Q . We note that we canaugment v with the height of u using Sort ( n ) I/Os in preprocessing,and thus we can forward F α i ( t ) to u during the sweep. Furthermore,for each negative saddle vertex v delimiting depressions α and β ,we forward F α ( t ) and F β ( t ) to the first vertices visited in α and β ,respectively. Computing fill-rates at negative saddles.
Let v be a negativesaddle vertex delimiting two depressions α and β . During executionof the sweep, we wish to compute the fill-rates of depressions α and β . We recall that the fill-rate of α can be computed as follows: F α ( t ) = R ( v , t ) + (cid:213) ( u , v )∈ E ( α ) ϕ ( u , v ) ( t ) . (9)We note that the flow-rates required to compute this sum and R ( v , t ) can be propagated using Q . However, that would in the worst-caselead to forwarding O ( nχ ) functions in total. We recall that F β v ( t ) isforwarded to v using Q . Furthermore, since F β v ( t ) = F α ( t ) + F β ( t ) ,it suffices to compute either F α ( t ) or F β ( t ) , whichever requiresthe fewest flow-rates to be forwarded. The number of flow-ratefunctions that need to be forwarded in order to compute F α ( t ) and F β ( t ) , respectively, can be precomputed by counting the number ofedges crossing the boundaries of α and β . This precomputation stepcan trivially be implemented by performing a scan of the verticesusing O ( Scan ( n )) I/Os and O ( ρ ) memory. We therefore preprocessfor which depressions we compute fill-rates and forward only theflow-rates required for computing those.We now bound the number of edges forwarded using a similarrecurrence as [13]; Let | α | denote the total number of edges with atleast one vertex contained in the depression α and note that | α | ≥| 𝐸 ( α )| . Letting T ( β v ) be the total number of flow-rates summed tocompute the fill-rates for all saddles contained in β v , we have that T ( β v ) = O (cid:0) min (| α | , | β |) (cid:1) + T ( α ) + T ( β ) . (10)Noting that the set of edges crossing into the two maximal depres-sions α and β are disjoint, it follows that | a | + | β | ≤ | β v | . Using this,the recurrence solves to T ( β v ) = O (| β v | log | β v |) , and we thus needto forward only a total of O ( n log n ) flow-rate functions. Since thecomplexity of each flow-rate function is bounded by O ( χ + k ) , wespend a total of O ( Sort (( χ + k ) n log n )) I/Os and O (( χ + k ) n log ( kn )) internal computation forwarding edges and computing fill-rates.Theorem 4.2. Given a triangulation of M with n vertices, a heightfunction h : M → R which is linear on each face of M and a raindistribution R , a terrain-flow query can be answered in O ( Sort (( χ + k ) n log n )) I/Os and O (( χ + k ) n log ( kn )) internal computation timeassuming ρ = O ( m ) , where χ is the height of the merge tree, k isthe number of times at which the rain distribution changes, ρ is thenumber of sinks in M , and m is the size of internal memory. The terrain-flow query can naturally be used to answer edge-flowqueries by returning ϕ e ( t ) for a query edge e = ( q , r ) ∈ M . Whilein practice the query time can be improved, the worst case runningtime under the MFD model remains the same as the terrain-flowquery. Under the SFD model, we can improve this running timesignificantly, building on the fast algorithm for the flood-time queryunder SFD given by Rav et al . [15], along with a linear-size datastructure supporting constant time reachability queries in planardirected graphs given by Holm et al . [10].The key idea of the algorithm is that under the SFD model whenwater falls on a vertex or spills from a negative saddle the waterflows along a single path to some sink in the terrain. Thus if wecan find which vertices and negative saddles from which water IGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA Aaron Lowe, Svend C. Svendsen, Pankaj K. Agarwal, and Lars Arge follows a path containing ( u , v ) , ϕ ( u , v ) will be the sum of waterfalling directly on or spilling from these sources. Before describingthe algorithm, we begin with some definitions.Let u be a negative saddle, let ( u , v ) and ( u , v ) be two downedges in T from u , and let ( w , u ) be the up edge from u . We call thedepression associated with ( u , v ) (resp. with ( w , u ) ) as the sibling (resp. parent ) (depression) of that associated with ( u , v ) . Any givenpoint q ∈ M is contained in a sequence of maximal depressions α ⊃ · · · ⊃ α k ∋ q . Each α i is delimited by a saddle v i and has acorresponding sibling depression β i . Note that these saddles forma path in T from q to the root. We refer to the maximal depressions β , . . . , β k − as the tributaries of q and denote them by b q .For any point q ∈ M we define the tributary tree G q as follows. G q is a directed graph with nodes corresponding to the tributariesof q plus β q . There is an edge ( β i , β j ) in G q if water spills fromthe saddle v i to a sink in β j when β i becomes full. Water spills toexactly one sink under the SFD model, so G q will be a tree rootedat β q .We present an O ( n log n ) -time algorithm for preprocessing Σ into a linear-size data structure that can answer an edge-flow queryfor a given rain distribution R ( t ) and query edge ( q , r ) . The querytakes O (| R | + b q k log n ) , where b q is the number of q -tributaries β with R ( β , t ) > ( q , r ) ∈ M and a rain distribution R ( t ) , webegin by assuming that water flows from q to r . Since water onlyflows from each vertex to one neighbor in the SFD model, if thiswere not the case then we would immediately have that ϕ ( q , r ) = ϕ ( q , r ) = ϕ q . For simplicity, we will instead compute theequivalent point-flow query.The algorithm begins by building G q . See Figure 2 for an example.As we have noted, water will reach q in one of two cases: rain fallson a vertex v of the terrain and follows a path which crosses q , orwater spills from a tributary of q and reaches q .Consider first the case when rain falls only on a single point p contained in some tributary β i . Take the path π in G q from β i to the root β q , ( β i , β i , · · · , β i k , β q ) . For each depression β i j let V i denote the depression volume of β i j and let τ j be the fill-time of β i j . The fill-time τ k , when βi k begins spilling into β q , will be whenthe volume of rain falling on p equals (cid:205) V i . Moreover we have F β q ( t ) = S β ik ( t ) = (cid:40) t < τ k F β i ( t ) t ≥ τ k . That is to say, instead of computing the fill and spill rates for eachtributary along the path, we can merge all the tributaries in thispath and treat it as if it were a single depression. Then to answer thequery, check whether there is a path from the saddle delimiting β i k to the query vertex q . If there is, we have ϕ q ( t ) = S β ik ( t ) , otherwise ϕ q ( t ) = ϕ q as a sum of the rainfall functions on verticesthat directly reach q , along with the spill-rates of parents of β q in G q that reach q . We begin by computing the initial fill-rate ofeach tributary in which rain falls directly. For each vertex v fromwhich rain flows to a sink contained in β q , check whether the paththe water takes crosses q . If so add their rainfall to the sum. If rainfalls in multiple tributaries that have disjoint paths in G q to β q β q π π π Figure 2:
The tributary tree G q rooted at β q with each vertexdenoting a tributary of q , rain initially falls in the tributaries markedas squares.(excluding the root β q ), we can simply perform the single-pointrain algorithm multiple times for each path and add each spill-ratesfrom depressions which reach q to the sum. However it might be thecase that two such paths intersect at a tributary γ before reaching β u (e.g. π and π in Figure 2.) Here, we can compute the spill-ratesof the two parent tributaries of γ as we did in the single-point rainalgorithm, and then recurse, treating γ as a single-point sourceof rain with fill-rate equal to the sum of spill-rates of its parenttributaries. Cases where more than two paths merge at a singlevertex can be handled in a similar manner.Now it remains to show how we can perform this algorithmefficiently. There are two main operations needed. First, we mustcompute the fill and spill-rates of the tributaries of q . Then we mustdetermine from which saddles and vertices water will reach q .To perform the first operation efficiently, build the linear-size fastdata structure for flood-time queries as described in [15]. With thisdata structure we can compute the spill-rates of parent tributariesof β q in O (| R | + b q k log n ) time, where b q as defined above, is thenumbe of q -tributaries in which rain is falling directly.To perform the second operation efficiently, we consider the sub-graph of the flow graph M taking only edges for which λ ( u , v , ) >
0, i.e., when v is the lowest neighbor of u . Then on this directed pla-nar graph, build the linear-size data structure reachability queriesas described in [10]. This data structure supports constant-timereachability queries. For a non-saddle vertex v , water falling at v reaches the query vertex q if q is reachable from v . For waterspilling from a negative-saddle v towards a node u we instead checkwhether q is reachable from this vertex u . Omitting all the details,we obtain the following:Theorem 5.1. Given a triangulation M with n vertices and aheight function h : M → R which is linear on each face of M , adata structure of size O ( n ) can be constructed in O ( n log n ) time sothat for a (time varying) rain distribution R ( t ) and an edge ( q , r ) anedge-flow query can be answered in O (| R | + b q k log n ) , where | R | isthe complexity of the rain distribution, b q is the number of tributariesin which rain is falling directly, and k is the number of times at whichthe rain distribution changes. So far we have assumed that water flows along the edges of Σ ,and computed the flow rate of water along these edges. But inreality the water is not constrained to edges and flows along 2Dchannels forming a 2D network of rivers. In this section we firstdescribe a model for determining the 2D channels given a 1D flow D and 2D Flow Routing on a Terrain SIGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA network. We then present an efficient algorithm for computingthese channels for a given path along the edges of M . We assume that we are given a path P along the edges of M . For eachedge e of P , let ϕ e ∈ R ≥ be the flow value along e . Unlike previoussections, we assume that ϕ e does not vary with time. But ϕ e mayvary with edges. The goal is to compute a 2D channel C : = C ( P , ϕ ) along which water flows. The channel C is defined by its leftand right banks and the height of water at every point on P . Moreprecisely, P is parameterized as P : I → R where I = [ x , x ] is aninterval. For every x ∈ I , we define lb ( x ) , rb ( x ) ∈ R as the left andright bank, respectively, of C at x , and ∆ ( x ) = h ( lb ( x )) = h ( rb ( x )) . The locus of points lb ( x ) (resp. rb ( x ) ), for x ∈ I , traces a curve lb(resp. rb), which is the left (resp. right ) bank of C . We overlay M with lb and rb. C is the portion of this overlaybetween lb and rb; see Figure 3. The complexity of C , denoted by | C | , is the number of vertices in C .To estimate lb ( x ) and rb ( x ) , we use Manning’s equation [14], awidely used empirical formula relating the channel geometry andflow rate as follows.Let x be a point on an edge e ∈ M with flow value ϕ e . Let ℓ x bethe line in the xy -plane passing through x and normal to e , and let Π = ℓ x × R be the vertical plane containing ℓ . Let Σ x = Σ ∩ Π bethe cross-section of Σ in Π , which we refer to as the profile of Σ at x . Σ x is a polygonal chain whose vertices (resp. edges) are theintersection points of edges (resp. faces) of Σ with Π . See Figure 3.Let ˆ x = ( x , h ( x ) be the vertex on Σ x corresponding to the point x ∈ e , i.e. ˆ x is vertically above x .We divide Σ x into two polygonal rays L x , R x at the vertex ˆ x ,with L x (resp. R x ) lying to the left (resp. right) of P . For a height z ≥ h ( x ) , let λ ( z ) (resp. ρ ( z ) ) be the first point on L x (resp. R x ) atheight z as we walk along L x (resp. R x .) Let A x ( z ) denote the areaof the polygon formed by the segment λ x ( z ) ρ x ( z ) and the portionof Σ x between λ x ( z ) and ρ x ( z ) , and let P x ( z ) denote the arc lengthof Σ x between λ x ( z ) and ρ x ( z ) . If the water has height z at x , thenManning’s equation [14] states that the flow rate at x is ϕ x ( z ) = A x ( z ) / σ / e µ e P x ( z ) / , (11)where σ e is the slope in the z -direction of the edge of Σ cor-responding to e , and µ e is Manning’s roughness coefficient. Weassume that we are given the value of µ e , which depends on thematerial of the terrain at e (e.g. concrete, dirt, light brush, etc.) Man-ning’s equation is typically used to compute the flow rate ϕ x ( z ) of rivers given a measurement of the river depth and approximatechannel geometry. Here instead we state an inverse problem: giventhe flow rate ϕ x at x , determine the depth and width of the riverat x . Let ∆ ( x ) be the value of z for which ϕ x ( z ) = ϕ e . We setlb ( x ) = λ ( ∆ ( x )) and rb ( x ) = ρ ( ∆ ( x )) , i.e., the river bank points on One method of generating P is computing the flow rate along all edges of M as inSections 3 and 4, fixing a time t = c , choosing a threshold ψ , extracting the edges e with ϕ e ( c ) ≥ ψ , post-processing these edges to construct a 1D flow network withsome desired properties, and finally decomposing this flow network into a set of paths.These steps are beyond the scope of this paper and an interesting direction for futureresearch. Here we assume that lb and rb are simple curves; if either of them is self-intersecting,then C has to be defined more carefully. Σ corresponding to x ; Let C x = Σ [ lb ( x ) , rb ( x )] be the profile of thechannel at x . See Figure 3.We first describe how we compute ∆ ( x ) for a fixed x and thendescribe how to track lb and rb as we vary x . For simplicity, wemake the following two assumptions:(A1) C x is unimodal for all x ∈ I .(A2) The point ˆ x is the unique minimum of C x .We discuss in Section 6.4 how to relax these assumptions. ∆ ( x ) Recall that we assume C x to be unimodal with ˆ x as its uniqueminimum. Without loss of generality, assume that the edge e con-taining x is parallel to the x -axis, so Π is parallel to the yz -plane. Weraise the value of z starting from h ( x ) and stopping at the heightof vertices of Σ x until we find a vertex ˆ v = ( v , h ( v )) such that ϕ x ( h ( v )) ≥ ϕ e . We now know the edges of L x and R x that containlb ( x ) and rb ( x ) . We then compute the points themselves.We now describe the procedure in more detail. Let f , x , f , x , · · · be the sequence of edges of R x , ordered from left to right, and let e i − , x , e i , x be the endpoints of f i , x ; e , x = x . Recall that each edge f i , x is the intersection of a face f i of Σ with the vertical plane Π , and each endpoint e j , x is e j ∩ Π for some edge e j of Σ . Foreach edge f i , x , let f ↑ i , x be the semi-infinite trapezoid formed bythe edge f i , x and the vertical rays in the ( + z ) -direction from theendpoints e i − , x , e i , x of f i , x . For a value z ∈ R , we define thetrapezoid f ↑ i , x ( z ) to be the intersection of f ↑ i , x with the halfspace z ≤ z ; f ↑ i , x ( z ) may be empty, or it may be a triangle. We define A i , x ( z ) to be the area of f ↑ i , x ( z ) and P i , x ( z ) to be the length of thebottom edge of f ↑ i , x ( z ) , which is a portion of the edge f i , x . Let e i , x = ( x , a i ( x ) , b i ( x )) denote the coordinates of e i , x as a functionof x , and set w i ( x ) = a i ( x ) − a i − ( x ) and h i ( x ) = b i ( x ) − b i − ( x ) .Then A i , x can be written as A i , x ( z ) = z < b i − ( x ) , ( z − b i − ) w i ( x ) h i ( x ) b i − ( x ) < z < b i ( x ) , w i ( x )( h i ( x ) + ( z − b i ( x )) b i ( x ) < z . (12)Similarly P i , x can be expressed as P i , x ( z ) = z < b i − ( x ) , ∥ f i ( x )∥ ( z − b i − ) h i ( x ) b i − ( x ) < z < b i ( x ) , ∥ f i ( x )∥ b i ( x ) < z . (13)We note that P i , x (resp. A i , x ) is a piecewise-linear (resp. piecewise-quadratic) function of z for a fixed x . We can define F j , x ( z ) , P j , x ( z ) and A j , x ( z ) for the edges f j , x of L x as well. We can now express P x and A x as: P x ( z ) = (cid:213) i P i , x ( z ) and A x ( z ) = (cid:213) i A i , x ( z ) , (14)where the summation is taken over all edges of Σ x that containa point of height at most z . P x and A x are also piecewise-linear andpiecewise-quadratic functions respectively. IGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA Aaron Lowe, Svend C. Svendsen, Pankaj K. Agarwal, and Lars Arge
Let z < z < z < · · · be the heights of vertices of Σ x . P x ( z ) = A x ( z ) =
0. Assuming P i , x ( z i − ) , A i , x ( z i − ) have been computedthen P i , x ( z i ) , A i , x ( z i ) can be computed in O ( ) time using (12)–(14).Let z k be the first value for which ϕ x ( z k ) ≥ ϕ e . Since P i , x ( z ) (resp. A i , x ( z ) ) is a linear (resp. quadratic) function for z ∈ ( z k − , z k ) ,the value of ∆ ( x ) ∈ ( z k − , z k ] can be computed in O ( ) time byplugging these functional forms in (11). Let f L , x (resp. f R , x ) be theedge of L x (resp. R x ) at which the vertical sweep stopped. Thenlb ( x ) (resp. rb ( x ) ) is the point on f L , x (resp. f R , x ) of height ∆ ( x ) and can be computed in O ( ) time. The total time spent by theprocedure is O (| C x |) . Hence, we obtain the following. Lemma 6.1.
For a given point x ∈ P , ∆ ( x ) , lb ( x ) and rb ( x ) can becomputed in O (| C x |) time. We now describe an algorithm for compute the channel C assuming(A1) and (A2) hold. Recall that for any x ∈ I , the vertices of C x are intersection points of the edges with the plane Π x . Let Γ x = ⟨ γ , · · · , γ u ⟩ denote the sequence of these edges, which implicitlydefine C x . We refer to Γ x as the combinatorial structure of C x . Wecompute C by varying x continuously from x to x and maintaining C x . As x varies, C x changes continuously, i.e., each vertex of C x traces a curve, but the combinatorial structure Γ x changes onlyat discrete values of x , called the events . The algorithm works bysweeping the line ℓ x along P , stopping at events as we traverse P .As long as x lies on the same edge of P , ℓ x simply translates. Atvertices of P , where the sweep line ℓ x shifts from one edge to thenext one in P , the algorithm continues by rotating ℓ x to make itnormal to the next edge. We first describe how we sweep along anedge of P and then describe how the sweep line rotates at a vertexof P . Edges.
Fix an edge e ∈ P . Without loss of generality, assume that e is parameterized as e : [ , ] → R .As we sweep along e and vary x , the left and right banks lb , rbtrace curves that lie inside fixed faces of Σ and the remaining ver-tices of C x trace the corresponding edges of Σ . The algorithmencounters the following two types of events at which Γ x changes:(1) the sweep line reaches an endpoint of an edge ( u ′ , v ′ ) of Γ x which is a vertex of Σ , or(2) lb or rb intersects an edge e ′ (bounding the face containingit) of Σ .The first event results in one or more edges in C x shrinking to thepoint v ′ , and one or more new edges starting at v ′ . The secondevent results in either the addition and/or removal of an extremaledge in C x (and thus insertion/deletion of an edge in Γ x ) dependingon whether the height of the channel is increasing or decreasing.The first type of events are easy to detect, as they correspond tothe vertices of Σ . The second type are more challenging, and wedetect them as follows. By maintaining functions representing thearea and wetted perimeter of C x as a bivariate function of x and z (using (12-14)) and using Manning’s equation we can express ∆ ( x ) as a function of x . We then compute when ∆ ( x ) reaches the top orbottom boundary of the face containing lb (resp. rb). These timeinstances correspond to the second type of event.Process the events in order, by popping the first event from thepriority queue Q . If it is the first type of event corresponding to a f f f f f f f f e e e e e e e e e x (cid:96)λ ( z ) ρ ( z ) P ( z )ˆ xL x R x f ↑ ,x ( z ) Figure 3:
Top: Σ with the edge e marked in purple and line ℓ . Theintersection of ℓ with edges of Σ marked with boxes. Bottom: Σ x ,with water at height z . The polygon P is marked in blue, and thewetted perimeter marked in red. The vertices correspond to theintersection point above in the top figure.vertex v of Σ , we remove from Γ the edges that end at v and insertinto Γ the edges that start at v . We add the other endpoints of thenew edges to Q as new first type of events. We also remove from A ( x , z ) and P ( x , z ) the terms corresponding to the old edges andadd terms corresponding to the new edges. If an edge of the face of Σ containing lb or rb changes, we also update the second type ofevent in Q .If the event corresponds to lb or rb reaching an edge of theterrain, either add the new face the it crosses into and/or removethe face it crosses out of. We update Γ as well as Q . We also A ( z ) and P ( z ) accordingly. A ( x , z ) and P ( x , z ) can be maintained using a height balancedtree so that its functional form can be updated in O ( log n ) time perinsertion/deletion of a term in A ( x , z ) and P ( x , z ) .We continue this process until we reach the event correspondingto x =
1, when the endpoint of e is reached. Let | C e | be the numberof total faces contained in the channel from C ( ) to C ( ) . Betweenany two events, lb ( x ) and rb ( x ) intersect a face only a constantnumber of times. Therefore the total number of events is O (| C e |) ,giving a total running time to sweep an edge of O (| C e | log | C e |) . Vertex.
When transitioning from one edge to another along thepath, we must join the two channels. We will assume the thereare no sharp bends between two edges of P , specifically the anglebetween the sweep line along the first and second channel is lessthan 90 degrees, Taking the two channels C ( u , v ) and C ( v , w ) wesee that on one side of v the two channels will intersect, while onthe other they will be disjoint. Let ℓ (resp. ℓ ) be the sweep lineperpendicular to ( u , v ) (resp. ( v , w ) ) at v , and p (resp. p ) be theintersection point between ℓ with the riverbank of C ( v , w ) (resp. ℓ with the riverbank of C ( u , v ) .) Then let ℓ ′ (resp. ℓ ′ ) be the lineparallel to ℓ at p (resp. ℓ at p ), and v ′ be the intersection point D and 2D Flow Routing on a Terrain SIGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA vu w(cid:96) (cid:96) p p (cid:96) (cid:48) (cid:96) (cid:48) v (cid:48) Figure 4:
The river profile is taken as the union of the profile ofthe two edges and that as we rotate around the vertex from ℓ to ℓ .between ℓ ′ and ℓ ′ . Now rotating about v ′ from ℓ ′ to ℓ ′ will sweepthe area where the two channels intersect. See Figure 4.In a similar manner in which we swept along each edge, let C v ′ , θ be the profile of v ′ with a line ℓ at angle θ . As we rotate the sweepline, we assume the water is flowing at the intersection of ℓ witheither ( u , v ) or ( v , w ) . We will similarly maintain the area function A ( θ , z ) (resp. the wetted perimeter function P ( θ , z ) ) as a sum offunctions on each face in the channel. Additionally when comput-ing Manning’s equation, as we rotate we will linearly interpolatebetween slope values as well as the roughness coefficients for thetwo edges. The main difference is now for each face f i ∈ F , h i ( θ , z ) and w i ( θ , z ) are not linear functions in θ as it was when we sweptalong an edge. However, we can still computed analytically theevents corresponding to the riverbank(s) crossing edges of faces.Letting | C v | denote the total number of faces in the channel ob-tained by sweeping at v . Then the total time spent is O | C v | log | C v |) Putting everything together we obtain the following:Theorem 6.1.
Given a triangulation M with n vertices, a heightfunction h : M → R which is linear on each face of M a path P in M , if the channel C x is unimodal for all x along the path, we cancompute the 2D flow network in time O (| C | log | C |) where | C | is thetotal number of faces in the channel obtained by sweeping along theedges and vertices of P . We will now show how to modify the algorithm when assump-tions (A1) and (A2) do not hold. When (A2) does not hold, i.e. ˆ x is not a local minima, the modification is straightforward. Simplywalk down the edges of Γ x until finding a local minima x ′ andthen run the algorithm using that vertex as the division point forthe polygonal rays L x ′ and R x ′ . It may be the case that the edgecontaining x ′ ends and a new minima x ′′ begins, but this can behandled like any other event in the algorithm. The only differencecomes in the interpretation and runtime analysis. When (A1) doesnot hold, it may be the case that the water level in the 2D channeldoes reach an edge along which water was flowing according tothe 1D flow network. However this is reasonable behavior, as the1D flow network assumes water only flows along edges of the ter-rain whereas in reality water will also flow along the faces. For the ˆ x f i e i e (cid:48) Figure 5:
A channel profile Σ x with local maxima y . If the flowin the primary channel on the right causes the channel height toexceed h ( y ) , excess spills to the secondary channel on the left.runtime analysis, we now include the edges searched along from x to x ′ to be included in the channel C e .When (A1) does not hold, i.e. the channel is not unimodal, whenlb or rb intersect an edge of Σ bounding the face containing it, thatedge may correspond to a local maximum in the profile. In this case,water flows over the ridge into a secondary channel. See Figure5 . We must now account for the decrease of the flow-rate in theprimary channel as well as determining the height of water in thesecondary channel. Updating flow rates.
Let f i be the first face containing a localmaximum encountered when sweeping along an edge maintaining Γ x , and let e i be the edge of Σ corresponding to the maximum. Us-ing Manning’s equation, we compute the flow rate correspondingto when this unimodal channel becomes full ϕ f i , x = ϕ x ( h ( e i , x ) .If ϕ e > ϕ f i , x a spill will occur, and the excess flow is sent tothe local minimum in the secondary channel. To account for thewater flowing out of the channel, consider the flow rate along e as a function of x , ϕ e ( x ) . Letting e ′ be the minimum edge inthe secondary channel, we set ϕ e ′ ( x ) = ϕ e − ϕ f i ( x ) . As we con-tinue the sweep ϕ e ( x ) will be a non-increasing function, we set ϕ e ( x ) = min ( ϕ e , ϕ f i , x , min x ′ < x ϕ e ( x ′ )) . Secondary channel.
When an overflow occurs we continue addingfaces to the profile until finding the next local minima correspond-ing to edge e ′ . Then in this secondary channel if the excess flowexceeds our threshold for the 1D channel, i.e. ϕ e ′ ( x ) ≥ ψ , we repeatthe unimodal channel algorithm in this channel using ϕ e ′ ( x ) . Itmay be the case that ϕ e ′ ( x ) is enough to fill the secondary channelabove the height of a ridge. If it is the same ridge that spilled intothe secondary channel we treat it as a single channel. If it is a ridgefurther out, we repeat the process in the tertiary channel. In this paper we presented algorithms for a number of flow routingproblems: First, we developed a fast internal-memory algorithmfor the terrain-flow query problem which given a rain distribution,computes the flow rate ϕ e for all edges e ∈ Σ . We also showed howthe algorithm could be made I/O-efficient. Next, we presented afaster algorithm if one is interested in computing the flow rate ofonly one edge, after some preprocessing. Finally, given a flow pathalong the edges of Σ , we proposed an algorithm to determine the2D channel along which water flows; our algorithm does not makeany assumption about the geometry of the channel. Here | C x | contains all faces between the leftmost and rightmost river banks, i.e., alledges marked red. IGSPATIAL ’20, November 3–6, 2020, Seattle, WA, USA Aaron Lowe, Svend C. Svendsen, Pankaj K. Agarwal, and Lars Arge
We conclude by mentioning a few directions for future work. Themodel of extracting 2D channels leaves a number of open questions.For instance, if the 1D flow network is a forest then channels alongdifferent paths will interact which leaves the question of how wemerge these channels to construct a 2D river network.While we consider the flow rate as a function of time, it onlychanges when the rain distribution changes or a spill event occurs.That is, the effects of such events are propagated to all reachablevertices instantaneously. While this assumption is reasonable forlocal effects and for flash floods when a large volume of rain fallsover a short duration, an interesting question is to make the modelmore general and account for the time it takes water to flow overthe terrain.
REFERENCES [1] A Aggarwal and JS Vitter. 1988. The input/output complexity of sorting andrelated problems.
Commun. ACM
31, 9 (1988), 1116–1127.[2] L Arge, M Rav, S Raza, and M Revsbæk. 2017. I/O-Efficient Event Based DepressionFlood Risk. In
Proc. 19th Workshop on Algorithm Engineering and Experiments .259–269.[3] L Arge and M Revsbæk. 2009. I/O-efficient Contour Tree Simplification. In
Intl.Sympos. on Algos. and Computation . 1155–1165.[4] L Arge, M Revsbæk, and N Zeh. 2010. I/O-efficient computation of water flowacross a terrain. In
Proc. 26th Annu. Sympos. on Comp. Geom.
Journal of Experimental Algorithmics
Journal of hydrology
Comp. Geom.
24, 2 (2003), 75–94.[8] Y Chiang, MT Goodrich, EF Grove, R Tamassia, DE Vengroff, and JS Vitter. 1995.External-memory graph algorithms. In
Proc. Sixth Annu. ACM-SIAM Sympos. onDiscrete Algos.
Proc. 17th Annu. Sympos. Comp. Geom. . IEEE, 370–389.[11] M Kreveld, R Oostrum, C Bajaj, V Pascucci, and D Schikore. 1997. Contour treesand small seed sets for isosurface traversal. In
Proc. 13th Annu. Sympos. on Comp.Geom.
Proc. 11th Intl.Sympos. on Spatial Data Handling . 137–148.[13] A Lowe and PK Agarwal. 2019. Flood-Risk Analysis on Terrains under theMultiflow-Direction Model.
ACM Trans. Spatial Algorithms Syst.
5, 4, Article 26(Sept. 2019), 27 pages. https://doi.org/10.1145/3340707[14] R Manning, JP Griffith, TF Pigot, and LF Vernon-Harcourt. 1890.
On the flow ofwater in open channels and pipes .[15] M Rav, A Lowe, and PK Agarwal. 2017. Flood Risk Analysis on Terrains. In
Proc.of the 25th ACM SIGSPATIAL Int. Conference on Advances in GIS . ACM, 36.[16] Peter Sanders. 2001. Fast Priority Queues for Cached Memory.
ACM J. Exp.Algorithmics