1) Let us first ignore the technical issues concerning a triangulated irregular network and just discuss the idealized smooth problem where there is given a smooth height profile $z=\phi(x,y)$ of the terrain. One can locally introduce a stream function $\psi=\psi(x,y)$ such that the curves
$$\psi(x,y)~=~{\rm constant}$$
represent the streamlines. (However, there are global obstructions, see item (4) below.) The gradients of $\psi$ and $\phi$ must be perpendicular,
$$ \nabla\psi \cdot \nabla\phi~=~ \frac{\partial \psi}{\partial x}\frac{\partial \phi}{\partial x}+\frac{\partial \psi}{\partial y}\frac{\partial \phi}{\partial y} ~=~0. \tag 1$$
This is a 1st-order linear PDE in $\psi$ in two variables $(x,y)$. Its solution $\psi$ formally solves the smooth problem locally.
2) If moreover the horizontal fluid velocity field $(u,v)$ is divergence-free
$$ \frac{\partial u}{\partial x}+ \frac{\partial v}{\partial y} ~=~0,$$
one can demand that
$$ u = \frac{\partial \psi}{\partial y}, \qquad v = -\frac{\partial \psi}{\partial x}. $$
3) If furthermore $(u,v)$ is also curl-free (=vortex-free), then the stream function $\psi$ becomes a harmonic function. See also flownets.
4) Global obstructions. The streamfunction $\psi$ is ill-defined in sources and sinks, i.e. in local extrema.
5) Now let us comment on the triangulated surface, with vertices, edges and faces. If one (instead of introducing a streamfunction $\psi$) just solve the problem one-streamline-at-the-time, one could end up with streamlines that un-physically cross each other because of numerical errors. Solving in terms of the streamfunction $\psi$ protects against such unphysical solutions locally.
6) Let $\phi$ and $\psi$ be defined on the vertices. The gradients $\nabla\psi$ and $\nabla\phi$ naturally live on the faces, i.e. the dual graph, so that the equation (1) can be made discrete.
