Detecting and Analyzing the Evolution of Subsidence Due to Coal Fires in Jharia Coalfield, India Using Sentinel-1 SAR Data

Abstract

Public safety and socio-economic development of the Jharia coalfield (JCF) in India is critically dependent on precise monitoring and comprehensive understanding of coal fires, which have been burning underneath for more than a century. This study utilizes New-Small BAseline Subset (N-SBAS) technique to compute surface deformation time series for 2017–2020 to characterize the spatiotemporal dynamics of coal fires in JCF. The line-of-sight (LOS) surface deformation estimated from ascending and descending Sentinel-1 SAR data are subsequently decomposed to derive precise vertical subsidence estimates. The most prominent subsidence (~22 cm) is observed in Kusunda colliery. The subsidence regions also correspond well with the Landsat-8 based thermal anomaly map and field evidence. Subsequently, the vertical surface deformation time-series is analyzed to characterize temporal variations within the 9.5 km² area of coal fires. Results reveal that nearly 10% of the coal fire area is newly formed, while 73% persisted throughout the study period. Vulnerability analyses performed in terms of the susceptibility of the population to land surface collapse demonstrate that Tisra, Chhatatanr, and Sijua are the most vulnerable towns. Our results provide critical information for developing early warning systems and remediation strategies.

1. Introduction

Subsurface coal fires constitute a significant socio-economic and environmental threat in at least 30 countries worldwide [1,2,3]. These coal fires are essentially self-smoldering coal seams triggered by exothermic reactions between coal and oxygen [4]. In addition to the loss of the largest non-renewable energy resource and associated economic impacts, coal fires have substantial adverse consequences on the local and global environment. Coal fire produces abundant toxic gases such as nitrogen oxides (NO_x), sulfur dioxide (SO₂), carbon oxides (CO_x), nitrous oxide (N₂O), and methane (CH₄) continuously, which also act as greenhouse gases [5,6,7,8]. While these gases cause severe respiratory and coronary diseases, thousands of people die by coming in direct contact with coal fires and coal fire-triggered explosions in mining environments [9]. Land surface collapse due to coal fire-triggered subsidence is also a major threat that took many human lives and remains a major cause of concern to the local administration and mine managers. Subsidence mainly initiates with vacant spaces being formed between subsurface rock layers due to the loss of coal volume as they burn to ash. As long as the coal fires persist, the extent and magnitude of these subsidence continue to increase and eventually reach hazardous proportions. Additionally, associated with land subsidence are the formation of deep cracks and ground fissures [10,11]. These cracks act as pathways for supplying oxygen to the subsurface coal fires to fuel them and also transfer toxic gases to the surface. Apart from the surface collapse hazards, subsidence causes severe damage to infrastructures such as road and railway networks, settlements, and buildings [7,12]. Besides, coal fires source various environmental and social catastrophes, including land degradation, vegetation loss, water pollution, loss of natural habitats, and challenges in mining activities [13]. Hence, proper monitoring of subsurface coal fire is key to effective coal fire control and efficient implementation of preventive measures. Unless these coal fires are controlled, major coal producing countries will continue to face economic losses in terms of wastage of coal resource, cessation of mining activities, and extensive damage to society and infrastructure.

If we consider the Jharia Coal Field (JCF) of India itself, until 2003, more than 37.6 million tonnes of coal have been lost to coal fires, which accounts for more than 3.87 billion US dollars according to current rates. Besides this, around 1.48 billion tonnes (accounts for more than 155.33 billion US dollars) of coal is inaccessible for mining due to the coal fires in JCF [11]. In addition, damaged infrastructures such as subsiding railways and cracked buildings have significantly impacted the economy and social lives. To date, more than 1300 people have died from the direct impacts of coal fires in India, and several thousands are exposed to toxic gases emanating continuously from these coal fires [9]. Thus, precise information on the occurrence and dynamics of coal fires are crucial for developing mitigations strategies and early warning systems in coal mining areas, particularly in India, where the population density is the highest amongst all the coal mining countries in the world [14,15].

Advancements in remote sensing techniques have contributed significantly to tackling these challenges by facilitating large-scale coal fire monitoring [16,17,18]. Diverse remote sensing methods have been used in coal fire monitoring, including analyses of subsidence related to coal fires [19,20], identification of surface temperature anomalies [21,22,23,24], and mapping of geo-environmental indicators of coal fires [25,26]. Remote sensing of coal fires is predominantly based on identifying thermally anomalous regions using thresholding techniques on thermal infrared (TIR) images [27,28,29,30,31]. The primary limitation of this method is the ambiguity in defining a temperature threshold that differentiates between coal fire pixels and background pixels [32,33,34]. Additionally, surface thermal anomalies due to coal fires are substantially impacted by surface cracks and vents [35]. Hence, the surface thermal anomaly may or may not be a direct indicator of the coal fire existing vertically below [12]. Similarly, mapping geo-environmental indicators like rock deformation and mineral deposition require extensive fieldwork, and not all coal fires produce noticeable indicators on the surface [11]. Consequently, mapping geo-environmental indicators are conducted for specific scientific purposes rather than large-scale coal fire mapping [36,37].

As an alternative, the potential of Interferometric Synthetic Aperture Radar (InSAR) was explored for coal fire mapping [38,39,40]. While overcoming the limitations of TIR-based studies and mapping geo-environmental indicators, InSAR also delivers additional advantages in coal fire mapping. All-weather capability, better seasonal stability, high-resolution measurements, and better scope for risk assessment make InSAR techniques unique from other coal fire mapping methods [11]. However, InSAR techniques are vulnerable to decorrelation and atmospheric artifacts, especially in dynamic regions like an active coalfield [12,41,42]. Coal fire regions are specifically prone to decorrelation because of continuous variations in surface topography due to surface collapses and opencast mining activities.

More recently, InSAR time-series techniques have been developed utilizing persistent scatterer (PS) points and distributed scatterer (DS) points. DS-based techniques, such as Small BAseline Subset (SBAS), usually yield better results in non-urban regions because of the availability of dense DS points [43,44]. On the other hand, PS-based studies are prone to overlook small deformation patches in non-urban regions, such as coalfields, due to sparse PS points [45,46]. In addition to limiting some of the constraints, InSAR time-series measurement also offers a platform for detailed temporal analyses [47]. However, despite the advantages, the potential of InSAR time-series for coal fire mapping has been rarely investigated and limited to Wuda Coalfields of China [19,20,48]. Furthermore, by default, InSAR techniques deliver surface deformation in a one-dimensional plane, which is in the line-of-sight (LOS) geometry of the sensor [49]. As a result, horizontal surface movements influence the LOS deformation measurements significantly. Even though the relevance of deriving the vertical surface deformation component has been explored in various geophysical studies [50,51,52,53], it is yet to be implemented for mapping coal fires.

In this study, we compute dense interferograms to generate a time-series of surface deformation using New-SBAS (N-SBAS) technique [54,55] for precise identification and characterization of coal fires in Jharia coalfield (JCF). N-SBAS is a modified SBAS processing chain optimized for monitoring surface deformation of small amplitudes over a longer time period in natural settings (especially non-urban regions such as coalfields) [54]. Sentinel-1A C-band SAR datasets acquired from April 2017 to April 2020 are used to develop the surface deformation time-series. We also compute vertical surface deformation by decomposing ascending and descending path LOS surface deformation measurements (hereafter, ascending and descending, respectively). Subsequently, the vertical surface deformation time-series is utilized to categorize coal fire regions according to time variations in subsidence rate. We also compare the N-SBAS based subsidence map with the more traditionally used coal fire map computed using Landsat-8 TIR image. The derived results are also validated using GPS locations of coal fires encountered during the field visits. Finally, village-level vulnerability to coal fire-triggered hazards is assessed using the derived subsidence map.

2. Study Area

JCF is the most productive coalfield and the only source of prime coking coal in India [11] (Figure 1). Simultaneously, JCF hosts one of the most extensive coal fires in the world and is also the most socially vulnerable coalfield due to the high population density [14,15,56]. JCF belongs to the Damodar valley basin of the Lower Gondwana group [57,58]. The Gondwana rocks in JCF formed between Lower Jurassic and Lower Permian periods are represented by Barakar Formation, Barren Measures Formation, Raniganj Formation, Talchir Formation, and Archean rocks [58]. The current characteristic of JCF, which include shallow coal seams, mining of multiple thickness coal seams (0.91 to 22.44 m), and the presence of various abandoned mines, favor the propagation of subsurface coal fires [59]. Moreover, the long duration and high density of coal fires in JCF have made an intense network of cracks that supply sufficient oxygen to sustain and fuel new coal fires [12]. Besides, mining operations in JCF are mostly opencast, which also provides conditions favorable for sustaining coal fires [11].

Coal mining in India started in 1774, and coal fires have been reported since 1916 [15]. In addition to the economic loss, coal fires in JCF have severely affected the local environment and livelihood through various pollution and catastrophic events, including fatal land surface collapses [13]. However, the coal industry got nationalized by 1973, which significantly improved the overall quality of coal mining in India and restricted public exploitations [1]. Since coal is the source of around 59% of India’s electricity and a necessity for major industries such as steel and cement, loss of coal reserves and challenges in coal extraction due to coal fires have significant ramifications to India’s economy [11]. The JCF covers a total area of 231 km², the northern parts belonging to the Barakar Formation consist of shallow coal seams. Because of the proximity of coal seams to the ground surface, dense opencast mining activity is found in these regions (visible in violet color gradient in Figure 1). In comparison, the coal seams are deeper in the southern portions of the coalfield. The only exception is Parbatpur, an isolated part of the Barakar Formation that hosts shallow coal seams.

3. Methods and Materials

The processes involved in forming the vertical surface deformation time-series are portrayed in a flow-chart (Figure 2). The vertical surface deformation is formed by decomposing ascending and descending LOS surface deformation time-series. Since the procedures involved in computing the surface deformation time-series of ascending and descending are the same, only ascending time-series formation is illustrated on the flow-chart. Detailed descriptions of each procedure are explained below.

3.1. Interferometry

We utilize Sentinel-1A SAR SLC datasets acquired in interferometric wide (IW) swath mode for studying the surface deformation in JCF. Three years of dense datasets (88 datasets for ascending and descending each) from April 2017 to April 2020 (hereafter 2017 to 2020) are utilized in this study. SAR data pairs having low temporal (maximum of 36 days) and perpendicular baselines (maximum of 177 m) are chosen to produce interferograms for better coherency and sensitivity. Since Sentinel-1 IW data are acquired using TOPS burst synchronization technique, it demands significantly high accuracy in co-registration of SAR pairs. A network-based, enhanced spectral diversity technique [60] is utilized for the co-registration process to overcome this. We use precise orbit ephemerides and SRTM 1 Arc-Second DEM for interferometry and topographic phase removal. Multi-looking (4 azimuth × 25 range) and Goldstein phase filtering are also applied on the interferograms for quality enhancement. Finally, the interferograms are unwrapped by Statistical-cost, Network-flow Algorithm for PHase Unwrapping (SNAPHU) algorithm [61]. ISCE 2.2 software [62] is used to produce 426 interferograms (212 in ascending and 214 descending).

3.2. Surface Deformation Time-Series

The time-series of surface deformation in ascending and descending LOS from 2017 to 2020 are computed using N-SBAS technique. N-SBAS consists of three steps: Interferogram stacking, pre-processing of the stack, and inversion. In the stacking process, unwrapped interferograms are stacked together as a data cube containing necessary auxiliary information. In the pre-processing stage, the phase delay due to atmospheric interactions and deramping caused by orbital errors are computed and corrected. European Centre for Medium-range Weather Forecasts (ECMWF) weather data are used to compute the atmospheric phase delays. This computation is executed using PyAPS module [63] integrated with GIAnT software. However, the actual time-series is computed in the inversion stage by converting the data cube to small subsets and connecting them temporally to form the entire time-series.

After producing the time-series of surface deformation, they are analyzed thoroughly to identify sudden jumps in the time-series. These jumps are characterized as a sudden hike or drop in a continuous time-series of surface deformation. They are caused mainly by the limitations in unwrapping and inversion procedures. After correcting jumps in the time-series, cumulative surface deformation maps in ascending and descending LOS are prepared, portraying the total surface deformation in JCF from 2017 to 2020.

3.3. Decomposition

By default, InSAR computes the surface deformation in the LOS direction, which is a composite of surface deformation in every direction at each point. The surface deformation at any point can be represented using a three-dimensional (3D) vector model. Each axis of the 3D model represents the surface deformation in one direction, i.e., vertical, East–West (EW), and North–South (NS). Most importantly, due to its polar orbit, Sentinel-1 images are not sensitive to surface deformation in the NS direction unless the surface deformation in NS is substantially high compared to the other two directions [49,51,53]. Thereby, unless any geophysical phenomenon causing substantial surface deformation oriented in NS direction, LOS measurements from Sentinel-1 can be assumed to consist of surface deformation in vertical and EW direction. From the field knowledge and literature survey, we do not find any phenomenon in JCF that causes substantial subsidence only in NS direction [38]. Hence, we have modified the 3D model to a 2D model by neglecting the surface deformation in NS direction. Further, the surface deformation in the vertical direction is derived by the decomposition of surface deformation measured in ascending and descending LOS directions [64,65], as mentioned below.

Let us consider ($\widehat{p}$), be the unit vector pointing from satellite to target on ground in local east (${\widehat{p}}_x$) and vertical (${\widehat{p}}_z$) direction as $\widehat{p}=[{\widehat{p}}_x\hspace{1em}{\widehat{p}}_z\left]=\right[\cos \varphi \sin \lambda -\cos \lambda ]$. Here $\varphi$ is satellite heading azimuth, and $\lambda$ is incidence angle at the ground target location. Let the displacement vector be represented as $d=[d_x\hspace{1em}{d}_z]$ in the same local reference frame. Hence, the observed surface deformation along LOS ($d_{los}$) for any given pixel is represented by Equations (1) and (2).

$$d_{los,a}={\widehat{p}}_a.d=[{\widehat{p}}_{x,a}\hspace{1em}{\widehat{p}}_{z,a}\left].\right[d_x\hspace{1em}{d}_z]$$

(1)

$$d_{los,d}={\widehat{p}}_d.d=[{\widehat{p}}_{x,d}\hspace{1em}{\widehat{p}}_{z,d}\left].\right[d_x\hspace{1em}{d}_z]$$

(2)

Here, the subscripts ‘a’ and ‘d’ denotes ascending and descending path respectively. These equations are resolved to estimate the surface deformation in vertical ($d_z$) direction for each pixel.

3.4. Thermal Anomaly Mapping

For comparison purposes, a coal fire map is also prepared by analyzing surface thermal anomalies in JCF using a wintertime Landsat-8 TIR image (Band 11). Before using this image for analysis, each pixel’s digital number (DN) value is converted to land surface temperature using the RS&GIS plugin in QGIS software. Since coal fires produce elevated temperatures on the surface, pixels covering those regions will exhibit excessive temperatures in TIR images. Therefore, thermally anomalous pixels in space-borne TIR images are considered as indicators of coal fire locations [66,67]. In this study, the thresholding method [68] is followed for coal fire mapping, in which a specific threshold temperature is identified from the TIR image histogram. The threshold temperature differentiates between coal fire pixels (thermally anomalous) and the background pixels (non-coal fire region). A significant challenge in the thresholding method is the conflict between coal fire pixels and high-temperature background pixels like those represent barren rocks. Since both the objects exhibit high-temperature in TIR images, high-temperature background pixels may be misinterpreted as coal fire. A TIR image acquired in the winter period (16 February 2020, Path:140, Row:44) is used to counteract this effect since the background pixels will not exhibit temperature comparable to those in coal fire regions.

4. Results

4.1. Surface Deformation Time-Series

The cumulative surface deformation in JCF from 2017 to 2020 is shown in Figure 3. Surface deformation in ascending (Figure 3a) and descending (Figure 3b) are computed separately using N-SBAS technique. Negative and positive values of surface deformation imply subsidence and uplift, respectively. However, surface deformation values between −30 mm and 30 mm (±10 mm/year) are considered background values. Regions with background values are typically those pixels that are not experiencing any surface deformation. The selection of this specific range for background pixels is based on the sensitivity/uncertainty measurement (Appendix A). Thereby, pixels with subsidence value less than 30 mm (for three years) are ignored in further analyses. Most importantly, significant sporadic changes on the surface will result in decorrelation, which the InSAR technique will not detect. Hence, surface collapses or surface deformations caused by opencast mining activities will not be construed as subsidence in the surface deformation measurements presented here. Moreover, even if some of the sudden changes are recorded in InSAR, they will be omitted while performing jump-corrections in the subsidence time-series.

While regions with significantly high subsidence are conspicuous and distributed sporadically throughout the JCF, considerable uplift is not observed in either deformation maps. The maximum cumulative subsidence in ascending and descending are 268 mm and 202 mm, respectively. Similarly, in most pixels, the magnitude of subsidence in ascending is considerably higher than in descending (Appendix A). However, few subsidence regions do not follow this pattern. For example, the highest subsidence in Parbatpur village in ascending mode is only 109 mm, while the same in descending mode has reached up to 155 mm. Additionally, some subsidence regions are found only in either of the maps. Even though slight differences between ascending and descending are typical due to the difference in LOS, significant differences between them indicate the role of horizontal movements in JCF. Subsidence in LOS direction depends on vertical and horizontal movements at that location, and the horizontal component varies according to LOS direction, resulting in variations in derived subsidence estimates in different LOS.

In both the maps (Figure 3), the highest subsidence is observed in the Kusunda colliery. This region exhibits significantly high subsidence (up to 266 mm) compared to other subsidence locations, and it is also the largest subsidence area (2.18 km²) in JCF. Apart from Kusunda colliery, major subsidence zones are observed in the Eastern JCF (up to 152 mm), Parbatpur (up to 156 mm), western JCF (up to 120 mm), and central JCF (up to 160 mm) (Figure 3).

4.2. Decomposition

In order to eliminate the influence of horizontal movements in the LOS surface deformation map, vertical surface deformation is derived by decomposing ascending and descending measurements (Figure 4). Map symbology used in Figure 3 are followed in this Figure 4 as well. However, the range of subsidence has considerably reduced compared to Figure 3. Further, the bias in the vertical deformation map is derived by analyzing 125 known geologically stable points. A mean bias of +8.82 mm with a standard deviation of ±2.43 mm is identified (Appendix B), and the bias value is subtracted from the vertical deformation. The estimated bias value of 11.25 mm (bias + standard deviation) is below the SBAS sensitivity of 30 mm chosen in this study. This bias is attributed to the errors involved in datasets used and procedures followed. By this estimate, we ensure that the subsidence regions considered in this study are not pseudo-subsidence regions that arise due to the errors in datasets used and processes involved. Most subsidence zones exhibit maximum cumulative subsidence of less than 100 mm in the vertical deformation map (Figure 4).

In contrast, vertical subsidence in Kusunda colliery (up to 219 mm) shows only a subtle difference compared to the subsidence estimated in LOS measurements (up to 240 mm) (Figure 3). Thus, it is inferred that, unlike other subsidence zones, the intensity of vertical subsidence in Kusunda colliery is significantly higher than the horizontal movement. In addition to Kusunda colliery, high subsidence is also observed in other collieries, including Kujama (east), Parbatpur (south), and collieries on the western parts of JCF. It is also noticed that the number of decorrelated pixels has increased in the vertical surface deformation map. It is an inherent limitation of the decomposition method since decorrelated pixels in either (ascending or descending) of the LOS measurements will result in the same in the products derived by decomposition. Thereby, the vertical deformation map contains decorrelated pixels from both ascending and descending surface deformation maps.

Additionally, we have identified that most of the subsidence regions are located in the vicinity of opencast mines using Google Earth. A primary exception to this is the subsidence area near Munidih. From field studies, we have identified that it represents the current working locations of Munidih underground (UG) mines. Munidih is the only UG mine in JCF that follows the longwall mining technique with a working depth of up to 650 m. Moreover, subsidence due to UG mining activities in Munidih is reported by the officials [70].

4.3. Temporal Variations

The contribution of groundwater over-drafting towards the observed subsidence is neglected since the derived subsidence regions are distributed discretely in small areas. Moreover, the outline of the subsidence zones does not follow the characteristic circular shapes [71,72]. Furthermore, the aquifers in the region are mostly unconfined, with shallow groundwater levels ranging from 1.4 to 11.1 m below ground level (mbgl) [73]. Since there are no other probable causes, subsidence regions identified in JCF are either due to coal fires or UG mining activities. Regions with subsidence due to UG mining are ascertained by field visits and by locating mine ventilation shafts in the vicinity of subsidence zones in Google Earth images. Furthermore, subsidence in UG mines are generally abrupt and of larger magnitudes. These instantaneous and inconsistent variations are identified as jumps in a surface deformation time-series, and they have been omitted from the analysis presented here. Therefore, it can be inferred that the rest of the subsidence regions identified in this study are caused by coal fires only. Consequently, temporal changes in the extent and magnitude of these subsidence zones are essentially representative of the spatial propagation and intensity variations of coal fires.

We have used N-SBAS technique to produce a time-series of vertical surface deformation to investigate spatio-temporal variations in subsidence due to coal fires in JCF. Before analyzing the temporal variations in subsidence, regions with cumulative subsidence of up to 30 mm (10 mm/year) are ignored according to the estimated background pixels range (Appendix A). Hence, only pixels having subsidence values more than 30 mm are considered subsiding regions. We have identified a total of 9.5 km² area subsiding due to coal fires. Further, the cumulative subsidence (comprising of three years) is split into three yearly subsidence datasets. Subsequently, based on the temporal variations in yearly vertical deformation, subsidence regions are categorized as ‘new’, ‘dormant’, ‘persistent’, ‘increased’, and ‘decreased’. The categorization is done by comparing the temporal variation from the initial condition (2017) to the final condition (2020). The pixels classified as ‘new’ are the newly formed subsidence spots that were not subsiding in 2017, while dormant represent stabilized regions that were subsiding in 2017. The ‘increased’ regions represent those areas registering continuous subsidence throughout the observation period with an increasing trend in subsidence rate. Similarly, ‘decreased’ regions are those registering continuous subsidence with a reduction in subsidence rate. On the other hand, the category ‘persistent’ characterizes regions where subsidence continued at a consistent rate throughout the observation period.

The Barakar Formation of JCF consists of shallow coal beds [74], making this formation more susceptible to coal fires. The occurrence of coal in shallow depths led to dense opencast mining and attracted unauthorized mining activities [75], which triggered new coal fires and enabled the continuous availability of oxygen to sustain coal fires [11]. This feature is also evident in Figure 5a, in which all the coal fire areas are located within the Barakar Formation. Most importantly, most of the subsidence regions are identified in opencast mine regions, especially near larger mines (e.g., Kusunda, Parbatpur). It is well-established that coal fires will persist for years after ignition, and the fire propagation is slow [15]. Agreeing with them, we discovered that 72.6% of the total coal fires in JCF persisted with similar intensity during the study period, while 11.7% of them have registered varying intensity (increased or decreased). Most importantly, despite intense coal fire prevention activities, new coal fire regions (9.8%) occupy nearly double the area of dormant coal fire regions. Areas with a reduction in coal fire intensity are concentrated in the Kusunda colliery region, while other categories of temporal variations are observed throughout JCF. It is also noticed that new coal fire spots are sparse in the Kusunda region and western parts of JCF. For better visualization, subsidence in the Kusunda colliery region is magnified and shown separately in Figure 5a.

Shown in Figure 5b–f are the time-series of vertical subsidence at five coal fire locations (marked B through F, respectively) in Figure 5a. Each of these locations represents a characteristic temporal behavior, as categorized earlier, and the category name is mentioned inside the subfigures. The best-fit line for temporal variation in subsidence for each year is shown with different colors. For example, Figure 5b represents a newly formed coal fire area while Figure 5d represents a coal fire region that was already present, but the rate of subsidence has increased. Even though these two subfigures look identical, the best-fit line of the first year (red line) shows a clear difference between the two categories. While the slope of the red line in Figure 5d exhibits a subsidence trend, that in Figure 5b represents a stable ground (parallel to the x-axis) that started subsiding at a later stage. Similarly, it is evident from Figure 5f that during the final year, the region has stabilized (horizontal magenta line), representing a dormant coal fire location. Whereas the same region in the similar-looking Figure 5e continued to subside, though at a reduced rate. In contrast, Figure 5c exhibits a continuous subsidence pattern over all the years, indicating a continuous coal fire region with a persistent subsidence rate. Figure 5c also showcases the intense coal fire activity in the Kusunda colliery region, characterized by significantly higher subsidence magnitudes than the other notable coal fire regions. Thus, dense input data has made the subsidence time-series significantly detailed to understand even short-term variations in subsidence models. It has also enabled to avoid corrupted subsidence values such as the few isolated blue point values in Figure 5b–f that are uncorrelated with other points.

4.4. Coal Fire Mapping Using TIR Image

A thermal anomaly map of JCF is prepared using a Landsat-8 TIR image following the thresholding method [68]. A TIR image acquired in the winter period is used to avoid conflicts between coal fire pixels and high-temperature background pixels. However, the final threshold value of 25.6 °C is chosen by the trial-and-error method, which involves multiple trials of varying threshold values to determine an optimal threshold. In the optimal threshold value, maximum coal fire areas are detected while the noise is negligible. Noise in the thresholding method results in small clusters of pixels exhibiting high-temperature being erroneously categorized as coal fire-affected regions. This noise is identified based on the small area of occurrences and their occurrences in a non-feasible location (example: inside a forest region). From the coal fire map derived from the TIR image, we have identified 9.07 km² of thermally anomalous areas in JCF. In order to investigate the correspondence between thermal anomaly and subsidence due to coal fires, outlines of thermally anomalous (>25.6 °C) areas are overlayed on vertical subsidence regions attributed to coal fires and UG mining separately (Figure 6a). The map also shows 13 GPS locations of coal fire sightings recorded during field visits to the JCF. These readings represent the safest possible access points from where coal fires were visible. Nearly all the points are lying within or along the boundaries of the subsidence zones. Additional photographs taken during the field-validation are shown in Appendix C.

It is evident from Figure 6a that subsidence zones due to coal fires are either within or close to thermally anomalous regions. The largest and the most conspicuous subsidence zone (Kusunda colliery) also correspond to the largest thermal anomaly. Other thermally anomalous regions are distributed throughout JCF with smaller spatial extents. It is also important to note that thermal anomalies are not found in the vicinity of any subsidence regions caused by UG mining. Most importantly, the areas where thermal anomaly does not coincide with subsidence are located mainly on the boundary of opencast mines. These boundaries are characterized by steep slopes that often cut across subsurface layers as part of the mining activities. As a result, UG channels through which hot gases from coal fires get transported to the surface are exposed, and the hot gases are emitted continuously to the surface through them. This results in surface temperature anomalies in those locations (Figure 6b). These inferences are also evident in the photographs taken during field surveys (Figure 6c,d).

Thermal anomalies are usually observed when underground channels transporting hot gases are exposed to the surface (Figure 6c). The authors of [12,35] had reported that surface thermal anomalies depend more on cracks and fissures on the ground rather than the occurrence of a coal fire burning below. The rationale for such observations is that the fast convective heat reaching the surface through cracks and fissures produces a more substantial impact on the surface temperature anomaly than the slow conductive heat passing through overlying bedrocks. As a result, hot gases passing through cracks and fissures are released to the surface and make that area thermally anomalous in TIR images. Therefore, even though coal fires are causing these surface thermal anomalies, the regions identified as thermally anomalous may not exactly correspond to the occurrence or extent of subsurface coal fires. This scenario is apparent from the map (Figure 6a) since most thermally anomalous regions are associated with a nearby subsidence zone, even though they may not overlap. Similar phenomena have been reported by the authors of [76] while integrating surface thermal anomalies and subsidence regions to derive the coal fire map. They have observed that most of the coal fire locations identified in field surveys are either nearby or on the boundaries of thermally anomalous regions identified in TIR images. Moreover, coal fires may take years to produce observable temperature differences on the surface [25,28]. Due to this delay in producing thermal anomalies on the surface, the spatial extent of the actual coal fire may significantly differ from the extent of the surface thermal anomaly. In some cases, coal fires will not cause any surface temperature anomaly [11].

It is important to mention here that there is no official database for the locations of subsurface coal fires in JCF that is publicly available, except for those reported by local newspapers about accidents and fatalities caused by coal fires. Conventionally, thermal anomaly maps derived using remote sensing techniques have been used as the best possible representation of subsurface coal fires in JCF. Hence, based on remote sensing and field evidence presented in this study, we can establish beyond reasonable doubt that the noted subsidence is due to subsurface coal fires. Moreover, no other geophysical phenomenon has ever been reported in the JCF that can explain such high rates of subsidence observed in this study.

4.5. Vulnerability Assessment

The large population and dense coal fires make JCF exceptionally vulnerable to various hazards compared to other coalfields in the world [77]. Village/town-level population in JCF is investigated to identify their vulnerability to coal fires. According to the Census of India, towns are defined as regions with a minimum population of 5000, having a minimum population density of 400 per km², and at least 75% of the male working population is engaged in non-agricultural jobs. Due to the non-availability of village-level administrative boundaries, the population distribution in JCF is represented by circles of varying diameters centered at specific point locations of the town or village. However, to accommodate for the significant difference in the population of towns and villages, they are represented by two different scales on the map. For towns, the diameter of representative circles (in meters) is determined by dividing the town’s population by 100. Whereas the diameter of circles representing village populations are determined by dividing the respective population by 10. This scaling factor also considers the population density, which is significantly higher in towns than villages. Therefore, a circle representing a town depicts a population that is ten times more than a village represented by a same-sized circle. This also implies that people living in towns close to coal fires are at greater risk than those living in villages in similar proximity to coal fires. After creating the circles representing the spatial distribution of the population, villages and towns within 1 km from the boundary of subsidence zones are identified and portrayed on a map (Figure 7).

While some subsidence zones are overlapping with populated regions, most others are close to populated areas. While the western parts of JCF are comprised of villages mostly, the eastern parts are dominated by towns having higher population density. Towns like Tisra, Chhatatanr, and Sijua are identified as the most vulnerable regions in JCF because of their large population and proximate coal fires. On the contrary, the subsidence zones in the western parts of JCF are surrounded by populated villages, in which some villages are already in the coal fire zones. While most rail and road networks are away from the subsidence zones, few railway stations such as Bansjora and Jharia are located close to subsidence regions. Most parts of the road networks are found safely away from the subsidence regions except in the western part of JCF. However, the map has shown only the major road networks while the rural road network is dense in JCF.

5. Discussion

Here we present the first study that utilizes vertical deformation estimates from InSAR to study subsurface coal fires instead of LOS measurements of surface deformation. It is also the first study in the Indian context, which utilizes any SBAS time-series technique to measure the subsidence due to subsurface coal fires. Furthermore, results from this study have demonstrated the advantage of using DS points to measure coal fire-induced subsidence. In PS-based studies, the detection of small-area subsidence regions will be limited due to the difficulty of obtaining sufficient PS points in a dynamic coalfield. The N-SBAS derived dense time-series of surface deformation has provided vital information to monitor continuous slow-rate subsidence. In addition to restricting atmospheric artifacts and decorrelation issues, N-SBAS also helped to understand the temporal changes in subsidence zones. Furthermore, N-SBAS overcomes the limitation of missing data in pixels due to decorrelation produced by sudden surface collapses and opencast mining activities. In N-SBAS technique, surface deformation before and after the event (causing the decorrelation) are connected temporally to estimate the total surface deformation in those pixels. Thereby it omits the effects (decorrelation and surface deformation) of the event causing decorrelation, which took place only during a specific period in the total study period.

This study also demonstrates the crucial role of subsurface geology in the occurrence and propagation of subsurface coal fires. The Barakar Formation holding shallow coal beds is more prone to spontaneous combustion (self-ignition) of subsurface coal by availing sufficient oxygen supply. Opencast mines stimulate coal fires and damage the local environment by subsidence, land and vegetation degradation, and emission of toxic gases [13,78]. Compared to opencast mines, UG mines cause less impact on the environment and less chance to trigger new coal fires. Additionally, the high population in the Barakar Formation imposes immense pressure on local livelihood. These facts make Barakar Formation at high risk, which is well established by past cave-ins and fatal accidents [79,80]. The major drawback that caused these accidents is the lack of any prior indications and early warning systems [81].

The vulnerability assessment map (Figure 7) shows areas prone to surface collapse hazards. Since the map has integrated villages, towns, major road networks, and railway stations, it provides vital information for identifying the high-risk zones, taking necessary precautions, and performing mitigation activities. However, the unavailability of administrative boundaries is a limitation. The administrative boundary could provide better population distribution rather than use a circle to represent the spatial extent of the population distribution. Nevertheless, the vulnerability map portrays critical locations that require necessary actions to avoid catastrophes while serving an early-warning assessment.

Vertical land subsidence in JCF investigated in this study, from 2017 to 2020, reveals that almost 10% of the subsidence regions are newly formed, while the dormant percentage is less than 6%. In effect, it denotes that the coal fire-affected area has increased by 4% (1.33% per year) within the study period. An annual coal fire growth of 1.33% amid extensive coal fire controlling actions is alarming in a coalfield burning for more than a century. However, only a long-term study can provide more insights into this matter. Moreover, long-term studies can also offer a better temporal assessment since the sensitivity will improve significantly. On the other hand, the sensitivity can also be reduced by using X-band SAR data, which has better sensitivity to subtle surface deformations. However, X-band-based InSAR measurements will be difficult to obtain in JCF since we have faced critical decorrelation challenges even in C-band InSAR analysis with SAR data pairs having small baselines.

In this study, a wide range of subsidence is identified within JCF due to coal fires. While coal fire, resulting in void spaces, is the reason for surface subsidence, the rate of subsidence is dependent on other factors, including subsurface geology, development of fractures, and depth of coal fire. Therefore, it is always challenging to obtain reliable spatial information about the coal fires. Additionally, the subsurface occurrence of coal fires over a wide area has only intensified the difficulty in monitoring temporal variations of coal fires, especially using field measurements. Despite the multitude of factors influencing surface deformation and challenges coal fire monitoring, the N-SBAS time-series has provided valuable information to execute coal fire controlling missions and risk assessment plans as a part of rehabilitation and remediation strategies.

New coal fire regions defined in this study may not necessarily be new coal fires. Instead, they could be dormant coal fire locations at the beginning of the study period, which got activated during the later stage. Accordingly, there will not be any sign of subsidence at the starting period but shows significant subsidence later on. Similarly, dormant coal fire regions can also be locations where the coal has burnt-out completely. These ambiguities can be resolved only by long-term studies. In addition to the active subsidence that co-occurs with the coal fire, there is a possibility of residual subsidence that can occur after the fire ceases to exist [82]. However, the active subsidence is significantly more prominent than the residual subsidence, and the duration over which residual subsidence can sustain is substantially less in shallow coal fire occurrences. Since most coal fires identified in this study are adjacent to shallow opencast mines in the Barakar Formation, we believe that these coal fires occur in shallow depths. Hence, the occurrence of residual subsidence can be ignored, and the link between the coal fire dynamics and spatio-temporal variations in subsidence is well established in JCF.

Discrepancies between subsidence regions and thermally anomalous regions are expected since they are fundamentally different since they are computed exploiting different characteristics of coal fires. TIR-based studies are less complicated, and they provide vital information about the presence of coal fires. However, they have limitations in detecting the actual subsurface occurrence and extent of coal fires. TIR-based studies often detect thermal anomalies due to hot gases coming out of vents and cracks via convection rather than the heat of coal fire that reaches the surface via conduction. Thereby, the thermally anomalous regions identified from TIR images depend on cracks and vents on the surface more than the actual coal fire extent. Additionally, there is a high chance for false alarms in TIR-based studies because of wind transporting hot gases to a different location, burning other materials, hot rocks, and industrial activities.

On the other hand, N-SBAS based subsidence mapping provides a weather/season resilient and better coal fire mapping in a time-series observation. Subsidence mapping is also vulnerable to false alarms like subsidence due to UG mining, but those can be isolated easily using ancillary information. The major disadvantage of N-SBAS processing is the intense and complex procedures requiring high computational resources. Additionally, the finer spatial resolution of SAR datasets enables the computation of significantly higher resolution coal fire maps compared to those based on TIR sensors. The accuracy and the detail by which coal fires can be monitored using InSAR techniques will be crucial in producing actionable information that would lay the foundation of early warning systems and hazard-management activities by the concerned authorities. Thus, proper monitoring of coal fires would have substantial economic and environmental impacts and prevent the loss of human lives.

6. Conclusions

A total of 9.5 km² coal fire area from different parts of JCF has been identified using N-SBAS technique. While proving the effectiveness of N-SBAS technique for coal fire mapping, it also correlates well with the 9.07 km² of coal fire regions derived from Landsat-8 TIR image using the thresholding technique. The GPS coordinates of coal fires located in the field also exhibit a good correlation with the subsidence attributed to coal fires. In both the maps, Kusunda colliery exhibits maximum coal fire activity in terms of the fire intensity and aerial coverage, producing maximum subsidence of 22 cm between 2017–2020. Even though severe decorrelation challenges arise in computing the surface deformation, the integration using dense SAR datasets has yielded a reliable surface deformation map of JCF. Vertical surface deformation time-series has provided vital information about temporal changes in the coal fire regions. There is a net 4% increase in the total coal fire in JCF, which results from the development of new coal fires (~10%) and existing coal fires becoming dormant (~6%) by 2020. In contrast, 73% of the coal fire regions remained unchanged in terms of coal fire extent and magnitude during the study period. While most transportation networks are not currently vulnerable to surface collapse due to coal fires, a few of the towns (e.g., Tisra, Chhatatanr, and Sijua) and railway stations (e.g., Bansjora and Jharia) are highly prone and needs immediate attention. The unawareness of the coal fire locations and subsidence-prone regions are the primary challenges of early warnings about surface collapse hazards in JCF. To conclude, this study has produced insights into the dynamics of coal fires and associated land–surface subsidence in JCF, which are vital for developing early-warning systems and future planning of fire-control actions.