# Calculating Topographic Wetness Index (choosing from different algorithms)

Topographic Wetness Index can be expressed as

`````` Ln(a/tanB) based on the idea of Beven and Kirkby (1979)
``````

where

``````  a is the specific catchment area (a=A/L, catchment area (A)divided by contour length(L))
``````

and

``````  tanB is the slope
``````

The basic idea here is simple, but as there are multiple ways to calculate both a and tanB, the results of a TWI can vary widely (Qin et al. 2011).

Flow accumulation and catchment area can be calculated, for example by:

`````` D8 (O'Callaghan, J.F. / Mark, D.M. (1984))
D-infinity  (Tarboton, D.G. (1997)
Triangular Multiple flow direction (Seibert, J. / McGlynn, B. (2007)
``````

algorithms, and there are many other algorithms available too.

Slope is usually calculated as local slope around the pixel (Sorensen et al. 2005). Local slope can also be calculated as minimum, mean and maximum slope around the pixel. Another way to calculate slope is presented by Hjerdt et al. 2004 where slope is calculated to a point d meters below cell center.

Slope is a basic tool in most of the GIS softwares, however the calculation can differ. Here are few examples: ESRI: http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Calculating_slope SAGA: http://sourceforge.net/apps/trac/saga-gis/wiki/Terrain%20Analysis%20-%20Morphometry%20module%20library

As you can see there are many options available to calculate both a and tanB. So, the question is, in practice, which is the proper (best) way to calculate TWI using different these algorithms? Or is there any?

I personally like working in SAGA, mainly because there are a large selection of open source hydrology tools.

P.s. I'm having a hard time to find out exactly how catchment slope is calculated in Saga GIS, and exactly what does it mean here. (Terrain analysis -hydrology: catchment area parallel).

EDITED: Answered by Volker Wichmann from SAGA Forums: "The catchment slope output grid of the Catchment Area (Parallel) module is computed like this: for each cell, the local slope is calculated using the approach of Zevenbergen & Thorne. These slope values are accumulated downslope. Finally, for each cell the accumulated slope values are divided by the derived catchment area of the cell. The unit of the grid are radians."

"The Topographic Wetness Index (TWI) module requires a normal slope grid as input. "

References:

Beven and Kirkby 1979. A physically based variable contributing area model of basin hydrology. Hydrological Sciences Bulletin, 24, pp. 43–69.

Hjerdt et al. 2004. A new topographic index to quantify downslope controls on local drainage. Water Resources Research, 40, W05602, doi:10.1029/2004WR003130.

O'Callaghan, J.F. and Mark, D.M. 1984. The extraction of drainage networks from digital elevation data.Computer Vision, Graphics and Image Processing, 28:323-344

Qin et al. 2011. An approach to computing topographic wetness index based on maximum downslope gradient. Precision Agric 12:32–43.

Seibert, J. and McGlynn, B. 2007. A new triangular multiple flow direction algorithm for computing upslope areas from gridded digital elevation models, Water Ressources Research, Vol. 43, W04501

Sorensen et al. 2005. On the calculation of the topographic wetness index: evaluation of different methods based on field observations. Hydrol. Earth Sys. Sci. Discuss., 2, 1807–1834

Tarboton, D.G. 1997. A new method for the determination of flow directions and upslope areas in grid digital elevation models, Water Ressources Research, Vol.33, No.2, p.309-319

I think this post in the SAGA GIS forum might prove useful in answering your question about how slope is calculated:

Also, based on my understanding of TWI (as a PhD hydrology student involved in hydrologic modeling), the D-Inf (Tarboton), MFD-md (Qin), DEMON (Costa-Cabral), and MFD (Quinn) with the exponent p=1.1 (Freemann) are the best options for determining the accumulation area 'a' in the TWI calculation.

I think Sorensen's work, and Qin's work, lends credence to my own semi-professional opinion. However, Qin's improved algorithm (MFD-md) hasn't been as throughly tested and used as much as the others have.

When I have used SAGA to calculate TWI, I first calculate the slope using the slope commend in the Terrain Analysis / Morphometry / Slope, Aspect, Curvature module, using the default algorithm mentioned in the forum post. Then I calculate catchment area using the Terrain Analysis / Hydrology / Catchment Area / Catchment Area (Parallel) choice, using either the MFD algorithm or the D-INF algorithm with 1.1 as the convergence factor (p=1.1, from Freeman).

Then I run the TWI option, in Terrain Analysis / Hydrology / Topographic Indices / Topographic Wetness Index (TWI), with the options to convert area to "1/cell size" and use standard calculation. I convert to specific catchment area because this is what the original formulation by Beven and Kirkby called for. As for the difference between the "Standard" and "TOPMODEL", I am not sure what it is - looking into that myself right now.

I forgot to add, this all assumes that the DEM has been preprocessed by filling sinks. That is another complicated topic, with two (main) separate options.

I hope this helps!

Tom

P.S. Edits for @reima:

This is something I have only recently dug into, and I can admit I don't think I've reached the bottom yet! I prefer the method of Lindsay and Creed, the minimum impact approach that chooses breeching or filling based on minimizing total topographical impact (officially named "Impact Reduction Algorithm" - IRA) which I thought was implemented in his Terrain Analysis Software tool (formerly TAS, now WhiteBox GAT - link: http://www.uoguelph.ca/~hydrogeo/Whitebox/).

However, even his tool seems to implement other filling schemes:

1. The sink/depression filling algorithm (basic, but insanely fast) by Wang and Liu (2006) - which I don't believe operates in the IRA manner, but similar to the way ArcMap fills sinks/depressions, just straight up without any breeching.

2. And the sink/depression filling of Planchon and Darboux (2001), which floods a DEM and then removes the water bit by bit - it can enforce a slope on the filed area, which I think might improve TI calculations.

ArcMap has a new "de-pitting" add-on (http://blogs.esri.com/esri/arcgis/2013/03/05/optimized-tool-for-dem-pit-removal-now-available/) that seems similar to Lindsay and Creeds IRA, but I haven't read the cited paper yet to determine how similar. This method might be worth a look.

I'm also interested in scrutinizing my assumption that TI calculations need filled DEMs. I have three different sized watershed DEMs (<100 sq km, 100-1000 sq km, >1000 sq km), clipped using a shape file from 10 m NED data. These are not filled, since the shape file already provided the watershed delineation. I am going to run the SAGA GIS TI calculation (MFD, p=1.1) on all three watersheds, on both filled and unfilled DEMs, using ArcMaps filling scheme (old and new), and the Wang and Liu algorithm (in Whitebox, maybe in SAGA), and the Planchon and Darboux algorithm (in Whitebox, maybe in SAGA). I will also be calculating the TI values using the TI calculation embedded in my hydrological model.

If you want, I can share these results with you. I might not have them for a month or so though, as I have other more pertinent research that my focus is currently on, but I need to refine my TI calculation process by mid May at the latest.

• Thanks for the entry Tom! It seems we have almost identical work flow. Preprocessing of DEM is indeed another topic, especially when there are roads, bridges etc. corrupting the natural flow direction and accumulation. May I ask you what is your default method for filling sinks and how do you handle lakes in the preprocessing of DEM? I ended up removing all the lakes because every fill sink option tended to raise these water baisins higher than they actually are, causing information loss near the lakesides. Apr 8, 2014 at 7:33
• @reima, see above. I don't do anything special for lakes, as of yet. All of my focus is on predicting discharge at the stream/river outlet - I haven't yet wrapped my mind around how lakes will interact with all this. Also, I feel like I should make clear, in the TI calculation you DON'T want "Catchment Slope". You want the normal slope - just want to make sure that was clear from the SAGA GIS forum link I provided. Apr 9, 2014 at 20:28
• Yes, (I tried italicizeing the part about local ("normal")slope but I'll guess it's left a little bit unclear(my bad)).I have ended up using the Wang Liu 2006 fill-method so far, but if I continue with this topic, I'll definitely have to dig deeper too. I'm very intrested about your results too! Apr 10, 2014 at 11:47