Evaluating Sinkhole Hazard Susceptibility Using Logistic Regression Model in Khlong I Pan Sub-Watershed, Surat Thani and Krabi Province, Thailand

INTRODUCTION

A sinkhole is a type of natural disasters that often occurs in limestone or karst topography (Cvijić 1925; Trofimova 2018; De Waele and Gutierrez 2022). Most of these topographical features are composed of carbonate bedrock such as limestone (CaCO₃) and dolomite (CaMg(CO₃)₂) (Zeng and Zhou 2019). When carbonate rocks contact with acidic rainfall and groundwater, the carbonate rocks are dissolved to form tall limestone pinnacles and subsurface caves, called karst topography, with large caverns containing stalactites and stalagmites and groundwater streams (La Rosa et al. 2018; De Castro et al. 2024). Accordingly, sinkholes generally form when the surface layer above the holes or caves is collapsed by the groundwater level dropping, landslides and subsidence of the upper surface, underground excavation, groundwater extraction, or earthquakes. Sinkholes cause substantial damage to life and property.

Researchers in the United States created karst topographic maps by collecting geological data and documenting the occurrence of valleys with historical evidence of sinkholes. However, a lack of data, including geological data, soil series, hydrogeology, and data on urban expansion in each region, has left the situation unclear (Veni 2002). Continuous improvements have been made to try to predict sinkhole formation, and it has been concluded that hydrogeological conditions are the main variable that makes sinkholes more likely to collapse (Nam et al. 2020; Wood et al. 2023). There are many reports of sinkholes around the world, such as in Tangshan, China (Hu et al. 2001), the Ebro Valley, Zaragoza, Spain (Galve et al. 2009a), Sango, Tennessee, USA (Siska et al. 2016) and major urban areas of Brazil (Galvão et al. 2015; de Queiroz Salles et al. 2018). All of these areas have subsidence phenomena due to their unique geological and geomorphological characterization (Stefanov et al. 2023). Geo-information technology is widely used in research studies to find sinkhole susceptibility because it is possible to define algorithms and conditions to evaluate and predict areas at risk of sinkhole formation. The researchers used spatial analysis to determine sinkhole susceptibility by using variables that affect sinkhole formation in multivariate analysis (Wu et al. 2018; Jia et al. 2019). There are many approached models of multivariate analysis, including spatial overlay analysis, probabilistic modeling, conditional probability, analytic hierarchy process (AHP), LR modeling, and machine learning. Research by Zhou et al. (2016) analyzed sinkholes that occur in Jili Village in Guangxi, China, using an LR model. The results show that the highest susceptibility area for sinkholes is in the foothills and Datou Hill. This information was beneficial for developing a mitigation plan for the communities living near the sinkhole area. Kim et al. (2018) applied an LR model to assess sinkhole susceptibility in urban areas caused by underground wastewater drainage. The model addresses key variables such as slope and pipe material that are sensitive to sinkholes in the area. Jia et al. (2019) used machine learning and a cloud model to analyze areas sensitive to sinkhole formation using topography and geology characteristics variables. This technique shows the sensitivity of the geological structure. Hu et al. (2021) used analytical hierarchy process and LR models to test the accuracy of the techniques to assess the susceptibility to sinkholes in the Wuhan city area of China. The results show that LR models have better performance than analytical hierarchy process models. In conditions of climate change, rainfall amounts change, causing rainfall patterns to become more variable. Xu et al. (2023) emphasized the significance of rainfall variables in causing soil erosion on karst terrain in Southwestern China, suggesting future terrain subsidence and erosion trends. Amin et al. (2023) utilized machine learning to analyze subsidence in Central Iran, identifying geological structure and underground pores as crucial factors in predicting sinkhole-prone areas. Maleki et al. (2023) investigated sinkhole susceptibility in Iran’s Bistoon-Parav karst region, considering 10 variables including precipitation, lithology, and vegetation. Their findings highlighted lithology as the most influential factor, accounting for 31.52% of sinkhole occurrences. Ramírez-Serrato et al. (2024) tried to find sinkholes in Mexico City by using linear regression models to import data on 13 variables that affect the occurrence of sinkholes, including population density, WD, distance to faults, fractures, roads, streams, elevation, slope, clay thickness, lithology, subsidence rate, geotechnical zones, and soil texture. The findings demonstrate the efficacy of regression models in predicting susceptibility to sinkhole occurrences. Utilizing advancements in geoinformatics technology enables the comprehensive analysis of diverse databases, facilitating the effective assessment of areas prone to sinkhole formation. Such analyses can generate risk maps, offering guidelines for managing vulnerable areas and implementing preventative measures to mitigate potential loss of life and property resulting from these disasters.

In Thailand, sinkholes are prevalent in karst landscapes, particularly in limestone formations dating back to the Permian period, around 286-245 million years ago. Notable rock groups susceptible to sinkholes include the Rachaburi, Saraburi, and Ngao groups. The Rachaburi group, originating from the Middle to Upper Permian period, spans the lower western and southern regions of Thailand (Sone et al. 2012). The Saraburi group extends across the lower Chao Phraya plains and western edge of the Korat plateau, formed during the Upper Carboniferous to Lower Permian period (Udchachon et al. 2022). Sinkholes are primarily found in the northern region of Thailand, where the Ngao group, formed during the Upper Permian period, is prevalent (Pondthai et al. 2023).

Human activities such as groundwater extraction, saltwater pumping, traffic vibrations, and construction exacerbate sinkhole formation, particularly in areas with limestone, dolomite, and marble bedrock prone to dissolution. The Department of Mineral Resources documented sinkhole occurrences from 1995 to 2005, with 66 areas experiencing large sinkholes, notably 25 following the 2004 earthquakes and tsunamis along the Andaman Sea coast in southern Thailand (Frost-Killian 2008; Szczuciński 2020). Urbanization and increased groundwater usage in southern cities contribute to heightened sinkhole risks due to land surface changes and associated land use alterations.

Today, climate change is a major challenge that greatly affects human life, the environment, and economic development. In particular, changes in seasonal rainfall patterns can make sinkholes more severe. This is particularly true for sinkholes that arise from fluctuations in the average annual rainfall and groundwater storage levels. To achieve preventive measures to reduce the loss of life and property caused by such disasters, and consistent with the principles of the United Nations Sustainable Development Goals (SDGs), Goal 13 addresses taking urgent action to combat climate change and its impacts. Arora and Mishra (2023) suggest that current climate change is the primary cause of increased frequency and severity of natural disasters, resulting in widespread damage to people and the economic system. However, the impact of climate change on sinkhole hazard has been barely explored. This issue is due to the lack of continuous recording of hydrological, hydrogeological, and meteorological data. The processes need to be studied and recorded for future research to more effectively assess sinkhole hazard areas. In Thailand, there is a policy for disaster risk reduction, which is the Disaster Prevention and Mitigation Act 2007 (Fakhruddin and Chivakidakarn 2014). Additionally, the Sendai Framework for Disaster Risk Reduction (2015-2030) also utilized in Thailand for creating the substantial reduction of disaster risk in local areas with national strategies (Kelman 2015).

Sinkholes are increasingly common in Krabi and Surat Thani Provinces, particularly in Plai Phraya and Khao Phanom Districts, as well as Phanom, Phrasaeng, and Chai Buri Districts. The upcoming Land Bridge project, part of the Southern Economic Corridor, poses additional risks due to planned infrastructure development. Despite this, the area has not conducted any sinkhole susceptibility studies. This research aims to assess sinkhole susceptibility in the Krabi and Surat Thani Provinces’ KIPs using GIS-based LR analysis, considering geological, geomorphological, and socio-economic factors. The findings will assist in managing vulnerable areas, offering insights for planning and monitoring potential hazards in regions lacking sinkhole risk maps at the local level.

MATERIALS AND METHODS

Study area

The KIPs lies within the Ta Pi watershed, a significant river basin in southern Thailand, spanning latitudes 8°10’ N to 8°50’ N and longitudes 98°40’ E to 99°20’ E, covering approximately 2,087.039 km². Situated between Krabi and Surat Thani provinces, it includes Plai Phraya, Khao Phanom, and Ao Luek Districts in Krabi, and Khian Sa, Phrasaeng, and Chai Buri Districts in Surat Thani. The area’s topography ranges from 7 to 1,346 meters above mean sea level, with the southern part being mountainous, including Khao Phanom peak within Khao Phanom Bencha National Park. The western side features scattered limestone mountains with karst topography.

The drainage pattern exhibits two key characteristics: the western upstream area displays a dendritic pattern with streams like Klong Ya and Khlong Bang Liao flowing into Khlong I Pan southwest to east. In the south, rivers like Khlong Sai Khao and Khlong Chang Tai flow south to north, forming a dendritic pattern upstream and a parallel pattern in the central basin, merging into Khlong Trom, which joins Khlong I Pan in the north.

Geologically, sedimentary rock predominates, including Quaternary alluvium along riverbanks, colluvium farther away, and old river terraces (Dheeradilok 1995; Leknettip et al. 2023). Tertiary sediments, notably the Krabi group, appear in the northeast, featuring shale, slate, calcareous shale, sandstone, siltstone, and minerals like limestone and lignite (Benammi et al. 2001). Limestone rocks are scattered throughout the upstream area, covering approximately 30-40% of the study area, from the Carboniferous-Permian-Triassic periods.

Data Preparation

This research utilized secondary data from various sources to evaluate sinkhole susceptibility in the KIPs area. The data comprised sinkhole occurrences between 2002-2022 and spatial data analyzing factors influencing their formation. Physical factors such as elevation, slope, soil type, geology, distance to streams, active faults, terrain wetness index, Total Hardness (TH), Total Dissolved Solids (TDS), and Chlorine (Cl) were collected, alongside socio-economic factors including land use, well density, and distances to wells and villages. These spatial datasets are crucial variables for assessing sinkhole risk areas, organized in a raster database format with a grid cell size of 30×30 m.

Method

The research process consists of the following steps: (1) sinkhole area analysis, (2) spatial database analysis of driving factors, and (3) statistical approach. The details of each step are briefly explained below (Fig. 2).

Sinkholes area analysis

To identify sinkhole occurrences in the research area from past to present, secondary data spanning 2002 to 2022 were sourced from the Department of Mineral Resources. Additionally, primary data obtained from interpreting sinkhole areas in 1:50,000 scale topographic maps from the Royal Thai Survey Department were utilized. Upon acquiring data from both sources, a conversion from analog to digital format was performed using ArcMap 10.4 software. Subsequently, the data underwent LR statistical analysis in the next stage of the process.

Spatial database analysis of driving factors

The selection of factors affecting sinkhole formation is important for susceptibility analysis (Wei et al. 2021; Hu et al. 2021). This research has selected the factors that directly affect and are related to the occurrence of sinkholes in this study area. The first important spatial data is elevation (digital elevation model-DEM) obtained from the Royal Thai Survey Department (RTSD) in shapefile format, including elevation point, contour line, and water source and water route information. Such data will be analyzed by spatial analysis using the Topo to Raster technique in ArcGIS 10.4 software. The result is DEM data with a 30×30 m grid cell size (Fig. 3a). This data can be analyzed for other variables such as slope (Fig. 3b) and TWI (Fig. 3g). TWI is an index indicating water accumulation and flow tendency to lower basin areas due to the Earth’s gravity (Chen and Yu 2011). High TWI values indicate prone areas, which may be the swamp areas, low slopes, or basins. This variable will be applied to search for areas where sinkholes occur in the study area. TWI analysis can be calculated from Eq. 1 (Hamid et al. 2020):

(1)

where, A_s indicates the definite catchment area and denotes the slope gradient.

Other physical factors that were applied in this research include soil (Fig. 3c), geology (Fig. 3d), distance to active fault (Fig. 3e), and distance to stream (Fig. 3f). The variables soil and geology are nominal data that are converted to raster data format. Data on variables distance to active fault and distance to stream were analyzed using the Euclidean distance technique in spatial analysis tools.

Furthermore, socio-economic factors are related to sinkhole formation, including LU (Fig. 4a), WD (Fig. 4b), distance to well (Fig. 4c), and distance to village (Fig. 4d). The LU data is nominal data received from the Land Development Department (LDD). The data was converted to raster data format. Data on variables distance to village and distance to well were analyzed using the Euclidean distance technique in spatial analysis tools. The WD variable can be analyzed using a mathematical function in ArcMap 10.2 with the Kernel Density (Eq. 2). The results will indicate the density of the artesian well. This variable is one of the variables that will be used for LR statistical analysis next. The Kernel Density (Nistor and Nicula 2021) calculated from the following equation (Eq. 2):

(2)

where ƒ is density ƒ; K is the kernel — a non-negative function; h > 0 is a smoothing parameter called the bandwidth; x₁, x₂, ..., x_n is a univariate independent and identically distributed sample.

The chemical property variables were Total Hardness (TH) (Fig. 4e), Total Dissolved Solids (TDS) (Fig. 4f), and Chlorine (Cl) (Fig. 4g). The data of all three variables were derived from the Department of Groundwater Resources. In this research, the data of all three variables were subjected to spatial interpolation using the Inverse Distance Weighting (IDW) method in ArcMap 10.2.

Statistical approach

From the sequence of steps to find sinkhole area analysis and spatial database analysis of driving factors, it is necessary to search for variables that affect the formation of sinkholes in order to know the context of the collapse that occurred in the KIPs. There were 14 variables used in the analysis, including elevation, slope, soil, geology, distance to stream, distance to active fault, TWI, Total Hardness (TH), Total Dissolved Solids (TDS), Chlorine (Cl), LU, WD, distance to well, and distance to village. Such variables will be analyzed using LR analysis.

LR is a technique for discovering the empirical relationships between a binary dependent and several independent categorical and continuous variables (Ozdemir 2016; Kim et al. 2020; Cao et al. 2020). LR analysis is calculated using the following Eq. 3:

(3)

where P is the flood-prone area, x_i are independent variables, and β is the coefficient value. The LR method was used to provide the variables that were used to analyze which variables influenced subsidence by considering the initial and dependent variables of every grid cell in the study area. The results from the LR analysis can be used to produce sinkhole susceptibility maps, which are classified into 5 classes: very high, high, moderate, low, and very low.

RESULTS

Sinkholes in KIPs area

The KIPs boast a unique karst geography characterized by sedimentary rocks interspersed with limestone formations, particularly evident in the northern, central, and western regions, where a north-south-oriented limestone mountain range prevails. Over the past two decades (2002-2022), 34 sinkholes collapsed in the area, covering a total of 3.52 km². Most sinkholes are concentrated in the northern and central-southern regions of the KIPs, with a notable cluster in Khao Phanom District, totaling 9 sinkholes. The largest sinkhole, spanning 2.049 km², appeared near Bankhoasamyot School in Khian Sa District, Surat Thani province, characterized by undulating plain topography and rubber plantations. This sinkhole formed on September 26, 2016, during heavy rain brought by a monsoon trough, resulting in a chasm nearly 100 meters deep. Another significant sinkhole, covering 0.47 km², occurred in the east of the Phanom Bencha mountain range in Khlong Noi Sub-district, Chai Buri District, Surat Thani, within a forested area. In Khao Phanom District, Krabi Province, the largest concentration of sinkholes (12) was found, primarily in community and agricultural areas, including rubber and palm oil plantations. Details are outlined in Table 2.

Affecting Variables in sinkhole formation in KIPs

The analysis focused on 14 variables crucial for investigating sinkhole formation in the KIPs, employing the LR statistical process. Table 3 illustrates the variables affecting sinkhole formation, with their impact expressed through β values. A positive β value signifies increased sinkhole susceptibility with higher variable values, while a negative value suggests the opposite. The relative operating characteristic (ROC) demonstrates the regression equation’s ability to predict sinkhole risk areas based on probability. The ROC value obtained for sinkhole area probability was 0.947 (Fig. 6), indicating high effectiveness, as values closer to 1.00 signify comprehensive analysis of sinkhole risk areas using all 14 variables.

All variables were significant at the p < 0.01 entry and p > 0.02 removal levels. ROC relative operating.

The study identified 8 variables with negative β values influencing sinkhole occurrence: geology, distance to well, distance to active fault, soil, terrain wetness index (TWI), distance to village, elevation, and distance to stream. Conversely, 3 variables exhibited positive β values: well density (WD), land use (LU), and slope. The Exp β values indicate the variables’ impact on sinkhole formation, with higher values indicating greater influence. WD emerged as the most influential variable, followed by geology, LU, TH, TDS, slope, Cl, distance to stream, elevation, TWI, distance to village, soil, distance to active fault, and distance to well, respectively, with Exp β values of 2.385, 1.880, 1.195, 1.119, 1.048, 1.007, 1.000, 0.998, 0.991, 0.957, 0.950, 0.687, 0.333, and 0.159 (Table 3).

Well density (WD) displayed the highest Exp β value, indicating its significant influence on sinkhole formation, particularly in areas with high artesian well density levels, ranging from 0.4-1.0 unit/km² in the eastern and central study areas. These regions, situated downstream in the KIPs, boast numerous artesian wells due to the natural groundwater flow from high to low areas. This phenomenon is prominent in densely populated districts like Chai Buri, Phrasaeng, and Khian Sa.

Geological variables are one of the important variables affecting sinkhole hazard and show strong negative β values. Since the geological variables are on a nominal scale, it is necessary to convert them to an ordinal scale and start from sedimentary rocks to igneous and metamorphic rocks. The study area is mostly covered by sedimentary rocks, including sedimentary rocks of the Carboniferous, Permian, Triassic, Tertiary, and Quaternary periods, with sedimentary rocks covering most of the area, and igneous rocks, including igneous rocks of the Cretaceous and Quaternary periods. From the LR analysis, the results of the geological variables show that the ancient sedimentary rocks affect the formation of sinkholes, especially the rocks in the limestone group that are in the Carboniferous and Permian periods.

Geological variables exhibit a notable Exp β value, albeit to a lesser extent, and a negative β value. The study revealed that grid cells with lower values in the study area predominantly consist of sedimentary rocks dating back to the Carboniferous, Permian, Triassic, and Tertiary periods. These rocks commonly contain limestone, dolomite limestone, dolomite, and gypsum formations, which are particularly prevalent in Permian-era formations (Pr). The study area features the Ratchaburi group (Pr), characterized by these rock types. Additionally, Carboniferous-Permian (CP) rocks include the Kaeng Krachan group (CPk), supporting the Pr series rocks. Some CP and CPk rock layers contain argillaceous limestone interspersed with other sedimentary rocks like siltstone and tuffaceous sandstone. Triassic-period sedimentary rocks include the Lampang group (Trl), featuring gray-black limestone. Tertiary sedimentary rocks, found in the Krabi group (Tkb), also comprise limestone, dolomite limestone, dolomite, and gypsum. Limestone predominantly covers the geological features in the study area, resulting in karst topography spanning from CPk, CP, and Pr formations continuously to the southern Phanom Bencha mountain range. The Trl and Pr sections cover the western mountain range with a northeast-southwest orientation, featuring scattered tower karst formations like Khao Lom Fang. The Tkb rock group encompasses low hills on the eastern side of the study area.

The LU variable ranks third in Exp β value after geology. It exhibits a positive β value, with higher values indicating urban, built-up, and agricultural areas. In the KIPs, agricultural land predominates, covering 1,805 km², accounting for 86.51% of the area. Much of this agricultural land is dedicated to oil palm and rubber plantations, covering 898.37 km² and 834.97 km², respectively. These agricultural practices involve land stripping and leveling, with some areas featuring drainage trenches. Such activities can render the soil surface fragile, and prone to subsidence, particularly in areas with underground burrows or karst terrain. Urban areas and buildings have a comparatively lesser impact on land subsidence, given the limited presence of buildings and infrastructure, affecting only specific areas.

After LU variables, the subsequent influential variables affecting sinkhole occurrence include slope, distance to stream, elevation, terrain wetness index (TWI), distance to village, soil, distance to active fault, and distance to well, respectively. Slope variable: Areas with moderate to high slopes, approximately 10-40 degrees, are prone to land collapse. Distance to stream variable: Proximity to waterways increases susceptibility to sinkhole formation. Elevation variable: Sinkholes are likely to occur at elevations ranging from approximately 10-100 meters above mean sea level. TWI variable: Low TWI values signify low terrain moisture levels, contributing to sinkhole formation. Distance from village variable: Areas not far from villages are particularly sensitive to collapse. Soil variable: Shallow soil-covered areas are at risk of subsidence. Distance from well variable: Proximity to artesian wells increases susceptibility to collapse.

Additionally, the distance to the active fault variable indicates the presence of the Khlong Ma Rui Active Fault in the northwest of the KIPs. Oriented northeast-southwest, this fault poses earthquake and tsunami risks, with potential tremors affecting karst topographic areas, which are also prone to sinkhole disasters.

The TH, TDS, and Cl variables showed positive β values. The results showed that higher concentrations of the chemical status of groundwater would result in an increased sinkhole probability.

Sinkhole susceptibility map in KIPs

According to the sinkhole susceptibility analysis in KIPs using the LR method by using the β value to create a database, it has been prepared as a map for disaster management at the sub-watershed level to represent the sinkhole susceptibility map (Fig. 7). The results of this study were to analyze spatial data in GIS as shown in Eq. 4.

(4)

The results of the β value of various variables make the findings a highlight of this research. That can express the level of risk as appropriate spatial data according to related variables and affecting the occurrence of sinkholes in the KIPs area in particular. Results of the study of KS susceptibility show that sinkhole susceptibility areas in the watershed can be classified into five levels: very high risk areas, high risk areas, moderate risk areas, low risk areas, and areas at very low risk of sinkhole, respectively (Table 4).

In the KIPs area, regions with a very high sinkhole susceptibility level cover approximately 399.86 km², constituting 19.16% of the total area. These high-risk zones are predominantly situated in the eastern part of KIPs, particularly near the convergence of Khlong I Pan and Khlong Trom streams. Additionally, mountainous areas, such as the central and western ranges, exhibit elevated sinkhole risk due to their karst geological formations. Surrounding the very high susceptibility zones, both high- and moderate-risk areas span over 421.73 km² (20.21%) and 422.18 km² (20.23%), respectively. Conversely, low- and very low-risk zones are primarily located in the southern region, notably in the upstream areas of KIPs along streams like Khlong Ya, Khlong Pho Thak, Khlong Chang Tai, and Khlong Sai Khao. These areas cover 430.41 km² (20.62%) and 412.85 km² (19.78%), respectively, with minimal sinkhole risk attributed to fewer artesian wells and the prevalence of modern sedimentary rocks, distant from active fault lines.

DISCUSSION

Sinkholes, influenced by complex factors, pose prediction challenges. Mitigation is aided by gathering diverse variables and using consistent analysis. Sinkhole susceptibility maps provide crucial spatial data for community readiness. This study selected 14 factors, integrated them into a GIS database for statistical analysis, and produced a map for effective local risk management.

Over the past two decades (2002-2022), sinkhole areas were studied using secondary data from sinkhole events, the Department of Mineral Resources of Thailand, and visual image interpretation. The study identified 34 sinkholes covering 3.52 km², concentrated along the Phanom Bencha mountain range. That area is an important upstream area of the Khlong I Pan tributary stream. Notably, a large subduction basin was discovered near Bankhoasamyot School in Surat Thani province. The area has an undulating plain topography.

As a result of this research, it can be seen that the technique for finding sinkhole areas using visual image interpretation is still a suitable method for finding actual sinkholes efficiently. Consistent with the research of Orhan et al. (2020) that used the principle of finding sinkhole areas with visual analysis of the aerial photographs of the district and on-site fieldwork. As a result, accurate sinkhole data was obtained and used to create a sinkhole susceptibility map using remote sensing and GIS processes in the Konya sub-watershed area, Türkiye. Even in regions with a tropical climate, such as the Metropolitan Region of Belo Horizonte, Minas Gerais, Brazil, sinkhole landforms can be searched using this technique by analyzing hydrograph control, topography, cross-section elevation profiles, and high-resolution satellite imagery (de Castro Tayer and Rodrigues 2021). As for the geologically sensitive region, the Zagros Mountains, Kermanshah Province, Iran, the area is the largest active fault in the Zagros overthrust. Techniques for finding sinkholes and subsidence areas using visual image interpretation of World Imagery and Google Earth can accurately display the actual sinkhole (Maleki et al. 2023). There are some weaknesses that must be noted with this analysis technique, especially human-assisted image processing which requires the individual’s experience (Dou et al. 2015). Certain regions may require reliance on visual analysis of remote sensing imagery and on-site research, despite the drawbacks of subjectivity, time intensiveness, limited reproducibility, and inadequacy for examining extensive territories. (Chen et al. 2018). Future research processes may require technology that can process data quickly and with high efficiency, using machine learning to create models that can discover and detect actual KS for collapse areas. Artificial intelligence should be trained to learn from the processing of visual image interpretation from experienced people frequently to obtain data on the collapse of the terrain caused by sinkholes.

Sinkhole occurrence factors in the KIPs area were analyzed using LR statistical methods. The study results found that socio-economic factors are sensitive to landscape collapse. In particular, prominent variables such as WD and LU show strong Exp β values at levels of 2.385 and 1.195, respectively. As for the physical factor variables that are sensitive to the collapse of the terrain, it was found that the geology and slope variables showed strong Exp β values at the levels of 1.880 and 1.007, respectively. It can be seen that the study area has a high level of artesian WD ranging from 0.4-1.0 unit/km² that is distributed in the eastern area continuing into the central part of the study area. Most of these areas have LU conditions where the area has been opened for agriculture to grow crops such as rubber, palm oil, and rice fields, further accelerating the collapse of the land more easily. It can be seen from the KS susceptibility map that the area is classified as the area with the highest risk of collapse. It covers most of the eastern area, continuing into the central part of the study area. This is different from the findings of Cahalan and Milewski (2018) which indicated high aquifer fluctuations and shallow overburden thickness. This is caused by activities on the topographic surface that are the primary variable affecting the subsidence of the land. In addition, the original geological landform with karst topography covers the central and western mountain ranges in the study area. It creates greater susceptibility to the collapse of the terrain in areas with very steep slopes. There are different variables that affect the occurrence of sinkholes in each area. These findings can be seen from the research of Wood et al. (2023) that sought to find sinkhole vulnerability in karst and pseudokarst regions of the contiguous United States, both present and anticipated. They identified important variables such as the thickness of the soil layer, soil type, and humidity that are the main variables that affect collapse until a sinkhole occurs, which is different from the results of this research.

The geochemical variables used in the study area were TH, TDS, and Cl. It was found that if the concentration level of these variables is high, it easily affects the formation of sinkholes. This is different from the results of the study by Ozdemir (2016) who stated that if the concentration level of Calcium (Ca) and Magnesium (Mg) increases, the probability of sinkholes decreases. However, the results of this study have shown that the TH variable is an important variable that affects the formation of sinkholes. The concentration levels of Ca and Mg were found in these variables. This conclusion is consistent with the study approach of Amin et al. (2023) who stated that Ca and Mg are likely to be dissolving carbonate rocks from limestone and dolomite. In the study area, the amount of this element was found to be quite high. It is possible that high levels of TH found in groundwater indicate that the underground caves have been severely eroded, resulting in high levels of Ca and Mg in the groundwater. The variables TDS and Cl were found to have little influence on sinkhole formation.

However, LR analysis has the capability to manage both continuous and categorical variables concurrently, which is not required in standard distributions. Additionally, the LR analysis technique can accurately identify regions prone to various geological hazards like KS (Zhou et al. 2016; Ozdemir 2016; Subedi et al. 2019; Kim et al. 2019). In summary, the overall picture from the discussion of the results reveals the differences in variables that affect the occurrence of sinkholes that occur in different contexts in each area. However, with the variables used in the LR analysis process, the results of the β value coefficient of the variables used in the study will vary according to the physical characteristics of the area. Therefore, if studying KS susceptibility in other areas, one should be aware of the factors that will be studied first in areas at risk of sinkhole formation.

KS susceptibility mapping is vital for disaster planning, yet the 2005 sinkhole risk map from the Department of Mineral Resources of Thailand lacks transparency and updates, making it ineffective for local disaster management. This study aimed to improve data resolution, revealing that nearly half of the study area has high susceptibility to sinkholes. Contrary to the outdated map’s three risk levels, this research identified five levels, reflecting nuanced susceptibility. Discrepancies between actual sinkhole occurrences and the department’s map highlight its inaccuracy. Local-level mapping provides crucial risk insight for communities and authorities, particularly in high-risk areas, aiding disaster preparedness and mitigating potential losses.

CONCLUSIONS

Sinkholes are unpredictable phenomena, often occurring in karst landscapes. This study aims to mitigate their impact by analyzing various factors contributing to sinkhole formation in the local area. Statistical methods, including LR analysis combined with GIS, were employed to create susceptibility maps for sinkholes. The research identified several key variables affecting sinkhole formation, such as well density, geology, land use, slope, distance to streams, elevation, terrain wetness index, proximity to villages, soil type, distance to active faults, and distance to wells. These variables were quantified using coefficient values and incorporated into the susceptibility mapping process. The findings revealed that over half of the study area exhibited a very high or high susceptibility to sinkholes. This research provides empirical evidence supporting the effectiveness of detailed, community-level sinkhole risk assessment, potentially replacing existing sinkhole risk maps maintained by the Department of Mineral Resources in Thailand. Future studies should expedite the exploration of sinkhole susceptibility across Thailand, considering the unique variables influencing sinkhole formation in different regions. Government agencies involved in disaster management should prioritize comprehensive research efforts to prepare for and mitigate sinkhole risks nationwide, ensuring the safety and sustainability of local communities in the future.