Chimaps in a few lines of code (3)

Posted on Updated on

In my previous posts (here and here), I introduced chimaps and how they are calculated in TopoToolbox. I showed that deriving a suitable stream network requires some work. Today, I will cover the dry part: the math behind.

Chi-analysis is predicated on the well-known stream power law (SPL). I will not detail the SPL but refer you to the paper of Dimitri Lague (Lague 2014). In its simplest form, the SPL shows that the energy expenditure required to incise into the stream bed is a function of stream gradient (S) and upslope area (A) and some parameters k, m and n. Combining the SPL with a simple vertical uplift rate U allows us to derive a very simple landscape evolution model for the fluvial domain:

Any change in elevation z with time t is caused by an imbalance of uplift and incision. If both processes are perfectly balanced, there is no change in elevation. We are in a state of a dynamic equilibrium or steady state:

We can rearrange this equation and solve for S=dz/dx where x is the stream distance measured from the outlet. At this time I denote that upslope area is a function of x. U and k are assumed to be spatially uniform.

The trick is then to take the integral of (dz/dx) with respect to x from a base level z(x0) (the elevation at the outlet). For convenience of units, we also introduce a reference area A0.

Now that looks complicated. Yet, by replacing the left-hand integral with z and the right-hand integral with

we can simply rewrite this equation as a function of a straight line:

A straight line? Yes… but only if we choose the right m/n ratio. I will demonstrate in the next post how to do this.

P.S.: This was really a short derivation of chi. I can recommend the paper of Perron and Royden (2013) for a more comprehensive account.


Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Process. Landforms, 39(1), 38–61, doi:10.1002/esp.3462, 2014.

Perron, J. T. and Royden, L.: An integral approach to bedrock river profile analysis, Earth Surf. Process. Landforms, 38(6), 570–576, doi:10.1002/esp.3302, 2013.

Chimaps in a few lines of code (2)

Posted on Updated on

In my last post, I started my series on chimaps. In this and the following posts I will go through the different steps of producing chimaps using TopoToolbox. These steps are:

  1. Load DEM and derive flow directions (FLOWobj).
  2. Choose a suitable threshold for stream initiation and derive a stream network (STREAMobj).
  3. Make sure that all outlets have the same elevation.
  4. Make sure that all drainagebasins are complete, i.e. that the DEM covers their entire extent. Remove incomplete drainage basins.
  5. Choose a suitable value for the mn-ratio.
  6. Calculate the transformed variable chi (Hint: Don’t use chiplot, but check the low-level function STREAMobj/chitransform!).
  7. Plot the results.

In this post I will cover step 1-4. I will use a DEM of the Mendocino region (northwestern California) provided by Amani Al Abri from the University of Boston. Her inquiry has actually prompted me to write about chimaps. The Mendocino region hosts the Mendocino Triple Junction where the Gorda plate, the North American plate, and the Pacific plate meet. The geometry of the plates and their movements have likely favored the development of a slab window and upwelling of astenospheric material, leading to differential patterns of uplift, seismic activity and extrusion of volcanic rock.

1. Load DEM and derive flow directions

DEM = GRIDobj('mendocino.tif')

The first step is to load the DEM. While this is a trivial task, let me briefly mention a few pitfalls one may encounter. While the constructor function GRIDobj tries to identify missing values (actually, it does so quite reliably with the hack of Simon (thanks!)), I usually test this using the function


which provides several information including the minimum and maximum value. If these values are far beyond expected values, they should be replaced by NaNs. If this value was -9999, following code snippet will set these pixels to NaN:


Moreover, the info function shows the coordinate system of the DEM. TopoToolbox requires projected coordinate systems (preferably WGS84 UTM) and equal horizontal and vertical units (e.g. meters). If your DEM comes in geographic coordinates, reproject the DEM with the function reproject2utm. Note also, that you should set missing values before reprojecting the DEM since the bilinear interpolation would generate some blurred regions between non-corrected missing-value regions and regions with valid values.

Missing values usually occur at the boundaries of the DEM or offshore. These are truly missing values that you may want discard. However, data gaps may also exist as individual blobs within the valid DEM domain either as single pixels or large areas of missing values. These areas should be filled because they might otherwise generate discontinuous flow paths. Unless they are too large, I usually use the function inpaintnans to fill these gaps.

DEM = inpaintnans(DEM);

Having prepared the DEM, it is now time to derive flow directions. TopoToolbox stores flow directions in an object called FLOWobj.

FD = FLOWobj(DEM,'preprocess','carve');

FLOWobjs contain all necessary information to calculate hydrological terrain attributes (see the function help for details and some more detailed description on carving and filling here, here and here).

2. Choose a suitable threshold for stream initiation and derive a stream network

Now that you have a DEM and flow directions, you can derive a stream network. TopoToolbox stores stream networks in a STREAMobj. In essence, a STREAMobj is a reduced version of the FLOWobj being limited to the channelized domain. Unless you have the locations of mapped channelheads, the standard approach to determine the channelized domain is to assume a minimum upslope area (MUA). You should neither choose a too large nor too small MUA. A large MUA may hide chimap details close to catchment divides; a small MUA may extend from the channelized domain into the debris flow or hillslope domain. Here we simply take an upslope area of 0.05 sqkm or 500 pixels:

S = STREAMobj(FD,'minarea',500);

The resulting stream network is shown in the Fig. 1:

Fig. 1: Stream network

3. Make sure that all outlets have the same elevation

Chimaps require the stream network to have outlets with the same elevation. As DEMs often fail to cover the entire altitude range of all individual streams it may be necessary to trim the stream network in a way so that all relevant streams have the same outlet elevation. The function GRIDobj/griddedcontour and STREAMobj/modify come in handy here (see also this post). The modify function is extremely versatile. Here we use its option ‘upstreamto’ that requires a GRIDobj of zeros and ones. The modify function will then modify the stream network so that it contains streams only upstream of regions with ones. We use the output of griddedcontour and modify it slightly using the function bwmorph to produce gridded contours that are four-connected. A four-connected contour line ensures that all streams crossing the contours are reliably detected.

C = griddedcontour(DEM,[50 50]);
C.Z = bwmorph(C.Z,'diag');
S = modify(S,'upstreamto',C);

The trimmed stream network is shown in Fig. 2.

Fig. 2: The stream network reduced to those streams having an outlet at an altitude of 50 m a.s.l.

4. Make sure that all drainage basins are complete

Some of the drainage basins may be located along the DEM borders and have part of their area outside the DEM domain. To avoid edge effects on the resulting chimap, these incomplete drainage basins and their associated streams should be removed from further analysis. Testing whether a drainage basin is complete, however, is not trivial and the approach used here may not be useful in all situations.

Here I choose to identify incomplete basins by testing if their boundaries touch NaN areas or DEM edges. This requires some hard-coding since there is no ready-made function in TopoToolbox so far.

% Get drainage basins
D = drainagebasins(FD,S);

% Get NaN mask and dilate it by one pixel.
I = isnan(DEM);
I = dilate(I,ones(3));

% Add border pixels to the mask
I.Z([1 end],:) = true;
I.Z(:,[1 end]) = true;

% Get outlines for each basin
OUTLINES = false(DEM.size);
for r = 1:max(D);
OUTLINES = OUTLINES | bwperim(D.Z == r);

% Calculate the fraction that each outline shares with the NaN mask
frac = accumarray(D.Z(OUTLINES),I.Z(OUTLINES),[],@mean);

% Grid the fractions
FRAC.Z(:,:) = nan;
FRAC.Z(D.Z>0) = frac(D.Z(D.Z>0));

The resulting GRIDobj FRAC is shown in Fig. 3. Several basins share a large proportion of their outlines with the DEM boundary and we don’t want these to be part of the analysis.

Fig. 3: Drainage basins colored by the fraction that their boundaries touch missing values in the DEM.

Based on the GRIDobj FRAC we now modify the STREAMobj so that it contains only streams where drainage basins share less than 1% with the NaN mask. Again, we use the function modify:

S = modify(S,'upstreamto',FRAC<=0.01);

The trimmed stream network is shown below. Unfortunately, a large portion of the network was removed.

Fig. 4: Trimmed stream network (white) that excludes incomplete drainage basins.

Wrap up until here

As you see, preparing the DEM and stream network for chi analysis takes some effort. I may have been somewhat over-optimistic regarding the title of this series of posts. We may have many lines of code. However, having now a clean network will make the following steps much easier and straightforward. I will cover these in my next posts.

Floods from breached landslide dams and the usefulness of crowd-sourced internet videos

Posted on

Frames from video of Beni Bridge taken during the peak of the Kali Gandaki flood (Source: Bricker et al. 2017)

On 24. May 2015, a landslide near the village Baisari, Myagdi District Nepal, dammed the Kali Gandaki. 15 hours later, the dam breached and released a flood wave. Fortunately, people downstream were warned of the imminent flood and no casualties were reported.

In a recent paper by Jeremy Bricker et al. (2017), we hindcasted the event using numerous 1D and 2D hydrodynamic models. Moreover, we used raw and smoothed topographies of the valley floor to investigate the effect of DEM uncertainties on flood wave propagation modelling. We show that using unsmoothed valley thalwegs result in delayed modelled flood wave arrival times that may be critical for the effectiveness early warning systems. A 2D model produced results most in line with field observations.

One of the most striking aspects of our study is the use of crowd-sourced video material available on youtube. We used video material recorded at two bridges crossing the Kali Gandaki to estimate flow depth and speed. If hydrological gauges are unavailable or destroyed, these videos provide an important source of information to assess the magnitude of these extreme events.


Bricker, J.D., Schwanghart, W., Adhikari, B.R., Moriguchi, S., Roeber, V., Giri, S. (2017): Performance of models for flash flood warning and hazard assessment: the 2015 Kali Gandaki landslide breach in Nepal. Mountain Research and Development, 37, 5-15. [DOI: 10.1659/MRD-JOURNAL-D-16-00043.1]

Chimaps in a few lines of code (1)

Posted on Updated on

Since Willett’s et al. (2014) paper on the growth and decay of drainage basins, chimaps have become a popular tool to illustrate the direction of divide migration. Whether a drainage basin is a victim or aggressor, i.e., will loose or gain catchment area in the future, is indicated by the variable chi in the vicinity of the divides. chi (Χ) is a transformed horizontal coordinate measured from the outlet of a drainage basin. The transformation uses upslope area to linearize the concave upward profile for a well-chosen value of the mn-ratio. If all drainage basin outlets have the same elevation, then chi values should be roughly the same on both sides of the divides in a steady-state situation. In this case, no drainage basin reorganization should be expected. Yet, if chi values differ, then one should expect that divides migrate from the basins with lower chi values (i.e. the aggressors) towards those with higher chi values (i.e. the victims) (see Figure below).

Drainage basin reorganization and chi values (modified after Willett et al. 2014, Fig. 1).

Now how is it possible to calculate chimaps using TopoToolbox? It is actually pretty easy, but there is no ready-made function. And there is a reason for this. The description above shows that plotting chimaps involves several modular steps, each of which should be carefully done. These steps are:

  1. Load DEM and derive flow directions (FLOWobj).
  2. Choose a suitable threshold for stream initiation and derive a stream network (STREAMobj).
  3. Make sure that all outlets have the same elevation.
  4. Make sure that all drainagebasins are complete, i.e. that the DEM covers their entire extent. Remove incomplete drainage basins.
  5. Choose a suitable value for the mn-ratio.
  6. Calculate the transformed variable chi (Hint: Don’t use chiplot, but check the low-level function STREAMobj/chitransform!).
  7. Plot the results.

I am on field work in Nepal right now and I don’t have access to a suitable DEM. I will thus go through these steps in my next blog entries.


Willett, S. D., McCoy, S. W., Perron, J. T., Goren, L. and Chen, C.-Y. (2014): Dynamic Reorganization of River Basins, Science, 343(6175), 1248765. [DOI:10.1126/science.1248765].

Blog on MATLAB Recipes for Earth Sciences

Posted on


Martin Trauth, apl. professor of paleoclimate dynamics and colleague at the University of Potsdam has just published his new blog. Martin authored two books that are standard works for both beginners and advanced MATLAB users in the Earth sciences. MATLAB Recipes for Earth Sciences (Springer 2015, 4th edition) covers a diverse range of topics from statistics, time series analysis to image analysis and directional data. MATLAB and Design Recipes for Earth Sciences (Springer 2012) co-authored by Elisabeth Sillmann focuses on how to visualize Earth science data and includes numerous examples on the use of internet resources, on the visualization of data with MATLAB, and on preparing scientific presentations. I can highly recommend these books!

Martin has now started his blog as a platform for regular updates and additions to his books, and to keep you updated on ongoing teaching and research activities. It will be on my blog reading list from now on.

TTLEM paper out

Posted on

This is a quick note that our paper on TTLEM in Earth Surface Dynamics has been revised and is published now. Benjamin Campforts added some major updates that improve performance and resolved some bugs in a previous version. You can find the software together with TopoToolbox here. We will certainly report here on examples, applications and new releases soon.


Campforts, B., Schwanghart, W., Govers, G. (2017): Accurate simulation of transient landscape evolution by eliminating numerical diffusion: the TTLEM 1.0 model. Earth Surface Dynamics, 5, 47-66. [DOI: 10.5194/esurf-5-47-2017]

Issyk Kul’s ups and downs: water-filled troughs and peak flows around one of the world’s largest mountain lakes

Posted on Updated on

dsc_0060Text and photos by Angela Landgraf and Swenja Rosenwinkel

Lake Issyk Kul in the northern Kyrgyz Tien Shan is the second-largest mountain lake worldwide and a proposed site for an ICDP drilling project to examine the climatic and erosional history of Central Asia over the past few million years. The lake is endorheic today, but different proxies indicate changing open- and closed-basin conditions through time. Shorelines and lacustrine deposits furnish ample evidence for major lake-level fluctuations during the Quaternary, and suggest the lake was connected to neighboring basins during lake high-stands. Better understanding the dynamics of spill and fill are paramount for correctly assessing changing environmental conditions in the lake basin. A special geomorphic puzzle is the appearance of massive lakebeds west of the present-day sill in the neighboring Kok Moinok basin. The western extent of these lakebeds coincides with a major knickpoint and an abrupt transition from a wide alluviated valley to a narrow bedrock gorge (the Boam Gorge).

In our recently accepted paper in Earth surface Processes and Landforms, we use a variety of field-, laboratory-, and modeling-based methods to assess potential linkages between natural damming to lake highstands and related outburst flood scenarios. These methods include geochronometry of lacustrine and related delta and fluvial deposits around the lake and at the gorge outlet, measurements of entrained mega-clasts for hydraulic paleoflood analysis, and paleoflood modeling using highstand-scenarios based on the elevation of mapped paleoshorelines. We refer to the full paper for a detailed discussion of these methods, and their specific links and assumptions.

We find that one or several catastrophic floods occurred through the Boam gorge in the late Pleistocene. A temporal succession of Issyk Kul’s lake-level drop and boulder deposition at the outlet supports a link between both. Paleoflood modeling, however, shows that catastrophic lake outbursts unconnected to Issyk Kul, i.e., from a separated Kok Moinok lake, could have sustained the necessary peak discharges for moving the boulders at the gorge exit. Although the overall geomorphic and sedimentary evidence around Issyk Kul records some of the largest catastrophic outburst floods in the Tien Shan mountains, if not Central Asia, direct links to documented lake-level changes of Issyk Kul remain elusive.


Rosenwinkel, S., Landgraf, A., Korup, O., Schwanghart, W., Volkmer, F., Dzhumabaeva, A., Merchel, S., Rugel, G., Preusser, F. (2017). Late Pleistocene outburst floods from Issyk Kul, Kyrgyzstan? Earth Surface Processes and Landforms, accepted. [DOI: 10.1002/esp.4109]