AUTOMATIC DETECTION AND DELINEATION OF KARST TERRAIN DEPRESSIONS AND ITS APPLICATION IN GEOMORPHOLOGICAL MAPPING AND MORPHOMETRIC ANALYSIS AVTOMATSKA ZAZNAVA IN OMEJITEV KRAŠKIH KOTANJ IN NJENA UPORABA V GEOMORFOLOŠKEM KARTIRANJU IN MORFOMERTRIČNI ANALIZI E. PARDO-IGÜZQUIZA1, J. J. DURÄN1 and P. A. DOWD2 Abstract UDC 551.435.8:528 E. Pardo-Iguzquiza, J. J. Duran and P. A. Dowd: Automatic detection and delineation of karst terrain depressions and its application in geomorphological mapping and morphometric analysis Digital elevation models (DEM) are digital representations of topography that are especially suitable for numerical terrain analysis in earth sciences and engineering. One of the main quantitative uses of DEM is the automatic delineation of flow networks and watersheds in hydrology and geomorphology. In these applications (using both low-resolution and precision DEM) depressions hinder the inference of pathways and a lot of work has been done in designing algorithms that remove them so as to generate depression-free digital elevation models with no interruptions to flow. There are, however, geomorphological environments, such as karst terrains, in which depressions are singular elements, on scales ranging from centimetres to kilometres, which are of intrinsic interest. The detection of these depressions is of significant interest in geomorphologic mapping because the development of large depressions is normal in karst terrains: potholes, blind valleys, dolines, uvalas and poljes. The smallest depressions that can be detected depend on the spatial resolution (pixel size) of the DEM. For example, depressions from centimetres to a few metres, such as some types of karren, cannot be detected if the raster digital elevation model has a spatial resolution greater than, say, 5 m (i.e., square 5m pixel). In this work we describe a method for the automatic detection and delineation of terrain depressions. First, we apply a very efficient algorithm to remove pits from the DEM. The terrain depressions are then obtained by subtracting the depression-free DEM from the original DEM. The final product is a digital map of depressions that facilitates the cal- Izvleček UDK 551.435.8:528 E. Pardo-Iguzquiza, J. J. Duran and P. A. Dowd: Avtomatska zaznava in omejitev kraških kotanj in njena uporaba v geomorfološkem kartiranju in morfomertrični analizi Digitalni model reliefa (DMR) je oblika predstavitve topografije, ki je še posebej primerna za numerično analizo površja. DMR je odlična podlaga za analizo rečnih mrež in njihovih zaledij. Pri takih analizah so (zaprte) kotanje problem, ker lahko predstavljajo lažno točko odtoka, zato je eden prvih korakov izdelava DMR brez zaprtih oblik. Na kraških območjih pa so kotanje različnih dimenzij in oblik ena glavnih značilnosti površja in predmet številnih študij. Računalniška zaznava brezen, slepih dolin, vrtač, uval, udornic in polij, je torej lahko velika pomoč pri analizi kraških površij. Velikost, ki jo lahko zaznamo, je navzdol omejena z ločljivostjo DMR; tako na primer kotanje velikosti od nekaj centimetrov do nekaj metrov (npr. različne vrste škrapelj) ne moremo določiti z DMR, ki ima prostorsko ločljivost 5 m. V tem delu predstavimo metodo za avtomatsko zaznavanje in omejevanje kraških kotanj. Najprej predstavimo učinkovit algoritem, ki kotanje odstrani iz DMR. Tak, "brezkotanjski" DMR odštejemo od izvornega in dobimo digitalno karto kotanj, ki je osnova za njihovo morfometrično analizo. Iz nje lahko izračunamo različne parametre, kot so površinska gostota, globine, oblike in te povežemo z drugimi parametri, npr. nadmorsko višino. Opisano metodo smo uporabili na območju Sierra de las Nievas v provinci Malaga, Španija. Gre za karbonatno območje, ki ga drenirajo trije glavni izviri, kjer kraške globeli igrajo pomembno vlogo pri napajanju. Na osnovi karte gostote dolin, ki temelji na 324 zaznanih vrtačah in uvalah, smo določili tri napajalna območja vodonosnika. Predstavimo tudi druge morfometrične rezultate povezane z velikostjo in usmeritvijo 1 Instituto Geológico y Minero de España (IGME). Ríos Rosas, 23, 28003 Madrid (Spain), tel. 00 34 91 349 5914, fax. 00 34 91 349 5742, e.pardo@igme.es, jj.duran@igme.es 2 Faculty of Engineering, Computer and Mathematical Sciences, University of Adelaide, Adelaide, SA 5000, Australia, peter.dowd@adelaide.edu.au Received/Prejeto: 15.10.2012 ACTA CARSOLOGICA 42/1, 17-24, POSTOJNA 2013 culation of morphometric features such as the geometry of the depressions, the mean depth of the depressions, the density of depressions across the study area and the relationship between depressions and other variables such as altitude. The method is illustrated by applying it to data from the Sierra de las Nieves karst massif in the province of Málaga in Southern Spain. This is a carbonate aquifer that is drained by three main springs and in which the depressions play an important role in the recharge of the aquifer. A doline density map, produced from a map of 324 detected dolines/uvalas, identifies three main recharge areas of the three springs. Other morphometric results related to the size and direction of the dolines are also presented. Finally the dolines can be incorporated into a geomorphology map. Key words: DEM, endorheic depressions, poljes, dolines, spurious depressions, cartography. kraških globeli. Metoda omogoča tudi enostavno umestitev globeli v geomorfološke karte. Ključne besede: digitalni model reliefa, zaprti bazeni, polja, vrtače, lažne globeli, kartografija. INTRODUCTION Digital elevation models (DEM) are numerical representations of the topography in the form of the altitude (e.g., metres above sea level) of the terrain on a regular grid of points. This is the raster format of the DEM, which is the most widely used because of its suitability for numerical analysis (Burrough and McDonnell, 1998). It has found a large number of applications in earth sciences and engineering especially in the context of geographical information systems (GIS). For example, DEM are extensively used in quantitative terrain analysis (Wilson and Gallant, 2000; Deng, 2007; Zhou et al., 2008). Apart from the calculation of terrain derivatives (slopes, aspect, curvature, etc) a major application of DEM is in hydrological terrain analysis in which the purpose is the automatic delineation of the drainage network and watersheds together with the calculation of parameters related to accumulated flow in the drainage channels. In these applications, depressions hinder the automatic delineation and a lot of work has been done in designing algorithms for the generation of depression-free DEM (Arnold, 2010), that is, DEM in which the depressions have been removed. However, the automatic detection of depressions and the automatic delineation of their geometry are of significant interest in various geomorphological environments, particularly in karst terrains in which maps of depressions are of interest for the evaluation of geological risk associated with potholes (USGS, 2011), the development of erosion models (Lopez-Vicente et al., 2009) and for geomorphology mapping (Siart et al., 2009). Case studies covering the use of GIS techniques for doline delineation and their use in morphometric studies of karst terrains may be found, among others, in Lyew-Ayee et al. (2007), Telbisz (2009, 2010) and Ravbar and Kovacic (2010). In this context, the automatic detection of depressions offers an advantage over manual delineation in various situations such as large study areas, inaccessible areas or areas covered by vegetation. In this work we propose a methodology for the automatic detection and delineation of terrain depressions. The methodology is illustrated by application to a karst terrain in Southern Spain where the endorheic depressions play an important role in concentrating the recharge of the karst system. We show various morphometric results that can be obtained from the map of depressions together with depression density maps that can be incorporated into a geomorphological map. We conclude with a discussion and conclusions. METHODOLOGY In hydrological terrain analysis, a lot of research has been done on the design of algorithms to generate DEM with the depressions filled (i.e., removed). In this context, the procedure proposed by Jenson and Domingue (1988) is the most widely used implementation in GIS software; e.g., the implementation in IDRISI (Clark Labs, 2001), which was used to generate the results that we present here. More recently, the basic Jenson and Domingue algorithm has been extended, improved and modified (Srivastava, 2000; Planchon and Darboux, 2001; Grimal-di et al., 2007). The basic Jenson and Domingue (1988) algorithm for removing depressions can be summarised 1) Fill single-cell depressions. 2) Compute flow directions. 3) Delineate watersheds of undefined flow and define a drainage points table. 4) Find the threshold value for watersheds with undefined flow and raise all watershed cells that are below the threshold. Definitions of flow direction, undefined flow, drainage point and threshold value of the watershed as well as a detailed description of the algorithm can be found in Jenson and Domingue (1988). The algorithm uses iterative spatial techniques (region-growing procedures) that extend the algorithm beyond the neighbourhood operation. Adjacency is defined as eight-way (d8) connectedness, which implies that two cells are deemed to be adjacent if they share a side or a corner. In brief, the method that we propose for the automatic detection and delineation of depressions is based on a simple idea: a digital map of depressions can be easily obtained by the map algebra operation of subtracting the depression-free DEM from the original DEM. As the method can be implemented using procedures already included in commercial GIS software, we will not ex- pand further here on the technical details. For example, in IDRISI the method requires the following two sequential commands: 1) Pit removal (the input is the original DEM and the output is a depression-free DEM), 2) Subtract the pit-free DEM from the original DEM. As will be shown in the case study section this is a simple, but powerful, idea that provides very good results. As with any automatic procedure, our method may yield false positives and spurious depressions. False positives are true depressions created by man-made structures - such as dams, quarries, opencast mines or water ponds - and are thus not relevant in the work described here. Spurious depressions, on the other hand, are artefacts created by the algorithm in entrenched gullies, in flat areas such as the flood plain of rivers or by the presence of a bridge. These are not real depressions as they have a drainage point but the precision of the DEM is insufficient to detect it; they are typically small and follow streams. Thus the number of spurious depressions decreases as the cell size (resolution) of the DEM decreases because the drainage of a given area is delineated with greater precision. This can be checked by coarsening the 5 m DEM to a 10 m DEM and then to a 20 m DEM. In all aspects the results from a 5 m DEM should be much better than those from the 10 m DEM and the 20 m DEM: the number of detected dolines should increase, the delineation of the rims of the dolines should improve and the number of spurious depressions should decrease. CASE STUDY The study area is the Sierra de las Nieves karst aquifer in the province of Málaga in Southern Spain, a few kilometres South-East of the city of Ronda (Fig. 1). This is a very rugged area with altitudes ranging from 400 m to 1919 m (the Torrecilla peak is the highest altitude in the area). It is a typical Mediterranean mountainous karst with high rainfall (average of 1000 mm annual rainfall). This karst massif is in the western-most part of the Betic Cordillera. References for studies of relevant aspects of the Sierra de las Nieves include Matín-Algarra (1985) for geology, Delannoy and Guendon (1985) for geomorphology and Liñan-Baena (2005) for hydrogeology. The methodology was applied to the highest available resolution DEM (Fig. 2a): a 5 m DEM from the Instituto Geográfico Nacional of Spain, freely available for non-commercial purposes at http://centrodedescargas. cnig.es/CentroDescargas/. The 5 m DEM was derived by automatic stereo image correlation of aerial photographs from the Spanish National Aerial Orthophotography Project (Plan Nacional de Ortofotografía Aérea - PNOA) with a resolution of 25 to 50 cm/pixel. The precision of the image is given as the root mean square error of 2 m with a 90 % confidence interval of ± 4 m. The Instituto Geográfico Nacional of Spain is responsible for maintaining the nation's geodetic system using the highest quality standards. Although the planimetric and altimet-ric precisions meet international cartographic standards, in practice we are dealing with pixels, which represent the mean altitude of an area equal to the pixel size and located at the position of the pixel. Fig. 2b shows the detected and delineated depressions (dolines and uvalas in green) for the area of interest inside the limits of the aquifer (white solid line in Fig. 2b). From inspection it was clear that the spurious depressions related to streams (see below) correspond very well to the smallest detected dolines. By trial and error, a threshold of 250 m2 (ten pixels) was selected as optimal in the sense of eliminating the spurious depressions. Figure 1. Location of the Sierra de las Nieves karst aquifer (red rectangle) in Southern Spain. Fig. 2b shows the results of applying the threshold; the total number of detected dolines is 324. Fig. 3a is a detail of Fig. 2b showing a doline field (with some uvalas) and Fig. 3b is a detail of Fig. 3a showing depths of the dolines measured from the rim of each doline (i.e., not an absolute depth with respect to a datum) as given by the algorithm. The Figures show very clear delineations of the dolines/uvalas that can be used (see below) for calculating the preferential geographic directions of the dolines. B 1 MWI ) J -, 'V- Figure 2. A: Raster DEM with a pixel size of 5 m (square pixel side length). b: Automatic detection of depressions (324 dolines/uvalas in blue). The solid white line is the boundary of the Sierra de las Nieves karst aquifer. A Mtlfrt 500.00 V* ■ . * i .. - . ; * " n V B ES figure 3. field of dolines and uvalas. A: Detail of figure 2B; the solid white line is the border of the Sierra de las Nieves karst aquifer. B: detail of figure 3A showing depths as measured from the rim of each doline. The point of greatest depth (black in the figure) is the location of the ponor. Figure 4. A: Detail of small zone near the centre of Figure 3A; b: topographic map of the same area. A comparison of the two Figures shows that significant dolines (A) have been detected, a captured doline (b) that has not been detected, spurious small depressions aligned along streams (C) and other smaller dolines (D and E) not evident from the topography map. A MAGNA geological map faults N = 375 0 0 figure 5. A: direction rose of 375 faults and fractures from the geological map. Each fault and fracture is weighted by its length. b: direction rose of233 directions of development of dolines and uvalas. Each direction has a weight of 1.0, that is, it is considered to be independent of the size of the doline, although in this work only dolines with an area greater than 250 m2 have been considered. Fig. 4a shows the detail of another area of the doline field in Fig. 3a prior to the application of the threshold so as to show the spurious depressions (C in Fig. 4a). Fig. 4b is a topographic map of the same area as that in Fig. 4a. A comparison of Figs. 4a and 4b highlights the strengths and weaknesses of the automatic algorithm: dolines that are not evident in the topographic map may be detected and delineated but captured dolines may be missed. The elongation of dolines and uvalas coincides with the principal directions of mapped faults and fractures, which reflects the strong structural control in their development. This structural control can be checked by comparing the directions of elongations of the geometries and directions of alignment of the dolines with the direction rose shown in Fig. 5. The histogram of directions in Fig. 5a was calculated from 375 discontinuities (faults and fractures) from the geological map of the study area. These Figure 6. Morphometric results. Distribution of the dolines (as accumulated area) with respect to altitude (metres above sea level). There are two important modes at 1050 m and 1700 m and a smaller mode at 1350 m. of N45°E, N65°E, N95°E and N135°E, which coincide with some of the principal fracture directions of the Betic Cordilleras (Benavente and Sanz de Galdeano, 1985). The directional histogram in Fig. 5b shows the preferential development directions of do-lines and was calculated from 233 different measurements of the dolines in Fig. 2b. The correspondence with the directions in Fig. 5a, especially for directions N45°E and N135°E, is clearly obvious. More morphometric results are shown in Figs 6, 7 and 8. Fig. 6 shows the distribution of dolines (accumu- A IK 10 E * - *v * ï : ' 1 c 3 5 O.l 10 100 1000 1QO0O Dofcw area (rn*> 100000 1000000 B „ 1S maximum depth = 2.5 mean depth 1S R' = 0.9527 * E u - f is y ' C i h . ». • I . 0 1 2 S 4 5 Muhi iliiGna ckipli (m) 6 7 Figure 7. Morphometric results. A: Relationship of doline area and mean doline depth. b: Relationship between mean doline depth and maximum doline depth. The depth of each doline is measured from its rim. discontinuities are of regional dimensions (from 300 m to several kilometres) and show clear modes with azimuths lated area) with respect to altitude (m above sea level). The most significant mode occurs at 1050 m, a secondary mode can be seen at 1700 m and there is a third mode, much smaller that the other two, at 1350 m. These three modes correspond to the three principal dolinized areas as shown below. Fig. 7a shows the relationship between doline area and mean doline depth; the dispersion of the positive figure 8. density of dolines in m2.km-2 on logarithmic scale in order to highlight the range of values. The white dots represent the three main springs of the aquifer and the white arrows represent the connections between the areas with higher density of dolines and the springs. The connections on the east has been proved by tracer tests while the other two (on the centre and at the left) are inferred from the known hydrogeology of the area. Figure 9. Geomorphology map of the Sierra de las Nieves aquifer with dolines, potholes, karst springs, karst scarps, streams and the boundary of the aquifer. This map combines the dolines in figure 2B with geomorphology features. correlation decreases as the size (area) of the dolines increases. Fig. 7b shows the expected strong linear relationship between mean doline depth and maximum doline depth. Figure 8 is a map of doline density (m2.km-2) plotted on logarithmic scale to highlight local variations. There are three main areas (dark green) that drain recharge water towards the three main springs (white dots). The white arrows show the connections between dolinized areas and karst springs. Finally, Fig. 9 is a geomorphology map of the Sierra de las Nieves aquifer that integrates the field of detected and delineated dolines with the principal geomorphology features. In conclusion, 324 dolines/uvalas were detected, more than half of which were not detected during previous field visits or from analysing aerial photographs. DISCUSSION AND CONCLUSIONS We have proposed a simple method for detecting and delineating depressions: dolines, uvalas and any other depression that occurs in the terrain as a false positive (e.g., quarries, ponds). Spurious depressions (depressions along streams, gullies or entrenched rivers) have been su-pressed by using an area threshold for the retained depressions, i.e., only dolines with an area greater than the threshold are retained as genuine karst depressions. A threshold of 250 m2, which corresponds to 10 pixels, enabled 324 dolines and uvalas to be detected and delineated and various morphometric parameters to be calculated. The proposed method for the automatic detection of dolines is an efficient tool for the geomorphological mapping of depressions. We have demonstrated its use for karst areas but it is equally valid for detecting seafloor depressions or craters on the Moon, provided that a corresponding DEM is available. To achieve the best re- sults the proposed method should be used in conjunction with all other relevant forms of information such as topographic maps, aerial photographs, satellite images and fieldwork. The method may detect depressions that otherwise would not be detected because of access problems to the area, large dimensions of the study area, the condition of the geomorphology of the depression (depressions almost filled) or by being hidden by vegetation cover. In addition, even if the doline is obvious in the field or in an aerial photograph, the method delineates the geometry of the depression with great accuracy. This method could be added to the geomorphologist's toolbox and, in addition to geomorphological mapping, it may be used for morphometric analysis (mean altitude of dolines, depression density and other morphometric parameters can be readily calculated). ACKNOWLEDGEMENTS This work was supported by Research Project CGL2010-15498, from the Ministerio de Economía y Competi-tividad of Spain. The authors acknowledge the support of Australian Research Council Discovery Grant DP110104766. We thank the reviewers for their comments and suggestions, which have significantly improved the final version of the paper. REFERENCES Arnold, N., 2010: A new approach for dealing with depressions in digital elevation models when calculating flow accumulation values. Progress in Physical Geography, 34, 6, 781-809. Benavente, J. & C. Sanz de Galdeano, 1985: Relación de las direcciones de karstificación y del termalismo con la fracturación en las cordilleras Béticas, Estudios Geológicos, 41, 177-188. Burrough, P.A. & R. A. McDonnell, 1998: principles of Geographical Information Systems. Oxford University Press, pp. 333, Oxford (UK). Clark Labs, 2001. The Idrisi Project. Worcester, MA, USA; idrisi@clarku.edu Delannoy, J.J. & J.L. Guendon, 1986: La Sierra de las Nieves (Málaga). La Sima G.E.S.M. Etude géomorphologique et spéléologique. Karstologia Mémoires, I, 71-85. Deng, Y. 2007: New trends in digital terrain analysis: landform definition, representation, and classification. Progress in Physical Geography, 31,4, 405419. Grimaldi, S., Nardi, F., Di Benedetto, F., Istanbulluoglu, E. & R.L. Bras, 2007: A physically-based method for removing pits in digital elevation models. Advances in Water Resources, 30, 10, 2151-2158. Jenson, S.K. & J.O. Domingue, 1988: Extracting Topographic Structure from Digital Elevation Data for Geographic Information System Analysis. Photo-grammetric Engineering and Remote Sensing, 54, II, 1593-1600. Liñan-Baena, C. 2005: Hidrogeología de acuíferos carbonatados en la unidad Yunquera-Nieves (Málaga). Publicaciones del Instituto Geológico y Minero de España. Serie: Hidrogeología y Aguas Subterráneas n° 16., pp. 322, Madrid. López-Vicente, M., Navas, A. & J. Machín, 2009: Geo-morphic mapping in endorheic catchments in the Spanish Pyrenees: An integrated GIS analysis of karstic features. Geomorphology, 111, 1-2, 38-47 Lyew-Ayee P., Viles H. A. & G.E. Tucker, 2007: The use of GIS-based digital morphometric techniques in the study of cockpit karst. Earth Surface Processes and Landforms. 32, 2, 165-179. Martín-Algarra, A. 1987: Evolución geológica alpina del contacto entre las zonas Internas y las zonas Externas de la Cordillera Bética. PhD thesis, University of Granada, pp. 1171. Planchon O. & F. Darboux, 2001: A fast, simple and versatile algorithm to fill the depressions of digital elevation models. Catena 46, 59-176 Ravbar N. & G. Kovačič, 2010: Characterisation of karst areas using multiple geo-science techniques, a case study from SW Slovenia. Acta Carsologica, 39, 1, 51-60. Siart, C., Bubenzer, O. & B. Eitel, 2009: Combining digital elevation data (SRTM/ASTER), high resolution satellite imagery (Quickbird) and GIS for geomor-phological mapping: A multi-component case study on Mediterranean karst in Central Crete. Geomor-phology, 112, 1-2, 106-121. Srivastava A., 2000: Comparison of two algorithms for removing depressions and delineating flow networks from grid digital elevation models. Master of Science, Virginia Polytechnic Institute and State University, pp. 119, Blacksburg, VA. Telbisz T., Dragušica, H. & B. Nagy, 2009: Doline Morphometric Analysis and Karst Morphology of Biok-ovo Mt (Croatia) Based on Field Observations and Digital Terrain Analysis. HRVATSKI GEOGRAFSKI GLASNIK, 71, 2, 5-22. Telbisz T., 2010: Morphology and GIS analysis of closed depressions in Sinjajevina mts (Montenegro). Karst development, 1, 1, 41-47. USGS (United States Geological Survey) 2011: Using GIS Techniques to Identify and Delineate Karst Features in Tennessee. Geohazards Impacting Transportation in the Appalachian Region- 11th Annual Technical Forum, pp. 36, Chattanooga, TN. Wilson J.P. & J.C. Gallant (eds.), 2000: Terrain Analysis. principles and Applications. John Wiley & Sons, pp. 485, New York. Zhou, Q., Lees, B. and Tang, G. (eds.), 2008: Advances in Digital terrain Analysis. Springer, pp. 462, Berlin.