Tutorial

This tutorial walks through a complete flow-routing analysis on a synthetic DEM.

Package loading

And seeding the random number generator for deterministic runs.

using WhereTheWaterFlows, CairoMakie
using Random; Random.seed!(42);

Random.TaskLocalRNG()

Build a synthetic DEM

We use a function that produces a landscape with a few hills, valleys, and a deliberate closed depression (pit):

function peaks2(n=100, randfac=0.05)
    coords = range(-π, π, length=n)
    return coords, coords,
           sin.(coords) .* cos.(coords') .-
           0.7 .* (sin.(coords .+ 1) .* cos.(coords')).^8 .+
           randfac .* randn(n, n)
end

x, y, dem = peaks2(200)

heatmap(x, y, dem; axis=(title="DEM",))

Run the flow router

(; area, slen, dir, nout, nin, sinks, pits, c, bnds) = waterflows(dem)

waterflows returns a NamedTuple with the following fields:

Terminology and special cells

waterflows distinguishes between several special cell types:

• NaN cells (dir == BARRIER): the DEM value is NaN. No routing is performed on these cells and their cellarea contribution is ignored.
• Sinks (dir == SINK): cells where flow leaves the active domain. By default this includes all domain boundary cells (bnd_as_sink=true) and cells adjacent to NaN values (nan_as_sink=true). Additional sinks can be specified via extra_sinks.
• Pits (dir == PIT): interior local minima that have no lower neighbouring cell under the D8 rule. With drain_pits=true (the default) the algorithm routes flow over the lowest spillway of each pit, so the returned pits vector is typically empty.

Positional arguments of waterflows

The simplest usage is:

waterflows(dem)

An additional positional argument can be used to route physical fluxes:

waterflows(dem, cellarea)

Keyword arguments of waterflows

The two last listed keyword-args are advanced features: See Feedback Functionality for process-coupled examples using feedback_fn.

Upslope area

plt_area plots the log₁₀ upslope area. Sink locations are marked in red (all of the domain boundary).

plt_area(x, y, area; sinks)

Catchments

plt_catchments colours each catchment uniquely. Catchments smaller than minsize cells not plotted.

plt_catchments(x, y, c; minsize=50)

Note that setting minsize=0 (the default) may run a lot faster on large DEMs.

Delineate a single catchment

Pick the cell with the largest upslope area along row 50 and trace everything draining into it:

i = 50
j = findmax(area[i, :])[2]

cc = catchment(dir, CartesianIndex(i, j))

fig, ax, _ = heatmap(x, y, Float64.(cc))
scatter!(ax, [x[i]], [y[j]]; color=:red, markersize=15)
fig

Delineate the catchment upstream of a box

Pass a CartesianIndices rectangle to catchment to collect every cell draining into that region:

box = CartesianIndices((90:110, 90:110))
cc2 = Float64.(catchment(dir, box))
cc2[box] .= NaN

heatmap(x, y, cc2)

Note that any collection of cells, passed in as a list of CartesianIndex, can serve to delineate a catchment. See Subglacially: Defining outlet groups documentation for an example.

Fill depressions

After waterflows with drain_pits=true (the default), flow crosses depressions (via a breach-type approach). fill_dem raises each pit cell to its spillway elevation, useful for visualising lake depth:

demf = fill_dem(dem, sinks, dir)

heatmap(x, y, demf .- dem; colormap=:blues, axis=(title="Lake depth",))

Routing physical fluxes

Set cellarea to a physical value (e.g. precipitation in m³/s per cell) to accumulate real fluxes rather than cell counts:

precip = 1e-3 .* ones(size(dem))   # uniform 1e-3 m3/s/cell
discharge = waterflows(dem, precip).area

plt_area(x, y, discharge; colorbar_label="log10(discharge [m³/s])",)

Note that because waterflows does not know the cell area, cellarea has to be set to the input into a cell, e.g. m³/s, and not an input per unit area, e.g. m/s. I.e. cellarea needs to be an extensive quantity.

Routing multiple quantities simultaneously

Pass a tuple of arrays as cellarea to accumulate several quantities at once. All quantities must be extensive (additive), e.g. energy not temperature. Here we route uniformly input water alongside a tracer injected into a single cell:

water  = ones(size(dem))
tracer = zeros(size(dem))
tracer[40, 40] = 1.0   # single-cell source

(water_area, tracer_area), = waterflows(dem, (water, tracer))

fig = Figure()
heatmap(fig[1, 1], x, y, log10.(water_area);  colormap=:blues, axis=(title="Water (log₁₀)",))
heatmap(fig[1, 2], x, y, tracer_area;          colormap=:reds,  axis=(title="Tracer",))
fig