arcgis desktop - 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

Answer

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

https://sourceforge.net/p/saga-gis/discussion/354013/thread/27ecdc6b/

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.

Reference page 106 onward of the linked pdf for more specific help: http://sourceforge.net/projects/saga-gis/files/SAGA%20-%20Documentation/SAGA% 20Documents/SagaManual.pdf/download

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.