Calculating basin-averaged ksn values

Posted on

The normalized river steepness (Ksn) is one of the most frequently used topogrpahic metrics in tectonic geomorphology. TopoToolbox has the function ksn that enables calculating this metric for each node in the river network. Often, however, researchers are rather interested in calculating basin average values of ksn rather than a ksn value for each river node. This is more tricky. Hence, here is a quick solution.

For our example, I am using the Taiwan DEM which you can download using the function readexample.

DEM = readexample('taiwan');

The DEM needs a little preprocessing including filling voids with nan and minima imposition which carves through artefact topographic sinks.

DEM = inpaintnans(DEM);
FD  = FLOWobj(DEM);
DEM = imposemin(FD,DEM);

Ksn requires upstream area which is calculated using flowacc. We use upstream area also to delineate the stream network which requires setting a minimum contributing area. Here we use 1000 pixels.

A   = flowacc(FD);
S   = STREAMobj(FD,A > 1000);

Ksn requires a common concavity index (theta) which is often set to 0.45. We then calculate Ksn using the ksn function.

theta = 0.45;
k   = ksn(S,DEM,A,theta);

Now things become a bit more tricky because we want to calculate the Ksn value for each drainage basin. Here is how:

First, we calculate the drainage basins using the drainagebasins function with the stream network as input. drainagebasins computes the outlets of the stream network and derives the catchments of these outlets.

D   = drainagebasins(FD,S);

Second, I use getnal to extract the drainage basin affiliation of each river node.

d   = getnal(S,D);

We use the drainage basin indices in d to calculate the average Ksn value in each basin.

km  = accumarray(d,k,[],@mean);

Finally, we need to map these average values back to a grid that spatially aligns with the DEM.

K   = GRIDobj(DEM)*nan;
K.Z(D.Z>0) = km(D.Z(D.Z>0));
imageschs(DEM,K,'colorbarylabel','K_{sn} [m^0.9]')


And that’s it.

Introduction to @DIVIDEobj

Posted on Updated on


This blog post was written by Dirk Scherler

In this blog entry, we demonstrate the funcionalities of the new class DIVIDEobj.

Getting started

We start by loading the well-known DEM of the Big Tujunga catchment into our workspace and derive flow directions using the D8 flow routing algorithm.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD  = FLOWobj(DEM);
ST = STREAMobj(FD,'minarea',1000);

The resulting DEM and stream network look familiar:

hold on
hold off


Drainage divide network

Let’s now obtain the drainage divide network from the stream network we just created. The flow accumulation threshold we chose controls the extent and branching of the stream network and thus also the divide network. We will look at this again.

axis image


The figure shows a plot of the (unsorted) drainage divide network, which consists of divide segments, junctions, and endpoints. Endpoints are generally located close to streams, where drainage basin boundaries can be thought of to start or end. Junctions are places where three or more divide segments meet. Let’s zoom in a little to get a better idea of what this means.

hold on
axis equal
xlim([3.9804 3.9900].*1e5)
ylim([3.7968 3.7975].*1e6)
hold off


The divides are shown in red and the rivers as black lines. All the divides that you see correspond to the boundaries of drainage basins that start at junctions. Even the divides that delineate the basin in the northeastern corner of the figure, which appear to touch each other are consistent with the outline of the associated drainage basin.

Connectivity of streams and divides

To understand how the position of divide segments is related to the gridded DEM, let’s zoom into the nodes of the divide segment and plot them on the grid structure of the DEM.

A = GRIDobj(DEM);
A.Z(1:2:end) = 1;
[x,y] = ind2coord(D,D.IX);
imagesc(A,[0 10])
hold on
xlim([3.8646 3.8687].*1e5)
ylim([3.8007 3.8009].*1e6)
hold off


Because divide segments follow drainage basin boundaries, their nodes are positioned on pixel corners and the divide edges that connect the nodes and constitute a divide segment have either a vertical or horizontal orientation. For each edge that connects two nodes, there exist two pixels from neighboring basins. This aspect allows us to assign attributes to each divide edge that can be related to various other grids derived from the DEM, for example. In doing so, we can either examine the two neighboring pixels directly, or use them to access other downstream pixel values. We will show how to do this when introducing the methods associated with a DIVIDEobj.

Contents of the DIVIDEobj

To get a better feeling for the DIVIDEobj, let’s have a look into the contents of the DIVIDEobj:



A DIVIDEobj is a structure with several fields, some of which are similar to those of a GRIDobj. These include ‘size’, ‘cellsize’, and ‘refmat’. The other fields are unique to the DIVIDEobj. The field IX contains linear indices to the divide nodes that constitute the divide segments and the divide network. The linear indices point into a GRIDobj that has an extent that equals that of the DEM with an additional row and column and that is shifted by half a cell size so that the cell centers coincide with cell corners in the DEM.

The fields ‘order’ and ‘distance’, if set, have the same length as the field ‘IX’ and arer thus properties of the nodes. The fields ‘ep’, ‘jct’ and ‘jctedg’ are related to endpoints, junctions, and the number of edges linked to junctions, respectively. The values of these fields are also linear indices into the grid as described above. Finally, the fields ‘issorted’ and ‘ordertype’ indicate if the divide network has been sorted and if it has been ordered. When creating an instance of DIVIDEobj, it will get automatically sorted unless you specify not to do so. However, the network is not yet ordered.

Edge effects

Before we proceed, let’s have a closer look at the edges of the DEM to see how the truncated topography affects the divide network.

axis image
D2 = cleanedges(D,FD);
axis image


When you take a closer look at the divide network on the edge of the DEM, you will see that most divides that lead towards the DEM edge have their final course along the DEM edge. This is clearly a spurious divide pattern that is related to drainage basins that touch the DEM edge. To avoid such spurious divides, it may thus be useful to eliminate all divides that follow the edge of the DEM.

Drainage divide order and distances

Higher up, we mentioned already that drainage divide segments can be assigned a distance and an order. Both propereties are related to the sorted nature of the drainage divide network. With sorted, we mean that the network can be considered to be directed. Just like stream networks, the divide network has a tree-like structure, in which the endpoints would correspond to the leaves of the tree and junctions correspond to branch forks. Let’s visualize this before we explain how order and distance are computed.

D = divorder(D2,'topo');
plot(D,'limit',[1000 inf],'color','k')
axis image
box on


This figure shows the drainage divide network, with the linewidth being proportional to the divide order and all divides starting at a divide distance of 1000 m. Note that the default settings of the plot function change according to whether divide orders have been calculated or not.

The sorting step in the derivation of the divide network can be regarded as the core of the algorithm. We provide a detailed explanation in the recently published paper (Scherler and Schwanghart, 2020a), but in brief: we iteratively compile the sorted drainage divide network by adding divide segments that are connected to endpoints and remove these divide segments from the unsorted collection of divide segments. After each iteration, junctions with only one segment remaining are identified and turned into endpoints. In some cases a junction may be encountered early in the search process, but it can take many more iterations before it turns into an endpoint. This is true for junctions located on the far eastern end of the drainage divides that delimit the Big Tujunga catchment. The itartion ends if all divide segments are transferred into the sorted drainage divide network, or if no more endpoints exist. At this point we should emphasize that the latter condition can be true if the DEM contains internally-drained basin, even if not all divide segments have been sorted. Therefore, expect issues when your DEM contains internally-drained basins.

Once the divide segments are sorted and the tree-like divide network has a direction, we can compute the divide distance along the divide network. We defined this distance to be the maximum directed distance from an endpoint. To illustrate this, consider to start walking up on the divide network from an endpoint at a stream junction. At each divide junction, the direction will be given by the subsequent divide segment that was added to the sorted divide network. Effectively, we are going down the tree from the leaves, to first smaller, than bigger branches and ultimately down the stem of the tree. The root of the tree and thus the greatest divide distance will be the junction that was the last one added to the sorted divide network. Let’s mark this point in the divide network shown in the figure.

[x,y] = ind2coord(D,D.IX(end-1));
hold on
hold off


The drainage divide network in the figure has thicker lines in places with greater distance, although the line thickness is proportional to the divide order. The divide order is computed in the same fashion it is done for stream networks and thus includes the Strahler and Shreve ordering schemes. In the TopoToolbox, we added another ordering scheme that we called ‘Topo’, in which divide orders increase by one at each junction. Because drainage divide segment lengths are approximately normally distributed the Topo orders are approximately linearly related to the divide distance. See more details on the ordering scheme in Scherler and Schwanghart (2020a).

Colored divide network plots

So far, we used the function plot to show the divide network, once to illustrate endpoints and junctions, and once to illustrate the divide order/distance. There exists another plotting function that allows illustrating the divide distance more precisely. The function plotc generates a colored plot of the divide network according to properties of the divide edges. The following command creates a figure of the divide network in which the divide edges are colored by the divide distance.

plotc(D,D.distance./1e3,'limit',[1000 inf])
box on
axis image
hc = colorbar;
hc.Label.String = 'Divide distance, d_d (km)';


In this call to the function, we provided a vector with attributes (divide distance) that has the same length as the field D.IX. Instead of such a vector, we can also derive attributes from a GRIDobj, following the logic outlined in section “Connectivity of streams and divides”. The following figure shows the divide network colored by divide elevation.

plotc(D2,DEM,'limit',[1000 inf])
box on
axis image
hc = colorbar;
hc.Label.String = 'Divide elevation (m)';


Any GRIDobj can be used an an input to color-code the divide network edges. For more complex calculations, we can also obtain the grid values adjacent to the divide edges using the function getvalue.

DZ = vertdistance2stream(FD,ST,DEM);
DZ.Z(isinf(DZ.Z)) = nan;
[p,q] = getvalue(D,DZ);
dz = abs(diff(abs([p,q]),1,2));

plotc(D,dz,'caxis',[0 300],'limit',[1000 inf])
box on
axis image
hc = colorbar;
hc.Label.String = 'Across-divide differences in hillslope relief(m)';




This example shows the divide network colored by the across-divide difference in hillslope relief. We may wish to normalize the differences by the sum of hillslope relief on either side of the divides to better depict the degree of asymmetry. We call this quantity the divide asymmetry index, which varies between 0 (symmetric) and 1 (most asymmetric).

dai = dz./sum([p,q],2);
dai(dai<0) = 0;
dai(isinf(dai)) = nan;

plotc(D,dai,'caxis',[0 0.5],'limit',[1000 inf])
box on
axis image
hc = colorbar;
hc.Label.String = 'Divide asymmetry index';


We discussed the spatial variations of this metric in the Big Tujunga Basin in Scherler and Schwanghartt (2020a) and tested it’s sensitivity to drainage divicde migration using experiments with a landscape evolution model in Scherler and Schwanghart (2020b).

We conclude this section with the function asymmetry, which computes the divide asymmetry index and yields a mapping structure with the divide network that can be exported as a shapefile. The optional structure S in the example below containes the sorted divide network with a number of additional properties such as x,y coordinates, the order and distance, as well as data on the divide orientation and asymmetry. We use these two outputs to create a figure that shows the divide network on top of a hillshade image, with the divides colored by the divide asymmetry index (DAI), and with arrows on the divide segments that point in the direction of the asymmetry. The length of the arrows corresponds to the average DAI.

[MS,S] = asymmetry(D,DZ);
for i = 1 : length(S)
  S(i).length = max(getdistance(S(i).x,S(i).y));

imageschs(DEM,[],'colormap',[.9 .9 .9],'colorbar',false);
hold on
plotc(D,vertcat(S.rho),'caxis',[0 0.5],'limit',[1000 inf])
axis image
hc = colorbar;
hc.Label.String = 'Divide asymmetry index';
ix = [MS.dist]>1000;
f = [S.length]./1e3;
title('Drainage divide asymmetry and direction of lower hillslope relief')


Mapping structures of divide networks with attribute data can also be obtained with the function DIVIDEobj2mapstruct. In the following example, we generate a mapping structure with several attribute fields. We refer to the help and example of this function for further information.

DX = flowdistance(FD,ST);
MS = DIVIDEobj2mapstruct(D,DEM,1000,...
  {'hr_mean' DZ 'mean'},{'hr_diff' DZ 'diff'},...
  {'fdist_mean' DX 'mean'},{'fdist_diff' DX 'diff'});
for i = 1 : length(MS)
  MS(i).dai = MS(i).hr_diff./MS(i).hr_mean./2;

% visualize
msdo = [MS.order];
msdai = [MS.dai];
symbolspec = makesymbolspec('Line',...
  {'order',[2 max(msdo)],'Linewidth',[0.5 6]},...
  {'dai',[0 1],'Color',flip(hot)});

hc = colorbar;
hc.Label.String = 'Elevation (m)';
hold on
ix = msdo>1 & not(isnan(msdai));
title('Divide asymmetry index: white=low -> red=high')


A benefit of mapping structures is that each divide segment has scalar-valued attributes that lend themselves for color-coding the divide network in different ways.

Divide properties by divide distance

Instead of visualizing the divide network in mapview, we can also think of ways to visualize it in profile view. We have not yet implemented this with a stand-alone function, like plotdz for STREAMobj, but it is straight forward to generat such a plot.

d = D.distance;
dz = getvalue(D,DZ,'min');



This figure shows the minimum height of the drainage divide edges above the adjacent rivers by divide distance, colored by the DAI. Divides at high divide distances (>10 km) typically hover around values of 300-500 m, but in some sections, which are also often quite asymmetric, the divides are very close to the rivers in elevation. We interpret this signature to reflect mobile divides that tap into existing drainage networks.


The new object class DIVIDEobj allows extracting and sorting of drainage divide networks from digital elevation models. We have shown how the divide network is derived and how it is structured. Several functions included as methods for the DIVIDEobj allow plotting and turning the divide network along with it’s attributes into a shapefile. Because of the ease of the extraction and display of the divide network, these functions also lend themselves for appying operations on a large amount of digital elevation models, like outputs from a landscape evolution model. In conjunction with the landscape evolution model study in Scherler and Schwanghart (2020b), we provided movies of all simulations, in which the divide network is colored by different attributes. If interested, you can find out more here:

We hope you enjoy the new features of the DIVIDEobj and find many useful applications. Please be aware that this is the first release of the new class and there may be bugs. Please report issues that you encounter through the comment option.



Scherler, D., Schwanghart, W., 2020a. Drainage divide networks – Part 1: Identification and ordering in digital elevation models. Earth Surface Dynamics, 8, 245–259. [DOI: 10.5194/esurf-8-245-2020]

Scherler, D., Schwanghart, W., 2020b. Drainage divide networks – Part 2: Response to perturbations. Earth Surface Dynamics, 8, 261-274. [DOI: 10.5194/esurf-8-261-2020]


geoglobe now in MATLAB R2020a

Posted on

In a previous post (here) I have shown how TopoToolbox is able to export data to kml-files which can be directly opened in Google Earth (or other digital globes such as ArcGIS Earth). Indeed, kml is a great way to share geographic data. However, if you want to visually explore your data, the step via kml and Google Earth may not be necessary any longer. Since its latest release, MATLAB’s Mapping Toolbox includes geoglobe, a geographic globe that allows navigating on the Earth surface and to add data. To this end, geoglobe will enable numerous applications that enhance the way we can explore data analyzed and generated using TopoToolbox.

geoglobe opens in a figure created with uifigure.

h = uifigure;
g = geoglobe(h);


Once initiated, you can start with the GeographicGlobe object g. For example, you can change the basemap, which is by default ‘satellite’.

g.Basemap = 'topographic';


There are numerous ESRI basemaps available, as well as some from Natural Earth. Terrain is based on the GMTED2010 model by the U.S. Geological Survey (USGS) and National Geospatial-Intelligence Agency (NGA) and hosted by MathWorks. The DEM has a spatial resolution of 250 m, which is sufficient for many large-scale applications. If you want to add higher resolution data, you can add DTED files with terrain data.

Finally, you can add your own geographic data. Here I show, how you can add stream network data generated with TopoToolbox. I will download some example data from Taiwan using the function readexample:

DEM = readexample('taiwan');
DEM = inpaintnans(DEM);
FD  = FLOWobj(DEM);
S   = STREAMobj(FD,'minarea',5000);
[lat,lon,z] = STREAMobj2latlon(S,DEM);


There are still some slight problems in visualizing the network because the lines remain partly hidden beneath the terrain. I reduced but not fully resolved this effect by adding 10 m to the elevations of the stream network. Notwithstanding, geoglobe is a great utility to visualize topographic data and I am sure there are more functions to come in the future that enable us to interactively work with geographic globes. So far, I particularly like how The Mathworks implemented navigation on the globe which works super smoothly.

Landscapes Live

Posted on

#shareEGU20 is currently setting an example. Virtual conferences are an effective means to bring scientists together and to spur vivid discussions.

To keep this scientific exchange alive during COVID-19, Philippe Steer, Vivi Pedersen, Stefanie Tofelde, Pierre Valla, Charlie Shobe, and me (my role was actually very minor) have initiated Landscape Live, a new remote seminar series focused on sharing exciting geomorphology research throughout the international scientific community.

The remote format allows for free, planet-conscious, and pandemic-proof attendance by anyone who is interested. Talks will take place using the Zoom meeting software.

Talks will happen during the academic year in blocks of 5-6 weekly talks, with long breaks in between. The first block will run from 28th May through 25th June 2020, with the next block taking place during the fall 2020 semester.

Each seminar in the first block will be held on a Thursday at 2 pm GMT (4 pm Central European Time). The first block features a fantastic lineup of speakers as follows:

  • Thursday 28th May at 2 pm GMT: Georgina Bennett, University of Exeter
  • Thursday 4th June at 2 pm GMT: Anneleen Geurts, University of Bergen
  • Thursday 11th June at 2 pm GMT: Liran Goren, Ben Gurion University of the Negev
  • Thursday 18th June at 2 pm GMT: Robert Hilton, Durham University
  • Thursday 25th June at 2 pm GMT: Fiona Clubb, Durham University

Please visit for the most up-to-date information, including the links to each Zoom meeting.

Suggestions for future speakers are welcome; please feel free to send names to any member of the organizing committee. We look forward to seeing you (virtually) at Landscapes Live!

This text is slightly modified from the version written by Charlie Shobe and previously published on the EGU geomorphology blog.


Posted on


Only two days left before #shareEGU20 opens its digital gates. I surely won’t be able to be online during the whole event, but I’ve my personal schedule of displays which I’ll try to view and discuss online inbetween homeoffice, homeparenting, homeschooling, …

If you haven’t a fixed schedule yet, consider getting involved in following displays which I authored or coauthored:

EGU2020-11177 – Divide mobility controls knickpoint migration on the Roan Plateau (Wolfgang Schwanghart and Dirk Scherler), Wed, 06 May, 14:00–15:45 | D1315, DISPLAY

EGU2020-19260 – The TopoToolbox v2.4: new tools for topographic analysis and modelling (Dirk Scherler and Wolfgang Schwanghart), Thu, 07 May, 10:45–12:30 | D818, DISPLAY

EGU2020-5609 – Uncertainties in Chi analysis: implications for drainage network and divide stability (Jens Turowski et al.), Tue, 05 May, 10:45–12:30 | D1132, DISPLAY

EGU2020-5900 – Evaluating the effect of variable lithologies on rates of knickpoint migration in the Wutach catchment, southern Germany (Andreas Ludwig et al.), Tue, 05 May, 08:30–10:15 | D1114, DISPLAY

EGU2020-8811 – Why do shelf-incising submarine canyons form? – Insights from global topographic analyses and regression trees (Anne Bernhardt and Wolfgang Schwanghart), Fri, 08 May, 08:30–10:15 | D1102, DISPLAY

EGU2020-3737 – Illuminating the speed of sand – quantifying sediment transport using optically stimulated luminescence (Jürgen Mey et al.), Fri, 08 May, 08:30–10:15 | D899, DISPLAY

New paper out: Divide mobility controls knickpoint migration

Posted on Updated on

In recent years, there has been a quite fierce debate about how landscapes evolve in response to lateral dynamics of river networks. These dynamics include laterally shifting rivers, their expansion or contraction in upstream and downstream direction, and the mobility of catchment divides. In fact, as we seek to gain insight into changes in climate and tectonics from the analysis of river networks, we often make the assumption that theses river networks are static, that their spatial configuration remains stable.

In our paper published now in Geology (Schwanghart and Scherler, 2020), we argue that this assumption must be made cautiously. We revisited the Parachute Creek basin, Colorado, US, which has been studied by Berlin and Anderson (2007). The site is a real natural laboratory, because it has a uniform climate as well as (sub)horizontally uniform bedrock which makes it possible to analyze the phenomenon of knickpoint migration into a relict landscape in isolation of other factors that govern knickpoint retreat.

View towards the west on the Parachute Creek basin. Imagery from Google Earth.

In this study, we applied the stream-power incision model to infer present day locations of knickpoints in river profiles of the Parachute Creek basin (see Figure above). The knickpoints emanated from a base level drop at the Upper Colorado River about ~8 Mio years ago. What we realized when looking at the misfit of predicted locations and actual locations, was that there is one subbasin where knickpoints consistently travelled further than our model would predict. Now, you would expect some randomness in knickpoint locations, but the systematic spatial pattern grasped our attention. Looking more closely at the DEM revealed that this subbasin shares some lengths of its divide with the plateau margin, a steep cliff that goes down to the Colorado River. That this margin is actively retreating became clear from numerous beheaded valleys visible in the DEM and Google Earth imagery. A hypothesis was quickly formulated: Knickpoints in this subbasin had migrated further than we expect, because part of the subbasin’s area had been lost to cliff retreat.

In the next few blog posts, I will talk a bit more about our approach to calculate the area lost, and how we actually could even gain some constraints on the timing of area loss. Until then, you may also check Dirk’s and my papers on divide mobility just published in ESURF (Scherler and Schwanghart, 2020a,b). We’ll write here about this two-part paper in due time.


Berlin, M. M. and Anderson, R. S.: Modeling of knickpoint retreat on the Roan Plateau, western Colorado, Journal of Geophysical Research: Earth Surface, 112(F3), doi:10.1029/2006JF000553, 2007.

Scherler, D. and Schwanghart, W.: Drainage divide networks – Part 1: Identification and ordering in digital elevation models, Earth Surface Dynamics, 8(2), 245–259, doi:, 2020a.

Scherler, D. and Schwanghart, W.: Drainage divide networks – Part 2: Response to perturbations, Earth Surface Dynamics, 8(2), 261–274, doi:, 2020.

Schwanghart, W. and Scherler, D.: Divide mobility controls knickpoint migration on the Roan Plateau (Colorado, USA), Geology, doi:10.1130/G47054.1, 2020. Supplementary material can be found here.

Calculating Hack’s Law using TopoToolbox

Posted on

Hack’s Law describes an empirical relationship between river length and drainage area (Hack 1957). The functional relationship is a power function with the equation L = c A^h where L is the length of the longest stream from the outlet to the divide, A is the drainage area above a particular locality, c is a constant, and h is the scaling exponent (see figure).

Measuring stream length and drainage area (Hack 1957).

In general, h is slightly below 0.6 (Rigon et al. 1996). Today, I show how to calculate the parameters of Hack’s Law. First, I will demonstrate how the parameters of Hack’s Law are derived for a single catchment. In the second part, I will then apply the technique to calculate the parameters for catchments draining Taiwan.

Big Tujunga catchment

In the first example, I will use the DEM of the Big Tujunga catchment, which is part of the TopoToolbox distribution. While Hack’s Law appears to hold for any point inside a basin (Rigon 1996), I am showing here how to derive it from set of pixels that comprise the longest path from the divide to the catchment outlet. As usual, I derive flow directions (FLOWobj) and the stream networks (STREAMobj). Then, I extract the largest basin and subsequently the longest stream in this basin that extends from the divide to the outlet.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD  = FLOWobj(DEM);
DIST   = flowdistance(FD);
S   = STREAMobj(FD,'minarea',1000);
S   = klargestconncomps(S);
D   = drainagebasins(FD,S) > 0;
DIST = clip(DIST,D);
[~,IX] = max(DIST.Z(:));
S   = STREAMobj(FD,'channelheads',IX);
imageschs(DEM,D,'falsecolor',[.0 0 0],'truecolor',[1 1 1], ...
hold on
hold off


In a second step, I calculate drainage area and distance. The function flowacc computes the flow accumulation (in pixels) for the entire grid, and the function getnal extracts the grid values for each node in the river network. In addition, I calculate flow distance using the function distance. By default, the function calculates the distance from each node to the outlet, but there are also other options. Here, we use the maximum distance from the channelhead.

A   = flowacc(FD)*DEM.cellsize^2;
a   = getnal(S,A);
d   = distance(S,'max_from_ch');

Now we can plot the two node-attribute lists versus each other.

xlabel('Area [m^2]')
ylabel('Stream length [m]')


Now lets use these data to fit Hacks Law. I am using nlinfit from the Statistics and Machine Learning Toolbox for this purpose which requires initial estimates of the parameters. The choice of these initial values can be challenging and way-off values may render nlinfit unable to converge to a good solution.

b = nlinfit(a,d,@(b,X) b(1)*X.^b(2),[0.002 0.6]);

b is a two-element vector and contains the parameters of Hack’s Law, where b(1) is the constant and b(2) is the scaling exponent. Here, b(1) is 0.0035 and b(2) is 0.8346. Having these parameters, we can now obtain estimates for distances for any value of drainage area.

hold on
fplot(@(X) b(1)*X.^b(2),xlim)
hold off
xlabel('Area [m^2]')
ylabel('Stream length [m]')


You may note that nlinfit returns an exponent different from Hack’s value of 0.6 and there may be different reasons for this. First of all, the basin may have a particular form because it extends on the relict landscapes of the Chilao Flats and features a series of knickpoints that suggest that it experiences a transient response to accelerated uplift. Second, the river contains both mixed bedrock-alluvial sections in the upper part, and a purely alluvial lower part. Finally, our data violates some of the assumptions that underly the inference of parameters, i.e. a identically and independent distribution (i.i.d.). There may be ways to at least partly account for these violations by random sampling of our data, but I don’t want to go into too much detail here.


In the second example, we will use the SRTM-3 DEM of Taiwan, which you can download using the readexample function. In this example, we calculate the parameters in Hack’s Law from multiple drainage basins.

DEM = readexample('taiwan');

We first inpaint missing values in the DEM and use FLOWobj to derive flow directions. In a second second step, we then calculate the trunk stream for each basin, again up to the divide. This is a bit tricky, if you have multiple drainage basins. Here, I show how to do this using FLOWobj2cell and cellfun. I omit many small drainage basins that occur along the coast of Taiwan by setting a minimum threshold of drainage basin are to an arbitrary size of 1000 pixels.

DEM = inpaintnans(DEM);
FD  = FLOWobj(DEM);
[CFD,~,a] = FLOWobj2cell(FD);
CFD = CFD(a>1000);
[~,IX]  = cellfun(@(fd) max(flowdistance(fd)),CFD);

S   = STREAMobj(FD,'channelheads',IX);
imageschs(DEM,[],'colormap',[1 1 1],'colorbar',false,'ticklabels','nice');
hold on


Again, I calculate upslope area and distance and plot both in semilogarithmic axes.

A  = flowacc(FD)*DEM.cellsize^2;
a  = getnal(S,A);
d  = distance(S,'max_from_ch');
OUTLETS = streampoi(S,'outlets','logical');
a  = a(OUTLETS);
d  = d(OUTLETS);
xlabel('Area [m^2]')
ylabel('Stream length [m]')
hold off


Then, I again use nlinfit to determine the parameters of Hack’s Law.

b = nlinfit(a,d,@(b,X) b(1)*X.^b(2),[0.002 0.6]);

hold on
fplot(@(X) b(1)*X.^b(2),xlim)
hold off
xlabel('Area [m^2]')
ylabel('Stream length [m]')

This time, nlinfit yields 0.5315 as exponent.


Finally, we might want to look at the spatial patterns of residuals from our predicted lengths.

rerr = (d - (b(1)*a.^b(2)))./d;
IXoutlet   = S.IXgrid(OUTLETS);
D = drainagebasins(FD,IXoutlet);
D.Z = double(D.Z);
D.Z(D.Z ~= 0) = rerr(D.Z(D.Z~=0));
D.Z(isnan(DEM.Z)) = nan;
imageschs(DEM,D,'caxis',[-1 1],'colorbarylabel','Rel. error')


The map shows that the residuals (or relative errors) from Hack’s Law. We can see that some catchments, particularly those draining the western slopes tend to have longer distances of the main river than predicted by Hack’s Law. These catchments are often rather elongated which might partly due to their coverage of large alluvial areas and lowlands.


Hack, J. T.: Studies of longitudinal stream profiles in Virginia and Maryland, USGS Professional Paper, 295, 45–97, 1957.

Rigon, R., Rodriguez-Iturbe, I., Maritan, A., Giacometti, A., Tarboton, D. G. and Rinaldo, A.: On Hack’s law, Water Resources Research, 32, 3367–3374, 1996.


Open PhD position: ‘Geochronology of the first Eurasian Ice Sheets’

Posted on

Hardanger Vidda (Photo credit: John Jansen).

Want to work with a great advisory team on the geochronology of the first Eurasian ice sheets? Well, here is your chance. My friend and colleague John Jansen (Czech Academy of Sciences) and his colleagues from Charles University and Aarhus University offer a highly exiting PhD opportunity. You will be working with geochronological techniques (cosmogenic nuclide dating) and numerical modelling to unravel the elusive traces of early Eurasian glaciations. There is definitely a lot to discover here.

Check the job advertisement here and apply!

Handling closed basins

Posted on Updated on

Much of the world’s terrestrial area is comprised of endorheic basins. These basins have no drainage to the oceans. Instead, water in these basins seeps into the ground or evaporates (wikipedia).

Internally drained basins come in different sizes. The Caspian Sea, the largest inland body of the world, has a catchment area of 3,626,000 sqkm. In contrast, dolines – small hollows in karstic terrain – have just a few tens of square meters.

In digital elevation model analysis, internally drained basins can be quite challenging, largely because we often do not know whether a closed basin in our DEM is a true closed basin or an artefact of the data.

When calculating flow directions in TopoToolbox, the default setting is that all closed basins are artefacts. However, when working in dryland areas, we must be cautious with this assumption.  In this post, I will show how to deal with closed basins. In addition, I’ll show how to use MATLAB Live Scripts that help us to interactively explore parameters that control the identification of closed basins.

To get started, please download the Tibetan Plateau example using the readexample function:

DEM = readexample('tibet');

The DEM features numerous closed basins that we identify using the fillsinks function. We then visualize the difference between the filled DEM and the actual DEM, which gives us an idea about the location and depth of internally drained basins.

DEMf = fillsinks(DEM);
'colormap',flipud(hot),'colorbarylabel','Sink depth [m]')


In a second step, we will investigate the properties of the closed basins in more detail. Here we look at two properties. The area of the sinks, and their depth. The function regionprops which is part of the Image Processing Toolbox comes in handy here.

SINK  = DIFF > 0;
stats = regionprops(SINK.Z,DIFF.Z,'Area','PixelIdxList','MaxIntensity');
sinkarea = [stats.Area]*DEM.cellsize^2;
sinkdepth = [stats.MaxIntensity];
xlabel('Area of sink [m^2]')
ylabel('Maximum depth of sink [m]')


There are more than 11,000 sinks in the DEM. Clearly, not all of them are true sinks. I propose that true sinks should have a large areal extent and be rather deep. Hence, we can apply some thresholds to classify sinks into true sinks and artificial sinks. As a first guess, lets assume a combination of an areal extent of 10 km2 and a depth of 10 m as suitable parameters of our sink identification model.

areathres = 10^7;
depththres = 10;
sinkstats = stats(sinkarea > areathres & sinkdepth > depththres);
sinkpixels = {sinkstats.PixelIdxList};

These thresholds returns six closed basins. In a third step, we calculate flow directions in a way that acknowledges the internal drainage of these basins. There are different ways to do this. Here, we use the trick that sets the pixels with the minimum elevation in each sink to nan. There are other ways to do this (see the function FLOWobj), but this approach is actually the fastest.

% A: Here we identify the linear indices of pixels 
% with minimum elevation for each sink
minelevix  = cellfun(@(ix) ix(find(DEM.Z(ix) == ...
% B: Make a copy of the DEM
% C: Set elevations of minimas to nan
DEM2.Z(minelevix) = nan;
% D: calculate flow directions
FD = FLOWobj(DEM2);
% E: calculate flow accumulations
A  = flowacc(FD);
% F: and visualize the results.
hold on
[x,y] = ind2coord(DEM,minelevix);
hold off


Finally, we may not be satisfied with the results. Based on our parameters, there may be to many or to few closed basins. Below is a gif movie that shows how to use the interactive tools in MATLAB Live Scripts to vary the sink-area and sink-depth parameters and instantly see their effects on the resulting flow network.


I hope you enjoyed this post. Applying TopoToolbox in regions with internally drained basins should now be no problem anymore.

Crossing divides

Posted on Updated on

Divide mobility and changes in flow network configuration are hot topics these winter days. How mobile are drainage divides, what controls their mobility, which are the involved time scales, and how do shifting divides affect other metrics derived from digital elevation models?

One of the issues addressed by some studies is which metrics are useful to characterize divide movements (Forte and Whipple, 2018, Scherler and Schwanghart, 2019). For example, there are chi maps which allow us to map spatial imbalances in drainage network configuration. However, whether these imbalances actually translate into divide movements, or whether these movements will occur in some distant future remains unclear. Gilbert metrics in turn quantify cross-divide differences in hillslope gradient (or other similar metrics) and thus may provide a more suitable proxy for processes that actually act along or close to divides. To this end, it may often be useful to look at both types of metrics.

Now, you know that I like visualisations. And today, I want to briefly show how to combine hillslope and divide geometry with longitudinal river profiles of two adjacent rivers. I am hoping to provide the computational means to automatically derive such profiles which might prove very insightful.

The example uses the Big Tujunga data which we first load. We derive flow directions and the stream network. We will only use the trunk rivers.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD  = FLOWobj(DEM,'preprocess','c');
S   = STREAMobj(FD,'minarea',1000);
S   = trunk(S);

In a next step, I am extracting two rivers which are sourced in the so-called Chilao Flats and drain this low-relief relict landscape. Note that extractconncomps is an interactive function and you will need to manually pick the trunk rivers yourself.


Plotting the two trunk rivers as profiles (plotdz) generates two profiles that have their outlets set to the distance value 0. That’s not so intuitive and I actually want that the channelheads should be opposing each other.

So, here is a little trick that I can use to connect the two rivers and thereby generate a path across the divide. First, I convert the STREAMobj of the two rivers to a GRIDobj. Then, I run a gray-weighted distance transform to calculate the least-cost path emanating from the outlet of one of the rivers to all other pixels in the DEM. The costs of moving from one pixel to another is elevation, but I decrease the costs at the rivers themselves, which forces the algorithm to run along the river network.

GS = STREAMobj2GRIDobj(S);
DEMaux = DEM;
DEMaux.Z(GS.Z) = 1;
IX = streampoi(S,'outlet','ix');
DEMaux.Z = graydist(DEMaux.Z,IX(1));


I now have a new GRIDobj DEMaux. I chose this variable name because DEMaux represents some auxiliary topography (the same approach is actually used to route through flat sections or closed pits in the DEM).

Now I am using this auxiliary topgraphy to calculate a FLOWobj, and a STREAMobj sourced at the outlet of the second channel. The resulting stream starts at the outlet of the second channel, runs upstream, crosses the divide, and then follows the course of the first channel until its outlet.

FDaux = FLOWobj(DEMaux);
Saux = STREAMobj(FDaux,'channelhead',IX(2));

% plot it
imageschs(DEM,[],'colormap',[1 1 1],'colorbar',false)
hold on


Ok, that’s quite neat. Now, in addition, I might want to plot the river and cross-divide profile together with chi-values.

A = flowacc(FD);
c = chitransform(S,A);
[~,zb] = zerobaselevel(S,DEM);
c = c + zb;

C = GRIDobj(DEM);
C.Z(S.IXgrid) = c;
caux = getnal(Saux,C);

hold on
h = colorbar;
h.Label.String = '\chi [m]';
hold off


Ok, that’s what I wanted. Have fun coding!


Forte AM, Whipple KX. 2018. Criteria and tools for determining drainage divide stability. Earth and Planetary Science Letters 493 : 102–117. DOI: 10.1016/j.epsl.2018.04.026
Scherler D, Schwanghart W. 2019. Identification and ordering of drainage divides in digital elevation models. Earth Surface Dynamics Discussions : in review (open discussion esurf-2019-51).