Combined Study of a Significant Mine Collapse Based on Seismological and Geodetic Data—29 January 2019, Rudna Mine, Poland

Abstract

On 29 January 2019, the collapse of a mine roof resulted in a significant surface deformation and generated a tremor with a magnitude of 4.6 in Rudna Mine, Poland. This study combines the seismological and geodetic monitoring of the event. Data from local and regional seismological networks were used to estimate the mechanism of the source and the ground motion caused by the earthquake. Global Navigation Satellite System data, collected at 10 Hz, and processed as a long-term time-series of daily coordinates solutions and short-term high frequency oscillations, are in good agreement with the seismological outputs, having detected several more tremors. The range and dynamics of the deformed surface area were monitored using satellite radar techniques for slow and fast motion detection. The radar data revealed that a 2-km² area was affected in the six days after the collapse and that there was an increase in the post-event rate of subsidence.

1. Introduction

Rock bursts and collapses are the most dangerous phenomena associated with seismicity in underground mines. The collapse of a mine roof can follow strong and shallow seismic events of magnitudes M > 3 [1,2,3,4]. They pose a serious threat to life for miners working in the vicinity of the collapsing rock. The strongest underground mining-induced seismic events are Volkershausen [5] and Newcastle [6], both M5.6 Although the magnitudes of such shallow seismic events are usually small in comparison with natural earthquakes, the ground shaking might be felt by local citizens and can damage surface infrastructure. Rock-burst issues have occurred in the Polish copper region known as the Legnica Głogów Copper District (LGCD). Thousands of events of M > 1.0 are recorded every year in the LGCD. Among these, several are followed by rock bursts or tunnel collapses. Both ground shaking (i.e., a seismic event) and local ground failure occur as consequences of mining. While the shaking effects are comparable to what can be felt during small to moderate natural earthquakes, the surface deformation is specific to mining areas. The subsidence is either a ground response to the closure of excavated mining panels, e.g., [7,8,9], or is related to strong tremors [10,11]. In this paper, we present a joint seismological and geodetic study that expands our knowledge concerning the connection between mine collapses and the ground response to seismic sources manifested by permanent ground deformations.

The analysis focused on the strongest seismic event recorded by LGCD in the last several years—namely, the 29 January 2019 M4.6 collapse in the Rudna Mine, Poland. The phenomenon was widely felt and was followed by massive collapse. While the seismological analyses were carried out using records from local and regional seismic networks, the geodetic data are comprised of Global Navigation Satellite System (GNSS) observations collected at station LES1 (Komorniki; part of the European Plate Observing System project [EPOS-PL] network) and the European Space Agency’s (ESA’s) Sentinel-1 radar observations. The LES1 station was aligned with a strong motion sensor (KOMR), part of the local seismic system, the Legnica–Głogow Underground Mining INduced Earthquake Observing System (LUMINEOS), located above the mine. Interferometric Synthetic Aperture Radar (InSAR) techniques were applied to analyse the co-seismic surface deformation.

2. The Collapse and Seismological Analyses

On 29 January 2019, Rudna Mine in western Poland was struck by a seismic event followed by a rock burst. It was one of a few such events in the area in recent years [4,12,13]. Using seismic signals recorded by the LUMINEOS surface local seismic network (Figure 1a) and a local velocity model [14], the epicentre was located at 51.5110°N, 16.1197°E. For our analysis, we assumed that the event depth was located at the level of the excavation; that is, 1 km below the surface. LUMINEOS consisted of 10 strong motion sensors (AC-73) and 17 five-second VE-53/BB instruments (both types manufactured by GeoSIG, Schlieren, Switzerland).

The source mechanism of the event was estimated using seismic signals recorded by broad-band seismometers (STS2 or Güralp CMG3ESP; Figure 1b) at regional distances of 70 to 300 km. The signals and corresponding station responses were obtained from the ORFEUS-EIDA project page [15]. The inversion was performed using the real displacement, derived from the original signals, for 11 seismic stations located around the epicentre (Figure 1b). The mechanism (Figure 2a) was estimated in terms of the full moment tensor, which can be directly used to interpret the mechanisms of mining tremors [16]. The Pyrocko toolbox [17] and KIWI tool software [18] were used for the calculations. The software had already been tested as a tool for determining a stable point source model in mines [19], including mine collapses [4,11]. In this work, we followed the procedure previously used in [4]. The inversion was based on the full waveform displacement and spectra in the specific frequency range. The regional velocity model proposed by [20] was used. The lowest misfit was obtained for data filtered between 0.07 and 0.11 Hz. The final moment tensor solution was decomposed according to [21] into its isotropic and deviatoric parts, and ultimately into the double-couple and compensated linear vector dipole. The resulting moment tensor (m11: −4.03 × 10¹⁵ m22: −3.91 × 10¹⁵ m33: −11.9 × 10¹⁵ m12: −0.07 × 10¹⁵ m13: −0.11 × 10¹⁵ m23: −0.22 × 10¹⁵ Nm) was dominated by the implosive component (isotropic part: −55%). This kind of mechanism can be explained by a tabular cavity collapse [4,22], and it was strongly supported by the observations recorded on the strong motion sensors belonging to the LUMINEOS network(see Supporting Material: Table S1). Figure 2 shows the distribution of the peak ground acceleration around the epicentre. The clearly visible peak ground acceleration pattern corresponds very well with the collapse source model.

3. GNSS Short- and Long-Term Monitoring

The deformations of mining-induced areas might be monitored by GNSS sensors. In the literature, some research concerning the ground deformation monitoring can be found (e.g., [23,24,25]). However, for short-term oscillations, the papers usually are related to natural earthquake analysis (e.g., [26,27,28,29]). For the current study, the GNSS observations at station LES1, collected at 10 Hz from all visible GPS, GLONASS, Galileo and Bediou satellites, were processed in two different modes: (1) as a long-term network daily solution, producing high-quality, daily-coordinates time-series for a period of six months (three before and three after the event); and (2) short-term, high-frequency oscillations in the horizontal and vertical coordinates, observed during the event, to produce a high-frequency picture of the seismic wave passage.

3.1. Long-Term Monitoring

The daily solution was performed for station LES1 using Bernese GPS Software v.5.2 [30]. Selection of the processing parameters was based on the standard EUREF Permanent GNSS Network (EPN) processing strategy in [31], and they included the following: reference frame: ITRF 2014; orbits and clocks: CODE Final; ionosphere: CODE GLOBAL IONOSPHERE; troposphere: ZTD a priori and mapping function: Vienna Mapping Function 1; reference stations selected from EPN, independent vectors; coordinates and velocities for reference sites: EPN_A_IGS14.SNX; solution type: daily, minimum constrain on translation; estimated parameters: daily coordinates transformed into east, north, up (ENU), troposphere delays, and gradients; time span: 1 November 2018–30 April 2019; quality of estimation measured by mean Helmert transformation root-mean-square error on the daily boundaries: 2.6 mm.

The solution, based on double differences and an independent network of vectors, showed a high level of consistency between consecutive days, as confirmed by a low Helmert transformation RMS. The daily solutions were not combined into a full-period solution, so as to allow the detection of jumps related to post-seismic deformations; therefore, the impact of the plate velocities on the reference stations was removed using the EPN solution [31]. The constant velocities in each component were estimated. The analysis of six months of coordinate time-series (Figure 3a) showed that, for the north component (top panel), five significant (larger than 2σ) jump events were recorded. For the east component (middle panel), it was five events, and for the up component, it was nine events. For the investigated case (red line in Figure 3a), none of the estimated coordinates exceeded 2σ threshold. However, other events (marked as a grey dashed line in Figure 3a) coincide with significant outliers; further studies should investigate whether there is a connection between deformations and induced seismic events.

3.2. Short-Term Monitoring

The high-rate GNSS data was recorded by station LES1, which was collocated with the KOMR accelerometer, 2.5 km north-east of the epicentre. Two hours of 10-Hz GPS data were processed using GAMIT/GLOBK v.10.9 software [32] and the double-differencing method in kinematic mode. The WROE reference station was located 79 km south-east of the epicentre, outside the deformation zone. For processing, CODE final orbits and differential code biases were used. The phase and code observables were combined to remove ionospheric refraction by the Melbourne–Wübbena combination [33,34]. Additionally, the Universitat Politècnica de Catalunya’s rapid Global Ionosphere Maps (UQRG) was used, which significantly reduced the noise in the resultant position time-series. The UQRG model had a 5.0 × 2.5 spatial resolution and a 15-minute temporal resolution [35].

The geocentric coordinates (XYZ) obtained were transformed into local topocentric (ENU) coordinates and displacements. To compare the GNSS-displacement time-series with the positions calculated from the accelerometric data, both time-series were reduced to the period starting 30 s before the event and finishing 90 s after (Figure 3b). The low-frequency noise of the GPS position time-series was removed using the second-order Butterworth filter and the high-pass frequency set to 0.15 Hz, based on the Fourier spectra of the seismological position time-series [27]. This reduced the standard deviations of the GPS position time-series from several millimetres to less than 1 mm. Pearson’s correlation coefficient reached 0.93, 0.89 and 0.85 for the E, N and U displacement time-series, respectively. This indicates good agreement between both time series.

4. InSAR Data and Analyses

Two InSAR techniques were used to obtain the areal pattern of the surface deformation: (1) a small baseline subset (SBAS) [36] for slow-rate displacements; and (2) a conventional differential InSAR (DInSAR) [37] for fast-motion monitoring around the epicentre. The selected data had the following parameters: orbits: ascending (track 73) and descending (track 22); type of images: Sentinel-1A/B (C-band), IW mode, VV polarization, SLC format; resolution: ~2.3 m in range, ~13.9 m in azimuth; time span: 28 October 2018 to 20 April 2019; revisiting time: six days, with the exception of one ascending 12-day pair (23 January–4 February 2019). In both techniques, SBAS and DInSAR, the differential interferograms were filtered with a Goldstein filter and were unwrapped using the minimum cost flow (MCF) phase unwrapping procedure [38]. The received phase differences present the surface deformation in the line-of-sight (LOS) direction to the satellite in grid images with 20-m resolution. All the interferograms were vertically translated to a common reference point, namely station LES1 (Figure 1a), and were roughly converted to vertical displacements by the division to cosine of the incidence angle. The 1 arc-second (30 m) SRTM-DEM [39] was used to remove the topographic phase component.

4.1. Slow-Rate Motion Monitoring with SBAS

This technique is typically applied for monitoring slow motions and strain accumulation. It faces some limitations regarding the capability to detect surface deformation rates above certain values depending on the wavelength of the radar signal [40]. Additionally, a bigger part of the studied area is covered by vegetation, thus it is very challenging to investigate it by using stacking-based methods (like SBAS), which depend on a long time coverage and are therefore sensitive to temporal decorrelation. Considering that SBAS has an advantage to estimate atmospheric components, and DInSAR is fully capable of detecting co-seismic displacement (see Section 4.2), these two techniques are complementary to each other and can also be evaluated with each other internally since their results are processed separately using two different softwares—SARscape [41] for SBAS and SNAP [42] for DInSAR. Thus, in SBAS processing, a stack of multi-looked interferograms was applied according to the typical SBAS method [43,44], which is optimised for resolution cells containing distributed scatterer (DS), which are more effective in vegetated areas. Additionally, to maximize the coherence and to avoid incorrectness in the series due to a short rapid movement of the collapse, we divided our dataset for the SBAS processing into two independent datasets from descending images, from three months before (15 images) and three months after the event (14 images). The maximum critical baseline percentage was set to 2%, and the maximum temporal baseline was set to 30 days. The low-quality interferograms were removed from further processing. Finally, 36 and 21 interferograms were used for interferometric component estimation for the first (Figure 4a) and second (Figure 4b) datasets, before and after the event, respectively.

It was not possible to capture the velocity in the centre of the subsidence bowl due to the forest land cover in that area (Figure 4c), which presented an obstacle for the C-band signal. Nevertheless, the results showed a slight acceleration of the movements around the main zone of deformation in the second period, and more prominent vertical displacement, in the range of 80 cm/year, in the most affected area after the earthquake.

4.2. Earthquake Surface Deformation Detected by DInSAR

For the DInSAR, the Sentinel-1 data were processed following the standard consecutive approach using the smallest temporal six-day baselines, with the exception of the one 12-day ascending pair mentioned earlier. Differential interferograms of adjacent SAR acquisitions were calculated and accumulated to provide complete time-series interferometric results. This technique uses small temporal baselines which minimize temporal decorrelation [44,45]. Vertical displacements at each pixel acquired by DInSAR estimates are presented in Figure 5.

4.3. InSAR Time Series

Four test points (Figure 4 and Figure 5) were selected in order to compare the SBAS and DInSAR vertical subsidence velocities before and after the seismic event (Figure 6a). Three of the test points (1, 2 and 4) coincided with ground objects with solid construction, such as roads and buildings, lying within highly coherent (≥0.9) pixels. Point 3 was chosen in the area of maximal detected SBAS subsidence. It was captured only in the second SBAS set.

The displacement values from the ascending and descending DInSAR deformation maps were extracted and compared by time-series (Figure 6b). The mean values of the vertical displacement in the time-series graphs were calculated by linear interpolation.

5. Discussion and Conclusions

Mining-induced seismicity and the associated ground deformation, measured by advanced geodetic observations, have not been studied by many authors, except in the cases of controlled salt mine cavity collapse [46,47,48], examples in the Upper Silesia Coal Basin [7,8,11,49,50], and in the LGCD [51,52]. Previous works concerning ground deformation studies in mining areas, however, were focused on geodetic data only (i.e., remote sensing or GNSS) ([46,47,48,51,53] and references therein) or used just very basic seismological parameters, usually epicentre locations provided by various seismological centres [13,50]. Some specific models of ground deformations have also been considered previously [54]. Since a mining collapse is actually an isotropic seismic source, these models seem to be in very good agreement with the seismological model of mining collapse. On the other hand, the impact of mining tremor deals with a very minor number of existing subsidence basins in mining areas. Hence, it is difficult to judge, using the geodetic data only, whether the subsidence is directly connected with a tremor or not. In all these studies, it was postulated that the increased seismicity manifested itself on the ground as a higher subsidence rate, with the measurable effects of post-seismic vertical deformation reaching up to 10 cm [52]. We applied several different techniques in order to gain a better understanding of the characteristics of the strong anthropogenic tremors induced by mining activities. These results not only support similar previous studies, but also contribute new findings.

First of all, we determined that shallow-induced seismic events with magnitudes of M4 and above can manifest not only as a relatively strong ground motion but can also provoke significant surface deformation in a very short time. Both effects have been observed in mining areas [6,11], but rarely in such a detailed way as presented here. The combination of different long- and short-time monitoring techniques provided an opportunity to track the process of evolving surface deformation in a complex regime of mining subsidence and anthropogenic seismicity in detail. Thanks to the GNSS and InSAR analysis, the rate and delineation of the subsidence was determined quite well, and we confirmed the rapid subsidence of the ground in relation to the mine roof collapse events.

High-rate GNSS observations have been applied in earthquake analysis and structural health monitoring. Previous studies of seismogeodesy have been concerned with large natural earthquakes with magnitudes exceeding M4.5 and amplitudes reaching decimetres and even meters (e.g., [27,55,56,57]). The event from 29th January is the first example presented of HR-GNSS application in short-term deformations caused by induced earthquakes. Studies related to induced shocks have involved long-term displacements and standard sampling frequencies of up to 1 Hz only [58]. In this study, the greater potential of the HR-GNSS for small events with centimetre-scale amplitudes was confirmed by a comparison with the seismological data. The short-term GNSS observation revealed the possibility of a wider interpretation of the source mechanism in the future. Data from the collocated GNSS and seismic sensors showed very good agreement in the displacement time-series, which indicates that the GNSS stations can be used for additional, alternative data, recorded around an epicentre in an area with sparse seismic sensors. On the other hand, the long-term GNSS study revealed the possible detection of other high-energy seismic events in the area. However, this is highly dependent on the distance between the GNSS station and the seismic sources and requires extended investigation.

The spatial range of the surface deformation derived from the InSAR studies corresponded well with the estimated location and focal mechanism of the event. The DInSAR results showed an undoubted capability to detect and delineate the scale of the surface deformations caused by such an event, as proved by the comparison between the ascending and descending data, even though the ascending co-seismic interferogram (12-day) had a lower quality. The combination of the SBAS and DInSAR results demonstrated the intensive subsidence with varying velocity associated with mining activities. Continuous subsidence occurred in the area at different velocities, depending on the dynamics of the material being extracted and also the local seismicity.

In principle, the source parameters can be derived directly or indirectly from seismograms recorded around the hypocenter. The most important and interesting in the case studied in this work is the source mechanism, which in fact is quite different in comparison with natural earthquakes that occurred on faults. In such regions, the dominance of the double couple mechanisms is observed, and no isotropic sources are expected. Concerning the mechanism, since the subsidence was completed in a very short time after the quake happened, the satellite observations can be used as an indicator for the source type determination. In our case, the geodetic studies strongly supported the full moment tensor estimated with seismological analysis. We also gained better knowledge about the size of the rupture thanks to the additional study of the seismic and geodetic data. This is quite unique in the case of seismicity related to underground mining. Our results also help in better understanding the influence of the collapse like events on surface vibration measured in terms of PGA/PGV. By geodetic techniques, we investigated the surface for possible deformations before and after the event in order to verify the proper detection of the affected area. Last, but not least, it was previously shown that seismic sources with visible rock bursts (e.g., tunnel collapses) in this area are complex (see, e.g., [4,12]). This is not obvious seismological observation and is crucial to understand the phenomena in the future. Since the subsidence is directly connected with the collapse, using geodetic analysis we would be ready to say more about the source physics.

The outcomes from the current study confirm the potential of the applied techniques for investigating ground vibrations and the surface deformation related to them. The studied event highlighted the need for a more detailed study, focused on several areas, starting with the highly complex nature of the ground deformations in areas where mining activities occur. The combination of factors provoking such ground changes, such as the very fast motions caused by moderate tremors with different magnitudes, the fast and slow nonlinear subsidence of the surface triggered by material extraction, and also the change in subsidence rates between the pre- and post-earthquake periods, the local geological structures, underground water-level changes, etc., produce a complex dynamic regime at the mine surface and in the behaviour of the cavity. The combination of different techniques provides a better understanding of the mining and its induced seismicity influence on the surface.

To increase the impact of such joint investigations, several recommendations can be made. A high density of GNSS stations, in areas with significant seismic risk also provoked by human activities, is recommended. A high density will allow a more detailed, long-term study of the horizontal and vertical trends of the surface in the area. On the other hand, the results from short-term GNSS measurements could be used to densify and verify the seismic recordings. Combined observation of long time series of seismological and geodetic data (GNSS, InSAR) could reinforce the possibility of detection of pre-seismic markers, enabling it to become a method for developing an early warning system in mining areas.

Undoubtedly, the investigation of this event has provided a new way for analysing hazardous rock bursts in underground mines. The current study has proved the high potential for obtaining a broader understanding of anthropogenic hazards by combining different Earth monitoring techniques. Joint geodetic and seismological observations can complement each other and help to overcome the individual disadvantages of the separate methods. This approach can be recommended as a complex tool, which can be used to gain a better understanding of the seismicity associated with all kinds of anthropogenic activities. The method would be especially useful in areas with sparse seismic monitoring and/or regions with already known isotropic source mechanisms, such as mining or volcanic.