Soil moisture modelling using TWI and satellite imagery in the Stockholm region

(1)

Soil moisture modelling using TWI and satellite imagery in the

Stockholm region

Jan Haas

Master’s of Science Thesis in Geoinformatics TRITA-GIT EX 10-001

School of Architecture and the Built Environment Royal Institute of Technology (KTH)

100 44 Stockholm, Sweden

March 2010

(2)

TRITA-GIT EX 10-001 ISSN 1653-5227

ISRN KTH/GIT/EX--10/001-SE

(3)

i

Abstract

Soil moisture is an important element in hydrological land-surface processes as well as land- atmosphere interactions and has proven useful in numerous agronomical, climatological and meteorological studies. Since hydrological soil moisture estimates are usually point-based measurements at a specific site and time, spatial and temporal dynamics of soil moisture are difficult to capture. Soil moisture retrieval techniques in remote sensing present possibilities to overcome the abovementioned limitations by continuously providing distributed soil moisture data at different scales and varying temporal resolutions. The main purpose of this study is to derive soil moisture estimates for the Stockholm region by means of two different approaches from a hydrological and a remote sensing point of view and the comparison of both methods.

Soil moisture is both modelled with the Topographic Wetness Index (TWI) based on digital elevation data and with the Temperature‐Vegetation Dryness Index (TVDI) as a representation of land surface temperature and Normalized Difference Vegetation Index (NDVI) ratio.

Correlations of both index distributions are investigated. Possible index dependencies on vegetation cover and underlying soil types are explored. Field measurements of soil moisture are related to the derived indices.

The results indicate that according to a very low Pearson correlation coefficient of 0.023, no linear dependency between the two indices existed. Index classification in low, medium and high value categories did not result in higher correlations. Neither index distribution is found to be related to soil types and only the TVDI correlates alongside changes in vegetation cover distribution. In situ measured values correlate better with TVDIs, although neither index is considered to give superior results in the area due to low correlation coefficients. The decision which index to apply is dependent on available data, intent of usage and scale. The TWI surface is considered to be a more suitable soil moisture representation for analyses on smaller scales whereas the TVDI should prove more valuable on a larger, regional scale. The lack of correlation between the indices is attributed to the fact that they differ greatly in their underlying theories.

However, the synthesis of hydrologic modelling and remote sensing is a promising field of research. The establishment of combined effective models for soil moisture determination over large areas requires more extensive in situ measurements and methods to fully assess the models’ capabilities, limitations and value for hydrological predictions.

Keywords: TWI, TVDI, soil moisture

(4)

ii

Acknowledgement

I would like to express my sincere gratitude to my supervisors Ulla Mörtberg and David Gustafsson, Department of Land and Water Resources Engineering, School of Architecture and the Built Environment, KTH ‐ Royal Institute of Technology for their continuous help and guidance, data and literature provision as well as conducting fieldwork together and data post‐

processing.

I also gratefully acknowledge the support and advice given to me by Dr. Yifang Ban, Head of Department of Urban Planning and Environment, School of Architecture and Built Environment, KTH ‐ Royal Institute of Technology. I would also like to express my gratitude towards Tuong Thuy Vu, Dorothy Furberg and Maria Irene Rangel Luna of the same department for their help, trust and technical support.

Furthermore I would like to thank Else‐Marie Wingqvist from Sveriges Meteorologiska och Hydrologiska Institut (SMHI) for the provision of detailed temperature and precipitation data for several gauging stations in and around Stockholm.

(5)

iii

List of Tables

Table 1 Technical specifications of Landsat missions 1 to 7 ... 29

Table 2 Summary of Landsat 7 ETM+ satellite sensor characteristics ... 29

Table 3 Landsat 7 ETM+ band characteristics ... 30

Table 4 Landsat scene information ... 36

Table 5 Land cover classification table ... 37

Table 6 Vegetation class definition table ... 38

Table 7 Summary of soil types and distribution ... 38

Table 8 Summary of day mean temperatures and precipitation ... 39

Table 9 DEM surface statistics before and after pit removal ... 40

Table 10 Band constants for radiance calculation from DNs ... 46

Table 11 Classes and index surfaces for combination ... 49

Table 12 Overview of all surface combinations ... 49

Table 13 Surface classification summary ... 50

Table 14 TWI surface statistics summary ... 55

Table 15 Correlation between TVDI/TWI for vegetation classes ... 59

Table 16 Correlation between TVDI/TWI for soil classes ... 59

Table 17 Correlation between TVDI/TWI for vegetation and soil classes ... 60

Table 18 TVDI and TWI surface comparison statistics ... 62

Table 19 Numerical comparison between ground truth and calculated TWI and TVDI values ... 63

(9)

vii

List of Figures

Figure 1 Simplified representation of the Ts/NDVI space from Lambin and Ehrlich (1996) ... 33

Figure 2 Ts/NDVI space as TVDI taken from Sandholt et al. (2002) ... 34

Figure 3 Geometrically corrected Landsat 7 ETM+ satellite image ... 36

Figure 4 Flow direction angle definition in the Dinf approach taken from Tarboton (1997) ... 41

Figure 5 Final TWI surface ... 42

Figure 6 Cloud mask ... 45

Figure 7 Haze mask ... 45

Figure 8 Final TVDI surface ... 48

Figure 9 In situ data collection site Törnskogen ... 51

Figure 10 Sample point overview ... 51

Figure 11 Unfilled DEM (pits and sinks visible) ... 52

Figure 12 Filled DEM (pits and sinks removed) ... 52

Figure 13 TWI surface without prior pit removal ... 52

Figure 14 Slope for TWI calculation in degrees ... 54

Figure 15 Slope for TWI calculation in percent ... 54

Figure 16 Slope calculated in TauDEM ... 54

Figure 17 Scatterplot of the complete TWI and TVDI surfaces with regression line ... 56

Figure 18 Visual comparison between ground truth and calculated TWI and TVDI values ... 63

(10)

viii

List of Appendices

Appendix A Original 50 metre resolution DEM as acquired from Lantmäteriet ... 78

Appendix B Precipitation over Stockholm in mm during July 2002 ... 79

Appendix C Average temperature around Stockholm in °C during July 2002 ... 80

Appendix D D∞ flow direction grid surface for TWI calculation within TauDEM ... 81

Appendix E D∞ slope grid surface for TWI calculation within TauDEM ... 82

Appendix F D∞ contributing area surface for TWI calculation within TauDEM ... 83

Appendix G NDVI surface calculated by Landsat 7 ETM+ bands 3 and 4 ... 84

Appendix H Temperature in °C derived from the ETM+ thermal band on August 4th, 2002 ... 85

Appendix I TWI surface calculation variants for 'NoData' ... 86

Appendix J TWI surface calculation variants for 'unfilled' ... 87

Appendix K TWI surface calculation variants for 'filled' ... 88

Appendix L TWI/TVDI scatterplots for vegetation classes 0 to 7 ... 89

Appendix M TWI/TVDI scatterplots for soil classes ... 90

Appendix N TWI/TVDI scatterplots for soil class bog and vegetation classes 0 to 7 ... 91

Appendix O TWI/TVDI scatterplots for soil class clay and vegetation classes 0 to 7 ... 92

Appendix P TWI/TVDI scatterplots for soil class rock and vegetation classes 0 to 7 ... 93

Appendix Q TWI/TVDI scatterplots for soil class sand and vegetation classes 0 to 7 ... 94

Appendix R TWI/TVDI scatterplots for soil class till and vegetation classes 0 to 7 ... 95

Appendix S TVDI surface combination statistics... 96

Appendix T TWI surface combination statistics ... 97

Appendix U TWI/TVDI and TVDI/TWI quantile classification scatterplots ... 98

Appendix V TWI/TVDI and TVDI/TWI natural breaks classification scatterplots ... 99

(11)

ix

List of Acronyms

ASCII ‐ American Standard Code for Information Interchange ATCOR ‐ ATmospheric CORrection

BFD ‐ BiFlow Direction CMP ‐ Common MidPoint CSV ‐ Comma‐Separated Values

DEMON ‐ Digital Elevation MOdel Networks DN ‐ Digital Number

EOS ‐ Earth Observation Satellite ETM ‐ Enhanced Thematic Mapper

FAPAR ‐ Fraction of Absorbed Photosynthetically Active Radiation FOV ‐ Field‐Of‐View

GA ‐ Genetic Algorithm GCP ‐ Ground Control Points GLCF ‐ Global Land Cover Facility

IRS‐WiFS ‐ Indian Remote Sensing Satellites Wide Field Sensor LAI ‐ Leaf Area Index

LiDAR ‐ Light Detection And Ranging MFD ‐ Multiple Flow Direction MIR ‐ Mid‐InfraRed

MSS ‐ MultiSpectral Scanner

NASA ‐ National Aeronautics and Space Administration NOAA ‐ National Oceanic and Atmospheric Administration NDVI ‐ Normalized Difference Vegetation Index

NIR ‐ Near InfraRed

PBMR ‐ Push Broom Microwave Radiometer RBV ‐ Return Beam Vidicon

(12)

x RMSE ‐ Root‐Mean Square Error

SCA ‐ Specific Catchment Area SFD ‐ Single Flow Direction

SGU ‐ Sveriges Geologiska Undersökning

SMHI ‐ Sveriges Meteorologiska och Hydrologiska Institut SPOT ‐ Satellite Pour l'Observation de la Terre

SWEREF99 ‐ SWedish REerence Frame 1999

TauDEM ‐ Terrain analysis using Digital Elevation Models TDRSS ‐ Tracking and Data Relay Satellite System

TFD ‐ Tracking Flow Direction TI ‐ Topographic Index

TIN ‐ Triangulated Irregular Network TIR ‐ Thermal InfraRed

TOPMODEL ‐ TOPography‐based hydrological MODEL TVDI ‐ Temperature‐Vegetation Dryness Index

TWI ‐ Topographic Wetness Index USGS ‐ United States Geological Survey VWC ‐ Vegetation Water Content WI ‐ Wetness Index

(13)

11

1 Introduction

Knowledge about soil moisture status and the processes defining moisture distribution is, and has always been essential in many respects. Starting to be of importance with the advent of agriculture in human history, information about soil moisture plays nowadays an important role in a much broader field, e. g. hydrology, meteorology and climatology, ecology, land surface modelling and most recently studies on global environmental changes (Verstraeten et al., 2006;

Kerr et al., 2001; Kerr, 2007 and Gruhier et al., 2008). As environmentally related issues increase, the demand for information on environmental parameters grows simultaneously. To gather and derive these is therefore a crucial task.

The understanding of hydrological processes and accurate determination of soil moisture has been subject of numerous studies in the field of hydrology. Soil moisture is usually expressed as indices in hydrological modelling. One broadly used index is the Topographic Wetness Index (TWI) originally developed by Beven and Kirkby (1979). It is based on surrounding topography and describes the tendency of a point to become saturated. Air‐ and spaceborne soil moisture retrieval has been a promising research field in remote sensing since the early 1970s. Relations between measured land surface temperature, vegetation and soil moisture could be found and expressed by the ‘triangle’ method. The Temperature‐Vegetation Dryness Index (TVDI) by Sandholt et al. (2002) is an index representing the abovementioned relations. Both indices have shown to be reasonable estimators of surface soil moisture, although some questions of validity remain. However, there seems to be a lack of knowledge about index correlations and if one particular index can be considered to give superior results in the region of interest.

One purpose of this study is to compare both index distributions to each other with respect to different soil types, varying vegetation cover and to explore possible correlations and dependencies among these. Furthermore, the results and methods are related to in situ measured soil moisture. Another purpose is to gain insight into current soil moisture distribution for the region of Stockholm that can be used in further analyses, research and might help in decision making and future planning.

(14)

12

2 Background

2.1 Topographybased Hydrological Model (TOPMODEL)

The TOPMODEL framework initially proposed by Beven and Kirkby (1979) is a compilation of different concepts for distributed hydrological modelling. It is based on the two fundamental assumptions that downslope subsurface flows can be adequately represented by a succession of steady state water table positions and that there is an exponential relationship between local storage (or water table level) and downslope flow rate. Furthermore it is assumed that the hydraulic gradients of the saturated zone can be approximated by local surface topographic slopes tan β , thus requiring a detailed analysis of catchment topography (Quinn, 1991 and Beven, 2004).

2.2 Topographic Wetness Index (TWI)

The TWI is a mathematically simple parameterisation of soil moisture status that has been applied in numerous studies. The index is based and calculated on slopes and depends therefore on digital terrain data. The TWI is firstly suggested within the TOPMODEL framework by Beven and Kirkby (1979) and can be calculated as:

⁄ (1)

where α = the cumulative upslope area draining through a point (per unit contour length)

tan β = the slope angle at that point

The index describes the tendency of water to accumulate at any point in a catchment (in terms of α) and the tendency for gravitational forces to move that water downslope (expressed in terms of tan β as an approximate hydraulic gradient) (Quinn, 1991). Thus, wet areas can arise either from large contributing drainage areas or from very flat slopes, whereas areas with low index values are drier, resulting from either steep slopes or small contributing drainage areas.

TOPMODEL is able to predict both distributed water tables and hydrographs. Modification of the local water table depth can give an estimate of local soil moisture status. Calculated local water table depths are dependent on the value of the index at a point (Quinn et al., 1995).

(15)

13

Statistical distributions and spatial patterns of the TWI vary greatly depending on calculation parameters (Quinn et al., 1995). The two most critical issues here are the resolution of the underlying elevation model and the choice of the appropriate flow routing algorithm. For successful TWI calculations it is assumed that any problems involved in DEM creation, i.e. sink and dam removal are resolved. Quinn et al. (1995) also noted that there is no standard solution for calculating the ln α tan β⁄ index but that its parametres need to be optimized according to the corresponding catchment.

Limitations and issues of validation and predictive uncertainty of the concept were investigated, detailed analyses of the model’s capabilities, underlying theories and basic assumptions were performed and improvements were suggested in several works (Franchini et al., 1996; Ambroise et al., 1996; Beven, 1997; Mendicino and Sole, 1997; Seibert et al., 1997; Brasington and Richards, 1998; Güntner et al., 1999 and Sørensen et al., 2006).

2.3 Digital Elevation Data

Digital elevation data is traditionally represented in three forms; gridded data, triangular irregular networks and contour based elevation models. All of them can be used in hydrological modelling.

2.3.1 Gridded data (DEM)

Digital elevation models are readily available and commonly used as representation of elevation information. Grid DEMs store height information in matrix form. Each matrix node of equal resolution is assigned one particular elevation value. Grid based DEMs are widely used in analyses of hydrologic problems but are disadvantageous in some respects. They lack e. g. the possibility of sufficiently representing elevation discontinuities. Another disadvantage is that mesh resolution affects the results and computational efficiency. Furthermore it is stated that grid spacing needs to be based on the roughest terrain in the catchment which results in redundancy in smoother areas and finally that the computed flow paths tend to zigzag instead of to follow drainage lines, thus becoming too long. Holmgren (1994) noted that grid based DEMs were relatively unsuitable for hydrological modelling before the development of more sophisticated flow direction algorithms, as runoff could only be modelled in a simplified way (limited to eight discrete directions).

(16)

14

2.3.2 Triangular irregular networks (TIN)

The representation of elevation data in TINs is based on the connection of nodes with stored 3‐

dimensiondal coordinates (x, y and z). The nodes are connected with edges to form triangles.

Additional information such as slope, aspect, surface length and area can then be derived for each of the resulting triangular facets. The TIN has become a rather unpopular format in recent years (Holmgren, 1994).

2.3.3 Contour Data

Contours are lines that connect points of equal (elevation) values. As the spacing between the contours widens, height differences become smaller and vice versa. Contour data is usually derived from digitizing of maps and represents height distribution continuously. Contour based models are superior to other forms of elevation data because they explicitly state in which direction the runoff flows (always perpendicular to the contours) (Holmgren, 1994). Aryal and Bates (2008) found that contour DEMs are less sensitive to changes in resolution than grid DEMs and provide more consistent estimates of TWIs at moderate to high DEM resolutions.

The quality of elevation data used in the analyses is of great importance. Yet, according to Quinn et al. (1995), Wolock and McCabe (1995), Beven (1997) and Mendicino and Sole (1997), there is no standard and superior set of digital elevation data as basis for TWI calculation.

2.3.4 DEM resolution

Selecting the optimal spatial resolution is a central issue in environmentally related analyses including digital elevation data (Rodhe, 1999; Aryal and Bates, 2008). The appropriate resolution should be chosen according to the modelling goals but is unfortunately often determined by data availability and economic aspects. The spatial resolution should always be modelled to fit the scale of the area under consideration (Vaze and Teng, 2007; Wu et al., 2007).

For example smaller areas with more detailed information require higher spatial resolutions. If the resolution becomes too low and hillslope lengths are comparatively small, no meaningful TWI can be derived (Beven, 1997). Studies show that highest resolutions do not always produce the most accurate prediction model (Wu et al., 2008). The studies of Refsgaard (1997), Molnár and Julien (2000) and Vázquez et al. (2002) highlighted the variability in hydrologic response to grid size and suggested an appropriate resolution as a trade‐off between minimized computational effort and retention of realistic model performance. Endreny and Wood (2003) found that runoff flows differently with varying spatial resolutions. Many studies found the TWI

(17)

15

concept to respond sensitively to grid cell resolution (Wolock and Price, 1994; Gallant and Hutchinson, 1996; Saulnier et al., 1997; Brasington and Richards, 1998; Higy and Musy, 2000;

Vivoni et al., 2005; Vaze and Teng, 2007; Wu et al., 2007; Wu et al., 2008), even in resolutions higher than 10 m (Hancock, 2005). Rather coarse data may not represent some important slope features whereas a too high resolution can introduce unwanted perturbations to flow directions and slope angles. This coincides with the findings of Wolock and Price (1994), Zhang and Montgomery (1994), Bruneau et al. (1995) and Thieken et al. (1999), who stated that low resolution DEMs generally show smoother terrains and shorter flow paths by statistics of slopes and topographic indices than high resolution data. Vaze and Teng (2007) found that maximum and mean slope values drastically decrease with a decrease in DEM resolution. TWI distribution is dependent on the degree of spatial resolution of elevation data. The findings in several studies show that a lower grid cell resolution results in higher TWI mean values and a broader range of smaller slope percentages (Hutchinson and Dowling, 1991; Jenson, 1991; Quinn et al., 1991 and 1995; Wolock and Price, 1994; Zhang and Montgomery, 1994; Saulnier et al., 1997; Wolock and McGabe, 1995; Wolock and McGabe, 2000; Hjerdt et al., 2004; Sørensen and Seibert, 2007; Wu et al., 2007 and 2008 and Aryal and Bates, 2008). Wu et al. (2008) additionally noticed that the standard deviation of plan curvature decreased with increasing DEM cell size.

Not only directly derived elevation attributes are dependent on resolution, but also the form of the index distribution (Iorgulescu and Jordan, 1994; Wolock and Price, 1994; Zhang and Montgomery, 1994; Bruneau et al., 1995; Quinn et al., 1995; Franchini et al., 1996; Saulnier et al., 1997 and Mendicino and Sole, 1997). Wolock and Price (1994) studied the effects of topography on watershed hydrology. They showed that DEM resolution has an effect on water table depth prediction, the ratio of overland flow to total flow, peak flow and variance and skew of predicted stream flows. It was found that predicted mean depths to water tables decreased and maximum daily flow increased with coarser DEMs. Wolock and McGabe (2000) compared 100 m and 1000 m resolution DEMs and found that slopes calculated from 1000 m resolution DEMs are generally smaller, whereas specific catchment area (SCA) and the TWI are larger compared to 100 m resolution DEMs. The reasons for these changes are attributed mainly to terrain discretization effects (dividing the terrain into a different amount of grid cells than initially present) and are found to be independent of underlying terrain types. Wu et al. (2007) examined the effects of elevation data resolution on the performance of topography based watershed runoff simulations.

Twelve DEMs of different grid sizes were compared in terms of their influence upon TWI distribution. It was found that the total watershed discharge increases continuously as resolutions decrease. Sørensen and Seibert (2007) also found decreasing variations in TWI between neighbouring cells with an increase in cell size. By increasing the resolution, TWI distributions turned out to become rather similar, but quite different patterns were to be

(18)

16

expected as well. Brasington and Richards (1998) examined frequency distributions of slope, tan , upslope contributiong area and TWI on DEMs with varying resolutions from 20 m to 500 m. They found a significant distribution dependency on grid size with most prominent changes in DEM resolutions from 100 m to 200 m. Saulnier et al. (1997) based their study on digital terrain maps with resolutions ranging from 20 m to 120 m with 20 m increments from initially digitized 10 m contour lines. Further scale dependence in the calibrated value of the saturated conductivity parameter K0 could be identified.

Horizontal Resolution

Different DEM resolutions as basis for calculation of topographic parameters were explored and used in numerous studies with varying goals. The resolutions under consideration range from meter level up to several kilometres.

Murphy et al. (2009) used two different data sets: A 1 m resolution Light Detection And Ranging (LiDAR) DEM and a 10 m conventional photogrammetric DEM. Both DEMs showed good results in comparison to field‐mapped conditions. The LiDAR DEM produced slightly better results.

Warren et al. (2004) compared slopes on different DEM resolutions and found the highest resolution of 1 m to be most adequate. Freer et al. (1997) calculated the TWI for subsurface storm flow analysis on two small catchments with 1 m and 2 m DEM resolutions. Zhang and Montgomery (1994) and Wu et al. (2008) based their analyses on 2 m resolution DEMs. Wu and Archer (2005) used DEM grids of 4 m resolution in their study. Sulebak et al. (2000) related soil moisture to primary and secondary topographic attributes on a 5 m resolution DEM. Sørensen and Seibert (2007) examined the scale dependency of a decreasing DEM resolution upon the TWI by resampling a 5 m resolution LiDAR DEM to 10, 25, and 50 m respectively and observed TWI values to be considerably sensitive to different grid resolutions whereas the estimation of upslope areas is even more affected than calculated slopes. Thompson et al. (2001) compared terrain attributes and quantitative soil‐landscape models on different horizontal resolutions (10 m – 30 m) and vertical precisions (0.1 m – 1 m). It could be observed that a decrease in horizontal resolution resulted in lower slope gradients on steeper slopes, steeper slope gradients on flatter slopes, narrower ranges in curvature, larger SCAs in upper landscape positions and lower SCAs in lower landscape positions. Zinko et al. (2006) calculated ground water flow based on a 20 m resolution DEM. Kim et al. (2008) investigated the strength of correlations between different soil properties and found a TWI at a 25 m resolution as best. Vaze and Teng (2007) also suggested a grid cell resolution of 25 m or higher (up to 1 m) to capture the scale of surface processes. The resolution is however determined by landscape variations.

Tombul (2007) used a 30 m resolution DEM for a comparison of wetness index (WI) and topographic index (TI). Prasad et al. (2005) explored the relationship of hydrologic parameters

(19)

17

and nutrient loads using a 30 m resolution DEM with moderate success. Garbrecht and Martz (1993) showed that most drainage features can be produced with a 30 m resolution DEM.

Endreny and Wood (2003) examined the spatial congruence of observed and DEM‐delineated overland flow networks with a 30 m resolution DEM. Cai and Wang (2006) observed that there are no significant changes in TWI values when comparing 30 m and 90 m resolution DEMs in landscapes with high elevation variation. Kim and Jung (2003) found a pixel size threshold between 30 m and 50 m for stable and consistent computation of spatial flow distribution patterns. Lineback Gritzner et al. (2001) calculated the TWI with DYNWET on a 30 m resolution DEM. Soil moisture could not be related to landslide risk. The authors suggest that a too coarse spatial resolution might account for the lack of correlation. This would be consistent with the findings of Zhang and Montgomery (1994) stating that a resolution of 30 to 90 m could be too large for modelling complex topography. Wu et al. (2007) examined the influence of DEM resolution upon TWI distributions by means of comparisons between twelve different DEMs ranging from 30 m to 3 km. Conoscenti et al. (2008) based their studies on a 40 m resolution DEM. Lin et al. (2007) successfully explored the relationships of topography and spatial variations in groundwater and soil morphology using a 40 m resolution DEM. Bruneau et al.

(1995) considered even a resolution as low as 50 m to be sufficient if the model calibration is able to compensate aggregation effects. Quinn et al. (1995) considered grid sizes of 100 m or more as too coarse for meaningful TWI calculations. The depiction of the topographic form of individual hillslopes required higher resolutions. A direct effect on the calculation of accumulated contributing areas could be observed for different grid resolutions. Hutchinson and Dowling (1991) examined a 2.5 km resolution DEM. Jain et al. (1992) conducted their study on distributed modelling even on a 4 km resolution.

Overall, a 10 m horizontal resolution was found most suitable in many studies and is therefore suggested (Wolock and Price, 1994; Zhang and Montgomery, 1994; Irvin et al., 1997; Saulnier et al., 1997; Thieken et al., 1999; Chaplot and Walter, 2003 and Fei et al., 2007).

Vertical Precision

A sufficiently high vertical precision is needed in order to determine accurate spatial locations of channel segments and thus the channel network (Thieken et al., 1999). Chaplot and Walter (2003) used a 0.3 m vertical resolution and a 10 m horizontal resolution DEM in their studies and achieved good correlations between the soil wetness and topographic attributes at all depths. Nelson and Jones (1995) point out the necessity of roundoff error removal in low vertical resolution DEMs before delineation of stream networks and basins of a watershed. They presented a method of smoothing already rounded digital terrain models with the effect of restoring the smoothness of the natural terrain as well as the removal of flat regions. The results

(20)

18

show that their algorithm is able to restore the natural slope of large flat areas where the lack of elevation can be attributed to roundoff errors in the elevation data. Different studies show that a decrease in vertical precision does not have any significant influence on slope gradients or SCAs (Quinn et al., 1991; Wolock and Price, 1994; Zhang and Montgomery, 1994 and Gyasi‐Agyei et al., 1995). However, Thompson et al. (2001) identified various points with zero slope gradients and zero slope curvature alongside a decrease in vertical precision. Pan et al. (2004) found that errors in all compared flow direction algorithms increase as the vertical resolution decreases.

This suggests calculating the TWI on the highest vertical resolution possible.

2.3.5 Depression removal in DEMs

Digital elevation models often contain unrealistic local depressions that affect calculations of flow directions and flow accumulation. These depressions sometimes reflect the real situation but are often introduced during the generation of the DEM through interpolation errors. Müller‐

Wohlfeil et al. (1996) and Irvin et al. (1997) emphasized the need of getting rid of DEM depressions which are described as spurious low spots that are generally relics of the DEM generation. This confirms the findings of O’Callaghan and Mark (1984) who stated that pits or depressions are usually rare topographic features (e. g. natural holes or quarries) or absent in most terrain types, but not uncommon in digital elevation models. They are therefore considered as errors. Nelson and Jones (1995) also identified depressions and pits as spurious and attribute them to errors in elevation data. Beven (1997) also noted the common problem of sinks in flow pathway determination in digital terrain analysis. Qin et al. (2007) pointed out the necessity of DEM pre‐processing for removal of depressions and flat areas preceding TWI calculation. Present flow direction algorithms need depressionless DEMs as input lest the flow be stopped in a sink. O’Callaghan and Mark (1984) introduced the techniques of elevation smoothing (reducing the number of artificial pits by a weighted sum of a point and its eight neighbours) and iterative interior pit removal. Elevation smoothing only identifies and removes shallow pits whereas deeper sinks remain. Pit removal involves identifying the point of a basin where water would overflow from the pit basin. This point lies on the boundary of the pit with minimum elevation difference to the pit. Another approach is to fill the identified pits, i.e. e. to match one particularly low elevation value to the lowest value of the surrounding pixels. Jenson and Domingue (1988) refined pit removal implementing neighbourhood techniques and iterative spatial techniques. Planchon and Darboux (2001) present a new method for filling depressions in DEMs. This new approach inundates the whole DEM surface and removes the excessive water afterwards instead of gradually identifying pits and filling them. The process is

(21)

19

simple to implement with low computational complexity. Furthermore, it is possible to add slope to the filled sinks to overcome the limitations of flow direction algorithms in flat areas.

2.4 Flow direction algorithms

Flow direction algorithms decide about how flow propagates from one grid cell to other grid cells. Thus, they determine the total upslope area entering one grid cell and the effective contour length orthogonal to the flow direction. The cumulative upslope area α needed in the TWI can then be calculated from the relationship:

⁄ (2)

where α = cumulative upslope area

A = total upslope area

L = effective contour length

Hillslope geometry and flow partitioning greatly affect the spatial patterns of the TWI (Aryal and Bates, 2008). The mixture of convex and concave features (Kim and Lee, 2004) as well as errors in elevation data (Walker and Willgoose, 1999) complicate the essential choice of the appropriate flow direction algorithm, which greatly influences the calculation of upslope contributing areas, flow accumulation and other processes, e. g. sediment transport capacity (Wilson et al., 2007). Furthermore it is impractical to prefer one specific flow direction algorithm to another. The algorithm should rather be chosen according to the underlying topography. One particularly critical issue related to the determination of flow paths that has not been resolved satisfactorily yet is flow direction routing in flat areas with no or only little height difference. Existing flow direction algorithms have resolved this problem more or less successful but it remains a limitation to most approaches.

2.4.1 Single Flow Direction (SFD) algorithms

SFD algorithms restrict the surface and subsurface runoff from a single grid cell to only one other cell without considering any other surrounding cells. No flow bifurcation and dispersive flow is possible. Water only flows to the cell with the highest gradient (the steepest slope). SFDs have been criticized for the lack of not being able to take flow dispersion into account, e. g. by Aryal and Bates (2008). The SFD is furthermore incapable of returning the right flow path if the

(22)

20

steepest gradient might coincidentally fall between two of the eight cardinal and diagonal directions (Wolock and McCabe, 1995).

D8

The D8 (deterministic eight‐node) algorithm is the first method of calculating flow directions from one cell to one of its eight neighbours by O’Callaghan and Mark (1984). It is rather simple and does not account for dispersive flow. Because flow is only assigned to the surrounding cell with the steepest downwardslope gradient, grid bias is introduced (Tarboton, 1997). The D8 method works well in valleys but produces many parallel flow lines and problems near catchment boundaries.

Rho8

Fairfield and Leymarie (1991) suggested a new method to overcome the limitations of the D8 algorithm by introducing a probability distribution proportional to altitude difference (slope) between the centre and neighbouring cells. In their method, parallel flow paths are broken up and a mean flow direction equal to the aspect is generated (Wilson et al., 2000). The flow distribution pattern results in more cells with no upslope connections but is uniquely generated each time. The randomness of the probability function has been criticised by Tarboton (1997) because SCAs are deterministic quantities and should therefore be calculated in a stable, repeatable way.

2.4.2 Biflow Direction (BFD) algorithms

Two Directional EdgeCentred Routing

The two directional edge‐centred routing algorithm (2D‐Lea) by Lea (1992) creates a planar surface prior to flow direction determination. Elevation information is interpolated and stored on the edges of four surrounding DEM cells, creating a plane. Flow is then expressed as a flow vector following the route a ball would take when put on the centre of a cell on the plane. This method allows flow to be determined precisely and expressed continuously with an angle between 0 and 2π. Flow is then directed and split towards the two cells closest to the angle. The method lacks the possibility to disperse flow into more than two cardinal directions.

Two Directional BlockCentred Routing

Jensen’s two directional block‐centred routing algorithm (2D‐Jensen) is based on the same planar concept as in 2D‐Lea with the difference that not only four but nine cells are used to determine the plane which is defined by averaging cell centre values rather than edges (Jensen,

(23)

21

2004). Both the SFD and BFD algorithms are problematic in finding the right flow directions in flat areas (Pan et al., 2004).

2.4.3 Multiple Flow Direction (MFD) algorithms

The MFD algorithm was initially proposed by Quinn et al. (1991) to improve the representation of convergence and divergence of flow. Unlike the SFD and BFD algorithms, the MFD algorithm allows dispersive flow from one grid cell proportionately to multiple surrounding cells. As the D8 algorithm, the MFD algorithm is based upon finding the steepest slope to the surrounding eight cells. Downward elevation gradients are determined for each possible direction and multiplied with the corresponding contour lengths, respectively. This implies that flow is always divided proportionally to all surrounding cells. Tarboton (1997) noted that this could introduce unwanted errors because dispersive flow from one pixel is not necessarily directed towards all pixels in the vicinity with lower elevation.

Quinn et al. (1991) used a slightly different formula for the contour length calculation than Wolock and McGabe (1995) who first described the implementation of the MFD algorithm in the ARC/INFO framework. According to Wolock and McCabe (1995), MFDs tend to give more realistic‐looking spatial patterns of accumulating area than single‐directional flow algorithms where flow is expressed in distinct lines. Although the MFD algorithms can not overcome the limitation of flow determination in flat areas where slope values become rather small, they still provide better results than the SFD and BFD algorithms (Pan et al., 2004). Kim and Lee (2004) pointed out another limitation of the initial MFD algorithm. Dispersive flow is only possible towards eight discrete directions. This drawback has been resolved by the concept of a flow tube network within the Digital Elevation Model Networks (DEMON) algorithm by Costa‐Cabral and Burges (1994) which enables flow in continuous directions from 0 to2π. Another disadvantage of the MFD is flow overdissipation and braiding near channel pixels especially in convex‐shaped topography as has been pointed out by Quinn et al. (1995).

FD8

The limitation of D8 algorithm only being able to direct flow to one of eight surrounding cells has been overcome by the FD8 algorithm that allows catchment spreading and flow dispersion to be represented (Moore et al., 1993). The algorithm allows dispersive flow in upland areas (above defined channels) but uses the initial SFD D8 algorithm below channel initiation points.

Although the problems of parallel flow and single flow could be resolved, some issues still remain. Mendicino and Sole (1997) noted that the FD8 algorithm simulates flow poorly in certain topographic conditions (e. g. alluvial plains). In these areas a pronounced expansion of

(24)

22

surface runoff along the alluvial plains occurs instead of a well delineated stream channel network.

FRho8

As the FD8 is for the D8 algorithm, the FRho8 is an advancement of the Rho8 algorithm. FRho8 can be described as the combination of the MFD algorithm by Quinn et al. (1991) and the SFD algorithm Rho8. Firstly, the algorithm identifies all downslope neighbouring points to a channel initiation point and respectively calculates the slope gradients. Secondly, flow is distributed among the surrounding cells according to the Rho8 procedure. The upslope contributing and SCAs are then determined using the flow width and flow accumulation approaches as in the D8 algorithm. The improvements made result in more realistic contributing areas and solve the problem of parallel flow paths. Wilson and Gallant (1996) noted that the flow dispersion algorithm should best be disabled in high contributing areas, because stream lines are usually quite well defined. According to Mendicino and Sole (1997), the algorithm leads to problems in channel network definition.

MFDmd

Qin et al. (2007) developed a modified MFD as proposed by Quinn et al. (1991). The MFD‐md algorithm is adaptive to local terrain conditions and uses as a linear function of maximum downslope gradient. Thus, the flow accumulation can be modelled more reasonably. In comparison to the original MFD algorithm, MFD‐md achieves fewer errors in the calculated cumulative upslope area α.

The MFD‐md algorithm results in smoother spatial TWI distributions which gives a more realistic state of flow, especially in flat areas. Additionally, TWI values tend to be slightly smaller than the ones calculated with the approach based on Quinn et al. (1991). As other MFD algorithms, the MDF‐md algorithm still can not get rid of the problems occurring in flat areas where soil moisture distribution should be homogeneous. Initial removal of sinks and pits is also necessary.

DEMON

The DEMON algorithm was developed by Costa‐Cabral and Burges (1994) and is based on the planar representation as initially described by Lea (1992). DEMON was created to overcome existing restrictions of the D8 algorithm. Flow is represented in two dimensions and directed by aspect in a flow tube network. Flow is generated at each source pixel and routed down a stream tube until a sink is encountered or the DEM boundaries are reached. The stream tubes are constructed from points of intersections of lines drawn in gradient direction and grid cell edges.

(25)

23

Contributing and dispersal areas can be determined. Main advantages as pointed out by Costa‐

Cabral and Burges (1994) are those of contour based models according to Moore et al. (1988).

Furthermore, DEMON is able to represent varying flow width over non planar topography and can be calculated on rectangular grid DEMs. Additionally, the DEMON algorithm represents the hillslope aspect better than the D8 algorithm. Curved flow paths are also represented in contrast to the D8 algorithm (Costa‐Cabral and Burges, 1994). Another advantage as pointed out by Wilson and Gallant (1996) is that DEMON shows the differences between convergent and divergent areas more accurately than the FD8 and FRho8 algorithms. However, the DEMON algorithm has problems with unwanted indentations in isolines which occur by the approximation of conical surfaces by a mosaic of planes. These could only be resolved by a non‐

planar surface fitting. However, this would yield large computational expenses and a high risk of error introduction and is therefore unsuitable. This limitation could be fixed by the D1 algorithm by Tarboton (1997) who found the DEMON algorithm not working correctly in flat areas and when exceptions in data sets occur.

D∞

The D∞ algorithm was invented by Tarboton (1997). It is based on the concept of MFD algorithms but incorporates the idea of a block‐centred surface fitting scheme as the 2‐D Lea, 2‐

D Jensen and DEMON algorithms. The new approach of flow dispersion is based on the calculation of a single angle between 0 and 2π being the steepest downward slope on the eight triangular facets centered at each pixel. Then, upslope area is determined by dividing the flow between two downslope pixels according to how close this flow direction is to the direct angle to the downslope pixel. Flow in flat areas is treated as proposed by Garbrecht and Martz (1997).

The algorithm is considered advantageous because of its abilities to minimise unrealistic dispersion in the calculation of the cumulative upslope area α and to handle ridges, pits and flat areas in natural catchments.

Tracking Flow Direction (TFD) algorithm

The TFD algorithm was created by Pan et al. (2004) to overcome the limitations of TWI calculation in flat areas. It treats minimum slope values as in Wolock and McGabe (1995) while the flow direction is determined according to the ARC/INFO flow direction module. The TFD algorithm is found to give superior results in flat areas in comparison to the other flow direction algorithms.

Spatially distributed flow apportioning algorithm (SDFAA)

Kim and Lee (2004) propose a new integrated flow determination algorithm. A spatially varying flow apportioning factor is introduced in order to distribute the contributing area from upslope

(26)

24

cells to downslope cells. An iterative process of parameter optimization using a Genetic Algorithm (GA) is implemented.

The SDFAA was found advantageous to existing approaches. A field example showed less overdissipation problems near channel cells, a robust parameter determination and the connectivity feature of river cells. However, and in contrast to other algorithms, the SDFAA requires the channel location as an additional input. Similar to other flow direction models, the removal of depressions and flat areas in elevation data is necessary prior to algorithm application. The SDFAA further resolves valley bottom overdissipation problems of the MFD and provides a more flexible allowance of flow distribution to downslope cells. Furthermore it reduces some problems concerning flow determination in complicated field topography.

However, the application of the spatially distributed flow apportioning algorithm to other watersheds still needs to be investigated (Kim and Lee, 2004).

2.4.4 Algorithm comparison

Wolock and McCabe (1995) found that mean TWI values increase and that smoother TWI patterns across the DEM surface result when using MFD algorithms in comparison to SFD algorithms. Desmet and Govers (1996) examined the influence of different flow routing algorithms on upslope contributing areas which are considerably smoother for multiple flow algorithms and the prediction of ephemeral gullies. The latter could be modelled more accurately with multiple flow models because single flow routing algorithms are highly sensitive to small errors in elevation data. On the whole, Desmet and Govers (1996) found the D8 and 2D‐

Lea algorithms superior for they visually seem to correlate better with the main drainage lines.

Mendicino and Sole (1997) noted that the D8 and Rho8 procedures, although conceptually different, produce similar distribution functions of the TI. The FRho8 and DEMON procedures were additionally identified to behave similarly as the Rho8 algorithm. For this reason it is suggested to calculate flow directions with a mixed procedure, based on the method of Quinn et al. (1995).

Gallant and Wilson (1996) stated that the FD8 and FRho8 algorithms are superior to the D8 algorithm because they result in more realistic distributions of contributing areas in upslope areas and their ability to resolve the problem of parallel unrealistic flow paths.

The DEMON and 2D‐Lea algorithms have the advantage that they are able to determine flow directions precisely and that they are deterministic. Errors may however be introduced in creating the plane through grid cells (Tarboton, 1997). Flow according to MFD as proposed by

(27)

25

Quinn et al. (1991) follows the topographic slope but introduces significant dispersion. Tests on the basis of the evaluation of test statistics and examination of influence and dependence maps on DEMs with different resolutions show that the D∞ algorithm performs better than D8, 2D‐

Lea and MFD as proposed by Quinn et al. (1991) in terms of upslope area calculation. Tarboton (1997) found the D∞ and DEMON algorithms to generally work best. They both produce smallest errors in the comparison of calculated and theoretical upslope. The D∞ algorithm additionally overcomes the problems of loops and inconsistencies of the DEMON and other plane‐fitting models.

Wilson et al. (2000) examined the effects of DEM source, grid resolution and choice of flow direction routing algorithms on five topographic attributes. SFD algorithms result in more cells with rather small upslope contributing areas. FRho8 and DEMON agreed best with each other in surface comparison.

According to Zhou and Liu (2002) the infinite flow direction algorithm proposed by Tarboton (1997) achieves better drainage configurations than other algorithms when using rasterized topographic data to calculate catchment areas.

Endreny and Wood (2003) explored the congruence which is explained as the overlap between field‐sketched flow path networks observed during a storm and the calculated moisture indices of observed and DEM‐delineated overland flow networks on a 30 m horizontal resolution DEM.

Five common flow direction algorithms were compared in terms of their spatial congruence. The D8 and MFD algorithms were assigned the lowest congruence ratings whereas the 2D‐Lea and D∞ algorithms, which direct flow to a maximum of two neighbours, produced the highest congruence ratings. The modifications to the D8 and MFD methods allowing limited dispersive flow could increase congruence. Dispersive flow to two or three surrounding cells is suggested as an optimum.

Pan et al. (2004) evaluated the performance of six different floating algorithms (combinations of flow direction and slope algorithms) by Root‐Mean‐Square Error (RSME) comparison on three different idealized hillslopes (planar, convergent and divergent (concave and convex)) and two combinations of vertical and horizontal terrain data resolutions, respectively. Among the tested algorithms, the MFD algorithm was found to be most accurate in all examined cases, followed by the BFD and SFD algorithms which especially produced significant errors when contour lines were not parallel to one of the eight cardinal flow directions. In all cases it could be observed that TWI errors decreased with an increasing vertical resolution independent on the underlying flow direction algorithm. The differences in algorithm‐based determination of contour lengths did not greatly affect the results. Additionally it could be found that the RMSE of the calculated

(28)

26

topographic indices compared to initial slopes and contour lengths increases with increasing slope values. The TFD worked best in flat areas where the SFD and BFD have problems resolving the right flow direction. Since the MFD algorithm calculates dispersive flow it propagates to a broader region than the SFD and BFD algorithms.

Güntner et al. (2004) recommended an intermediate approach between SFDs and MFDs.

Results from different studies (Wolock and McCabe, 1995; Desmet and Govers, 1996; Tarboton, 1997; Wilson et al., 2000 and Endreny and Wood, 2003) show that existing algorithms behave differently according to the underlying landscape (e. g. planar, conical, hilly and steep or flat).

The current state of research is not yet able to recommend a specific flow routing model suitable for all of the varying landforms. Wilson et al. (2007) compared the performance of several flow direction algorithms on a 10 m resolution DEM across six different fuzzy‐defined landscape classes and tested the hypotheses that the performance of five flow routing algorithms does not change with the flow descending from steeper to flatter terrain (1) and that the performance of those algorithms does not vary across different fuzzy‐classified landscapes (2). The lowest mean SCA occurred on hilltops and ridgelines becoming higher as the terrain became flatter. Overall, landscape can be divided into two groups (high elevation and low elevation) based on the SCA calculation results. All in all, the results of the comparison show which flow routing algorithms give similar or different results and where differences lie, but they can not indicate which algorithm works best for what landscape type in regard of measured soil moisture or hydrologic theory. All in all, the results suggest that the D∞, MFD, and DEMON algorithm performed better than the D8 and Rho8 algorithms.

Sørensen et al. (2006) compared TWIs calculated with different modifications and evaluated them in terms of their correlations with five soil moisture related parameters. Upslope area, creek cell representation and slope were determined according to different theories. Although correlations could be observed, no TWI combination was found superior for all measured variables. It is therefore assumed that the best method varies and is site specific.

2.5 Soil moisture retrieval in remote sensing

With readily available high resolution and consistent, long‐term remote sensing records and progress in digital image processing techniques as well as the tendency towards macro scale modelling, the combination of remote sensing and hydrological modelling becomes more and more promising. Remote sensing can provide large scale distributed data sets where ground based measurements are unavailable. However, the determination of small scale hydrologic

(29)

27

features is not sufficient. Determination of fluxes and parameterization of hydrologic processes, e. g. evaporation, infiltration, runoff or drainage still require detailed hydrological models (Dubayah et al., 1995). There have already been some attempts to combine remote sensing and hydrologic modelling (Milly and Kabala, 1985; Houser et al., 1998; Bauer et al. 2006). Lin et al.

(1994) compared remotely sensed and simulated moil moisture over a heterogeneous watershed and found both techniques to be consistent with ground measurements. Kite and Pietroniro (1996) presented a comprehensive overview of remote sensing applications in hydrologic modelling.

2.5.1 Soil moisture retrieval techniques

Soil moisture retrieval with remote sensing techniques can be achieved in all regions of the electromagnetic spectrum. Comprehensive comparisons between different retrieval techniques can be found in Bryant et al. (2003) and Moran et al. (2004). Engman and Gurney (1991) present an overview of soil moisture retrieval techniques and analyse their capabilities, advantages and disadvantages. The following methods mentioned have been successfully used to model soil moisture and are briefly described:

Gamma radiation techniques

The detection of soil moisture is based on the different natural gamma radiation flux emitted from dry and wet soils and captured by airborne sensors at a height of 100 m to 200 m, because the atmosphere attenuates the radiation flux. Background soil moisture and gamma count rate need to be obtained by antecedent calibration flights as additional parameters.

Hyperspectal techniques

Reflected solar radiation is measured. Soil moisture can be estimated by unequal albedos (wet soils have lower albedos than dry soils). Many confounding factors such as surface roughness, angle of incidence, vegetation cover and soil colour and texture complicate moisture determination.

Thermal techniques

Thermal techniques make use of the relation of diurnal surface temperature change and soil moisture. Surface temperature is primarily dependent on thermal inertia of soil, which in turn is related to thermal conductivity and heat capacity. With increasing soil moisture, both the thermal conductivity and heat capacity rise. Using thermal techniques requires anterior knowledge of the soil type, is limited to bare, unvegetated soils and retrieves soil moisture only

Soil moisture modelling using TWI and satellite imagery in the Stockholm region