Mapping the karstification potential of central Cebu, Philippines using GIS

Karstification is influenced by numerous factors, such as climate, geomorphology, and geology. Karstified areas are faced with numerous problems such as ground subsidence and groundwater vulnerability. In rapidly urbanizing areas like Cebu in the Philippines, it is imperative to know the karstification potential, especially if karst landforms are already present. We evaluated the karstification potential of central Cebu using Raster Overlay Analysis (ROA) in a geographic information system (GIS) platform. ROA uses multiple maps which represent different factors affecting karstification, such as precipitation, temperature, vegetation, elevation, slope, drainage, lineaments, and geology. These factors were then correlated with the karstification potential map using the Cell Statistics tool. ROA reveals that 17% of central Cebu has very low karstification potential, while < 65% has low to moderate karstification potential and 18% has high to very high karstification potential. High karstification potential is generally associated with high precipitation (> 2100 mm/yr), low temperature (< 24.5 °C), very dense vegetation, high elevation (> 787 masl), gentle slope (< 10°), very high lineament density (> 0.45), very low drainage density (0), and older limestone bedrock. Areas with moderate to very high karstification potential coincide with cave locations. Cell Statistics (Range) reveals that the karstification potential map is strongly correlated with precipitation and temperature, and weakly correlated with drainage density and geology. While assumed to greatly affect karstification, the weak correlation of local geology with the karstification potential is attributed to lithological heterogeneities in central Cebu. Information on an area’s karstification is beneficial to effective land use planning and reduction of karst-related hazards.


Introduction
Karstification is the process by which highly soluble bedrock such as carbonates and evaporites are dissolved through the chemical and mechanical action of water (Monroe 1970). This dissolution process leads to the formation of a distinct topography known as karst terrain, which is characterized by geomorphological landforms such as caves, subterranean rivers, pinnacle karst, and sinkholes. Karstification can occur at any latitude and elevation, allowing karst topography to be found anywhere in the world (Cahalan and Milewski 2018). Around 20% of the Earth's landscapes are karstic (Ford and Williams 1989). In the Philippines, around 10% of the land area is karst (Wagner 2013). This is manifested through a multitude of karst features that can be found in the archipelago, some of which have been developed into tourist destinations and protected areas (e.g., Puerto Princesa Subterranean River National Park in Palawan, Chocolate Hills in Bohol, and the karst towers in El Nido, Palawan). Hazards associated with karst features such as sinkhole collapse and land subsidence have also occurred in the country. A notable instance is the appearance of close to 100 sinkholes in Bohol following the 2013 M w 7.2 Bohol earthquake, some of which have engulfed whole houses (Matus 2013). Aside from these, groundwater in karst areas is also prone to pollution and contamination (Brinkmann and Parise 2012).
Several factors influence the degree and potential of karstification in an area. Ford and Williams (1989) used the term "karst denudation" to pertain to the dissolution of carbonates or evaporites through chemical and mechanical means resulting in karst topography. Karst denudation rates are generally affected by climate, geology, and geomorphology (Ford and Williams 1989). Differences in the relative influence of each factor can also result in variations in the denudation rate. The relative influence of each factor differs spatially, making karstification studies sitespecific in general (Wilson and Beck 1992). To evaluate the relative degree of influence of each factor, geospatial methods are employed through geographic information systems (GIS) (Hubbard 2001;Shofner et al. 2001;Farrant and Cooper 2008;Doctor and Doctor 2012). Multi-criteria analysis of geological and geomorphological parameters using GIS has been previously utilized for potential suitability and hazard assessment modeling, which eventually contributed to the selection of areas that are potentially susceptible and not susceptible to natural hazards (Bathrellos et al. 2012;Skilodimou et al. 2019).
Numerous karst studies have already been conducted by several authors around the world, most of which are concentrated on cave networks, sinkhole occurrence, and groundwater systems (Whitman et al. 1999;Hung et al. 2002;Bakalowicz 2005). With the advent of modern methods in geomorphological studies, karst studies commonly utilize remote sensing and GIS tools that are sometimes accompanied by little to no field validation. Papadopoulou-Vrynioti et al. (2013) used a multivariate statistical method and GIS to map the karst collapse susceptibility of the northern suburbs of Athens. Seif and Ebrahimi (2014) used GIS modeling and introduced a karstification index to evaluate the degree of karstification of carbonate rocks in Iran; this study concentrated on outcropping carbonate rocks and used higher resolution data (e.g., LiDAR) than previous studies on the same area Moradi et al. (2016) used fuzzy logic and analytical hierarchy process (AHP) models to map the karstification potential of northeastern Khuzestan Province in Iran. All studies are concentrated on a semi-arid to semihumid environment, with landscapes that have barren to little vegetation. The two latter studies also aimed towards understanding karst development as well as looking into the groundwater potential of the area. An easier approach, however, can also be used in dealing with raster data, especially if rankings cannot be given. Simple raster addition or other mathematical operations applied on a local scale have been previously proven as sufficient (Dumperth et al. 2016).
Additionally, several approaches to vulnerability in karst areas have been explored by previous authors. The Epikarst, Protective Cover, Infiltration conditions, and Karst (EPIK) network development is a method that assesses groundwater sensitivity of karst areas through multi-attribute weighting and rating (Doerfliger et al. 1999). Rocks, Epikarst, Karstification, and Soil Cover (REKS) is another method following EPIK that presents an easy alternative for karst groundwater vulnerability estimation (Malík and Svasta 1999). An improved version of SINTACS, a method used to assess intrinsic groundwater vulnerability, has also been used to assess the karstic region of Apulia, southern Italy (Marsico et al. 2004).
In this study, the karstification potential of central Cebu (Fig. 1), a tropical karst landscape, is analyzed through a similar multi-criteria approach using Raster Overlay Analysis (ROA) in GIS. Thematic layers corresponding to various factors affecting karstification are used (i.e., elevation, slope, lineament density, drainage density, precipitation, temperature, geology, and vegetation density). The resulting karstification potential map is one of the first attempts to provide insights as to where most of the karst features can be found or developed in the area. It also shows which factors have a stronger influence on karst denudation. Karstification potential maps can be used as a valuable tool in land use and urban planning to avoid or mitigate the hazards that result from a karstic landscape.

Study area
Cebu is one of the Visayan Group of Islands located in central Philippines (Santos-Yñigo 1951). It is an island that trends northeast-southwest and was regarded as a large anticline whose stratigraphy and structure are representative of the Visayan Basin (Rangin et al. 1989). Consequently, central Cebu is characterized by a mountainous middle region (termed as central highlands) with younger strata found closer and dipping towards the coast (Santos-Yñigo 1951) (Fig. 1). Central Cebu is the most lithologically diverse part of Cebu Island. Basement rocks of Cretaceous age can be found juxtaposed with younger Cenozoic formational units (Santos-Yñigo 1951;Müller et al. 1989). This distinct configuration of formations with varying ages (i.e., very old units with very young units) is the product of the complex tectonic history of the Visayan Basin and the Negros Volcanic Arc (Rangin et al. 1989). Due to uplift caused by anticline formation and subsequent erosion, the older basement units were exposed in central Cebu. Sea-level changes, especially during the Middle Oligocene to Pleistocene, allowed the deposition of the younger sedimentary units ).
More than 70% of Cebu Island is underlain by sedimentary formations, most of which are carbonate-bearing ; only the central part of the island and nearby municipalities are underlain by non-sedimentary rocks such as schists, diorites, and volcanic rocks (Bureau of Mines andGeo-sciences 1983a, b, c, d, 1985) (Fig. 1). Unlike most areas of karstification potential studies, karst landscapes in Cebu are commonly covered with dense vegetation. Of the 13 sedimentary formations in central Cebu, only the Luka and Pandan Formations are not carbonate-bearing (Fig. 1).
The resulting drainage pattern of the area is dendritic to trellis, which is due to the fault and fold structures that are manifestations of the formation of the anticlinorium (Fig. 2). The major fault system in the study area is the Central Cebu Fault System (Philippine Institute of Volcanology and Seismology 2015; Mendoza et al. 2022). Major streams include Abucayan River, Maingit River, Talavera River, and Sapangdaku River, which drain into the Tañon Strait in the west; and Cotcot River, Butuanon River, and Mananga River which drain into the Camotes Sea in the east (Fig. 2). The easternmost areas of central Cebu are vast reclaimed lands which include parts of Cebu City, Mandaue City, and Talisay City.
Due to its vast karst landscapes, karst features such as caves and sinkholes are common on the island. Cebu currently faces karst-related phenomena such as sinkhole collapse and water scarcity (Leyson 2017;Mondragon and Erram 2018). Most of Cebu Island's population, furthermore, is in its central area, since it also includes its capital, Cebu City (Philippine Statistics Authority (2016)).

Analysis of factors promoting karstification
Raster Overlay Analysis involves the manipulation of raster data through different GIS toolsets to extract information on elevation, slope, and lineament density, among others (Environmental Systems Research Institute [ESRI], (2016)). The Raster Calculator, under the Spatial Analyst toolbox in Arc-GIS, executes different mathematical operations (Map Algebra) given multiple raster datasets (Environmental Systems Research Institute 2016. In this study, several factors from available and published sources (Table 1) were reclassified into five categories (very high to very low; Tables 2 and 3) based on their general influence on karstification and added through Raster Calculator (Fig. 3) to generate a karstification potential map of central Cebu. Vector data layers were converted into raster through Rasterization so they can be inputted into the Raster Calculator expression. The different thematic layers corresponding to the various factors are subsets of the following major factors: climate, vegetation, topography, geomorphology, and geology. Each layer is reclassified from 1 to 5, wherein higher values reflect conditions that positively induce karst denudation following the classification of Farrant and Cooper (2008), Green (2015), and Moradi et al. (2016). After reclassification, the layers were added using Raster Calculator (Fig. 3). Higher values reflect areas with higher karstification potential. The karstification potential map was then compared with previous studies and locations of existing karst features to check its validity. Cell Statistics (Range) was then used to compare each factor layer with the output karstification potential values to determine relative correlation. Authority (NAMRIA)). Blue lines represent the major streams: 1-Abucayan River, 2-Maingit River, 3-Talavera River, 4-Sapangdaku River, 5-Cotcot River, 6-Butuanon River, 7-Mananga River (NAMRIA)  BMG (1983aBMG ( , 1983cBMG ( , 1983dBMG ( , 1985PHIVOLCS 2015;Mendoza et al. 2022); previous field surveys Drainage density IFSAR (streams, drainage basins) 5 m a NAMRIA Geology Geology Geologic Map of selected Cebu Quadrangles 5 m a Mines and Geosciences Bureau (BMG 1983b(BMG , 1983a(BMG , 1983c(BMG , 1983d(BMG , 1985; previous field surveys

Climate
Climate data such as precipitation and temperature were used in this study. Higher rainfall, which is related to higher elevations and lower temperatures, has been identified as conducive to karst development (Santo et al. 2007;Green 2015) since water action is needed for dissolution. Cebu Island belongs to the largest climate group in the Philippines, which has a relatively flatter seasonal cycle and an average rainfall peak of 230 mm from July to October (Corporal-Lodangco and Leslie 2017). Precipitation (Fig. 4A) and temperature (Fig. 4B) data were averaged and clipped from the average monthly data from 1970 to 2000 by Fick and Hijmans (2017). We added the 12 average monthly precipitation datasets to obtain a mean annual precipitation layer (in mm/yr). The same procedure was applied to 12 average mean temperature datasets to yield a mean annual temperature layer (in °C). Both layers were then reclassified into five categories (very high, high, moderate, low, and very low) based on manually adjusted equal interval classification of the data (Table 2). Equal interval classification divides the data into equal subranges and was used for data with familiar ranges, such as precipitation, temperature, and elevation. Equal interval is also recommended to emphasize value differences if the values in the data are evenly distributed, such as lineament density (Environmental Systems Research Institute 2016). The precipitation and temperature maps of central Cebu appear very similar in distribution but exhibit an inverse relationship. Precipitation is highest in the central areas and decreases towards the coast. Temperature is lowest in the central areas and increases towards the coast. Precipitation and temperature play opposite roles when it comes to karstification potential. Higher precipitation leads to higher infiltration, yielding higher karstification potential. Heavy rainfall may also lead to sinkhole collapse, further favoring karstification (Green 2015). Lower temperature has been previously linked to higher karstification potential (Moradi et al. 2016).

Vegetation
Vegetation cover (Fig. 4C) is obtained through the reclassification of the 2010 land cover vector file from the National Mapping and Resource Information Authority (NAMRIA) into five categories (very sparse, sparse, moderate, dense, very dense; Table 3). Very dense areas correspond to areas classified as closed forests and mangrove forests. Dense areas correspond to areas classified as open forests. Moderate areas correspond to areas classified as grasslands or shrubs. Sparse areas correspond to areas classified as cultivated. Very sparse areas correspond to barren land, inland water, built-up areas, or fishponds. Vegetation cover has been used by Moradi et al. (2016) as a proxy for infiltration, wherein denser vegetation begets higher infiltration, leading to high karstification potential. In tropical areas, higher vegetation density is also linked with karst terrain development, since denser vegetation is conducive to near-surface percolation, allowing groundwater to be stored in karst bedrock that is nearer to the surface (Chenoweth 2003;Marsico et al. 2004;Moradi et al. 2016).

Topography
Topographic data such as elevation and slope were extracted from a 5-m resolution Interferometric Synthetic Aperture Radar (IFSAR)-derived digital elevation Table 2 Range of values per category generated from raster data using manually adjusted equal interval classification. Slope intervals were modified from the classification of Green (2015) Reclassification value model (DEM) from NAMRIA. Elevation does not directly influence karstification, but has been used as a proxy for rainfall, recharge, and hydraulic gradient; karstification potential is increased with increased elevation in a region since rainfall and groundwater recharge are also increased (Moradi et al. 2016;Berhanu and Hatiye 2020). Flat slopes are also conducive for karst denudation since a lower slope gradient allows concentration, and therefore dissolution, in certain areas (Farrant and Cooper 2008;Green 2015). However, this can also hinder karst erosion which is more prevalent on steep slopes (Santo et al. 2007). The elevation map (Fig. 4D) was reclassified into five categories based on the manually adjusted equal interval classification of values, as there is no basis on how to classify elevation data with respect to elevation, but the elevation values in the area are well distributed ( Table 2). The slope map (Fig. 4E) was also reclassified into five categories (Table 2), where very low/flat slope corresponds to areas with less than 5° slope while very high/steep slope corresponds to areas with greater than 35° slope (modified from Green 2015).

Geomorphology
Geomorphology data include lineament density and drainage density. The presence of lineaments, which contribute to secondary porosity and higher permeability, is also associated with karst formation and development (Brook and Allison 1986;Doctor et al. 2008;Faivre and Pahernik 2008). Lineament density was calculated from the average frequency of lineaments in the area, which was generated using the Kernel Density Tool, under the Spatial Analyst toolset, in ArcGIS 10.3 (Fig. 4F). These lineaments correspond to the different faults from existing maps (Bureau of Mines andGeo-sciences 1983a, b, c, d, 1985; Philippine Institute of Volcanology and Seismology 2015; Mendoza et al. 2022) which were validated during field surveys. Drainage density was obtained by dividing the total length of streams by the Fig. 3 Flowchart of methods used in this study. Raster files are highlighted in orange while vector files are highlighted in green. GISbased processes are italicized, and thematic layers are highlighted in gray. A IFSAR-derived data were subjected to Clipping, Terrain Analysis, Field Calculator, and Kernel Density tools to yield the elevation, slope, drainage density, and lineament density maps, respec-tively. B Both mean annual precipitation and mean annual temperature were calculated from the global model raster file by Fick and Hijmans (2017). C Vector files geology and land cover were rasterized first. D All the thematic maps generated from the previous processes were subjected to Reclassification into five categories (very high to very low) total area of each drainage basin in central Cebu (Fig. 4G). Both streams and drainage basins were extracted and delineated from the IFSAR DEM through automated methods (i.e., Hydrology tools) in ArcGIS 10.3. Both the lineament density and drainage density maps were reclassified into five categories based on equal interval classification (Table 2).

Geology
The main factor leading to karstification is the presence of highly soluble bedrock (i.e., carbonate or evaporite) (Monroe 1970;Farrant and Cooper 2008). Limestone, in particular, dissolves more easily if it occurs on its own ("pure") as compared to if it was mixed with other insoluble materials (Sweeting 1979). This porosity factor, adding to high surficial permeability, is also related to the characteristically poor drainage resulting in low drainage density in karst landscapes (Ford and Williams 2007;Farrant and Cooper 2008). Previous geologic maps of the Mines and Geosciences Bureau (Bureau of Mines andGeo-sciences 1983a, b, c, d, 1985) were adopted in this study. Rock units were reclassified into five categories: non-sedimentary, clastics only, clastics with limestone, young limestone, and old limestone ( Fig. 4G; Table 3). Non-sedimentary areas correspond to those underlain by either igneous or metamorphic units, and these were given a value of one (1). Clastic areas pertain to areas that are underlain by sedimentary formations that do not contain limestone, and these were given a value of two (2). Clastics with limestone pertain to areas with both clastic and limestone members, and these were given a value of (3). Young limestone areas pertain to areas underlain by limestone younger than Middle Miocene, and these were given a value of four (4). Old limestone areas pertain to areas underlain by limestone older than Middle Miocene, and these were given a value of 5. Older limestone areas have more developed karst features such as wider, coalesced sinkholes (Lumongsod et al. 2020).

Karstification potential map of central Cebu
The karstification potential map generated from the raster addition (i.e. overlay) of the eight thematic layers (Fig. 4) yielded values ranging from 13 to 33. The map was then reclassified into five categories (very high to very low) (Fig. 5) using Natural Breaks classification. The resulting spatial resolution of the map is 30 arc-seconds (approximately 1 km), following the lowest resolution datasets. Majority of central Cebu is under low (33.70%) and moderate (32.45%) karstification potential while a very small region falls under very high karstification potential (5.95%). The rest belongs to either very low (16.93%) or high (10.97%) karstification potential. The areas of moderate to very high karstification potential coincide with the presence of known caves, such as Cantipla Cave and Satuhan Cave ( Fig. 6A-C). Ground-validated sinkholes in Cebu City by Manzano et al. (2017) were also plotted against the resulting karstification map but did not show a significant relationship in terms of frequency. However, the karstification potential is very high towards the west of Cebu City and becomes lower towards the coast (Fig. 6A). This is comparable with the results of Lumongsod et al. (2020), wherein more mature sinkholes are observed farther from the coast, as this is where the relatively older karstic formations are. In terms of yielded data values, it is observed that the precipitation and temperature layers share the most similarity (i.e., least deviation of cell values) with the karstification potential map. More than 40% of the precipitation dataset matches the cell values in the output karstification potential (Fig. 7). The layers which weakly correlate with the karstification potential map are the drainage density (20.41%) and geology (20.28%) layers.

Discussion
Karstification maps highlight either karstification potential or the degree of karstification. The resultant map in this study reflects karstification potential since the input values resemble factors promoting dissolution and karst development. The degree of karstification usually refers to the occurrence of persistent karst features such as sinkholes (sinkhole/doline index), caves, or even a network of various features characteristic of a karst landscape. While portions of the study area are underlain by non-sedimentary or non-carbonate rocks, the karstification potential Fig. 7 A Cell Statistics (Range) multiple line graph plots showing the differences between the cells from each input (thematic) layer compared with the cells of the output (karstification potential) layer. The precipitation B and temperature C layers show a good correlation with the karstification potential map, while the drainage density D and geology layers E appear to show the least correlation analysis was still performed across the whole study area to account for the possible presence of buried karst or if carbonate units underlie non-carbonate formations, since older rocks in central Cebu have been reported to be in thrust contact with younger rocks (Santos-Yñigo 1951).
Central Cebu is dominantly underlain by sedimentary formations, most of which are carbonate-bearing ( Fig. 1;  Fig. 4H). Majority of central Cebu receives moderate to low rainfall (Fig. 4A) and is subjected to moderate to high average temperature (Fig. 4B). Only a small portion of the study area is densely covered with vegetation (Fig. 4C). Most of central Cebu is low-lying, especially the coastal areas in the west and east, while steep highlands are observed in the central portion ( Fig. 4D-E). Numerous lineaments and streams can be found in the area as well (Fig. 4F-G).
Cell Statistics is a simple tool in ArcGIS that uses the cell values of multiple input rasters to calculate different statistics and produce a corresponding output raster (Environmental Systems Research Institute 2016). In this study, we used the Cell Statistics-Range function to find the differences in the maximum and minimum values of each cell between each factor raster with the generated karstification potential map. The normalized values of the differences, minimum of 0 (no difference) and maximum of 4 (5-1; minimum in one raster and maximum in another), for each input layer were plotted to observe their correlation with the generated karstification potential map (Fig. 7A).
Based on cell value deviations, the karstification potential map strongly correlates with the precipitation (Fig. 7B) and temperature layers (Fig. 7C), and moderately with the elevation layer. Aside from their direct correlation with karstification potential itself, this strong correlation (Fig. 7) can also be attributed to the high similarities between their cell value distributions. Areas with high precipitation coincide with low temperature and high elevation; all these three factors are conducive to karst denudation. Conversely, the generated map least resembles the drainage density (Fig. 7D) and geology (Fig. 7E) layers and reflects little resemblance with the slope layer, contrary to previous studies that highlight geology as the main factor with respect to karstification potential (Farrant and Cooper 2008;Green 2015). Such evaluation may be due to the fact that the employed analysis assigned equal weights to each factor, whereas varied weights are more likely in the natural setting (i.e., karst features cannot be formed in igneous and metamorphic units). The disparity between the resultant map from this study with the results of previous studies can then be resolved if areas underlain by non-sedimentary rocks are removed from the layers altogether, thus requiring pre-raster calculator filtering of each layer. This then removes the possibility for karstification to occur in non-sedimentary areas and makes geology the more dominant factor in promoting karstification.
The karstification potential map in this study suggests that the central Cebu highlands are the most prone to karstification. This implies that there is a higher possibility of karstification-related phenomena such as sinkhole formation, ground subsidence, and groundwater vulnerability occurring in these areas. Usually, typical karst mitigation practices such as grouting would lead to negative environmental impacts such as aquifer pollution which only aggravates the problem (Bonacci et al. 2009).

Conclusions
Karstification potential is the degree of susceptibility of an area to karst denudation and can be affected by multiple factors. In this study, Raster Overlay Analysis in GIS was performed to evaluate the different climatic, topographic, geomorphological, and geologic factors contributing to the karstification potential of central Cebu. The following conditions are conducive to karstification: high mean annual precipitation, low mean annual temperature, very dense vegetation, very high elevation, gentle slope, very high lineament density, very low drainage density, and older limestone bedrock. Most of the central Cebu area exhibits low to moderate karstification potential. Areas with moderate to very high karstification coincide with nine caves of known location. Precipitation, temperature, and elevation show a strong relative correlation with karstification potential; vegetation cover and lineament density show moderate correlation, while drainage density, geology, and slope show weak correlation. Using raster analysis with equal-weighted overlay was used in this study as the first step in analyzing the karstification potential of a karstic area such as Central Cebu. However, karstification potential maps such as these are static and relatively simplified, since digital elevation models (DEMs) only reflect a single or statistically summarized time frame. DEMs or satellite images of varied temporal scales need to be analyzed to reflect any seasonal changes in precipitation and temperature. Much of karstification also occurs below the surface so the karstification potential can be better understood if subsurface data can be added to the analysis.
It is further recommended to explore more statistical or geospatial analysis to further quantify the influence of each of the driving factors. Other predictive tools such as Analytical Hierarchy Process and Geographically Weighted Regression can be used and compared with the equal-weighted overlay approach.
This study utilizes a straightforward method in assessing karstification potential as well as the relative correlation of the karstification potential map with different climate, geomorphological, and geological factors in a manner that is easy to understand and communicate with lay audiences. It is also one of the few studies that have looked into karstification processes and evaluated the karstification potential of a rapidly urbanizing karst terrain in the Philippines. Aside from paving the way for more advanced karstification studies in the Philippines in the future, this study also encourages the participation of the community in understanding and relaying information regarding karst and karstification processes.
Acknowledgements This study was supported by the research project entitled Coastal Tectonics of Central Philippine Islands and Its Implications for Seismic Hazards (under the program "Integrative studies on the Geodynamics of the Philippines and Taiwan  Funding The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: