dem - Getting Negative Topographic Wetness Index (TWI) values in SAGA GIS?

I am getting negative topographic index wetness values in SAGA GIS - anyone else every get this?

Not all negative values, but the minimum of the TWI range is negative for some steep watersheds we are examining. Ranges get as negative as (approx.) -1.5.

For this to happen, the value resulting from (specific catchment area / tan slope ) must be less than 1, since the ln of values between 0 than 1 = a negative value.

It is a very steep watershed, so we have slopes that approach 1.57 radians (almost 90 degrees). In fact, the maximum slope for the watershed in question is 1.509 radians (or 86.459 degrees). And the minimum Specific catchment Area value is 2.2175m. So we can see that if we calculate ln(2.2175/tan 1.509), which further simplified is ln(2.2175/16.16159) and then ln(0.137208), we get a negative value (-1.8 or something) that shows the theoretical absolute minimum TI value for our watershed. Its theoretical cause the steepest slope and smallest SCA don't combine, but some set of close values do, to get a real minimum of approx. -1.5.

To get more specific, in our steep watershed, the DEM has a resolution of 8.86 m and the Total Catchment Area minimum calculated by SAGA is this cell size squared (78.4996 sq m). Which makes sense, as there are cells in the catchment that don't accept flow from any other cells and therefore have a Total Catchment Area equal to the cell area. Now, for these cells, at localized peaks in the watershed (or even maximum peaks at the edges), the Specific Catchment Area should be equal to the Total Catchment Area for these cells divided by the contour width of downslope flow paths. These peaks flow to all neighboring cells, this 78.4996 is divided by (as currently coded in SAGA) 35.4 (which is 8 * 0.5 * 8.86). This gives the Specific Catchment Area minimum of 2.2175m.

My main problem with this negative TI result, is that I don't think it is possible, or at least have never really seen negative TI values talked about/graphed in a paper.

Any thoughts form anyone?

Answer

This is a good question, and one that I tend to get asked from time to time. First, as you've pointed out, the equation for TWI = ln(a / tan(B)), where a is the 'specific' catchment area (i.e. the upslope inflowing area normalized for a measure of contour length) and B is the slope gradient, in radians, at the grid cell. As you correctly pointed out TWI will only take on negative values where a < tan(B) since the natural log of any value less than 1 will be negative. So, negative TWI values should be surprising. Remember, TWI is just an index. All it is telling you is that two grid cells that have similar TWI values will likely become saturated at the surface under similar basin moisture conditions given their combined contributing areas and local slope gradients. It's also giving you a relative ranking of when grid cells will become saturated, such that cells with higher values are likely to become saturated before cells with lower values. So, those negative TWI values are really among the driest locations in the basin. The real genius of Beven and Kirkby's TWI is that if you hook this up to a simple hydrological model (e.g. TOPMODEL), which then keeps track of the basin's current moisture conditions, you can relate the current extent of the surface saturated area to a threshold TWI value. All grid cells with a TWI greater than this threshold value define the current surface saturated area which will be predicted to experience overland flow and possibly return flow.

For a particular DEM, you can actually figure out what conditions you'll likely find these negative TWI values. For a DEM of resolution g, the lowest possible a value will be equal to (or approximately equal to depending on how the TWI implementation determines contour length) g (i.e. g^2 / g). We can call this minimum value amin. Grid cells with this a value occur where there are no inflowing cells, which will only really happen along topographic ridge lines. Interestingly, these are the exact places that are likely to experience high local slope gradients as well. So, now, you can calculate the slope threshold at which point TWI starts to take on negative values as: Bthresh = atan(amin). Since you have a DEM resolution of 8.86 m, then any grid cells with this a value that have a slope greater than 83.56 degrees will take on negative TWI values. You can see why many people believe that TWI can't be negative; that is a fairly high slope and the conditions under which you would get a negative TWI are somewhat rare. You're lucky then to be working in such an interesting environment. (NOTE: TWI was never intended to be used across the very wide range of landscapes that it is commonly applied to today.)

In the TWI raster above, I've highlighted each of the grid cells with negative TWI. As you can see, it isn't all that common in the grid. In fact, in this case, it is completely confined to grid cells associated with off-terrain objects (buildings in this LiDAR DEM) because of their very steep sides and small a values. This particular case isn't very realistic from a process point of view (what does the wetness index of a building mean?), but nonetheless, it highlights my point well.

In his comment, Jeffery Evans raises the possibility that the reason for negative values is because of the way that zero slope values are handled. In fact, these flat areas will have very large TWI values rather than negative values. Flat areas, which also tend to be very low-lying valley bottom positions, are among the wettest and first to saturate locations in a basin. Now obviously if B is zero, tan(B) is also zero and you get a division by zero, which of course is infinity. This is like saying, if you have a perfectly flat location, it's always going to be saturated since it cannot shed water to a downslope position. Again, remember that TWI is meant to be a relative index of wetness based on topography, we're not saying that a cell with a TWI value of 2 is twice as saturated as a cell with a TWI of 1. However, obviously we can't have a grid with infinity values, and depending on the quality and source data of your DEM, zero gradients are quite likely. So it's true that developers of TWI tools need to have special rules in place to handle grid cells with zero local gradient. Generally you'll simply substitute in a 'really high value' to let the user know this site is quite wet. The problem with this approach is that, just as you've asked why you get negative TWI values, users will commonly say, 'why am I getting TWI values as high as 9999?'. So an alternative is simply to stick a NoData value in these cells.

I hope that helps to answer you questions.