Original Article

Geomechanical analysis of subsidence-induced sinkholes and drempels from longwall coal mining using the material point method: A case study from South-Limburg, the Netherlands

J.L. González Acosta, B.B.T. Wassing, J.E. Martins, E. van Linden and J. Hasselman

TNO Geological Survey, Utrecht, the Netherlands

Abstract

Subsidence-induced surface instabilities, such as vertical sharp steps in the topography (known as ‘drempels’ in Dutch), subsidence, and sinkholes, can severely impact infrastructure and public safety. These surface displacements features have been observed during and immediately after coal mining in the South Limburg province of the Netherlands. In this study, a numerical framework for modelling large deformation processes, namely the material point method (MPM), is applied to investigate the underlying mechanisms driving the formation of these complex surface features during coal seam excavation. The modelling of this study focuses on a case study from Heerlen, South Limburg, where a shopping centre in the ’t Loon area partially collapsed due to the development of a sinkhole. Different seam configurations and excavation procedures are tested to assess their influence on the magnitude and spatial distribution of ground deformation. The MPM simulations demonstrate a clear connection between the observed deformation patterns and the early longwall mining processes. Modelled stress and strain concentrations coincide with the observed locations of both drempels and the sinkhole at the ground surface. Furthermore, additional insights were obtained. For example, the direction of seam excavation was found to contribute significantly to the overall distribution of deformations, with larger deformations occurring near the starting position and decreasing towards the end excavation. In contrast, reducing the seam dip angle further amplified the deformations, regardless of the excavation direction, with horizontal seams producing more pronounced effects. These behaviours arise in the absence of geological heterogeneities, indicating that operational configuration alone can predispose sites to the development of surface anomalies. This set of results demonstrates that specific mining and excavation configurations can trigger distinct surface deformation features at predictable locations. The study highlights the potential of MPM to capture the complex mechanisms driving mine collapse and ground subsidence, offering both improved understanding and a means to identify vulnerable areas for post-mining geohazard assessment based on subsurface configurations.

Keywords: drempels; longwall mining; material point method; numerical modelling; subsidence

 

 

Introduction

Longwall mining is a coal extraction technique that originated in the late 17th century, when it was known as the Shropshire method (Mangal, 2013). Between 1910 and the late 1940s, numerous mechanical advancements occurred in both the United States and Europe, culminating in the development of modern mechanised longwall mining in the early 1950s (Berryann & Voelker, 2005; Peng et al., 2019). This technique involves excavating a coal seam along rectangular blocks, known as longwall panels, which typically range from several hundreds of metres to a few kilometres in length and are generally 200–300 metres wide. Coal is extracted from the working face and transported to the surface. The roof immediately behind the working face is temporarily supported to provide a safe working environment for miners. As the coal face and supports gradually advance along the length of the mining panel, the support for the overlying rocks behind the coal face is removed, leading to the collapse of the rock into the void, referred to as ‘goaf’ or ‘gob’ (Figure 1). Some methods, such as roof bolting (Peng, 2015) and backfilling techniques (Zha et al., 2011), are applied both during and after excavation to control and reduce the effects of roof collapse. Nevertheless, some degree of collapse and the subsequent compaction of the excavated void generally remain unavoidable. The collapse of the mine causes deformation of the overlying strata, which manifests over both short- and long-term periods (Dudek et al., 2022). Although mine collapse can be considered the principal cause of surface deformation, additional hydro-mechanical processes also play an important role. In particular, underground exploitation requires prior dewatering of the rock mass and connected aquifer systems, leading to a reduction in pore pressure, an increase in effective stress, and the consequent compaction of the aquifer system (Guzy & Witkowski, 2021). As a result, the surface deformation observed represents the superposition of two processes: direct mining-related subsidence and compaction induced by groundwater depletion. However, these processes may occur at markedly different magnitudes, with mining-induced subsidence typically reaching the order of metres, whereas compaction due to groundwater depletion is generally limited to centimetre-scale displacements.

Figure 1. Change of permeability due to mine collapse and ground deformation (e.g. Khanal et al., 2019).

During the deformation, the permeability of the underground strata also increases due to the formation of fractures and pathways, enabling easier water flow (Figure 1). In the short term, subsidence caused by the successive removal of hydraulic supports and the collapse of the mine roof can lead to the formation of drempels (vertical sharp steps in the topography). Over longer periods, even long after mining activities have ceased, uplift may occur due to the recovery of groundwater levels (Malinowska et al., 2020). Finally, in extreme cases, groundwater inflow into an already fractured and mechanically weakened rock mass may enhance suffusion processes, leading to the development of sinkholes (Pan et al., 2018; Walczak et al., 2025; Xiao et al., 2018). This process can take many years to fully manifest and have been observed all around the world (e.g. Bell et al., 2000; Vervoort, 2020; Zhao & Konietzky, 2021).

These local deformation processes can result in significant differential displacements over short distances, potentially causing severe damage to infrastructure and buildings. Therefore, anticipating the development of such deformation processes in connection to longwall mining activities is critical. This article aims to demonstrate that the observed drempels and ground subsidence are a direct consequence of early longwall mining processes. To investigate this, a large-deformation numerical framework is implemented, namely the Material Point Method (MPM). This study explores whether the observed spatial distribution of the sinkhole and associated drempels at ’t Loon can be reproduced numerically. In particular, the excavation direction and seam dip angle are examined to determine their role in the magnitude and distribution of the ground deformations.

This investigation focuses on a case study in Heerlen, South Limburg, where a shopping centre collapsed due to the formation of a sinkhole. Firstly, an overview of the historical and recent mining context in South Limburg is provided. This is followed by a concise explanation of the MPM framework. Finally, the results are presented, illustrating a strong correlation between the numerical findings and the observed ground anomalies.

Case history

Coal mining in the Netherlands has a rich history dating back to the 11th century in South Limburg (Figure 2a), and has evolved significantly over the centuries (Vis et al., 2020). Initially, coal was extracted using simple (and sometimes rudimentary) methods. However, by the end of the 19th century, technological advancements allowed the efficient extraction of coal at depths up to 1100 metres. Major development efforts in the Limburg coal field during the late 19th and early 20th centuries resulted in substantial increases in coal production. In 1974, all Dutch mines were closed primarily due to the large-scale adoption of natural gas as an alternative energy source (de Jong, 2004). Following the closure of the mines, distinct land surface changes have been observed. During the mining period, these changes manifested as subsidence, while in the post-mining period, uplift occurred due to the increase in water pressure following the flooding of the mines and the cessation of pumping activities (Caro Cuenca et al., 2013; Vis et al., 2020). These ground movements have caused significant damage to infrastructure and buildings in the surrounding areas. Moreover, such effects are expected to continue for decades to come, as mine water levels are still rising and have not yet reached equilibrium.

Figure 2. (a) Coal district located at the Limburg province, south of Netherlands, and (b) seam distribution and sinkhole location at the ’t Loon centrum (adapted from the map provided by Ministerie van Volkshuisvesting en Ruimtelijke Ordening, Kaarten van Nederland).

One of the most striking examples of post-mining ground instability occurred in late 2011, in the city of Heerlen, South Limburg. Between the end of August and early December 2011, cracks began to appear in several structural elements of the ’t Loon shopping centre. In December 2011, a sinkhole formed, ultimately causing the collapse of the structure.

Extensive investigations revealed that the shopping centre was located above an old mine excavated during the 50s (Figure 2b). It was theorised that the event was triggered by the washing away of fine particles due to the rising water table, combined with the formation of fractures and cavities resulting from early coal mining activities. More detailed information on the technology used to investigate the subsidence and a precursor signal to sinkhole formation can be found in a study by Chang and Hanssen (2014).

Numerical investigation

The MPM

The MPM is a numerical technique that originated in the 1990s and was developed to simulate large deformation problems considering solid materials (Sulsky et al., 1994, 1995). Over the years, the method has evolved into a robust tool and is now widely used to analyse a broad range of geotechnical problems (Abe et al., 2014; González Acosta et al., 2021; Sołowski & Sloan, 2015; Zhao et al., 2022). Furthermore, MPM has been successfully applied in the analyses of longwall mining, capturing the unique characteristics of mine behaviour during excavation (Xia et al., 2024; Zhang et al., 2024; Zhou et al., 2023b).

The principles of MPM are analogous to those of standard finite element method (FEM). Figure 3 illustrates these principles, where the red dots depict the material points (analogous to Gauss integration points in FEM) and the grey area represents the discretised body. Initially, material properties are mapped to the mesh nodes (Figure 3a), after which the equilibrium equations are solved at the nodes, resulting in mesh deformation (Figure 3b). A key distinction, however, is that the mesh used for solving the equilibrium equations is not essential for the geometric discretisation of the analysed body. This allows the mesh to be reset to its initial position after each time step, thereby avoiding mesh distortion and preventing the accumulation of deformation (Figure 3c).

Figure 3. MPM solution steps. (a) nodal solution, (b) mesh deformation, and (c) mesh reset (after González Acosta et al., 2024).

Since MPM is based on the standard FEM formulation, the equations of mass and momentum conservation in MPM are as follows:

where ρ is the mass density, v is the velocity vector, t is the time, a is the acceleration vector, σ is the Cauchy stress tensor, b is the body force (gravity), and ∇ · denotes the divergence vector. Furthermore, supplementing these previous equations, stress-strain laws are considered as follows:

where ε represents the strain increment vector, ∇ is the gradient operator, and D is the material elastic matrix. Then, considering FEM discretisation procedures, these equations can be formulated in matrix notation as follows:

where

where nmp is the number of material points in the computational domain, p is the material point index, M, K, F^ext and F^int are the mass, stiffness, external force, and internal forces matrices, respectively, ā and ū are the vectors of nodal acceleration and displacement, respectively, N and B are the shape function and strain-displacement matrices, respectively, W is the material point weight, and J is the Jacobian. Finally, by considering Newmark’s constant average acceleration time integration formulation (Newmark, 1959), Equation 5 can be expressed in its implicit scheme formulation as follows:

where

where is the vector of nodal velocities, m is analogous to M, and ∆t is the time increment. These governing equations are the same as those used in FEM and MPM. However, because material points (which serve a role analogous to Gauss points in FEM) move through the mesh, enhancements to the standard MPM formulation are required to preserve accuracy. Details of these enhancements, as well as the specific MPM implementation used in this investigation, can be found in González Acosta (2020) and González Acosta et al. (2020). Note that this implementation has been developed in an in-house framework coded in FORTRAN, with the constitutive integration following the implicit return mapping approach described in a study by de Souza Neto et al. (2008).

Case study settings and boundary conditions

To investigate the deformation processes that led to the formation of drempels and the sinkhole at the ’t Loon site, a numerical model was developed using historical mining records and subsurface exploration data. The model was refined to incorporate site-specific geological and structural conditions based on findings from underground investigations and archived documentation. Figure 4 illustrates the urban area surrounding the ’t Loon building, the location of the underground mine, and the approximated position of the observed superficial irregularities, including the drempels and the sinkhole. Figure 5 presents the soil stratification based on data from the borehole indicated in Figure 4. The stratigraphy comprises a sand stratum measuring 55 m, a limestone layer of 25 m, and a claystone layer of 31 m.

Figure 4. Distribution of the urban area around ’T Loon shopping mall, mine seam, sinkhole, drempels and borehole at the ’t Loon site.

 

Figure 5. Underground stratification at the ’t Loon site from borehole data.

Based on the approximated mine position and its depth shown in Figures 4 and 5, Figure 6 has been created, where the geometry used to perform the numerical simulations is delineated. The domain configuration is defined by a horizontal distance of 438 m, from which 138 m corresponds to the width of the excavated seam. The coal seam is modelled as 4 m high, with a dip angle of approximately 9.5° (m = –0.166), with the left and right edges located at a depth of 85 m and 108 m, respectively. Note that the original coal layer height is smaller, approximately 1 metre; however, to avoid small-scale problems, specifically those associated with computational limitations, and facilitate the collapse of the gallery, the height was increased to 4 metres. The approximated position of the ’t Loon centrum (A – A’) is at 40 m (to the left) and 33 m (to the right) of point B corresponding to the starting point of the excavated seam, and the sinkhole position is at 25 m (to the right) of point B. Regarding the material properties, Figure 7 shows the vertical distribution of the elastic modulus E derived from a cone penetration test (CPT), and unconfined compressive strength (UCS) tests. To define the material strength, the Mohr-Coulomb constitutive description has been considered, and the values of friction (ϕ) and cohesion are shown in Figure 8. Note that the full procedure for deriving E and ϕ can be found in Robertson and Cabal (2022) and Kulhawy and Mayne (1990), respectively.

Figure 6. Updated model considering borehole, mine dipping and urban map (not to scale). Sinkhole has been estimated to be (horizontally) 25 m away from the mine top edge.

 

Figure 7. Approximated elastic modulus E (in red) based on CPT exploration (Robertson & Cabal, 2022), and UCS tests, where I_c is the soil index, q_c are the CPT lectures and σ_v is the underground vertical stress.

 

Figure 8. Strength properties as function of depth. (a) friction angle estimated using ϕ = 17.6 + 11 log (q_c−σ_v) (Kulhawy & Mayne, 1990), and (b) cohesion estimated from UCS lectures.

Regarding the numerical boundary conditions of the model, the external boundaries (left, right, and bottom) are defined as rollers, allowing only tangential displacements along the boundaries while restraining normal displacements. The background mesh used for the analysis is constructed with square elements measuring 2 × 2 m, initially containing four material points per element. In all analysis, the excavation is simulated by removing portions of 4 m of material points from the computational domain at each step (Figure 9), where each step consist of a period of 5 seconds, and the time discretisation is of ∆t = 0.005 s.

Figure 9. Excavation procedure steps (not to scale).

Results

This section is divided into two parts. The first part presents the results of the problem described in the previous section, referred to as the base case. Additionally, two base case scenarios are analysed: the excavation from B to B’ (top-to-bottom) and from B’ to B (bottom-to-top) to evaluate the impact of the excavation direction on ground deformation. The second part, referred to as the sensitivity analyses, consists of simulating the excavation considering different dip angles of the seam, aiming to better understand the influence of the dip angle on the ground deformation. These results summarise the performed tests into the mining configuration factors influencing ground deformation, thereby enhancing our understanding of the process.

Base case: influence of mining direction

Figure 10 illustrates the vertical displacement after the mine has been fully excavated. As shown, the excavation direction significantly influences the results. Excavating from top-to-bottom produces substantial vertical deformation (Figure 10a), whereas excavation from bottom-to-top results in smaller vertical deformation by approximately 20% less vertical deformation (Figure 10b). The most notable aspect of this figure is that the estimated position of the sinkhole and observed drempels lie within the fracture zone (i.e. triangular zone prone to increased permeability) indicated in Figure 1. This conclusion is further supported by Figure 11, which shows the horizontal strains after excavation. Here, it is more evident that the maximum strains, developed at point B and the extent to the ground surface, encompass the area where the drempels and the sinkhole are located. Similar to Figure 10, when the excavation starts from the top (Figure 11a), a distinct concentration of strains is observed. In contrast, when excavation starts from the bottom, the strains appear more diffuse (Figure 11b).

Figure 10. Vertical displacement of the dipping mine considering an excavation from (a) top to bottom, and from (b) bottom to top.

 

Figure 11. Horizontal strains of the dipping mine considering an excavation from (a) top to bottom, and from (b) bottom to top.

Figure 12 summarises the various ground surface outcomes resulting from the mining process. Each subfigure highlights a distinct feature associated with the observed drempels and sinkhole. The subsidence and horizontal displacement figures (Figure 12a and b) indicate that the anomalies are located where ground movement begins rather than where the movements are more significant. This observation is more evident in Figure 12c and d, where the strains clearly depict significant amplifications at the location of the anomalies. Notably, the sinkhole is situated at the point of maximum vertical positive strains, indicating vertical stretching of the ground. Conversely, the drempels are located at the points of maximum horizontal positive strains, reflecting horizontal stretching of the ground. Additionally, the drempels and the sinkhole are enclosed within the area of high shear strains (Figure 12e).

Figure 12. Ground surface (a) subsidence, (b) horizontal deformation, (c) horizontal strains, (d) vertical strains, and (e) shear strains.

It is important to observe that, although the formation of the gob is not directly represented in the numerical model, the vertical displacement of the mine roof and its eventual contact with the floor serves as an approximation of the stabilising influence typically provided by the collapsed material.

Sensitivity analyses: influence of the dipping angle

A series of analyses were conducted to quantify the magnitude of ground vertical displacements as a function of the dipping angle. The slope values chosen for the analysis are m = –0.166 (9.5°) (Figure 13a), –0.10 (5.7°) (Figure 13b), –0.05 (2.8°) (Figure 13c), and 0.0 (Figure 13d). It should be noticed that these results, along with subsequent findings, are based solely on the top-to-bottom excavation framework, as this approach produces greater amplification of surface irregularities. The analyses show that the slope of the excavated mine has a significant effect on the maximum vertical displacement. As the slope m decreases, the maximum vertical displacements increase, a trend also reported by Luo (2015). In this specific case, the ratio between the maximum subsidence (δ_v-m) when m is = -0.166 and m is = 0 is approximately δ_v-0.0 = 0.66 δ_v-0.166. Similarly, the lateral extent (LD) of vertical deformation increases as the slope decreases. This is more evident in Figure 14, where the lateral extension of the vertical deformation is highlighted by a yellow triangle. In this case, the ratio between the lateral extension (δ_LD-m) when m is –0.166 and 0 is approximately δ_LD-0.0 = 0.63 δ_LD-0.166. The correlation between the mine slope and the draw angle α, which describes the vertical distribution of lateral deformation, is shown in Figure 15. This nearly linear correlation can be expressed as: α = –1.1781m + 30.551.

Figure 13. Vertical displacement considering a changing slope m of the excavated mine equal to (a) –0.166 (9.5°), (b) –0.10 (5.7°), (c) –0.05 (2.8°), and (d) 0.0.

 

Figure 14. Lateral extension (LD) of the vertical deformation and angle (α) created in function of the mine slope m, where the slope considered are (a) –0.166 (9.5°), (b) –0.10 (5.7°), (c) –0.05 (2.8°), and (d) 0.0.

 

Figure 15. Linear correlated between the slope m of the mine and the formed lateral extension angle α.

Figure 16 illustrates the horizontal strains as a function of the dipping angle. Consistent with the results shown in Figure 13, a reduction in m leads to an increase in the magnitude of the horizontal strains. In Figure 16a and b, the positive strains developed at the right-side mine edge exhibit a narrow vertical distribution. However, in Figure 16c and d, the reduction of the dipping angle broadens the width of the shear stress distribution. Additionally, the shear strains do not follow a purely vertical pattern but instead adopt a curved trajectory, increasing the area of the ground surface influenced by the strains.

Figure 16. Horizontal strains considering a slope m of the excavated mine equal to (a) –0.166, (b) –0.10, (c) –0.05, and (d) 0.0.

Figure 17 and Table 1 provide a summary of the results, specifically the ground surface features, as a function of the slope m. Figure 17a and b show that as m reduces, the amplitude and extent of the vertical and horizontal displacements increase substantially, with significant displacement amplitudes observed at the approximate sinkhole location. Although these two figures do not clearly indicate the presence of drempels, Figure 17c and d provide stronger evidence. This is demonstrated by the pronounced strains at the locations of the drempels. Since strains represent accumulated changes in dimension per unit length, pronounced strains over a short length may indicate the presence of characteristic vertical sharp steps, or drempels. Furthermore, this observation is reinforced in Figure 17e, where the shear strains display even sharper amplifications at the same location.

Figure 17. Ground surface (a) subsidence (δ_V), (b) horizontal displacements (δ_H), (c) horizontal strains (ε_h), (δ) vertical strains (ε_v), (e) shear strains (ε_hv).

Table 1. A summary of surface ground irregularities in function of the mine slope m.

Slope m	Max. δV (m)	Max. δH (m)	Max. εv (–)	Max. εh (–)	Max. εhv (–)	LD (m)	α (°)
–0.16	–0.83	0.53	0.016	0.014	0.0018	30.9	20
–0.10	–1.05	0.55	0.017	0.014	0.0016	37.9	23
–0.05	–1.23	0.69	0.021	0.018	0.0020	44.0	27
0.0	–1.24	0.77	0.022	0.019	0.0032	48.9	31

Finally, Figure 18 shows the distribution of the plastic strain invariant as a function of m, which is computed as where is the plastic strain tensor. This metric is used to identify the regions where plasticity, and therefore failure, is more prominent. This analysis focuses on the worst-case scenario, namely the top-to-bottom excavation. In all cases, plastic strains are concentrated at the corner of the seam where the excavation begins. Although a reduction in m decreases the spread of plastic strains across the domain, higher m values concentrate larger plastic strains at the seam corner, forming a vertical band that almost reaches the ground surface when m is = –0.166 (Figure 18a). Notably, it is only for m = –0.166 that the plastic strains nearly reaches the surface, coinciding with the location of the observed drempels and the subsequent sinkhole. These observations help to explain the asymmetry in subsidence profiles, particularly in the horizontal case. The extensive development of plastic strains at the start of the excavation results in more pronounced and widespread deformations above the initial excavation point. As the excavation progresses, the reduction in plastic strains along the direction of advance leads to smaller subsidence magnitudes towards the end of the excavation. This suggests that, in the presence of plastic deformations, subsidence is generally larger above the starting point than above the end point, unless modified by boundary conditions. In addition, plastic strains also concentrate at the interface between the sand and limestone layers (Figure 18b, c and d), which may lead to the formation of voids. These voids can enable material flow and potentially triggering the sinkhole.

Figure 18. Accumulated plastic strain invariant considering a slope m of the excavated mine equal to (a) –0.166, (b) –0.10, (c) –0.05, and (d) 0.0.

Discussion

Comparison with literature

The numerical results show patterns consistent with those reported in historical longwall mining cases and previous numerical studies. The simulated evolution of the ground settlement profile during excavation is consistent with observations from historical, field measured, case studies (National Coal Board, 1975; Sasaoka et al., 2015; Zhou et al., 2023a) and with results obtained from other numerical frameworks (Li et al., 2022; Zhang et al., 2022). The draw angle (α) obtained from different geometrical configurations also matches values reported in previous studies (e.g. Karmis et al., 1982; Strzałkowski, 2021). In addition, the stresses redistribution patterns near the seam boundaries during excavation and in the final configuration agree with those estimated analytically (Li et al., 2019; Suchowerska et al., 2013) and numerically (Du et al., 2020; Zhang et al., 2016).

Although fracturing and material migration were not explicitly modelled, as demonstrated in other MPM studies of longwall mining (e.g. Zhou et al., 2023b), the present implementation was nevertheless sufficient to reproduce the observed deformations at the ’t Loon location. In particular, the position of the observed drempels and the emergence of the sinkhole align spatially with the zones of maximum deformation and strain concentration. These findings support the conclusion that the MPM is robust enough to simulate and improve the understanding of longwall mining mechanisms.

New insights from the simulations

Beyond reproducing well-established subsidence mechanisms, the simulations provided new insights that have not been widely reported in subsidence research. The results showed that excavation direction controls both the magnitude and localisation of deformation, with asymmetric settlement troughs forming and higher strain accumulation near the starting side. Excavating from top-to-bottom generated approximately 20% greater vertical displacement than from bottom-to-top, with deformation zones concentrated above the initial levels of advance. These asymmetries correspond to higher strain accumulation near the starting panel edge, a mechanism not widely documented in empirical subsidence frameworks. One possible explanation for this limited documentation is the assumption that most mining operations progress from the deepest coal seams to the shallowest, taking advantage of gravity to facilitate coal extraction. However, several mining maps from the Limburg coal district indicate that this was not always the case. Historical mining maps suggest that the excavation direction varied among different seams.

In addition, reducing the dip angle consistently amplified both vertical and LD of deformation, with horizontal seams producing the strongest deformation response. These findings extend empirical observations of subsidence by demonstrating that excavation configurations can systematically influence the development of surface anomalies such as drempels and sinkholes.

Limitations of the investigation

Despite the strong correlation between the numerical results and the observed ground features, drempels and sinkholes, further refinements are necessary to fully reproduce their geometries and formation mechanisms in the simulations. To address this limitation, four key enhancements are proposed to improve the accuracy of MPM simulations in modelling longwall mining and the formation of drempels and sinkholes:

• Element discretisation: The size of the background mesh is a critical factor in simulating geomechanical processes such as mine collapse, drempels, and sinkholes. Coarse meshes may suppress shear localisation and underestimate deformation magnitudes leading to inaccurate results (Beuth et al., 2008; Ma et al., 2005). A mesh with sufficient geometric resolution is therefore essential for representing these features. However, the use of finer meshes is directly linked to computational cost, making them often infeasible in large-scale simulations. This problem was also encountered in this study, where a mesh size had to be chosen that was small enough to capture deformation patterns but large enough to ensure computational feasibility.
• Constitutive model: The choice of constitutive model strongly influences the representation of soil failure and the resulting distribution of shear stresses and strains. The standard Mohr-Coulomb criterion, used in this research, cannot fully capture soil behaviour such as hardening and unloading (Obrzud, 2010). Advanced constitutive formulations, such as those developed in González Acosta et al. (2024), can improve simulations by representing brittle fracture processes, that are common in mine collapses, and therefore scaling the amplitude of the deformations. These models also reduce mesh dependence effects (González Acosta & Mánica, 2024), leading to a more realistic representation of failure paths. Although not included in this study, adopting such models would enhance the predictive power of MPM in post-mining environments.
• Time-stepping and numerical phasing: This study showed that excavation direction greatly influences deformation during mine collapse. In a dynamic formulation, however, the phasing of excavation (i.e. the time between excavation steps) combined with the chosen time step also plays a critical role, as it can either suppress or accelerate the accumulation of plastic deformations. This aspect was not explored here due to computational limitations and has rarely been addressed in other studies, yet it is important to consider in future work to better understand the full range of simulation outcomes.
• Validation with levelling measurements at the time of the mining: Levelling measurements were carried out during the mining period. Comparing the magnitude of the modelled strains with such measurements is crucial, as it allows validation of the model against real data. However, these measurements are only available in paper archives, making it impossible to extract them within the timeframe of this study.

There are several additional aspects of the model that could lead to different results, including material properties, the location of the boundary conditions, alternative loading conditions, and the use of other frameworks such as a hydro-mechanical MPM. Exploring all aspects simultaneously was beyond the scope of this study due to computational complexity, primarily due to the extremely large numerical effort required to explore them in detail. Nevertheless, future research is encouraged to examine these conditions further in order to advance the understanding of numerical modelling in the context of longwall mining problems.

Conclusions

This study has demonstrated the effectiveness of the MPM in simulating longwall mining processes and the associated ground deformations. By replicating a set of excavation scenarios comparable to real longwall mining conditions, such as mining direction and dipping angle, the MPM applied model successfully reproduced realistic deformation patterns, including subsidence and strain distributions. Although the simulations did not explicitly reproduce sinkholes or drempels, the results showed strong spatial correlation with observed ground surface anomalies. Excavation direction was found to significantly influence the resulting deformations, with top-to-bottom excavation producing approximately 20% more vertical deformation than bottom-to-top excavation. In addition, seam dip was shown to further amplify deformations. Vertical deformation increased as the dip angle decreased, while the LD of deformation increases with steeper dips. A nearly linear correlation was found between the seam slope and the angle describing the distribution of lateral deformation. Notably, both the sinkhole and the observed drempels are located within the modeled fracture zone, where increased permeability is expected.

These findings suggest that the method captures the underlying mechanisms driving the development of such features. Furthermore, the MPM proved to be a robust tool for testing complex hypotheses related to ground deformation. In particular, the hypothesis that both drempel and sinkhole formation may be linked to the collapse of abandoned mine workings was supported by the simulation outcomes, highlighting the method’s potential as a predictive tool post-mining geohazards assessments.

References

Abe, K., Soga, K. & Bandara, S., 2014. Material point method for coupled hydromechanical problems. Journal of Geotechnical and Geoenvironmental Engineering 140(3): 04013033. DOI: https://doi.org/10.1061/(ASCE)GT.1943-5606.0001011

Bell, F.G., Stacey, T.R. & Genske, D.D., 2000. Mining subsidence and its effect on the environment: some differing examples. Environmental Geology 40(1): 135–152. DOI: https://doi.org/10.1007/s002540000140

Berryann, R.J. & Voelker, J.A., 2005. Evolution of longwall mining and control systems in the United States. Triadelphia, WV: Mine Safety and Health Administration, Division of Electric Safety.

Beuth, L., Benz, T., Vermeer, P.A. & Wieckowski, Z., 2008. Large deformation analysis using a quasi-static material point method. Journal of Theoretical and Applied Mechanics 38(1–2): 45–60.

Caro Cuenca, M., Hooper, A.J. & Hanssen, R.F., 2013. Surface deformation induced by water influx in the abandoned coal mines in Limburg, The Netherlands observed by satellite radar interferometry. Journal of Applied Geophysics 88: 1–11. DOI: https://doi.org/10.1016/j.jappgeo.2012.10.003

Chang, L. & Hanssen, R.F., 2014. Detection of cavity migration and sinkhole risk using radar interferometric time series. Remote Sensing of Environment 147: 56–64. DOI: https://doi.org/10.1016/j.rse.2014.03.002

de Jong, T.P.R., 2004. Coal mining in the Netherlands; the need for a proper assessment. Geologica Belgica 7(3–4): 231–243.

de Souza Neto, E.A., Peric, D. & Owen, D.R., 2008. Computational methods for plasticity: theory and applications. John Wiley & Sons, New York.

Du, B., Liu, C., Yang, J. & Wu, F., 2020. Abutment pressure distribution pattern and size optimization of coal pillar under repeated mining: a case study. Arabian Journal of Geosciences 13(23): 1261. DOI: https://doi.org/10.1007/s12517-020-06281-y

Dudek, M., Sroka, A., Tajduś, K., Misa, R. & Mrocheń, D., 2022. Assessment and duration of the surface subsidence after the end of mining operations. Energies 15(22): 8711. DOI: https://doi.org/10.3390/en15228711

González Acosta, J.L., 2020. Investigation of MPM inaccuracies, contact simulation and robust implementation for geotechnical problems. PhD thesis, Delft University of Technology, Delft.

González Acosta, J.L. & Mánica, M.A., 2024. Investigación preliminar sobre el deslizamiento de laderas utilizando el método del punto material y una regularización no local. In: N.P. López Acosta & A.L. Espinosa Santiago (Eds.), XXXII Reunión Nacional de Ingeniería Geotécnica: Volume 2, 3–7 September 2024 (pp. 665–670). Ciudad de México: Sociedad Mexicana de Ingeniería Geotécnica.

González Acosta, J.L., Mánica, M.A., Vardon, P.J., Hicks, M.A. & Gens, A., 2024. A nonlocal material point method for the simulation of large deformation problems in brittle soils. Computers and Geotechnics 172: 106424. DOI: https://doi.org/10.1016/j.compgeo.2024.106424

González Acosta, J.L., Vardon, P.J. & Hicks, M.A., 2021. Study of landslides and soil-structure interaction problems using the implicit material point method. Engineering Geology 285: 106043. DOI: https://doi.org/10.1016/j.enggeo.2021.106043

González Acosta, J.L., Vardon, P.J., Remmerswaal, G. & Hicks, M.A., 2020. An investigation of stress inaccuracies and proposed solution in the material point method. Computational Mechanics 65: 555–581. DOI: https://doi.org/10.1007/s00466-019-01783-3

Guzy, A. & Witkowski, W.T., 2021. Land subsidence estimation for aquifer drainage induced by underground mining. Energies 14(15): 4658. DOI: https://doi.org/10.3390/en14154658

Karmis, N., Chen, C.Y., Jones, D.E. & Triplett, T., 1982. Some aspects of mining subsidence and its control in the US coalfields. Minerals and the Environment 4(4): 116–130. DOI: https://doi.org/10.1007/BF02085973

Khanal, M., Adhikary, D., Poulsen, B. & Guo, H., 2019. Mine site specific longwall-induced permeability changes: a case study. Geotechnical and Geological Engineering 37(4): 3271–3282. DOI: https://doi.org/10.1007/s10706-019-00842-z

Kulhawy, F.H. & Mayne, P.W., 1990. Manual on estimating soil properties for foundation design (Research Project 1493-6, Report EL-6800). New York: Cornell University, Geotechnical Engineering Group.

Li, X., Wang, Y., Hu, Y., Zhou, C. & Zhang, H., 2022. Numerical investigation on stratum and surface deformation in underground phosphorite mining under different mining methods. Frontiers in Earth Science 10: 831856. DOI: https://doi.org/10.3389/feart.2022.831856

Li, Y., Lei, M., Wang, H., Li, C., Li, W., Tao, Y. & Wang, J., 2019. Abutment pressure distribution for longwall face mining through abandoned roadways. International Journal of Mining Science and Technology 29(1): 59–64. DOI: https://doi.org/10.1016/j.ijmst.2018.11.018

Luo, Y., 2015. An improved influence function method for predicting subsidence caused by longwall mining operations in inclined coal seams. International Journal of Coal Science and Technology 2: 163–169. DOI: https://doi.org/10.1007/s40789-015-0086-x

Ma, J., Lu, H., Wang, B., Roy, S., Hornung, R., Wissink, A. & Komanduri, R., 2005. Multiscale simulations using generalized interpolation material point (GIMP) method and SAMRAI parallel processing. Computer Modeling in Engineering & Sciences 8(2): 135.

Malinowska, A.A., Witkowski, W.T., Guzy, A. & Hejmanowski, R., 2020. Satellite-based monitoring and modeling of ground movements caused by water rebound. Remote Sensing 12(11): 1786. DOI: https://doi.org/10.3390/rs12111786

Mangal, A., 2013. Data interpretation of periodic strata loading characteristics in a mechanized powered support longwall caving face-a trend setter. Indian Journal of Engineering 4(10): 31–36.

Martins, J., van Linden, E., Fruijter, C. & Davids, B., 2025. Atlas of the Limburg coal mining region: a comprehensive geospatial perspective from the subsurface to satellite observations, Submitted, this issue.

Ministerie van Volkshuisvesting en Ruimtelijke Ordening, Kaarten van Nederland. https://www.kaartenvannederland.nl

National Coal Board, 1975. Subsidence engineers’ handbook. London: Great Britain.

Newmark, N.M., 1959. A method of computation for structural dynamics. Journal of the Engineering Mechanics Division 85(3): 67–94. DOI: https://doi.org/10.1061/JMCEA3.0000098

Obrzud, R.F., 2010. On the use of the Hardening Soil Small Strain model in geotechnical practice. Numerics in Geotechnics and Structures 16: 1–17.

Pan, Z., Jiang, X., Lei, M., Guan, Z., Wu, Y. & Gao, Y., 2018. Mechanism of sinkhole formation during groundwater-level recovery in karst mining area, Dachengqiao, Hunan province, China. Environmental Earth Sciences 77(24): 799. DOI: https://doi.org/10.1007/s12665-018-7987-0

Peng, S.S., 2015. Topical areas of research needs in ground control – a state of the art review on coal mine ground control. International Journal of Mining Science and Technology 25(1): 1–6. DOI: https://doi.org/10.1016/j.ijmst.2014.12.006

Peng, S.S., Du, F., Cheng, J. & Li, Y., 2019. Automation in US longwall coal mining: a state-of-the-art review. International Journal of Mining Science and Technology 29(2): 151–159. DOI: https://doi.org/10.1016/j.ijmst.2019.01.005

Robertson, P.K. & Cabal, K., 2022. Guide to cone penetration testing (7th ed.). California: Gregg Drilling LLC.

Sasaoka, T., Takamoto, H., Shimada, H., Oya, J., Hamanaka, A. & Matsui, K., 2015. Surface subsidence due to underground mining operation under weak geological condition in Indonesia. Journal of Rock Mechanics and Geotechnical Engineering 7(3): 337–344. DOI: https://doi.org/10.1016/j.jrmge.2015.01.007

Sołowski, W.T. & Sloan, S.W., 2015. Evaluation of material point method for use in geotechnics. International Journal for Numerical and Analytical Methods in Geomechanics 39(7): 685–701. DOI: https://doi.org/10.1002/nag.2321

Strzałkowski, P., 2021. Duration of the final phase of mining area deformation process in the conditions of Upper Silesia (Poland). Bulletin of Engineering Geology and the Environment 80(7): 5769–5780. DOI: https://doi.org/10.1007/s10064-021-02272-9

Suchowerska, A.M., Merifield, R.S. & Carter, J.P., 2013. Vertical stress changes in multi-seam mining under supercritical longwall panels. International Journal of Rock Mechanics and Mining Sciences 61: 306–320. DOI: https://doi.org/10.1007/s10064-021-02272-9

Sulsky, D., Chen, Z. & Schreyer, H.L., 1994. A particle method for history-dependent materials. Computer Methods in Applied Mechanics and Engineering 118(1–2): 179–196. DOI: https://doi.org/10.1016/0045-7825(94)90112-0

Sulsky, D., Zhou, S.J. & Schreyer, H.L., 1995. Application of a particle-in-cell method to solid mechanics. Computer Physics Communications 87(1–2): 236–252. DOI: https://doi.org/10.1016/0010-4655(94)00170-7

Vervoort, A., 2020. Long-term impact of coal mining on surface movement: residual subsidence versus uplift. Mining Report Glückauf 156(2): 136–141.

Vis, G-.J., van Linden, E., van Balen, R. & Cohen, K., 2020. Depressions caused by localized subsidence in the Netherlands, Belgium and Germany: a link with coal mining? In: P.A. Fokker & G. Erkens (Eds.), Proceedings of IAHS, 22 April 2020 (pp. 201–205). Delft.

Walczak, S., Witkowski, W.T., Stoch, T. & Guzy, A., 2025. Detecting sinkholes and land surface movement in post-mining regions using multi-source remote sensing data. Remote Sensing Applications: Society and Environment 38: 101560. DOI: https://doi.org/10.1016/j.rsase.2025.101560

Xia, B., Zhou, Y., Zhang, X., Zhou, L. & Ma, Z., 2024. Physical and numerical investigations of target stratum selection for ground hydraulic fracturing of multiple hard roofs. International Journal of Mining Science and Technology 34: 699–712. DOI: https://doi.org/10.1016/j.ijmst.2024.05.003

Xiao, H., Li, H. & Tang, Y., 2018. Assessing the effects of rainfall, groundwater downward leakage, and groundwater head differences on the development of cover-collapse and cover-suffosion sinkholes in central Florida (USA). Science of the Total Environment 644: 274–286. DOI: https://doi.org/10.1016/j.scitotenv.2018.06.273

Zha, J.-F., Guo, G.-L. & Wang, Q., 2011. Mining subsidence control by solid backfilling under buildings. Transactions of Nonferrous Metals Society of China 21: s670–s674. DOI: https://doi.org/10.1016/S1003-6326(12)61660-4

Zhang, C., Mitra, R., Oh, J. & Hebblewhite, B., 2016. Analysis of mining-induced valley closure movements. Rock Mechanics and Rock Engineering 49(5): 1923–1941. DOI: https://doi.org/10.1007/s00603-015-0880-1

Zhang, Y., He, F., Kong, J., Zhu, Y. & Wang, L., 2022. Relationship between surface subsidence range and geological mining conditions using numerical simulation and machine learning. Scientific Programming 2022(1): 8720831. DOI: https://doi.org/10.1155/2022/8720831

Zhang, Y., He, K., Hu, X., Liu, W., Zhang, S., Wu, J. & Xi, C., 2024. Mechanism of surface subsidence and sinkhole formation in mining areas: insights from MPM. Bulletin of Engineering Geology and the Environment 83(8): 330. DOI: https://doi.org/10.1007/s10064-024-03822-7

Zhao, J. & Konietzky, H., 2021. An overview on flooding induced uplift for abandoned coal mines. International Journal of Rock Mechanics and Mining Sciences 148: 104955. DOI: https://doi.org/10.1016/j.ijrmms.2021.104955

Zhao, L., Qiao, N., Huang, D., Zuo, S. & Zhang, Z., 2022. Numerical investigation of the failure mechanisms of soil–rock mixture slopes by material point method. Computers and Geotechnics 150: 104898. DOI: https://doi.org/10.1016/j.compgeo.2022.104898

Zhou, B., Yan, Y. & Kang, J., 2023a. Dynamic prediction model for progressive surface subsidence based on MMF time function. Applied Sciences 13(14): 8066. DOI: https://doi.org/10.3390/app13148066

Zhou, L., Li, X., Peng, Y., Xia, B. & Fang, L., 2023b. Material point method with a strain-softening model to simulate roof strata movement induced by progressive longwall mining. International Journal of Rock Mechanics and Mining Sciences 170: 105508. DOI: https://doi.org/10.1016/j.ijrmms.2023.105508