A high-resolution global dataset of topographic index values for use in large-scale hydrological modelling

Introduction Conclusions References


Introduction
Land Surface Models (LSMs) are widely used for predicting the effects of global climate change on vegetation and land surface temperature (Prentice et al., 2007;IPCC, 2013).However, the simulation of hydrological dynamics within LSMs remains relatively simplified because these models are run at grid-box scales (∼ 60-300 km resolution) and the physics they follow is based predominantly on approximations of processes that occur at much finer spatial scales (Ducharne, 2009;Wainwright and Mulligan, 2013).Correctly characterising hydrology is very important because meso-scale/landscapescale water movements and changes in the water cycle control many effects ranging Figures

Back Close
Full from local energy and carbon fluxes to land-atmosphere feedbacks to the climate system to potentially-catastrophic changes in vegetation distributions.When coupled to atmospheric models, most Land Surface Models (LSMs) can simulate a wide variety of natural and human-modified processes from soil moisture feedbacks on precipitation (Seneviratne et al., 2006(Seneviratne et al., , 2010;;Coe et al., 2009) and river flow (Gedney et al., 2004(Gedney et al., , 2006;;Clark and Gedney, 2008;Milly et al., 2008;Falloon and Betts, 2010;Sanderson et al., 2012) through to vegetation development and carbon productivity (Prentice et al., 2007;Marthews et al., 2012;IPCC, 2013).Although usually applied at mesoscale resolutions (gridcells of scales 10-100 km, e.g.Harding and Warnaars, 2011), LSMs are increasingly finding applicability at finer resolutions approaching 1-10 km, at which the physics they encapsulate begins to approach the more detailed scales (100-1000 m) typically required in process-based hydrological models or used in catchment-based water resources assessments (cf. Wood et al., 2011(cf. Wood et al., , 2012;;Beven and Cloke, 2012).A growing body of work has lately emerged using LSMs to produce large-area projections of current and future water resources for use in applications related to climate change impacts assessment (Gedney and Cox, 2003;Gerten et al., 2004;Falloon and Betts, 2010;Wood et al., 2012;Zulkafli et al., 2013;Harding et al., 2013).
Land Surface Models require a representation of surface and subsurface runoff.Models of runoff production used in regional and continental applications typically contain parameterised physics that is based on statistical representations of processes known to operate at finer scales (Ward and Robinson, 2000;Clark and Gedney, 2008).This can lead to inaccurate predictions in data-sparse regions and generally high uncertainty.The large quantity of detailed topographic information now widely available at sub-mesoscale resolutions offers an opportunity to improve the fidelity of large-area simulations of the hydrological cycle, for the benefit of both climate and hydrological models (Dharssi et al., 2009;Wainwright and Mulligan, 2013).
Currently, the most common approach is to use a statistically-generalised runoff production scheme such as TOPMODEL, which partitions runoff from the soil column Introduction

Conclusions References
Tables Figures

Back Close
Full into surface and subsurface components (Beven and Kirkby, 1979;Quinn et al., 1991Quinn et al., , 1995;;Beven, 1997Beven, , 2012)).One of the most important configurational parameters for TOPMODEL is the well-known topographic index (defined in Appendix A), which is widely used in hydrology and terrain-related applications (Ward and Robinson, 2000;Wilson and Gallant, 2000).
The HYDRO1k global values for Compound Topographic Index (CTI) were released by USGS in 2000(USGS, 2000) and they have since become the most commonly used global ancillary files for topographic index values.HYDRO1k was a great step forward in the development of global hydrological modelling applications: it allowed spatially-explicit hydrological routines to be incoporated in LSMs for the first time and large-scale applications of the TOPMODEL hydrological model to become standard (Beven, 2012).However, because of its relatively coarse resolution (30 arcsec, approximately 1 km at the equator) which limits precise slope calculations, and because it was based on mosaicked elevation data of differing quality over different geographical areas (USGS, 2000), CTI ancillary files are no longer considered ideal and there is a need for improvement.
The limitations of HYDRO1k CTI values become most apparent when considering wetland ares.Wetlands are critical nodes in the Earth System where land-atmosphere fluxes are strongly dependent on seasonal and inter-annual hydrological variability (Coe, 1998;Baker et al., 2009;O'Connor et al., 2010;Dadson et al., 2010).In wetlands, the availability of water introduces important feedbacks on climate via surface fluxes of energy and water and these areas form a key link between the hydrological and carbon cycles (Ward and Robinson, 2000;Gedney et al., 2004;Seneviratne et al., 2006Seneviratne et al., , 2010;;Coe et al., 2009;Dadson et al., 2010).Some analyses based on CTI values have persistently overestimated the extent and persistence of tropical wetlands of various types.Notably, simulations using the Earth System Model HadGEM2, which is parameterised using CTI (Collins et al., 2011), predict much larger and more persistent Amazonian wetlands than actually exist according to current surveys (e.g.Lehner

Conclusions References
Tables Figures

Back Close
Full  , 2004;Prigent et al., 2007;Junk et al., 2011), which may at least partly be caused by the quality of the HYDRO1k CTI.

Döll
In the context of LSMs, the need for high-resolution topographical data across wide spatial domains has recently been highlighted (Lehner et al., 2008;Wood et al., 2011;Lehner and Grill, 2013).With the advent of satellite-based global mapping, notably the Shuttle Radar Topography Mission (SRTM), there has been a significant improvement in the availability of high-resolution datasets with continental coverage, such as in the high-resolution global HydroSHEDS database (Lehner et al., 2008), but unfortunately such datasets have generally not yet been utilised to support large-scale hydrological modelling studies (Wood et al., 2011(Wood et al., , 2012)).
Responding to the needs for higher-resolution data for use in LSMs, in this study we have three main aims: (1) to calculate the topographic index using the GA2 algorithm based on high-resolution global HydroSHEDS data; (2) to compare our values to the current standard for values of this index (the CTI of HYDRO1k) and (3) to discuss current developments in large-scale hydrological modelling and how models can benefit from higher-resolution parameter maps such as these.
The algorithm required for calculating this index is relatively simple (Appendix A), but it has not previously been applied to generate a global dataset at very high resolution because (1) the index must be calculated from harmonised topographic information, which only became available in the 2000s and (2) Land Surface Models (LSMs) have only recently become sophisticated enough to make use of such a high-quality layer (Prentice et al., 2007).Introduction

Conclusions References
Tables Figures

Back Close
Full The HydroSHEDS "hydrologically-conditioned" layers Grid-based topographic index calculations require a Digital Elevation Model (DEM) and we have used the Hydrological data and maps based on SHuttle Elevation Derivatives at multiple Scales (HydroSHEDS) DEM (Lehner et al., 2008; http://www.hydrosheds.org/).The HydroSHEDS database was derived from raw SRTM data at 3 arcsec pixel resolution (approximately 90 m at the equator) through the application of hydrological conditioning in a sequence of correction steps (Lehner, 2013) which resulted in a globally consistent suite of grid layers which were subsequently upscaled to a resolution of 15 arcsec (approx.450 m at the equator).We acquired the HydroSHEDS DEM and also a layer of pre-calculated contributing upstream catchment areas for each 15 arcsec pixel (UPLAND in m 2 , B. Lehner unpubl.data 2013).As of April 2014, HydroSHEDS data has only been produced at its highest quality for all land areas south of 60 • N, so for areas at higher latitude we substituted the HYDRO1k DEM (disaggregated to 15 arcsec resolution) to provide seamless global grids.

Generating the ancillary files
Our calculations had to be carried out over domains composed of complete watersheds, so we mosaicked both the DEM and UPLAND tiles into a global data layer using ArcGIS 10.1 (Esri Inc., Redlands, California).These two input layers were then converted to NetCDF format using gdal (OSGF 2011).Topographic Index values were calculated using the GA2 algorithm, which is the widely-used GRIDATB algorithm with some modifications for use with HydroSHEDS data (see Appendix A for details).Resulting index values for the global land surface were then filtered to remove areas for which topographic index values are invalid or meaningless, including lakes and reservoirs (masked out using the Global Lakes and Wetlands Database, Lehner and Döll, 2004), mountain glaciers and ice caps (using the Randolph Glacier Inventory, Pfeffer et al., 2014) and the Greenland ice sheet (using Lewis, 2009).Introduction

Conclusions References
Tables Figures

Back Close
Full GA2 was run on the ARCUS server for all continental-scale calculations, a 1344core computer cluster at the Oxford e-Research Centre (OeRC).Zonal histograms were plotted using ArcGIS 10.1 and subsequent statistics calculated using R (R Development Core Team, 2013).

Results
We produced a layer of topographic index (TI) values following the GA2 algorithm for all ice-free land pixels worldwide.TI values calculated this way are not just relative measures but consistent and comparable between catchments (Appendix A), so we may compare global values: As expected, TI values are low at ridge-tops (minimal catchment area) and high in valleys (along drainage paths and in zones of water concentration in the landscape, Wilson and Gallant, 2000), yielding a global range of 0.00-25.00and average of 5.99 (Fig. 2).
Wetter areas of the globe generated generally higher TI values (Fig. 1), although there are many exceptions to this (e.g. in desert areas where high TI values do not correlate with high flow accumulation).Zonal statistics calculated for the various lake and wetland types of the world (as defined by the Global Lakes and Wetlands Database, see Table 2) show that pixels representing rivers had the highest TI values (global mean 8.81 over 0.42 × 10 6 km 2 ), but also the highest variance with some river pixels scoring below the global mean for ice-free land outside lakes, reservoirs, wetlands and wetland complexes (global mean 5.88 over 128.99 × 10 6 km 2 ).In terms of TI, wetland complexes in Asia (mostly occurring in India and Tibetan China, Table 2) and mires (mostly occurring in boreal Canada and the Russian Federation) were indistinguishable from dry land (Fig. 3), indicating that wetlands in these areas are maintained by factors other than topography.
In comparison to HYDRO1k, the new TI values from GA2 based on HydroSHEDS were higher for river pixels and slightly higher for intermittent wetlands and lakes.TI Introduction

Conclusions References
Tables Figures

Back Close
Full values were lower at pixels in tropical swamp forests and inundated forests and also slightly lower in coastal wetlands.
The new TI values from GA2/HydroSHEDS were in line with HYDRO1k values for Compound Topographic Index (CTI, USGS, 2000) at most global pixels, but in certain areas there were significant divergences.Considering all river catchments larger than 10 6 km 2 in particular (Table 1), CTI values were higher for many basins, most notably the Amazon, Congo, Paraná, Niger and St. Lawrence rivers, in the case of the Amazon as much as 20 % higher than the values from GA2 (Fig. 4).According to our calculations, the catchments with the highest spatially-averaged TI values were the Murray-Darling, Nelson-Saskatchewan, Nile and Niger (compared to the order Amazon, St. Lawrence, Niger and Nelson under the CTI calculations, although n.b.HYDRO1k's CTI included no estimates for the Murray-Darling, Table 1).Although it might be expected that the size of the Amazon floodplain would be enough to ensure it scored the highest TI, please note that (i) there is no globally consistent correlation between wetland area and TI (Fig. 3) and (ii) because these are spatial averages, the density of wetland within each catchment is more important than the absolute wetland size (and the Nelson-Saskatchewan, for example, is known for a high density of wetland terrain).

Discussion
Modelling soil water flow and runoff generation is of critical importance for simulating land-surface fluxes, predicting water table dynamics, wetland inundation and river routing and, at a regional scale, quantifying surface evaporation rates and the growth, transpiration and seasonality of vegetation (Ward and Robinson, 2000;Baker et al., 2009;Dadson et al., 2010;Marthews et al., 2014).Meso-scale or landscape-scale hydrological processes are therefore key elements in modelling land surface-atmosphere Introduction

Conclusions References
Tables Figures

Back Close
Full exchange processes and critical to the successful use of coupled LSMs to predict the effects of climate change at larger scales.The hydrological routines of LSMs have received much attention and undergone steady improvement in recent years (e.g.Gedney and Cox, 2003;Gerten et al., 2004;Dadson and Bell, 2010;Dadson et al., 2010Dadson et al., , 2011;;Zulkafli et al., 2013;Wood et al., 2011;Wainwright and Mulligan, 2013;MacKellar et al., 2013).However, these mesoscale processes remain generally less well-modelled than processes operating at the finer local-scale (e.g.photosynthesis models) or larger continental-scale (e.g.climate circulation models).Arguably, the development of meso-scale processes has been relatively slow not just because of a lack of complete understanding of the processes involved, but also, more simply, by the limited availability of high-resolution parameter maps for the models concerned (Wood et al., 2011;Wainwright and Mulligan, 2013;Marthews et al., 2014).Because LSMs are now being applied at increasingly high spatial resolution in order to analyse the distribution and movement of water resources, model development is gaining momentum.Large-scale gridded simulations based on high-resolution drivers are now becoming routine, and this has led to an increasingly recognised need for the high-resolution datasets required to drive those simulations (e.g.Wood et al., 2011Wood et al., , 2012;;Beven and Cloke, 2012;Castanho et al., 2013).
Since, 2000, there have been many advances in the base data from which topographic index values may be calculated (i.e.improved DEMs such as SRTM, and improved hydrological derivatives such as HydroSHEDS).A strong demand has arisen in the land surface modelling community for higher-resolution simulations and there are now even calls for "hyperresolution" global simulations at resolutions down to 100 m (Wood et al., 2011(Wood et al., , 2012;;Beven and Cloke, 2012).Now is an appropriate time to revisit and update the availability of global hydrological datasets and in this study we have revised and recalculated the widely-used CTI layer of the HYDRO1k database.We have presented new, high-resolution topographic index maps in formats appropriate for use with the latest LSMs.Introduction

Conclusions References
Tables Figures

Back Close
Full

High-resolution hydrological modelling
TOPMODEL was originally applied at the scale of small catchments, using pixels less than 50 m×50 m in extent (Quinn et al., 1991(Quinn et al., , 1995;;Ward and Robinson, 2000;Beven, 2012), with the index values understood to have relative significance only (i.e.similar values calculated in different catchments do not necessarily imply hydrological similarity, see Chappell et al., 2006).There have been many developments from this basic framework over the years (e.g.see Wilson and Gallant, 2000;Hjerdt et al., 2004;Beven, 2012) and this study has likewise taken a novel approach.Notably, we have applied our calculations at continental scales with larger pixels (approximately 450 m×450 m at the equator), using the resolution correction of Ducharne (2009; also see Moore et al., 1993;Wolock and McCabe, 1995;Clark and Gedney, 2008).Additionally, because our calculations are carried out over complete continental land masses, the index values derived may be considered to be consistent and comparable between catchments.
Although we accept the arguments of Beven and Cloke (2012) that moving to higherresolution data sets is not the only line of development that should be followed, ultimately we support the ideas of Wood et al. (2011Wood et al. ( , 2012) that increasing the resolution at which global hydrological simulations are carried out will have many beneficial knock-on effects.Methane production in wetlands, for example, is critically dependent on the level of the water table (Gedney et al., 2004;O'Connor et al., 2010;Pangala et al., 2013), models of which are in turn dependent on accurate representation of the topography, therefore higher resolution simulations involving improved topographic index values should of necessity improve the representation of wetland fluxes of heat, water and trace gases to the atmosphere (Gedney et al., 2004) and overall estimates of methane release.
In this study we have refined the standard topographic index calculations and greatly improved their spatial resolution.We have presented our new maps of topographic index values both by wetland type (using the Global Lakes and Wetlands Database, Lehner and Döll, 2004) and also in terms of the largest river catchments occurring Introduction

Conclusions References
Tables Figures

Back Close
Full on each continent, finding that in comparison to our revised values, HYDRO1k's CTI topographic index values were significantly higher in some catchments (Table 1).

Limitations of the GA2 algorithm
The topographic index is a measure of the relative propensity for soil to become saturated to the surface as a result of local topography (Beven, 2012).We have calculated it using a robust algorithm (GA2) based on the original implementation of these calculations (GRIDATB, Appendix A).Although topographic index values are comparable between different areas, it is important to remain careful when interpreting their meaning in different regions, such as arid vs. humid, or shallow vs. deep soils (i.e. when factors other than topography influence water accumulation in the landscape).
A well-known limitation of topographic index values is that they are not absolute because the maximum value in any particular catchment is dependent on the catchment's area and slope profile.Therefore, when different calculation methods only result in a change of index distribution shape leaving the minimum the same (at 0), as is the case when comparing our TI values from GA2/HydroSHEDS vs. the CTI from HYDRO1k (Fig. 4), only a partial validation is possible.Therefore, we cannot state conclusively that our revised values are more correct than those of the CTI from HYDRO1k, but the consistency and rigour of the algorithm used and our closeness to the original GRIDATB implementation as well as the improved HydroSHEDS base data used for the calculation lead us to believe that our values are indeed more robust.
A second limitation of our method is that we have used global base elevation data that is not on an equal-area projection.The HydroSHEDS data layers are projected using the World Geodetic System (WGS) 1984, i.e. a grid of unrotated cells that are increasingly stretched in the north-south direction as latitude increases.This implies that slopes will be underestimated in east-west directions at higher latitudes as true pixel distances are getting shorter (Appendix A).There is no appropriate method, however, to avoid uncertainty completely in the slope calculations as the underlying SRTM elevation measurements are already unequally spaced, and as there is no commonly agreed Introduction

Conclusions References
Tables Figures

Back Close
Full upon slope calculation method (e.g., FD8, MD8, D8).We assume that our calculations of steepest gradients with average pixel distances provide a reasonable compromise to approximate the real slope of each pixel (see Appendix A).

Conclusions
In this study we have calculated a new high-resolution, spatially consistent data layer of topographic index values for all ice-free land pixels worldwide based on the hydrologically-conditioned HydroSHEDS database (Lehner et al., 2008).These data layers are at four times the resolution of the HYDRO1k compound topographic index layers (USGS, 2000) -actually within the range of "hyperresolution" data as proposed by Wood et al. ( 2011) -and we believe represent the most robust global-scale calculation of topographic index values that exists to date.
LSMs have now been applied over many years to the problem of explaining and predicting global climate change (Prentice et al., 2007;IPCC, 2013).Recent developments in land-surface modelling and Earth Observation have attempted to incorporate better hydrological understanding into these applications, with a particular focus on a better characterisation of the physical processes that control the water cycle (Coe, 1998;Gedney and Cox, 2003;Coe et al., 2009;Dadson and Bell, 2010;Dadson et al., 2010Dadson et al., , 2011;;Zulkafli et al., 2013).LSMs are continuing to develop rapidly, but in some ways code development has progressed more quickly than the development of the parameterisations on which code simulations are based (Wainwright and Mulligan, 2013).We have made a significant step to redress this by improving the quality of the widely-used topographic index parameter maps.We hope this will lead to more robust simulations of hydrological dynamics and all the processes that depend on landscape-scale water movements and the water cycle in general.Introduction

Conclusions References
Tables Figures
The topographic index is a measure of the relative propensity for the soil at a point to become saturated to the surface, given the area that drains into it A and its local outflow slope β (Beven, 2012; increasing A will tend to increase the accumulation of water, but increasing β will tend to reduce it by increasing gravitational outflow, Quinn et al., 1991).The index is often calculated using an algorithm called GRIDATB, originally written in 1983 by K. Beven of the Hydrology Group, University of Lancaster (revised for distribution 1993-1995 by P. Quinn and J. Freer and described in Quinn et al., 1991Quinn et al., , 1995)).
We calculated topographic index values for each pixel using the GA2 algorithm, which is a slightly modified version of GRIDATB version 95.01 (FORTRAN program gridatb.f ) written specifically for this study based on the basic loop structure implemented in Buytaert (2011) with some modifications to allow for the use of HydroSHEDS data.GA2 calculates the outflow gradient of each pixel (Fig. A1) and uses precalculated UP-LAND values from HydroSHEDS for the catchment area A of each pixel (corrected for latitudinal projection distortions, B. Lehner unpubl.data 2013).
Because of the use of the HydroSHEDS DEM, we made three small modifications in GA2 to the standard GRIDATB calculations: Introduction

Conclusions References
Tables Figures

Back Close
Full -We applied the correction for DEM resolution suggested by Ducharne (2009) to allow calculations to be carried out at continental scales (see below).
-GRIDATB used the multiple flow direction algorithm of Quinn et al. (1991Quinn et al. ( , 1995)), also known as the FD8 or MD8 routing model (Wolock and McCabe, 1995;Zhao et al., 2009;Lang et al., 2013).However, in GA2 we instead used a direction-ofsteepest-descent model: the Deterministic Eight Node (D8) routing model (Moore et al., 1993;Wolock and McCabe, 1995;Wilson and Gallant, 2000;Zhao et al., 2009).This was for consistency with the HydroSHEDS drainage direction approach used to derive UPLAND areas in this study, which were calculated using D8.
-The HydroSHEDS DEM does not have uniformly-sized grid-cells because of its native geographic projection (WGS84) where pixel dimensions vary with latitude (i.e. the real height of a pixel increasingly exceeds its width towards the poles).Because slope directions are restricted to the eight cardinal and diagonal directions, we account for varying pixel dimensions in our slope calculations by taking an average distance between neighboring pixels (rather than direction-dependent): we approximated DX as the square root of the area of each cell (with latitudecorrected pixel areas calculated using the Met.Office Unified Model routine are-alat1.f90written by T. Oki in 1996; Dadson and Bell, 2010).When away from the equator, this implies that slopes will be slightly overestimated in north-south directions and underestimated in east-west directions.
-Finally, because the value of dfltsink is undefined on plains (i.e.areas of no outflow and no inflow, which occur more often when vertical resolution is lower) we followed USGS (2000) and Evans (2003) in applying a minimum of 0.001 to tan(β ).Introduction

Conclusions References
Tables Figures

Back Close
Full  Chappell et al., 2006): how to compare topographic index values calculated from DEMs at different resolutions?Although not part of the original topographic index calculations, it has become accepted that topographic index values as calculated above should be reduced to the "equivalent" value for a 1 m resolution DEM by subtracting ln(DX ) (and restricting the result to be ≥ 0).Although not universally implemented (e.g.neither Evans, 2003nor Buytaert, 2011 applied it), applying this scale-correction has become standard: e.g.see Ducharne (2009; also see Moore et al., 1993;Wolock and McCabe, 1995;Clark and Gedney, 2008).

A2 Slope calculations
There are many different methods for calculating slopes from DEM data and there is not yet a universally-agreed method (Wilson and Gallant, 2000;Zhao et al., 2009), so the use of D8 above is not an unreasonable modification.Additionally, slope values depend on the resolution of the DEM (being by default higher for smaller resolutions), therefore both the use of D8 (above) and the resolution correction (Ducharne, 2009) modify the slope values in our calculations.Introduction

Conclusions References
Tables Figures

Back Close
Full  1. Topographic index values from the GA2 algorithm applied to HydroSHEDS data (Appendix A) compared to CTI values from HYDRO1k for all global river basins larger than 10 6 km 2 .Note that some sources quote much higher index values, but these are often not scale-corrected values and are therefore not directly comparable (e.g.Yang et al., 2007).2000), both of which applied the scale-correction of Ducharne (2009).Boxes show mean ± SD index values for the global distribution of that wetland type.For reference, the mean topographic index value for ice-free land outside wetlands is shown by a broken line on all panels (Table 2).
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | A1 Correcting for DEM resolution A question arises when comparing catchments digitised at different resolutions (e.g.
Discussion Paper | Discussion Paper | Discussion Paper | Table

Figure 2 .Figure 3 :Figure 3 .
Figure 2. Histogram of global topographic index values (vertical line shows global mean of 5.99; global maximum is 25.0044 at a pixel within a river island at the confluence of the Amazon and Xingú rivers in Brazil).