Post Mining Ground Deformations Transition Related to Coal Mines Closure in the Campine Coal Basin, Belgium, Evidenced by Three Decades of MT-InSAR Data

Abstract

Spatio-temporal ground-movement measurements and mappings have been carried out in the Campine coalfield in Belgian Limburg since the closure of the mines to document post-mining effects. MT-InSAR measurements are compared to groundwater head changes in the overburden and to height data from the closest GNSS stations. Radar interferometry is used to estimate the extension and the velocity of ground movements. In particular, the MT-InSAR technique has been applied to SAR acquisitions of the satellites ERS-1/2 (1991–2005), ENVISAT (2003–2010), COSMO-SkyMed (2011–2014), and Sentinel-1A (2014–2022). The images were processed and used to highlight a switch from subsidence to uplift conditions in the western part of the coal basin, while the eastern part had already been affected by a rebound since the beginning of the ERS-1/2 acquisitions. Following the closure of the last active colliery of Zolder in 1992 and the subsequent cease of mine-water pumping, a recharge of mine-water aquifers occurred in the western part of the basin. This process provoked the change from subsidence to uplift conditions that was recorded during the ENVISAT period. In the center of the coal-mining area, measured uplift velocities reached a maximum of 18 mm/year during the ENVISAT period, while they subsided at −12 mm/year during the ERS-1/2 period. Mean velocities in the western and eastern parts of the coalfield area have decreased since the last MT-InSAR measurements were performed using Sentinel-1A, while the Zolder coal mine continues to rise at a faster-than-average rate of a maximum of 16 mm/year. The eastern part of the coalfield is still uplifting, while its rate has been reduced from 18 mm/year (ERS-1/2) to 9 mm/year (Sentinel-1A) since the beginning of the radar–satellite observations. Time-series data from the two GNSS stations present in the study area were used for a local comparison with the evolution of ground movements observed by MT-InSAR. Two leveling campaigns (2000, 2013) were also used to make comparisons with the MT-InSAR data. The station’s measurements and the leveling data were in line with the MT-InSAR data. Overall, major ground movements are obviously limited to an extension of the actual underground-mining works and rapidly diminish outside of them.

1. Introduction

In Belgium, the majority of natural, regional geohazards that can provoke ground movement are of the following types: earthquakes, landslides, karst, and unconsolidated or poorly consolidated sediments. On the other hand, human-induced geohazards are those related to the overexploitation of underground resources (groundwater and gas exploitation) and those associated with underground-mining or quarrying exploitation (coal, metallic minerals, natural stones, ore deposits, etc.). The study area of this research project comprises the Limburg coalfield area in Belgium and its surroundings. The aim is to create a spatio-temporal overview of the extent of ground movement due to the coal-mining and post-mining history of the area. All coal-mining activities in this coalfield ended a few decades ago, with the last colliery of Zolder closing in 1992 [1]. Mining exploitation was based on the longwall method, whereby subsidence started immediately after exploitation. This was due to the reallocation of rock masses above mining works that resulted to a new equilibrium after the excavation of coal from the longwall panel. This provoked a gradual propagation of deflections and fractures towards the ground surface [2]. A dewatering of the aquifers of coal measures and overlying strata during the exploitation period was the second active process of subsidence. The groundwater extracted by pumping inside underground facilities produced an increase in the effective stress of its rock mass that resulted in a compaction of dewatered layers and, therefore, subsidence on the ground surface [3,4,5]. This movement is still continuing at a lower rate a few decades after the end of mining exploitation. After the mines closed and pumping stopped, the excavated areas and the surrounding rock reservoir were flooded again. The subsidence movement was slowly reversed to an uplift caused by groundwater rebounding into the dewatered aquifers. The groundwater level in the mining aquifers is slowly rising to ultimately recover its hydrostatic equilibrium [6]. A flooding of the mine workings resulted in an increase in their hydrostatic pressure, followed by a reduction in their effective stress, which induced the uplift phenomenon. Since the accessible areas of the coalmines are reduced compared to the total original areas where coal and steriles were extracted due to roof collapse and compaction, the total recorded uplift will only be a fraction of the subsidence caused by coal extraction. Over the years, several methods have been developed to estimate surface movement and to predict its trends. Modeling prediction methods are based on either numerical simulations [7,8,9] or empirical approaches [10]. The first one requires a complicated approach, but, on the other hand, it can model surface movement based on physical laws in which mining and the geological factors are taken into account. The empirical method is easier to handle, but it is based on site-specific measurements. Radar interferometry results can be a source of information to feed this method. Nevertheless, this approach has weaknesses since it is not able to consider complicated geological conditions.

The uplift of former mining areas is a common feature that has also been observed and studied in Great Britain [11], South Limburg in both the Netherlands and Germany [12,13], the Ruhr district in Germany [14], France [15,16], the Czech Republic [17], and both coal [18,19,20,21,22] and copper [23,24] mines in Poland. In the majority of these studies, either only uplift or only subsidence was analyzed, and the analyses usually took place within shorter time periods than the one this research uses. In this paper, the transition from subsidence conditions to uplift conditions as a direct post-mining effect is presented. This required monitoring the phenomenon via a long and almost continuous measurement of time series of displacements. These time series were provided by multi-temporal radar interferometry. Ground-proof data, such as that of leveling and Global Navigation Satellite Systems (GNSS), provide an undeniable source of information for confronting the multi-temporal radar-interferometry data.

The monitoring of the subsidence or uplift of the Earth’s surface related to mining activities and their impact on surrounding built facilities is a challenging task. Techniques such as GNSS, LiDAR, and conventional leveling survey methods for measuring displacement are time consuming and, therefore, expensive for covering large areas with millimeter-to-centimeter precision. For instance, the selection of a reference point for leveling purposes that is not under the influence of ground movement may result in a demanding, iterative process. In addition, a major drawback is a lack of spatial coverage. Synthetic Aperture Radar Interferometry (InSAR), as a remote sensing technique [25,26], has evolved step by step over the last 30 years into a well-established tool possessing the capabilities required for the measurement and monitoring of ground-surface movement. InSAR utilizes radar images acquired by satellite sensors which are independent of daytime and weather conditions but have high spatial resolution and cover wide areas (250 Km for Sentinel-1A). InSAR uses phase information from the images of different SAR (Synthetic Aperture Radar) satellites to measure the land-surface displacement between successive satellite passages. However, conventional InSAR is usually affected by temporal, geometric, and atmospheric decorrelation that limits the uses of the technique [27]. The temporal effect is related to physical changes in the surface during the time period between the two acquisitions. Changes in the reflectivity of the surface due to vegetation, the melting of snow, soil humidity, or other phenomena decorrelate the measurements significantly. Spatial decorrelation is due to differences in the incidence angles between the two acquisitions. All the backscattered signals from a resolution cell on the ground are added up slightly differently, and the measurements are not reproduced exactly. Typical sources of decorrelation are vegetation (especially for X-band and C-band); newly built areas (houses, roads, etc.); erosion processes that occur because the land surface changes all the time; and rapid movements, such as landslides and earthquakes, that rapidly change the Earth’s surface due to heavy modifications and/or the destruction of the area. To overcome these limitations, several advanced techniques known collectively as Multi-Temporal InSAR (MT-InSAR) were developed in the late 1990’s, including the Persistent Scatterer Interferometry (PSI) [28,29,30] and Small Baseline Subset (SBAS) techniques [31,32]. More recently, new algorithms based on Distributed Scatterer (DS) have been developed; for example, SqueeSAR was proposed by [33] to overcome the limitations of the PSI technology and to avoid the processing of large quantities of interferograms needed for the SBAS approach.

In this study, the PSI algorithm was used to process C-band and X-band SAR images covering the period between 1992 and 2022 over the coalfield of Belgian Limburg with the aim of creating a complete spatio-temporal overview of the extent of deformation due to the coal-mining history of the area. The present study shows that the coal basin has been affected by long-term subsidence related to the mining activities of the past but has shifted towards an uplift more recently after the closure of the last collieries. A gradual rise in groundwater levels in the former workings is responsible for the observed uplift.

The Resumption of Coal Mining in the Campine Coal Basin

A drilling campaign conducted by André Dumont, Jr., a professor at the Catholic University of Louvain, led to the discovery of a concealed coal deposit in Limburg. The first coal seam was hit in August 1901 in As at a depth of 541 meters. In 1917, coal mining started at the Winterslag colliery. In total, seven collieries were active in the Campine Basin. To be profitable in spite of technical difficulties (great depth, unstable strata, aquifers, etc.), the mine sites were more extensive than those of the Walloon Basin in the South of the country, such as in the Liège area. The first mine to be closed was the Zwartberg Mine in 1966, after the discovery of natural gas in the North Sea Basin. In 1992, the last colliery, Zolder, closed. In total, these coal mines produced 440 million tons of coal, with the highest production per year of about 8–10 million tons. Peak production was in 1956 with 10.468 million tons; at that time, 44,000 employees were working in the mines. Coal extraction occurred at depths between 400 m and 1100 m. Longwall mining was the preferred extraction method once the mines became mechanized. Usually, the standard geometry of a coal panel was a rectangular shape of ±800*200 m. At the beginning of exploitation, the extraction panels were much smaller.

Geographical and Geological Setting

The Region of Interest (ROI) is situated in the northeastern part of Belgium in the province of Belgian Limburg and includes part of the South Limburg province in the Netherlands. The most important cities in the ROI are Hasselt, Beringen, and Genk (Figure 1A,B). The Albert Canal and the Demer and Meuse rivers, the latter forming the border with the Netherlands, are the major waterways. A major part of the ROI is in the Scheldt Basin, which flows westward in the direction of the North Sea through Antwerp, while the east part is in the Meuse Basin, flowing northward. The morphology of the study area is largely determined by the Campine Plateau. This partly fault-bound plateau originated during the early Quaternary period as a periglacial alluvial fan of the Meuse and later came into positive relief due to differential erosion and incision of the Meuse [34]. The topography is globally flat with an average height of 70 m asl on the Campine Plateau. The altitude values drop from 105 m in the Southeast to less than 40 m asl in the west and in the Meuse Valley (Figure 1A,B).

The ROI belongs to the tectono-sedimentary Campine Basin (Figure 2), which is located north of the London–Brabant Massif, dropping into the Roer Valley Graben (or Roermond Graben) in the Northeast. The basement consists of Lower Paleozoic rocks (Cambrian to Silurian), as an extension of the London–Brabant Massif. The sedimentary succession of the Campine Basin was formed during three main subsidence phases, which are detailed in Doornenbal and Stevenson (2010) and the references therein [35]. The sedimentary history of the Campine Basin started during the Middle Devonian period, after the final Acadian phase of the Caledonian orogeny. This basin was formed as a back-arc extensional basin (the onset of subsidence) and the area in extension was constantly shifting north under the influence of the northward progression of the Rhenohercynian deformation front [36].

During this first subsidence phase, mixed clastic–carbonate Devonian (Givetian to Famennian) rocks were succeeded by Lower Carboniferous (Tournaisian to Visean) dolostone and limestone, which differ locally in thickness due to active-block fault tectonics. During the Variscan orogenesis at the end of the Westphalian period, the deposits mainly underwent block faulting and were slightly tilted [36,37,38]. The coal seams that were mined in the ROI belong to this group.

The second subsidence phase encompassed the Permian to Lower Jurassic periods and started with the erosion of Variscan massifs during the Permian period, resulting in an only locally preserved conglomerate assigned to the Rotliegend area and succeeded by carbonate Zechstein deposits. The sedimentary environment changed from continental to marine between the Triassic and Jurassic periods. Sedimentation was interrupted after the Lower Jurassic period during a new tectonic phase associated with the Cimmerian orogeny, which caused a general tilting towards the north to northeast in direction of the center of the graben. Sedimentation resumed during the Upper Cretaceous (Santonian) to Campanian periods’ sea-level highstand [37] and resulted in a thick package of chalk, which evolved into calcarenite in the shallowing sea during the Maastrichtian and Lower Paleocene (Danian) periods [39,40]. Inversion movements during the Upper Cretaceous period led to a reduction in thickness in the direction of the Roer Valley Graben, succeeded by a relaxation event during the Danian period (lowest Paleogene). Danian carbonates terminated the largely carbonate deposits of the Chalk Group [41]. Eustatic sea-level fluctuations during the Paleogene–Neogene periods were responsible for the complex sedimentation of marine and continental sands, silts, clays, and marls [42,43]. During the Quaternary period, the area was permanently above sea level and the sediment deposits of aeolian or alluvial origin subject to erosion.

The Campine Basin is bound to the NE by the Roer Valley Graben. The main faults inside the Campine Basin have an NNW to SSE or NW to SE direction, which is parallel to the Roer Valley Graben boundary faults; their occurrence defines the neotectonic graben shoulder. The current configuration of the Roer Valley Graben originated during the Upper Oligocene (Chattian) period due to rapid subsidence, which continues into modern times because it is a seat of active intraplate tectonics [44,45]. Structurally, the Campine coalfield forms part of the graben shoulder, just as the adjoining South Limburg coalfield in the Netherlands forms the other part.

Moreover, the Campine coalfield is split by the striking NNE–SSW Donderslag transpressional fault into two sub-basins: the western Campine coalfield, encompassing the Beringen, Zolder, and Houthalen collieries; and the eastern Campine coalfield, encompassing the Winterslag, Zwartberg, Waterschei, and Eisden collieries. The eastern and western coalfields show different sedimentological and burial histories [1,38,46,47], leading to different mining histories as well, with significant repercussions on ground movement.

2. Materials and Methods

To calculate ground movement using PSI, more than 170 C-Band SAR images from ascending and descending tracks of different satellites were used, providing good coverage of both time and space in the ROI. Additionally, 29 X-Band COSMO-SkyMed SAR images, which were acquired in Stripmap (SM) mode (3*3 m), were processed (

Table 1. Characteristics of the SAR-image datasets used in the ROI.

Satellite	Track	Pass	Number of Scenes	Acquisition Period	Reference Point (° LON/LAT)	Master	Processing Software	Avg PS Density (PS/km2)
ERS-1/2	380	Descending	69	1992–2005	4.8583/50.4709	6 April 1997	ROI_PAC, Doris, StaMPS	56
ENVISAT	380	Descending	67	2003–2010	4.8583/50.4709	12 August 2007	ROI_PAC, Doris, StaMPS	61
COSMO-SkyMed	/	Descending	29	2011–2014	5.1183/51.1806	30 May 2012	Doris, StaMPS	164
Sentinel-1A	Path 88 track 163 Swath 1 and 2	Ascending	154	2016–2022	4.3575/50.7981	9 October 2019	ISCE StaMPS	31

). Figure 3 corresponds to the star graph of the baseline vs. acquisition times, with the master in the center of the plot. The available COSMO-SKYMED images were not acquired systematically though the ROI, which caused gaps and non-identical time differences. A considerable gap exists (244 days between September 2011 and May 2012) and the time between recurrent acquisitions varies between 3 and 32 days. These two factors have an influence on the quality of the results. However, as the focus in the ROI concerns a relatively strong ground displacement observed above former coal mines, the time-series behavior is globally less impacted. The average acquisition-sampling time is 35 days for ERS-1/2 and ENVISAT and 12 days for Sentinel-1A. A list of the processed data is given in the Table 1. The ERS-1/2, ENVISAT, and COSMO-SkyMed data were processed using the Doris-StaMPS suite [30,48], while the recent Sentinel-1A data were processed using the InSAR Scientific Computing Environment (ISCE) [49] and StaMPS [30]. For ERS-1/2 and ENVISAT, the data in our possession were in L0-RAW format, which needed to be focused before having the SLC format. This process was realized using ROI-PAC [50]. The COSMO-SkyMed and Sentinel-1A data were already in the SLC format. The interferometric processing of the ERS-1/2, ENVISAT, and COSMO-SkyMed data was accomplished using the Doris software [48]. During processing, the external Digital Elevation Model (DEM) data from the Shuttle Radar Topography Mission (SRTM 3-arc second), which had a 90 m horizontal resolution [51], were used to remove the topographic component of the interferometric phase. StaMPS was used afterward to finalize the PSI processing. The selected master images for ERS-1/2 descending, ENVISAT descending, COSMO-SkyMed descending, and Sentinel-1A ascending were acquired on 6 April 1997, 12 August 2007, 30 May 2012, and 27 July 2017, respectively. After processing, only Persistent Scatterer (PS) points characterized by a temporal coherence greater than 0.7, a common value in radar interferometry studies [52], were selected The study area benefits from many available images; however, there is a one-year gap between the last COSMO-SKYMED acquisitions (with partial coverage) and the first Sentinel-1A image. Furthermore, between the periods spanned by ENVISAT and Sentinel-1A, there are six years of data gaps for those areas which are not covered by the COSMO-SKYMED frames. The first two years (2014–2016) of available Sentinel-1A data were not used in this processing for the benefit of a longer period of acquisition and a larger ROI (processing was realized in the framework of covering the entire country of Belgium).

3. Results

During the ERS-1/2 period (1992–2005), the average annual LOS velocity map (Figure 4) highlights that the entire mined coalfield, contained within the area of the coal-mining concessions, is characterized by two opposite ground movements visible at the surface. The western part (covering the mine concessions of Beringen, Zolder, and Houthalen) shows negative LOS velocities ranging from −12 mm/year to −0.5 mm/year while the eastern part is characterized by positive LOS velocities up to 18.7 mm/year (in the coal concessions of Zwartberg, Winterslag, Waterschei, and Eisden). In this study, a range of LOS velocities between −0.5 mm/year and 0.5 mm/year are considered as stable. Negative LOS velocities with values under −0.5 mm/year concern subsidence, while LOS values above +0.5 mm/year concern uplift. The spatial extension of the highest LOS velocities in absolute value is mainly limited to the borders of the defined coal concessions. It must be noticed that positive velocities are also observed in the KS-Reserve C and Reserve A areas and also along the former coal mines in the Netherlands. The total area in Belgium affected by negative and positive values has a more-or-less elliptic shape with its main axis oriented east–west. The surface of the affected zone is tantamount to ~380 km². The large zone including the KS-Reserve C and Neeroeteren areas, located in the eastern part of the coal basin, barely contains any PS. These areas correspond to the Hoge Kempen National Park, which lacks objects capable of reflecting the radar signal. However, there was limited (KS-Reserve C) or no (Neeroeteren) mining activity underneath these areas. The observed uplift phenomenon continues in the Netherlands where the same coal deposits as in Belgium were exploited across the national border [10,12,53,54].

Over the ENVISAT period (2003–2010) (Figure 5), the LOS velocities of PS points in both the eastern and the western parts are positive and range from 0.5 mm/year to 10.8 mm/year. The spatial extension of the highest positive values is constrained to the coal-mining limits, just as was the case for ERS-1/2. A clear LOS-velocity gradient is observable between the center and periphery of the mining area. The highest positive velocity values are found in the center of coal concessions, while values decrease toward the borders of the concessions as a result of differences in the intensity and depth of the mining, with multiple overlying coal seams exploited in the concession centers.

For the COSMO-SkyMed period (2011–2014) (Figure 6), the map of LOS-velocity values of PS points does not differ much from that of ENVISAT, except for the fact that the SAR track does not completely cover the eastern part of the coal basin. The LOS-velocity values inside coal concessions are positive, but the maximum rate of 7.9 mm/year is lower than that of ENVISAT. The highest positive velocities are observed in the centers of the coal concessions of Waterschei, Houthalen, Zolder, and Beringen. The areas to the south, in the concession of Beringen, or to the southwest, in the Zolder concession, are characterized by negative PS velocity values ranging from −2.6 to −5.8 mm/year, are located south of the Beringen fault (marking the southern boundary of the graben shoulder), and were never mined.

For the Sentinel-1A period (2016–2022) (Figure 7), positive LOS-velocity values of PS points still characterize the entire coal basin as well as most of the observed zone. The average LOS-velocity values of the PS points located inside coal concessions is equal to 3.7 mm/year. The highest positive LOS-velocity values are observed in the coal concessions of Zolder and Houthalen in the western part of the basin. To the north of Reserve-A, there is a large zone characterized by positive LOS-velocity values ranging from 0.5 to 3.4 mm/year. This zone was previously affected by negative LOS velocities during the COSMO-SkyMed interval. Most of Reserve-A is now showing positive LOS velocities, as is the enclosed coalfield.

4. Discussion

4.1. Effects of Groundwater Recharge

The above results demonstrate that, during the studied time span, the abandoned Belgian Limburg coalfield was mainly affected by centimetric positive velocities in the entire coal basin with one exception. The western part, during the ERS-1/2 observation period (1992–2005), was characterized by negative LOS-velocity values. Throughout that time, the highest LOS-velocity values inside the coalfield were slowly decreasing. A comparison of two time series, one located in the vicinity of the Zolder coal concession (5.33068°E 51.02713°N, WGS84) and the second in the area of the Winterslag coal concession (5.53475°E 50.97011°N, WGS84), allows us to better understand the displacement evolution over three decades in the western and the eastern zones (Figure 8A,B). The figure clearly shows that the rate of uplift is reducing with time. This is particularly visible in the slope reduction of the COSMO-SkyMed (2011–2014) and Sentinel-1A (2016–2022) time series. Differences in the incidence angles and geometries of the different satellites has an impact that could reach up to 30% on the actual vertical displacement in comparison with the LOS-velocity values. Therefore, the rate reduction observed during the Sentinel-1A period may be slightly lower than observed in the plots. In the coal district of Liège (east of Belgium), positive velocities are also observed and the rate was on average 2 mm/year [55]. This has been occurring throughout the past 20 years since the closure of the last mine. Therefore, through a comparison of the two areas, we could expect the ongoing uplift conditions to continue for the next two decades by taking into consideration comparable geological contexts while several features (e.g., the number of the coal veins exploited, the thickness of the coal seams, methods of exploitation, the depth of the exploitation panels, structural tectonics, etc.) differ both spatially and between the two areas. However, since the overburden is thicker in Limburg, there will be a slower recharge rate and, therefore, an even longer time for the uplift velocities to slowly decrease.

The subsidence and uplift phenomena result from two different processes, both linked to coal mining. The first one consists of the underground excavation of longwall mine panels, followed by a controlled collapse of the roof of the mined-out coal panel (‘goaf’) and, to a minor extent, a volume reduction of the galleries and driveways (the latter create high permeability channels in a naturally low permeable environment resulting in the homogenization of groundwater responses for each coal mine). This process created additional induced fractures and increased vertical-fracture connectivity above and below the collapsed panels. It should be noted that older mining operations included the backfill of goafs, impeding subsidence and fracture growth, but this practice encompassed only a small portion of the mined-out volume. The second step reverses the infiltration process through the installation of groundwater pumps inside mine galleries, keeping the mines in operation dry. After the closure of the mines and the removal of the pumps, groundwater was able to rise and flood the mine workings. However, discharge data from underground workings do not provide reliable data on inflow rates from the surrounding ‘wet’ strata. This is because the pumps were set on centralized gathering points and most of the pumped-up water was fresh water that was previously pumped down [56]. During the entire period of exploitation (1917 to 1992) and until the closure of Zolder, the ‘mine-water’ aquifers, which consist of coal measures and sandy basal layers from the overlying Cretaceous aquifers, were drained. The subsidence phenomenon appearing shortly (weeks rather than months) after the mining of a coal panel is largely related to the collapse of mined-out coal panels, whereas the remaining subsidence (in terms of years rather than months) may be due to the rearrangement and compression of grain fabric after the removal of pore water. After the closure of all the mines (Zolder in 1992 was the last), a recharging and rising of water tables in the mine aquifers caused the mechanism responsible for uplift, a process gradually slowing down with diminishing differences between the realized hydrostatic rebound and the original hydrostatic pressure from prior to mining. The fact that subsidence initially occurred in the western coalfield after the end of exploitation suggests that the underground mine workings were not yet flooded. To recharge the fractured mining aquifers and overlying Cretaceous aquifers sufficiently to create enough pore water pressure is a gradual process; it takes several years to completely flood the 300 to 500 m high underground infrastructure network. Additionally, as long as the void spaces in coal measures are not completely saturated, there will be no hydrostatic rebound [57]. Thus, the switch from negative to positive LOS velocities occurs with a certain delay following the closure of mines and the stopping of pumps. This switch marks the complete flooding of underground workings. The recharging of coal measure aquifers is mostly a result of seepage through their Cretaceous chalk and into the voids created by coal seams and roof collapses. This delay is because the coal measures, which are predominantly composed of claystone, are nearly impervious in their natural, undeformed state. Additionally, the sandstone beds contained within this claystone environment and the basal sandy layers of the overlying Cretaceous aquifers (which form a single aquifer together with the coal measures) were likely completely drained during the exploitation. A hydrogeological model should be proposed in the future to assess, with more detail, the relationship between subsidence and uplift as they relate to groundwater extraction in the Limburg coal-mining area. However, similar studies [10,12] on Dutch Limburg have shown that uplift recovers 5 to 10% of the vertical displacement that occurred during the subsidence period, which is in line with the degree of compaction and the remaining residual volume in the underground workings. Figure 9 presents one virtual borehole located in Houthalen (5.44674°E 51.01202°N, WGS84) and a second one located in Winterslag (5.5334531°E 50.9706586°N, WGS84), showing the different geological layers.

PS results from the ERS-1/2 images (1992–2005) show that the western part of the Campine Coal Basin is continuing to subside while the eastern part is already uplifting. The main reason for this resides in the fact that coal mining and the associated groundwater pumping for dewatering underground workings stopped earlier in the eastern part than the western one (1988 vs. 1992, respectively, Figure 10). In a situation with connected coal mines, the mine water would have flowed from east to west as a result of gradient differences between the two coalfields. However, because of the presence of the Donderslag fault system, which alters the structure of the coal basin and vertically displaces individual coal seams by 1000 m (Figure 11), the coalfields and mine-water aquifers are not connected. A second reason for this is that the overburden along the southern edge of the eastern coalfield is thinner and sandier, allowing for a more rapid replenishment of the Lower Cretaceous–coal measure aquifer from above (Figure 9). A third reason results from differences in the diagenesis of the coalfields; the eastern one reaches higher thermal maturity in the coal seams [59,60], and its surrounding rocks are more brittle, allowing for open fractures. The surrounding rocks of the western coalfield, especially of the Beringen mine, are notoriously more plastic as a result of the expansion of claystones and the closure of fractures, leading to less stability for the underground workings and also less inflow from the surrounding aquifer.

Multiple piezometers are present inside the study area, but few measure the Cretaceous aquifer and even fewer of those have long time series. None reach the coal measures. Based on the chemical water composition of seepages into the dewatered underground workings, D’Hooge (1990) [56] postulated that there is no connection between the mine-water aquifer and the deepest freshwater aquifer of the Maastricht Formation (top Cretaceous). However, this cannot be verified anymore because of the closure of the mines and the sealing of the shafts [57]. Therefore, it is almost impossible to develop a quantitative approach based on the recharging of aquifers and on the subsidence/uplift that affect the study area. Taking into account these considerations, two piezometers have been selected from the Flemish online database: PZ 7_0546, installed in the Maastricht Calcarenite Formation of the confined Cretaceous aquifer at a 330 m depth (HCOV 1100); and PZ 7_0350 in the phreatic Miocene Bolderberg sand at a 100 m depth (HCOV 0253). Figure 10 shows their location, and Figure 12A,B show the evolution of the water table in the considered aquifers. The water table in PZ7_0350 is variable, but the greatest changes are limited to a maximum of +1.4 m. Obviously, the observed changes of the water-table levels are not related to the recharging of the mining aquifers. The aquitard of the Boom Formation plays the role of an impermeable horizon. For PZ 7_0546, the measurements did not start until 2007, but the water-table level has been rising slowly since then for a total of 4 m. This value is very small in comparison to the 77 m groundwater rise measured in 1994 in the Dutch coal mines [12] right after the stop of pumping activities. This difference is due to the sandier character and reduced thickness of the Cretaceous aquifer in South Limburg. In the Campine Coal Basin, at least in the western coalfield where the piezometer is located, the upper calcarenite aquifer is separated from the lower sandy reservoir by a thick layer (in the order of 150 m) of rather impervious chalk. A graph of the water levels in the Dutch shafts shows a rapid increase in water-table levels in the Wilhelmina shaft, reaching an asymptote followed by a slower linear increase after the year 2000. The quasi-linear recharge presented by PZ7_0546 (Cretaceous) is, thus, similar to the end of the measurement curve observed in the Netherlands.

4.2. Effects of the Actual Mining Works

PS data from the different time periods (see Section 3, Results) show that ground movements, both positive and negative, are located inside the Limburg Coal Basin as delimited by the coal concessions. These boundaries are important from the perspective of mining history but do not correspond with the actual extension of mined areas, which are more precisely represented by main galleries leading towards extracted coal-seam panels. Figure 13 shows all the galleries seen from the top, which means that there are vertical superpositions in this view that are also of extracted coal seams. This is particularly the case in the coal concessions of Waterschei and Winterslag. In Beringen, it is especially clear that the extension of positive LOS-velocity values fits perfectly with the frame formed by the galleries. This situation is also true for the other concessions but less conspicuous. Northward galleries were driven in advance and were followed by the longwall mining of coal panels within reach of the new gallery. However, some coal deposits remained unmined because the mine closures interfered with the normal exploitation program (Zolder colliery) or because they were driven underneath an industrial estate, which successfully prevented mining (Eisden colliery). Nevertheless, from these observations, it appears that uplift areas are spatially limited by the extension of galleries and, de facto, by exploited coal panels. The coal concession of Neeroeteren was never mined and only the extreme NE area of the Oostham concession was exploited. This situation is well reflected by the absence of ground movement in these two concessions. The uplift observed east of the Eisden colliery is also related to the exploitation of coal across the Belgian border in Dutch Limburg [12,53].

A series of cross sections have been extracted to present the evolution of the velocity values of PS points associated with each SAR dataset along the two paths, A and B, of Figure 13. The paths are drawn to have one cross-section in the western basin and the second in the eastern part. Figure 14 shows that the evolution of the velocity values of PS points averaged an over 200 m buffer across each path for the 4 SAR datasets (Table 1).

Cross-section A for the ERS-1/2 period (1992–2005) shows that negative velocities were limited to extraction zones, which were concentrated towards the northern part of the colliery and for which Reserve A constituted a border. During the ENVISAT period (2003–2010), the same limited zone was strongly affected by positive ground movements reaching nearly 10 mm/year. These velocities were averaged by considering all PS points included in a 200 m buffer along the cross section. This averaging flattens the velocity values given along the cross sections. During the time associated with COSMO-SkyMed (2011–2014), cross-section A first showed negative velocity values in the southern part right after crossing the Reserve A. After that zone, the remaining part of the cross-section was characterized by positive values. This negative zone was only present during the COSMO-SkyMed period and could not have been related to any pumping activities or works in the vicinity. A possible explanation is the presence of cohesive glauconiferous sand and ferruginous sandstone beds in the near-surface Miocene Diest Formation, which could delay the upward movement of subsidence until finally yielding to rupture and subsidence breakthrough up to the surface [61]. The shape of cross-section A for Sentinel-1A was quite flat, which could indicate that the uplift phenomenon was slowly reducing its speed.

During the ERS-1/2 period, for cross section B located in the eastern basin, PS velocities became strongly positive when entering the mined area. In the southern part of the curve for ERS-1/2 and ENVISAT, the velocities showed a steep increase when entering the Winterslag colliery. After Winterslag, the velocities slowly decreased toward the north over the Zwartberg colliery, where exploitation had already stopped in 1966 but which continued to be dewatered so as to not hamper mining in the adjoining (and partly connected) Winterslag colliery. The COSMO-SkyMed cross section was noisier, but clearly demonstrated an increase in velocities within the mined area. The Sentinel-1A curve for cross-section B was very flat and could confirm the attenuation of the uplift phenomenon as a result of the end of coal-mining activities. However, since the length of the time series in Sentinel-1A is shorter in comparison to ERS-1/2 or ENVISAT, it is possible that the trend is not well defined.

To analyze the impact on built environments, fieldwork was conducted for a project named the Geotechnical and Patrimonial Archives Toolbox for Architectural Conservation in Belgium (GEPATAR) [62,63]. It was highlighted that several buildings had cracks that were related to the subsidence phenomenon. On the other hand, it was not possible to designate any fracture as being only related to uplift.

In conclusion, the measured and mapped impacts on ground stability and identified ground motions practically coincide with the extent of the effective extraction area, and are, thus, clearly related to past mining activities.

4.3. GNSS vs. PSI Time Series

In this paper, PSI results are compared to the ellipsoidal height time series of two GNSS antennas located in the ROI. The GNSS time series are from reference antennas of the FLEPOS RTK-GNSS network. Specifically, the time series are from the antennas Maasmechelen (MAME) and Houthalen (HOUT) and cover the periods of 2004–2014 (GPS week 1279–1797) and 2009–2022 (GPS week 1217–2222), respectively (Figure 15A,B). The coordinates of the antennas were determined in post-processing using the Bernese v5.2 software [64]. The daily RINEX files of the stations were processed, together with all other Belgian GNSS reference stations, using precise ephemerides. The daily solutions used in this study were the result of a re-processing of all data from the Belgian GNSS reference stations covering the period of 2002 to 2022, performed beginning in 2022, in order to maintain consistency in the results. The re-processed daily solutions were afterwards combined using sequential adjustment techniques (normal equation stacking) [65] to determine the coordinate solution for each GPS week. The time series of weekly solutions contains less noise compared to the time series of daily solutions, which have a larger standard deviation because of higher influence from external factors. Finally, because the PSI data processing was performed in monthly intervals, the GNSS week solutions were also averaged over monthly intervals to enable comparison.

Figure 15A,B show a comparison of the vertical displacement detected by GNSS and PSI (in LOS) for the MAME and HOUT stations. We note that the vertical displacement of the GNSS time series was visualized by plotting the evolution of its ellipsoidal height, while the vertical displacement of PSI was calculated as the average vertical displacement of all PS points within a buffer of 250 m around the GNSS stations over the time span covered by the GNSS time series. For the MAME station, a total displacement of 52 mm was observed for the GNSS time series between 14/07/2004 and 18/06/2014, resulting in an average rate of 5 mm/year. Note that an antenna change in 2010 resulted in a shift in its absolute height value, but its displacement velocity remained consistent. For the station of Houthalen, a total displacement of 121 mm was observed for the GNSS time series between 07/05/2003 and 10/08/2022, resulting in an average rate of 6 mm/year.

The time series in Figure 15A,B indicates a clear agreement between both techniques, which is supported by both auto-correlation and cross-correlation coefficients. For the MAME station, auto-correlations for the GNSS-derived height changes between −26 and +26 weeks are all above 0.82, which is in line with the uplift signal and the absence of a clear seasonal signal in the time-series plot. To compare the GNSS-derived height changes with the PSI data, only the ENVISAT data could be used due to the decommissioning of the MAME GNSS reference station in 2014. The cross correlation with a lag of 0 between the monthly averages of uplift, determined by GNSS and ENVISAT, is 0.92. Using time-lagged cross-correlation coefficients, minor differences (a difference of 0.04 between the minimum and maximum cross-correlation coefficients) can be observed over a window of ±12 months. For the HOUT station, the auto correlations that were determined for the GNSS-derived height changes between −26 and +26 weeks are all above 0.88, confirming the presence of an uplift signal in the time series. Given the auto-correlation values, seasonal signals appear to not be relevant for this work. The cross correlations with lags of 0 between the monthly averages of uplift, determined by GNSS and ENVISAT, SENTINEL-1A, and COSMO-SkyMed are 0.99, 0.96, and 0.95 respectively. The slight differences between the coefficients might be explained by the number of reflecting elements in the target zone of 250 m around the GNSS station. For SENTINEL-1A, only one reflection point was found, while the movements of ENVISAT and COSMO-SkyMed resulted from an average of over five and three reflection points, respectively. When determining cross correlations with time lags, maximal correlation coefficients were found at −8, −7 and +11 months for ENVISAT, SENTINEL-1A and COSMO-SkyMed, respectively (with respective cross-correlation values exhibiting differences of 0.01, 0.03 and 0.07 between their minimum and maximum cross-correlation coefficients). The data density and sampling, however, do not allow for further research on whether or not the environment surrounding the GNSS antennas is effectively deforming in a similar manner, but with a lag between movements. In Figure 16, four other GNSS stations surrounding the coal-mining area are also presented: DIES (Diest), MEEU (Meeuwen), BREE (Bree), and MAAS (Maaseik). The numbers on the map that are associated with the locations of each GNSS antenna represent the vertical velocities of the GNSS antennas in a mm/year format, calculated for the period of 2004–2010. As expected, the GNSS antenna located outside of the mining activity is characterized by slower, but still positive, velocities than the one located near the coal mines. This is another insight showing that ground movements are restricted to occurring a short distance from the coal mines. Nevertheless, it is also worth mentioning that the surroundings of coal concessions are also affected by persistent movement after the end of exploitation. Unfortunately, there is no GNSS station installed in the study area during the acquisition time of ERS-1/2, which could have shown the subsidence pattern of the western part of the basin. In conclusion, the time series of displacement monitored by the GNSS antennas and the PS data are in good agreement and show the same trends.

4.4. Leveling vs. PSI Time Series

In Belgium, the National Geographic Institute has conducted several leveling campaigns in the area since its origin. Among them, two surveys were performed during the period covered by the PSI data. The first survey used in this study was conducted in 2000 in support of a leveling network adjustment and in preparation for the calculation of the Belgian geoid model hBG03. The second survey was performed in 2013, as the leveling network required an update in the mining zone as a result of ground deformations. Note that the height values of both datasets are expressed in the same frame (2003) of the Belgian orthometric second general leveling height system to avoid a comparison of height values that are influenced by different network adjustments. About 250 locations were measured in the vicinity of the coal concessions. A comparison with PSI data was performed by calculating the height differences for each leveling point between the measurements made in 2013 versus the ones made in 2000. The height differences were then converted into velocities by dividing them by the time between the two leveling campaigns. The velocities of the PSI data for ENVISAT were used to produce the maps because the ENVISAT time period is in line with the period covered by the leveling data (Figure 17A,B). The maps show clearly that the leveling measurements are in line with the PSI data. The maximum velocities of the uplift were well recorded in the center of the basin, while they were reduced in the zones that were not mined. One must pay attention to slight differences that could be attributed to the reference points used for the leveling that were intended for lines still under the influence of the mining-ground movements mapped by PSI. Therefore, the base reference could have been shifted. We also noticed that the velocities that were determined using the differences between the two leveling campaigns are the result of a simple averaging over 13 years, a technique that cannot grasp the varying velocities of ground deformation as well as the ones that can be observed using GNSS or MT-InSAR. Finally, because a leveling campaign takes multiple weeks to perform, ground deformation impacting the ROI during the time of the measurement campaigns might have influenced the final accuracy as well.

5. Conclusions

The processing and analyses of SAR data since 1992 from four successive satellites has demonstrated that ground movement occurs and is still being measured (up to 2022) in the coal-mining basin of Belgian Limburg. Regarding the western part of the coal basin in the transitional period between the ERS-1/2 (1992–2005) and ENVISAT (2003–2010) acquisitions, the area first subsided as a result of remaining mechanical compaction processes associated with underground coal mining and uplifted afterwards. This was due to the recharging of mine-water aquifers after the active dewatering of underground workings ceased, which resulted in a hydrostatic rebound. On the other hand, the eastern basin had been uplifting the entire time because mining and pumping activities ended earlier there, just before the recording time. The measured rates are the fastest ground-movement velocities recorded in Belgium, and they were measured at an average of 20 mm/year (positive or negative). An amplitude reduction of the recorded velocities of PS points inside the uplift area during the Sentinel-1A period tended to suggest that the phenomenon is becoming less and less intense. Future observations should confirm or deny this trend. A mapping of the galleries and driveways into extracted coal-seam panels has shown that ground movements are almost entirely limited to the extent of underground mine workings. This is also attested in Vervoort and Declercq (2018) [66]. Thanks to the presence of GNSS antennas in the study area, it was possible to compare and correlate the recorded ground movements using two different geodetic approaches. The time series of displacement of the GNSS stations and PS data show a high correlation, which confirms the quality of SAR processing and measurements. The leveling campaigns performed in 2000 and 2013 also confirm the observations that were made using the PSI data. The absence of a sufficient number of deep piezometers for assessing the variation of water tables in the Cretaceous and coal measure aquifers did not allow for a robust analysis of the recharging of mining aquifers. However, in the future, the development of a basic hydrogeological modeling based on the little data available should be envisaged to model and forecast, in both space and time, the impact of rebound on ground movements observed at the surface. This paper also highlights the importance of the geodetic technique to monitor ground deformations affecting rural and partly urbanized areas above old coal mines and past exploitations. This aspect is of prime importance for future developments and monitoring human infrastructures in such mined areas on a long-term basis.