The Münsterdorf sinkhole cluster: void origin and mechanical failure

Abstract

Since 2004, collapse sinkholes occur on the sports field of Münsterdorf, a village north of Hamburg in Germany. The sinkholes, around 2–5 m in diameter and 3–5 m deep, develop in peri-glacial sand, which at around 20 m depth is underlain by Cretaceous chalk. The chalk has been pushed up close to the surface by a salt diapir. The sinkhole formation initiated suddenly and occurs with a frequency of one every 2 years.

We use a variety of geophysical results (e.g. gravity, electrical resistivity imaging, ground-penetrating radar) from previous fieldwork campaigns together with a new data set from direct-push-based methods to infer mechanical and hydrological properties of the material beneath the sports field (peri-glacial sand, glacial marl, Cretaceous chalk).

Based on the derived material properties, we develop a mechanical model for the sinkhole collapse, starting from simple analytical considerations and then moving towards a three-dimensional distinct-element model explaining the sudden onset of collapse sinkholes for the sports field.

The mechanical model supports our hypothesis that the sudden onset of sinkholes is triggered by changes in groundwater level.

1 Introduction

Collapse sinkholes are cylindrical to elliptical surface structures with diameters ranging from 1 m to several hundred metres and similar depth extensions (e.g. ). Collapse sinkholes form suddenly and pose a substantial risk to infrastructure because of their sudden inception and almost no warning signs before (e.g. ). Often, collapse sinkholes form in karst terrains.

A karst terrain is defined as an area comprising a soluble rock, such as limestone, dolomite, anhydrite, gypsum, or salt. Common to all these soluble rock types is the ability of water flowing through rock fractures and bedding partings to dissolve material from the rock mass. While for the first two types of soluble rocks mentioned, limestone and dolomite, dissolution is substantial only when the water is acidic, e.g. by dissolving carbon dioxide from air and/or soil to form carbonic acid, for the latter rock types – anhydrite, gypsum, and salt – dissolution occurs in pure water.

The dissolution of soluble rock enhances the initially low primary permeability of the rocks, which is controlled by fractures in the sub-millimetre range. With time, the developing secondary permeability becomes orders of magnitude larger. The reason is the enlargement of fissures and bedding partings in the subsurface, in parts to the metre scale, to form subsurface voids and caves. This process occurs naturally over timescales of 10 000–100 000 years (limestone, dolomite), 1000–10 000 years (anhydrite and gypsum), and even faster in salt.

Once the subsurface voids reach a critical size, they can become mechanically unstable and start to collapse. If the void is deep below the surface, the collapse results in roof and wall breakdown, and breakdown can enlarge the void space by forming a mechanically stable shape. If the voids develop closer to the surface, breakdown might slowly migrate towards the surface, until only a thin ledge of rock remains. It is the final collapse of this thin ledge which causes the sudden appearance of the collapse sinkhole (e.g. ).

Depending on the overburden of the soluble rock, collapse sinkholes are classified into bedrock-collapse sinkholes, caprock-collapse sinkholes, and cover-collapse sinkholes (see ; and references therein). All of the three types of sinkholes listed above describe collapse as a process to form the sinkholes, but the surface layer affected by the collapse is different: soluble bedrock in the first case, unconsolidated deposits in the second case, and insoluble rock in the third case.

We discuss a group of collapse sinkholes formed on the sports field in the village of Münsterdorf in northern Germany. As the material above the soluble rock comprises both unconsolidated peri-glacial sand and more consolidated glacial marl, these sinkholes fall somewhere between the cover- and caprock-collapse sinkhole categories. The sinkholes started to appear in 2004 with a rate of one per year, being steep-walled, with a diameter of 2–3 and 2–5 m deep. We explore the mechanical stability (or instability) of the locality and try to answer the following questions:

• i.

  Why did the collapse sinkhole formation start suddenly in 2004?

• ii.

  What mechanical conditions are needed to initiate sinkhole collapse?

• iii.

  Where is the location of the initial subsurface void responsible for the sinkhole formation?

We have organised the paper as follows. In Sect. , we introduce the locality and its geological set-up, and we describe the sinkhole cluster. In Sect. , we briefly summarise the results of previous geophysical campaigns, and we derive a lithological model for the subsurface structure. We then discuss locations of the initial voids responsible for the sinkhole collapse in terms of depth and size, and we investigate the chance to spot these voids prior to sinkhole formation with geophysical methods, depending on their depth. In Sect. , we describe several direct-push-based methods performed on and around the sports field in Münsterdorf. We then discuss hydrological and mechanical properties derived from these direct-push-based methods. In Sect. , we first develop a simple mechanical model for sinkhole collapse, and then we apply a discrete-element model to simulate the collapse sinkholes occurring on the sports field. We discuss the mechanical and hydrological properties likely to trigger the collapse sinkholes. In Sect. , we come back to our hypothesis of the void origin and present geochemical measurements favouring the deep void origin in the chalk layer. In Sect. , we refine our previous hypothesis on the temporal evolution of the collapse sinkholes on the sports field in Münsterdorf. We finally summarise our findings in Sect. .

2 Site

The village of Münsterdorf is located around 50 km north of Hamburg in the northern part of Germany (Fig. ), just south of the river Stör, which runs in a westerly direction close to Münsterdorf, then turns south to flow into the river Elbe. The river Stör is the natural local base level for the region and is influenced by the tides of the North Sea.

Figure 1

Topographical map of the Münsterdorf–Lägerdorf area. Shown are the elevation (colour-coded); the river Stör as the local base level in the north; the villages of Münsterdorf and Lägerdorf (grey dots); the three open-pit chalk mines Saturn (abandoned), Schinkel, and Heidestrasse (black contour lines), including their actual extension (dashed blue line) not present in the SRTM (Shuttle Radar Topography Mission) model; the depth to the chalk layer (red contours); and the recent sinkholes since 2004 (blue dots). Insert: overview map with working area marked in red.

[Figure omitted. See PDF]

2.1 Geological set-up

As most of the landscape in northern Germany, the surface morphology is dominated by the Pleistocene ice-age cycles. The repeated advance and retreat of the Fennoscandian Ice Sheet carried glacial and peri-glacial sediments, which form the main features in the landscape: Münsterdorf sits on a geest ridge, a sandy deposit from deposition below the ice, about 20–30 m above sea level in the region. Beyond this geest island, the landscape is lower, between 5–15 m above sea level, and dominated by swampy marshland. The geest ridge was deposited during the Elsterian (455–320 ka) and Saalian (300–125 ka) glacial cycles, while during the last glacial cycle, the Weichselian, the ice did not reach the area around Münsterdorf.

The geest ridge consists of unconsolidated and consolidated peri-glacial sand, interbedded with glacial till, in parts composed of clay lenses.

The entire region is tectonically controlled by the Krempe–Lägerdorf salt ridge (e.g. ), a salt diapir, with salt sequences from the Permian Zechstein period. These Zechstein rocks can usually be found in a depth range between 3–5 km underneath northern Germany, but as the salt layers in the Zechstein evaporite they can flow plastically under stress; numerous salt structures have been pushed up by the large stress imposed by the overburden. The vertical uplift of these salt domes has pushed up the overburden by several kilometres, often reaching the surface. In the Münsterdorf region, the Krempe–Lägerdorf salt ridge has pushed up Cretaceous chalk from its original depth of 0.5–1 km to a position close to the surface (e.g. ).

Just 1.5 km south-west of Münsterdorf, the chalk is just below the surface (6 m below the ground), and several large open-pit mines exploit the material. The deepest open-pit mine reaches down to 90 m below sea level, which requires a substantial effort to dewater the mines (e.g. ). The pumped groundwater from the open-pit mines is routed via an artificial channel in the east towards the river Stör. In the mines, numerous karst features can be observed, such as dissolutionally enlarged fissures and bedding planes and small collapse sinkholes.

The local sports field is located on the southern rim of the village of Münsterdorf. The chalk can be found here at around 20 m depth below the sports field (see Fig. ), confirmed by several boreholes (see ; and references therein). Above the Cretaceous chalk, we find mainly peri-glacial sand, inter-bedded with glacial till.

From the geological and morphological set-up, we identify three key aspects responsible for the presence of karst structures in the Münsterdorf region:

• salt diapirism responsible for local tectonics and uplift of Permian and Cretaceous sequences;

• Cretaceous chalk around 20 m below the sports field;

• and landscape morphology dominated by ice-age cycles, with deposition of glacial till as a flow barrier.

The glacial till with its lower hydraulic permeability confines groundwater flow, whilst the salt diapirism caused chalk layers to uplift close to the surface (see also ).

2.2 Sinkhole cluster

Since 2004, collapse sinkholes have occurred regularly on the sports field (Fig. ). These collapse sinkholes are about 2–3 m in diameter and about 3–5 m deep, exposing the peri-glacial sand. The sinkholes occur suddenly, often underneath the grass of the sports field, and are initially steep-walled, with near-vertical sides. Most of the sinkholes are filled with artificial material immediately but, if left open (as in Fig. ), develop towards a mechanically more stable funnel shape. No groundwater can be found in the sinkholes, as the groundwater table nowadays is lower, at about 5–7 m depth, when compared to previous decades (1–3 m depth). Fluctuations due to rainfall are in the sub-metre range. The sinkholes are aligned along a narrow band about 50–70 m wide and 300 m long, striking in a west–easterly direction.

Figure 2

Sports field in Münsterdorf, view from the south towards north. Areal photo draped onto topography. The red cylinders mark the locations of the sinkholes, which occur along a narrow west–east trending line. The blue cylinders are direct-push-based probe locations (shown along a south–north profile in Fig. ). Digital orthophoto © GeoBasis-DE/LVermGeo SH.

[Figure omitted. See PDF]

Figure 3

Typical sinkhole (2015, east of sports field in the meadow), with a diameter of 2 m, about 3 m visible depth, and the originally vertical slopes already collapsing to form a more stable form (Photo: Georg Kaufmann).

[Figure omitted. See PDF]

The sudden onset of the sinkhole occurrence and the limitation to the narrow band suggest a relation to a lithological feature in the subsurface and a possible inception by a change in hydraulic conditions (e.g. ). Thus, we regard these two aspects as important constraints for the explanation of the collapse sinkholes along the sports field in Münsterdorf: (i) sudden onset in around 2004, (ii) restriction to the narrow west–east band.

3 Geophysical site exploration

Since the sudden onset of sinkholes on and around the sports field in 2004, several geophysical investigations have been carried out (e.g. ), with the aim to (i) find the cause of the sinkholes and (ii) characterise the subsurface beneath the sports field.

3.1 Previous work

From the numerous geophysical methods used (gravity, geoelectrics, electro-magnetics, georadar, seismics, borehole-based spectrally induced polarisation; see Table  for a summary), especially geoelectrics and georadar have been successful in delineating the glacial till from the peri-glacial sand. In almost all profiles (taken in a north–south direction), the high-resistivity peri-glacial sand is thin (1–2 m) in the north and underlain by more conductive glacial till, while the peri-glacial sand becomes thicker (5–8 m) towards the south, and the glacial till seems to be thinner and more isolated in the south. These results can be explained by the geest island, on which Münsterdorf rests, which peters out towards the south (Fig. ). Cored boreholes confirm this set-up locally (e.g. ).

Table 1

Selected geophysical fieldwork carried out around the Münsterdorf sports field.

Method	Year	Institute	Reference
Electrical resistivity sounding	2005	LLUR Kiel
P-wave reflection seismics	2006	geoFact Bonn
Airborne electromagnetics	2006	BGR Hannover
P-wave reflection seismics	2007	LIAG Hannover
S-wave reflection seismics	2007	LIAG Hannover
Electrical resistivity tomography	2009	LIAG Hannover
Electrical resistivity tomography	2014	CAU Kiel
Ground-penetrating radar	2014	CAU Kiel
S-wave reflection seismics	2014	CAU Kiel
Direct-push-based measurements	2014	UFZ Leipzig
Spectrally induced polarisation	2017	TU Berlin
Gravity	2018	FU Berlin
Ground-penetrating radar	2018	FU Berlin
Electrical resistivity tomography	2018	FU Berlin

Figure 4

Sports field in Münsterdorf, view from the west towards east. The red cylinders mark the locations of the sinkholes, which occur along a narrow west–east trending line. The colour-coded cylinders are boreholes from the LLUR. The blocky area underneath the sports field shows mapped high electrical resistivities above 1000 $\mathrm{\Omega }$m from ERT profiles. In the interpretation, the anomaly delineates the dry peri-glacial sand from the more conductive glacial till. Digital orthophoto © GeoBasis-DE/LVermGeo SH.

[Figure omitted. See PDF]

The chalk layer at about 20 m depth has not been identified in the geoelectric and georadar measurements because electrical conductivities between the chalk and wet till are too similar for the former method and because of depth restrictions for the latter method. However, in seismic p-wave profiles, the chalk layer can be seen.

From the geophysical and borehole measurements, we have derived a simplified lithological set-up underneath the sports field:

• three layers – peri-glacial sand, glacial till, soluble chalk;

• peri-glacial sand varying from 2 to 8 m thickness, becoming thicker in the south;

• massive glacial till layers in the north (several metres thick), petering out towards the south;

• chalk at an almost constant depth of 20 m below the surface underneath the sports field.

We discuss this simplified set-up in more detail later, with more sub-divisions of glacial sand and till, based on our results of the direct-push-based methods.

We argue that the soluble chalk at shallow depth is accessible to groundwater and thus has the potential to be dissolved, resulting in karstified structures. However, as stated below in our first hypothesis, there is also the possibility that voids are present in the sand layers.

3.2 Hypothesis

To start our discussion on the origin of the subsurface voids, we propose two different competing scenarios for the development of the collapse sinkholes.

• Non-karstified. The initial void causing the sinkhole collapse has developed above the lower till layer, at about 6–8 m depth in the peri-glacial sand.

  As we know the sinkhole sizes ($w=\mathrm{2}$–3 m diameter, $d=\mathrm{2}$–5 m depth), the initial void should have a volume around $V=\mathit{\pi }{\frac{w}{\mathrm{2}}}^{\mathrm{2}}d\in [\mathrm{6},\mathrm{35}]$ m${}^{\mathrm{3}}$, a substantial initial void, which at this depth should possibly be detectable by geophysical methods. We try to estimate the effect on gravity and ERT (electrical resistivity tomography) surveys below.

  We note that in this scenario the initial void would be located in insoluble unconsolidated rocks; thus the origin of the void cannot be dissolution of material, only either removal of material by erosion or melting of a block of ice. We then have to ask ourselves why such a void is only present along the narrow band on the sports ground, and not much more widespread, as expected in a landscape dominated by glacial features.

• Karstified. The initial void causing the sinkholes developed in the chalk at around 20 m depth by chemical dissolution of the rock. The initial void volume must be similar to that above.

  In we have argued that initial voids can have developed in the top part of the chalk because surface water undersaturated with respect to calcium is diverted into deeper zones. The reason for surface water to reach the chalk layer is the lower hydraulic conductivity of the glacial till, which is more prominent underneath the northern part of the sports field. This massive glacial till in the northern section of the sports field acts as a flow barrier, forcing surface water to greater depths. Here, the aggressive water with calcium concentrations below the calcium equilibrium concentration creates (on timescales of 10 000–100 000 years) dissolutionally enlarged voids, having metre-size dimensions. The voids are mechanically meta-stable because the water in the voids with its buoyancy counteracts gravity (which we show later by modelling).

  In this case we speculate that the overconsolidated lower peri-glacial sand layer, which is directly above the initial void, has been loosened and then slowly drops into the karst void, but the layers above the peri-glacial sand were kept in place by the finite strength of one of these layers. Thus the initial void started migrating upwards, into the saturated zone. A drop of the groundwater table through mining (as discussed in ) then removed the buoyant support for the overburden, and the upper layers finally collapse.

3.3 Geophysical forward modelling

In our summary of the geophysical work done over the years, we mentioned the identification of the lithological set-up of the region around the sports field. However, geophysical exploration was less successful in identifying the subsurface voids causing the collapse sinkholes. Only old sinkholes show a faint signal in georadar and (in parts) in geoelectrical measurements.

We therefore present results from forward modelling of gravity and geoelectrical signals induced by subsurface voids in a simplified set-up (Figs.  and ): a top layer of dry peri-glacial sand (0–10 m), followed by a wet layer of peri-glacial sand (10–20 m), and below the soluble chalk ($>\mathrm{20}$ m). We do not include a glacial till layer, as (i) the model depicts the situation close to the transition zone from thicker glacial till to thin, more isolated glacial till, and (ii) the material properties of the glacial till, especially electrical resistivity, are not too different from the wet chalk.

Figure 5

Prediction for Bouguer gravity (a) and electrical resistivity (b) resulting from a given lithological structure (c). The structural model (c) comprises dry and wet peri-glacial sand (brown and red), chalk (blue), and either air-filled (white) or water-filled (grey) voids. For the first hypothesis, we located air-filled voids at around 6 m depth with radii (from left to right) of 1.1, 1.7, and 2.0 m. Density and electrical resistivity values can be found in Table . Gravity predictions are based on a box geometry model, and electrical resistivity predictions are based on a Wenner set-up with 5 m electrode spacing.

[Figure omitted. See PDF]

Figure 6

As Fig. , but voids are water-filled (1000 kg m${}^{-\mathrm{3}}$ and 0.5 $\mathrm{\Omega }$m) and start at 20 m depth.

[Figure omitted. See PDF]

Depending on the hypothesis, voids are either (i) in the peri-glacial sand (air-filled; see Fig. ) or (ii) in the chalk (water-filled; see Fig. ). Material properties are listed in Table .

Table 2

Material properties for different lithological layers beneath the Münsterdorf sports field. Density and electrical resistivity from . Young modulus and Poisson ratio from (sand and till) and (chalk).

Material	Dry sand	Wet sand	Glacial till	Chalk	Air-filled void	Water-filled void
Density $\mathit{\rho }$ (kg m${}^{-\mathrm{3}}$)	1800	1900	2100	2600	1.2	1000
Electrical resistivity ${\mathit{\rho }}_{\mathrm{e}}$ ($\mathrm{\Omega }$m)	1000	200–300	5–200	100–700	$\sim \mathrm{800}$	$\sim \mathrm{0.5}$
Young modulus $E$ (MPa)	10–25	50–81	10–700
Poisson ratio $\mathit{\nu }$ (–)	0.3–0.4	0.4–0.4	0.2–0.5
Friction angle ${\mathit{\alpha }}_{\mathrm{f}}$ (${}^{\circ }$)	40–42	40–42	37
Shear strength $u_{\mathrm{c}}$ (Pa)	300	500–600	70–150

The gravity forward model is based on box models for the lithological layers and spheres for the voids (e.g. ). All density values are referenced against the chalk density; thus the calculated Bouguer gravity values will be negative, as peri-glacial sand and voids have lower densities when compared to the chalk.

For the electrical resistivity model, we first create a 2D forward model with res2Dmod (e.g. ) to calculate apparent resistivities for the given set-up. Then, these apparent resistivities are perturbed by random noise (2 %) and then loaded into res2Dinv (e.g. ) to find the resulting resistivity model from inverse modelling.

Following our two hypotheses for the void origin, the non-karstified and the karstified case, we define two forward models with voids at different depths: (i) shallow origin (6–8 m) in insoluble dry sand – air-filled voids, (ii) deeper origin (20–22 m) in soluble chalk – water-filled voids. We locate three voids in these two different depth ranges, with increasing radius from left to right (1.1, 1.7, and 2.0 m) to account for the uncertainty in initial void volume introduced before.

• Non-karstified. In the non-karstified set-up (Fig. ), the voids are located in the dry peri-glacial sand. The Bouguer gravity is negative because the density of the peri-glacial sand (1800–19 000 kg m${}^{-\mathrm{3}}$) is lower than the density of the chalk (2600 kg m${}^{-\mathrm{3}}$), which we took as a reference density. The air-filled voids are modelled with a density of 0 kg m${}^{-\mathrm{3}}$. All three initial voids are visible in the Bouguer signal, but the amplitude (less than 1 mGal) would be hard to detect within the accuracy of a typical relative gravimeter.

  In the ERT cross section, we used electrical resistivities of 1000 $\mathrm{\Omega }$m for the dry peri-glacial sand and 800 $\mathrm{\Omega }$m for the air-filled voids. Below, both the wet peri-glacial sand and the chalk are more conductive, with values around 200 $\mathrm{\Omega }$m. The resulting resistivity model can distinguish between dry sand and wet sand/chalk but does not reveal the void spaces, as their resistivity is too close to the one of the dry sand.

• Karstified. In the karstified set-up (Fig. ), we moved the voids down to the top of the chalk, with a water-filled void density of 1000 kg m${}^{-\mathrm{3}}$ and a void resistivity of $\mathrm{0.5}$ $\mathrm{\Omega }$m, resembling water with dissolved ions. The Bouguer gravity signal now is rather uniform, reflecting the mass deficit of the peri-glacial sand, but no clear indication of the voids. In the ERT cross section, resistivities again identify the interface between dry and wet sand, but the low-resistivity voids within the low-resistivity chalk cannot be seen.

From this simple forward-modelling exercise, we would expect to find hints for initial voids from geophysical signals (mainly gravity), if these initial voids would be in the shallow peri-glacial sand. For a deeper-seated origin of the voids, e.g. in the chalk, no clear geophysical signal in either gravity or geoelectrics can be expected. As outlined in , we have not found evidence for subsurface voids in the geophysical data.

4 Direct-push-based methods

In this section, we discuss results from direct-push-based methods applied to several locations on the sports field, which we use to characterise the hydrogeological regime and relevant mechanical parameter values for the site. The results from the direct-push-based methods provide high-resolution data for geomechanical and hydrogeological properties of the subsurface.

These new field data will complement the existing geophysical data set used to derive a lithological model for the Münsterdorf sports field and its vicinity. For the numerical modelling discussed later, the results from the direct-push-based methods are mandatory to calibrate the models.

The area around the sports field has been probed extensively by several direct-push-based methods (e.g. ). These versatile in situ measurements can be used in unconsolidated and weakly consolidated sediments to obtain vertical high-resolution profiles of a number of material properties, e.g. mechanical failure criteria, hydraulic parameter values, and electrical conductivities, as we outline below. From the measured parameter values listed above, several hydrological and mechanical material properties can be derived.

• i.

  HPT (hydraulic profiling tool). Water is injected with a flow rate $Q$ (m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$) from the direct-push-based probe into the surrounding unconsolidated rock, and the back pressure $P$ (Pa) is measured and used to delineate hydro-stratigraphic units within the subsurface. The ratio of flow rate and back pressure (corrected for air pressure and hydro-stratigraphic pressure) can be used as a proxy that can be correlated to hydraulic conductivity, $K$ (m s${}^{-\mathrm{1}}$), of the material probed (e.g. ). Relative HPT hydraulic conductivity, referring to the ratio of HPT injection rate $Q$ and back pressure $P$, can be correlated with hydraulic conductivity $\mathrm{\Delta }{K}_{\mathrm{h}}$ (m s${}^{-\mathrm{1}}$; e.g. ):1$$\mathrm{\Delta }{K}_{\mathrm{h}}=\mathrm{21.14}\ln \left({\displaystyle \frac{Q}{P}}\right)-\mathrm{41.71}.$$Note that $Q$ is in millilitres, and $P$ is in pounds per square inch in this ratio.

• ii.

  EC (electrical conductivity). By injecting an electric current $I$ (A) from the probe into the ground and measuring the resulting electrical potential $U$ (V), electrical conductivity ${\mathit{\sigma }}_{\mathrm{e}}$ (Sm) can be calculated (e.g. ).

• iii.

  CPT (cone penetration testing). A cone can be pushed at constant speed into the ground, and the force at the tip of the cone can be measured and normalised to the cone surface area to obtain the cone pressure $q_{\mathrm{c}}$ (Pa). Additionally, the force acting along the side of the cone can be measured and normalised by the so-called sleeve area to obtain the sleeve friction $f_{\mathrm{s}}$ (Pa) . To estimate the failure of the material under stress, both the friction angle ${\mathit{\alpha }}_{\mathrm{f}}$ and the shear strength ${\mathit{\mu }}_{\mathrm{0}}$ can be derived. The friction angle is related to the cone pressure $q_{\mathrm{c}}$ :2$${\mathit{\alpha }}_{\mathrm{f}}=f\left(q_{\mathrm{c}}\right),$$with the function $f$ tabulated in geotechnical manuals (e.g. ).

  The undrained shear strength $u_{\mathrm{c}}$ (Pa) can be calculated from the cone pressure $q_{\mathrm{c}}$ and the overburden stress $q_{v\mathrm{0}}$ (Pa) as 3$$u_{\mathrm{c}}={\displaystyle \frac{q_{\mathrm{c}}-q_{v\mathrm{0}}}{N_k}},$$with $N_k$ an empirical cone factor, around $N_k\sim \mathrm{17}$ for overconsolidated clays (e.g. ) and $N_k\sim \mathrm{11}$–19 for normally consolidated marine clays (e.g. ). The overburden stress can be approximated as $q_{v\mathrm{0}}=\mathit{\rho }gh$, with $\mathit{\rho }$ (kg m${}^{-\mathrm{3}}$) the average density, $g$ (m s${}^{-\mathrm{2}}$) gravitational acceleration, and $h$ (m) the thickness of the overburden.

  With the two quantities friction angle and undrained shear strength, the shear stress ${\mathit{\sigma }}_{\mathrm{s}}$ (Pa) as a function of normal stress ${\mathit{\sigma }}_{\mathrm{n}}$ (Pa) can be defined as4$${\mathit{\sigma }}_{\mathrm{s}}={\mathit{\mu }}_{\mathrm{0}}+{\mathit{\sigma }}_{\mathrm{n}}\tan \left({\mathit{\alpha }}_{\mathrm{f}}\right),$$which can be used to evaluate the likelihood that the material will break under certain applied stresses. Here, we assume ${\mathit{\mu }}_{\mathrm{0}}\simeq u_{\mathrm{c}}$. From cone pressure $q_{\mathrm{c}}$ and sleeve friction $f_{\mathrm{s}}$, the soil-behaviour type (SBT) index can be calculated following :5$$\mathrm{SBT}=\sqrt{{\left[\mathrm{3.47}-\log \left({\displaystyle \frac{q_{\mathrm{c}}}{p_{\mathrm{a}}}}\right)\right]}^{\mathrm{2}}+{\left[\log R_{\mathrm{f}}+\mathrm{1.22}\right]}^{\mathrm{2}}},$$with our cone pressure $q_{\mathrm{c}}$ defined earlier, $p_{\mathrm{a}}$ (Pa) the atmospheric pressure, $R_{\mathrm{f}}=\frac{f_{\mathrm{s}}}{q_{\mathrm{c}}}$ the friction ratio, and $f_{\mathrm{s}}$ (Pa) the sleeve friction defined earlier. The SBT index characterises soil types, as it maps the $q_{\mathrm{c}}$–$f_{\mathrm{s}}$ parameter space, and thus the position in the Mohr circle.

• iv.

  SMP (soil moisture probe). Electrical permittivity ${\mathit{ϵ}}_{\mathrm{r}}$ (–) can be measured in the frequency domain to derive the soil moisture content $\mathrm{\Theta }$ (–), which for the saturated zone is equal to the porosity $\mathrm{\Phi }$ (–; e.g. ). From the dielectric permittivity ${\mathit{ϵ}}_{\mathrm{r}}$ measured, we can derive the volumetric water content $\mathrm{\Theta }$ (–) from an empirically fitted relation (e.g. ):6$$\begin{array}{rl}\mathit{\theta }& =-\mathrm{5.3}\times {\mathrm{10}}^{-\mathrm{2}}+\mathrm{2.92}\times {\mathrm{10}}^{\mathrm{2}}{\mathit{ϵ}}_{\mathrm{r}}\\ & -\mathrm{5.50}\times {\mathrm{10}}^{\mathrm{4}}{\mathit{ϵ}}_{\mathrm{r}}^{\mathrm{2}}+\mathrm{4.3}\times {\mathrm{10}}^{-\mathrm{6}}{\mathit{ϵ}}_{\mathrm{r}}^{\mathrm{3}}.\end{array}$$In the saturated zone, the volumetric water content is equal to the porosity $\mathrm{\Phi }$ (–): $\mathrm{\Theta }=\mathrm{\Phi }$.

All direct-push-based methods listed above have been carried out along a south–north transect on the western edge of the sports field; thus the locations are aligned perpendicular to the sinkhole zone (for locations see Fig. ).

The following lithological units can be identified from the combined interpretation of the direct-push-based methods (Fig. ): (i) well-sorted sand and silty sand (peri-glacial sand), (ii) poorly sorted sand/till layers (glacial till), (iii) dissolved chalk along the maximum depth reached by the direct-push measurements. Thus, the direct-push-inferred set-up provides more details when compared to the simpler lithological structure derived earlier from other geophysical methods. The glacial till can be further sub-divided into an upper discontinuous layer (3–5 m depth) and a lower, more continuous till layer (10–12 m depth). The lower till layer can be characterised as aquitard based on the results from the HPT measurements.

Figure 7

Results from direct-push-based methods, plotted along a south–north transect. Direct-push-based locations are labelled on top (EC_x,CPT_x, with $x$ the core number). Borehole transects show the calculated soil-behaviour type (SBT; legend on top) and additionally electrical conductivity (EC; solid lines). The colours in the background represent the structural model interpolated from the direct-push-based results (orange: peri-glacial sand; green: glacial till; blue: chalk).

[Figure omitted. See PDF]

It is interesting to note that the direct-push-based profile performed at investigation point 5 (see Fig. ) reveals an isolated depth interval of increased relative hydraulic conductivity within the chalk. This isolated interval is also characterised by increased values of measured electrical conductivity. The combination of these two elevated parameter values can be interpreted as a karstified zone, having a high hydraulic conductivity and water enriched with dissolved species from the chalk, e.g. calcium and bicarbonate. This is also supported by a CPT log that was performed through an artificially filled sinkhole at approximately 10 m distance from the HPT probing location 5. This CPT log, which was in contrast to the other CPT measurements able to reach into the chalk layer, shows an area with very low cone resistance (indicating very weak and/or loose material) between 22–27 m below the ground surface.

Mechanical properties have been estimated from the CPT measurements (see also Table ). While the glacial till seems to have a friction angle around ${\mathit{\alpha }}_{\mathrm{f}}\simeq \mathrm{37}{}^{\circ }$ and an undrained shear strength around $u_{\mathrm{c}}\simeq \mathrm{70}$–150 Pa, the peri-glacial sand is stronger, with a friction angle around ${\mathit{\alpha }}_{\mathrm{f}}\simeq \mathrm{40}$–42${}^{\circ }$ and an undrained shear strength of $u_{\mathrm{c}}\simeq \mathrm{300}$ Pa in the upper part and up to $u_{\mathrm{c}}\simeq \mathrm{500}$–600 Pa in the lower part of the cross section. It seems that the lower package of peri-glacial sand is overconsolidated, a result from the high normal stress induced during the ice-covered phases (Elsterian and Saalian). The peri-glacial sand above the lower till layer is also strong, except in the zone where the sinkholes occur (friction angles around ${\mathit{\alpha }}_{\mathrm{f}}\simeq \mathrm{35}{}^{\circ }$).

From the results of the direct-push-based methods, with the almost continuous lower till layer being (at least in parts) an aquitard, we continue discussing mechanical models for the sinkhole collapse to elucidate which of our two hypotheses is more likely.

5 Mechanical models

In this section, we introduce and discuss mechanical concepts and models, which we use to describe the sudden occurrence of collapse sinkholes on the sports field. We start with simple analytical stability considerations, which serve as a base for more complicated discrete-element mechanical deformation models.

5.1 Simple mechanical models

We first discuss two simple end-members of a situation prone to collapse sinkhole generation (e.g. ): a first case considering buoyancy only and a second case assuming that some of the stresses are supported by a stronger elastic material.

5.1.1 Buoyancy-driven model

We define a two-layer set-up (Fig.  top), with consolidated sediments of thickness $H$ (m) on top and a density of ${\mathit{\rho }}_{\mathrm{s}}$ (kg m${}^{-\mathrm{3}}$), underlain by soluble rock of density ${\mathit{\rho }}_{\mathrm{c}}$. Below the depth $H$, a circular void of radius $a$ (m) is present, created by dissolution of the rock. The void is either water-filled, with the density of water ${\mathit{\rho }}_{\mathrm{w}}\simeq \mathrm{1000}$ kg m${}^{-\mathrm{3}}$, or air-filled, with a density of ${\mathit{\rho }}_{\mathrm{a}}\simeq \mathrm{0}$. Thus, in this set-up, a vertical cylinder (thin dashed line) with cross section $A=\mathit{\pi }{a}^{\mathrm{2}}$ (m${}^{\mathrm{2}}$) and volume $V=\mathit{\pi }{a}^{\mathrm{2}}H$ (m${}^{\mathrm{3}}$) hangs above the void and is kept in place by buoyancy and friction along its side walls.

Figure 8

Simple mechanical sinkhole set-up. (a) Buoyancy-driven sinkhole, void as thick dashed line, sediment pile on top as thin dotted line, forces as black arrows, water table as blue line. $H$ is the sediment thickness. (b) Flexure-driven sinkhole, parameter values as above; $\mathrm{2}h$ is the thickness of the elastic plate.

[Figure omitted. See PDF]

We then can derive a force balance: 7$${\mathit{F}}_{\mathrm{G}}-{\mathit{F}}_{\mathrm{B}}-{\mathit{F}}_{\mathrm{R}}=\mathrm{0},$$ with ${\mathit{F}}_{\mathrm{G}}$ and ${\mathit{F}}_{\mathrm{B}}$ (N) the gravitational and the buoyancy forces, acting in the normal direction, and ${\mathit{F}}_{\mathrm{R}}$ the frictional force, acting on the side of the column. Considering only the vertical component of the force balance (e.g. ${\mathit{F}}_{\mathrm{G}}=(\mathrm{0},\mathrm{0},F_{\mathrm{G}}{)}^T$) and solving for the frictional force, we obtain 8$$F_{\mathrm{R}}=F_{\mathrm{G}}-F_{\mathrm{B}}={\mathit{\rho }}_{\mathrm{s}}gV-{\mathit{\rho }}_{\text{w/a}}gV=({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\text{w/a}})g\mathit{\pi }{a}^{\mathrm{2}}H.$$ Note that the variable ${\mathit{\rho }}_{\text{w/a}}$ means either water or air density. Rewriting Eq. () to a shear stress by using ${\mathit{\tau }}_{\mathrm{R}}=\frac{F_{\mathrm{R}}}{A_{\mathrm{R}}}$ (Pa), with $A_{\mathrm{R}}=\mathrm{2}\mathit{\pi }aH$ the side hull of the cylinder, we obtain 9$${\mathit{\tau }}_{\mathrm{R}}={\displaystyle \frac{({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\text{w/a}})ga}{\mathrm{2}}}.$$

We now define a standard model with parameter values similar to the Münsterdorf sinkholes: sinkhole radius $a=\mathrm{1.5}$ m, depth of unconsolidated rock column above void $H=\mathrm{20}$ m, and gravitational acceleration $g=\mathrm{9.81}$ m s${}^{-\mathrm{2}}$. For a sediment density of ${\mathit{\rho }}_{\mathrm{s}}\simeq \mathrm{2000}$ kg m${}^{-\mathrm{3}}$, a water density of ${\mathit{\rho }}_{\mathrm{w}}=\mathrm{1000}$ kg m${}^{-\mathrm{3}}$, and air density of ${\mathit{\rho }}_{\mathrm{a}}=\mathrm{0}$, as shear stress we find for a water-filled and an air-filled void, respectively, $$\begin{array}{cc}{\displaystyle }{\mathit{\tau }}_{\mathrm{R}}\left(\text{water-filled}\right)\sim & {\displaystyle }\mathrm{7.4}\mathrm{kPa},\\ {\displaystyle }{\mathit{\tau }}_{\mathrm{R}}\left(\text{air-filled}\right)\sim & {\displaystyle }\mathrm{15}\mathrm{kPa}.\end{array}$$ This result exemplifies the minimum shear stress needed to support the hanging column in both cases, water-filled and air-filled cavity. If the measured shear stress is above the maximum value of 15 kPa, the hanging column remains stable, and the void is kept open. For shear stresses below this threshold, the column drops into the void, and a collapse sinkhole forms.

5.1.2 Flexure- and buoyancy-driven model

For the second example (Fig.  bottom), we consider the same set-up as above, but this time the vertical cylinder is held in place by a thin ($\mathrm{2}h$) elastic layer. From the force balance, we can derive the classical equation for a thin elastic plate under surface loading $q$ (e.g. ) with isotropic and homogeneous material properties:

10$$D{\mathrm{\nabla }}^{\mathrm{4}}w=-q,$$ with $D$ (Nm) the flexural rigidity, $w$ (m) the vertical deformation, and $q$ (Pa) the load acting on the plate (in the vertical direction). For the flexural rigidity, we use (e.g. ) 11$$D={\displaystyle \frac{(\mathrm{2}h{)}^{\mathrm{3}}E}{\mathrm{12}(\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}})}},$$ with $E$ (Pa) the Young modulus, $\mathit{\nu }$ (–) the Poisson ratio, and $\mathrm{2}h$ (m) the plate thickness. Considering a thin circular plate of radius $a$ (m), which is fixed around its perimeter, an analytical solution of Eq. () can be derived within the framework of Love–Kirchhoff plate theory (e.g. ): 12$$w\left(r\right)=-{\displaystyle \frac{q}{\mathrm{64}D}}{\left(a^{\mathrm{2}}-r^{\mathrm{2}}\right)}^{\mathrm{2}}.$$ Defining the load as the sum of the vertical sediment cylinder ${\mathit{\rho }}_{\mathrm{s}}gH$ and the buoyancy induced by the water-filled void ${\mathit{\rho }}_{\mathrm{w}}gH$ and assuming wall friction ${\mathit{\tau }}_{\mathrm{R}}$ along the side of the column, we arrive at 13$$\begin{array}{rl}{q}_{\text{full-load}}& ={\mathit{\rho }}_{\mathrm{s}}gH-{\mathit{\rho }}_{\text{w/a}}gH-{\mathit{\tau }}_{\mathrm{R}}\\ & =({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\text{w/a}})gH-{\mathit{\tau }}_{\mathrm{R}}.\end{array}$$ Inserting Eq. () into Eq. (), we can now calculate the deflection of the supporting layer in the centre ($r=\mathrm{0}$). For a Young modulus of $E\simeq \mathrm{10}$ GPa, a Poisson ratio of $\mathit{\nu }\simeq \mathrm{0.3}$, and a plate thickness of $\mathrm{2}h=\mathrm{0.6}$ m, the deflection for a water-filled and for an air-filled void, respectively, is $$\begin{array}{c}{\displaystyle }w\left(\text{water-filled}\right)\sim \mathrm{0.08}/\mathrm{0.06}\mathrm{m},\\ {\displaystyle }w\left(\text{air-filled}\right)\sim \mathrm{0.16}/\mathrm{0.14}\mathrm{m},\end{array}$$ with the left values assuming no wall friction (${\mathit{\tau }}_{\mathrm{R}}=\mathrm{0}$) and the right value assuming a wall friction of ${\mathit{\tau }}_{\mathrm{R}}=\mathrm{50}$ kPa (a higher value than the minimum shear stress derived before).

However, the more important parameter is the thickness of the elastic layer $h$ as, according to Eq. (), it scales with the power of 3.

We now want to estimate when the elastic layer will break. We therefore need to calculate the stresses in the thin elastic plate: $$\begin{array}{cc}\text{14}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{rr}=-{\displaystyle \frac{\mathrm{3}qz}{\mathrm{32}{h}^{\mathrm{3}}}}\left[(\mathrm{1}+\mathit{\nu })a^{\mathrm{2}}-(\mathrm{3}+\mathit{\nu })r^{\mathrm{2}}\right],{\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathit{\theta }\mathit{\theta }}=-{\displaystyle \frac{\mathrm{3}qz}{\mathrm{32}{h}^{\mathrm{3}}}}\left[(\mathrm{1}+\mathit{\nu })a^{\mathrm{2}}-(\mathrm{1}+\mathrm{3}\mathit{\nu })r^{\mathrm{2}}\right],\\ {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{r}\mathit{\theta }}=\mathrm{0},\end{array}$$ with $z$ (m) the vertical coordinate, $z\in [-h,h]$. As the shear stress ${\mathit{\tau }}_{\mathrm{r}\mathit{\theta }}$ is zero, the two normal stresses ${\mathit{\sigma }}_{rr}$ and ${\mathit{\sigma }}_{\mathit{\theta }\mathit{\theta }}$ are the two principal stresses, which we term ${\mathit{\sigma }}_{\mathrm{1}}={\mathit{\sigma }}_{rr}$ and ${\mathit{\sigma }}_{\mathrm{2}}={\mathit{\sigma }}_{\mathit{\theta }\mathit{\theta }}$. The maximum values for these two normal stresses are obtained for $r=a$ and $z=h$: $$\begin{array}{c}{\displaystyle }{\mathit{\sigma }}_{rr}^{\max }=+{\displaystyle \frac{\mathrm{3}q{a}^{\mathrm{2}}}{\mathrm{16}{h}^{\mathrm{2}}}},\\ \text{15}& {\displaystyle }{\mathit{\sigma }}_{\mathit{\theta }\mathit{\theta }}^{\max }=+{\displaystyle \frac{\mathrm{3}q{a}^{\mathrm{2}}}{\mathrm{16}{h}^{\mathrm{2}}}}\mathit{\nu }.\end{array}$$

Using the material parameter values from our numerical example above, the normal stresses are $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{rr}^{\max }\left(\text{water-filled}\right)\sim \mathrm{920}\mathrm{kPa},\\ {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathit{\theta }\mathit{\theta }}^{\max }\left(\text{water-filled}\right)\sim \mathrm{276}\mathrm{kPa},\\ {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{rr}^{\max }\left(\text{air-filled}\right)\sim \mathrm{1839}\mathrm{kPa},\\ {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathit{\theta }\mathit{\theta }}^{\max }\left(\text{air-filled}\right)\sim \mathrm{552}\mathrm{kPa}.\end{array}$$ Note that by reducing the plate thickness to 50 %, both normal stresses increase by a factor of 4!

For the two principal stresses, we can calculate the centre $M_{\mathrm{c}}$ and the radius $M_{\mathrm{r}}$ of the Mohr circle as (e.g. ) $$\begin{array}{c}{\displaystyle }{M}_{\mathrm{c}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}}{\mathrm{2}}}=\mathrm{598}/\mathrm{1196}\mathrm{kPa},\\ \text{16}& {\displaystyle }{M}_{\mathrm{r}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}}}{\mathrm{2}}}=\mathrm{322}/\mathrm{644}\mathrm{kPa},\end{array}$$ with the values given for a water-filled and an air-filled void, respectively.

The two Mohr circles are shown in Fig. , along with the lines of failure for the glacial till and the peri-glacial sand. For the failure lines, parameter values have been taken from the results of the direct-push-based methods discussed before. If the column above the void is supported by the buoyancy of water, the Mohr circle is below the line of failure. Removing the buoyant support moves the Mohr circle closer to the failure line. Here both the thickness of the elastic plate, $\mathrm{2}h$, and the Poisson ratio, $\mathit{\nu }$, strongly control the failure: reducing either one or both of these values results in a Mohr circle crossing the failure line.

Figure 9

Mohr–Coulomb circle for stresses induced by column resting on elastic plate, with buoyancy (dark blue) and without buoyancy (light blue). The line of failure is shown for the weaker glacial till (solid red) and the stronger peri-glacial sand (dashed red), with parameter values derived from results of the direct-push-based methods.

[Figure omitted. See PDF]

Thus we have shown with our simple analytical examples that a collapse sinkhole with size and depth typical of the Münsterdorf sinkholes can form if the buoyant support of the material above the void is reduced.

Next, we extend our simple set-up to a three-dimensional mechanical model, based on the discrete-element method.

5.2 Discrete-element mechanical models

The discrete-element method (DEM) is a computational method to simulate the interaction of a large number of particles (e.g. ). Often, the particles are spheres with pre-described radius and some material properties such as density $\mathit{\rho }$ (kg m${}^{-\mathrm{3}}$), Young modulus $E$ (Pa), and Poisson ratio $\mathit{\nu }$ (–). Spheres are packed into a three-dimensional region to represent a macroscopic material, with packing being either regular or irregular. Neighbouring spheres can interact, with a large choice of mechanical laws defined. Modelling of the mechanical deformation consists of several steps: (i) set-up of modelling domain, (ii) packing with spheres, (iii) definition of material properties and interaction laws, (iv) time-stepping procedure to track position and velocity of each sphere.

We use the open-source DEM model YADE (e.g. ). Our modelling domain is 50 $\times$ 50 $\times$ 50 m in size, thus 10 times larger than the collapse sinkholes in Münsterdorf (Fig. ). This large modelling domain ensures that the fixed boundary spheres (also shown in red as the chalk) have no influence on the deformation (see , for more details on the implementation). The material properties for the different layers (upper and lower peri-glacial sand, glacial till, chalk) have been taken from the literature (Young modulus, Poisson ratio) and from the direct-push-based methods discussed above (friction angle and cohesion).

Figure 10

Set-up of YADE mechanical model. The model is 50 $\times$ 50 $\times$ 50 m in size and filled with spheres irregularly packed. The colours of the sphere represent different layers: upper peri-glacial sand (light brown), glacial till (orange), lower peri-glacial sand (light brown), chalk (red), and initial void (blue). Note that the sides of the modelling domain (also shown in red) are fixed.

[Figure omitted. See PDF]

The spheres in the initial model (upper and lower peri-glacial sand: light brown; glacial till: orange; chalk: red) are packed, and the model is stable; thus no deformation is observed (Fig.  top). We then initiate a void in the chalk (red spheres) by removing spheres within a given radius (removed spheres marked in blue). The removed spheres represent either a water- or an air-filled initial void. This initial void in the chalk destabilises the model domain mechanically. The material above the void starts to break down into the void until a new mechanical equilibrium is reached.

Figure 11

Mechanical breakdown simulated with YADE. (a) Initial condition; the dissolved void in the chalk is shown in blue. (b, d) Two snapshots of collapse with full buoyant support. (c, e) Two snapshots of collapse with partial buoyant support. Note that we show a thin slice of the 3D model, crossing the void in the chalk.

[Figure omitted. See PDF]

Figure 12

Evolution of collapse sinkhole. Light-blue colours indicate chalk as soluble rock, yellow colours peri-glacial sand, and brownish colours glacial till, the latter two as insoluble rocks. The dark-blue line is the hypothetical water table, and the white area is a void in the chalk.

[Figure omitted. See PDF]

Two scenarios are shown in Fig.  (middle and bottom row).

• i.

  Buoyant support. In the left column, the subsurface is characterised by a high water table, and the initial void is filled with water. Thus the difference between gravitational force and buoyancy force is smaller, (${\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{w}})g$, and thus the buoyant force compensates the gravitational force significantly. The model experiences breakdown of sand layers above the chalk, but the elastic strength of the till supports the weight of the column above, and the void does not migrate towards the surface. The initial void enlarged by chemical dissolution becomes enlarged by additional breakdown but remains stable. No collapse sinkhole will form (unless the boundary conditions are changed).

• ii.

  No buoyant support. In the right column, we reduce the buoyant force, and the initial void is air-filled, simulating a drop in water table, as possibly initiated by groundwater withdrawal. The difference between gravitational force and buoyancy force becomes larger, (${\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{a}})g$ with ${\mathit{\rho }}_{\mathrm{a}}\simeq \mathrm{0}$, and thus the overlying column experiences a stronger downward force. In this model, the initial void also enlarges through roof collapse, but the overlying layers cannot support the weight anymore, the entire overburden slips into the void, and a collapse sinkhole appears suddenly on the surface.

By comparing these two models we have shown that a simple change in the hydraulic boundary condition, here a drop in water table, which reduces buoyancy and thus increases the weight of the column hanging above the void, can trigger a collapse sinkhole.

We stress that the timescales of dissolution and mechanical breakdown are completely decoupled: while the void in the chalk, created by chemical dissolution of the material on timescales of 1000–100 000 years, might initiate roof collapse into the insoluble formation, the surface remains stable because the overlying material is partially supported by the mechanical strength of the layers. However, a drop in water table through groundwater withdrawal on timescales of years can destabilise the system and initiate a sudden collapse of the surface, creating a collapse sinkhole.

Of course one can argue that a similar situation might occur for a void located in the peri-glacial sand. However, we have no evidence for a mechanically sound layer in the top part of the insoluble overburden (down to 6–8 m); thus a potential subsurface void in the peri-glacial sand is unlikely to be mechanically meta-stable but will collapse more or less immediately. Thus we argue that the void origin is more likely to be in the chalk layer. We support this argument further in the next section.

6 Chemical considerations

In this section, we pick up the short-term mechanical failure models developed in the last section in view of the long-term evolution of the sinkhole cluster on the sports field. With the mechanical models presented before we show that a water table drop can trigger the sinkhole collapse. The initial void, however, can be present in the subsurface for a long period before the collapse.

We want to continue the discussion on the origin of the voids in the surface. In Sect. , we discuss the lower glacial till layer, which seems to be a continuous structure at around 12–15 m depth, with possible interruptions in the sinkhole zone. Hints for these interruptions come from the higher hydraulic conductivities and the elevated electrical conductivities in boreholes crossing the sinkhole zone, which might indicate water from the chalk, enriched with dissolved matter (we come back to this argument later).

6.1 Chemistry

From the direct-push-based electrical conductivity measurements, we found values for the electrical conductivity around ${\mathit{\sigma }}_{\mathrm{e},\mathrm{f}}\simeq \mathrm{30}$–50 mS m${}^{-\mathrm{1}}$ where direct-push profiling reached the hydraulically conductive part of the chalk, in parts even higher. We convert this electrical conductivity representative of the fluid, ${\mathit{\sigma }}_{\mathrm{e},\mathrm{f}}$ ($\mathrm{\mu }$S cm${}^{-\mathrm{1}}$), to the total amount of dissolved solids, TDS (mg L${}^{-\mathrm{1}}$), using an experimentally determined linear relation (e.g. )

17$$\mathrm{TDS}=k_{\mathrm{e}}{\mathit{\sigma }}_{\mathrm{e},\mathrm{f}},$$ with $k_{\mathrm{e}}\in [\mathrm{0.55},\mathrm{0.80}]$ (cm mg L${}^{-\mathrm{1}}$ $\mathrm{\mu }$S${}^{-\mathrm{1}}$) a correlation factor and the electrical conductivity ${\mathit{\sigma }}_{\mathrm{e},\mathrm{f}}$ ($\mathrm{\mu }$S cm${}^{-\mathrm{1}}$) in this equation referenced to 25 ${}^{\circ }$C. Using $k_{\mathrm{e}}=\mathrm{0.6}$, we obtain a TDS value around 180–300 mg L${}^{-\mathrm{1}}$.

When we argue that a large amount of the TDS value can be explained with dissolution of the chalk, we derive the concentration of calcium $c$ (mol m${}^{-\mathrm{3}}$) from the TDS value, using 18$$c={\displaystyle \frac{\mathrm{TDS}}{m_{\mathrm{r}}}},$$ with $m_{\mathrm{r}}=\mathrm{0.101}$ kg mol${}^{-\mathrm{1}}$ the atomic mass of calcite, resulting in $c\sim \mathrm{2}$–3 mol m${}^{-\mathrm{3}}$.

If we compare this calcium concentration $c$ derived from TDS to the calcium equilibrium concentration $c_{\mathrm{eq}}$ (mol m${}^{-\mathrm{3}}$) for calcite dissolved under closed-system conditions (decoupled from the atmosphere), which is in the range of $c_{\mathrm{eq}}\simeq \mathrm{2}$–5 mol m${}^{-\mathrm{3}}$ (e.g. ), we find that the dissolution of chalk is close to or at its equilibrium value, thus at its maximum.

Arguing that the high calcium concentration in the probed water in the chalk has been attained locally, the EC results from direct-push-based methods support the idea of that enlargement of voids in chalk by chemical dissolution has created the initial void space at around 20 m depth. The isolated area of high permeability within the chalk that was identified using the direct-push-based HPT as well as the depth interval with strongly reduced cone pressure (CPT measurement) supports the theory that the high calcium concentration has been generated by the local dissolution, creating a large secondary permeability in the chalk.

7 Interpretation

In , we have developed a simple hypothetical model for the formation of sinkholes on the sports field in Münsterdorf. With the additional information gained in this work from the direct-push-based methods, the short-term mechanical collapse models, and estimates of the chemical composition of the groundwater, we are now able to refine our previous hypothesis.

The cross sections in Fig.  follow a temporal evolution (from top left via bottom towards top right). Shown are three lithological units: chalk (light blue), peri-glacial sand (yellow), and glacial till (brown). The thick blue line depicts the water table.

From our evolution models, we speculate that groundwater is forced to flow vertically down in the middle of the section because the glacial till present in the northern part is less hydraulically conductive. The groundwater, undersaturated with respect to calcium, reaches the chalk at around 20 m depth and dissolves it, creating voids up to the metre range. This dissolution process, a long-term process, occurs on the 10 000–100 000-year timescale.

The voids developing in the chalk start to migrate upwards through roof breakdown, once they reach a certain size. However, thin glacial till layers in the south, mechanically stronger than the peri-glacial sand, stabilise the situation, as the entire void and collapsed roof part is located in the saturated zone; thus buoyancy supports the strength of the thin glacial till layers in the southern part. This situation is meta-stable, as long as there are no significant changes in groundwater level.

A drop in groundwater level, which can be initiated through pumping in the nearby open-pit mines, removes the buoyant support to the glacial till above a void, and the weight of the overburden becomes too much and triggers a collapse, which results in a sinkhole in the overlying peri-glacial sand.

8 Conclusions

The collapse sinkhole cluster on the sports field of Münsterdorf started in 2004, with an occurrence rate of about one per year. The sinkholes, typical cover-collapse sinkholes, are usually 2–3 m in diameter and 3–5 m deep. We have investigated the mechanical stability of collapse sinkholes with similar dimensions with a discrete-element model to simulate the stability conditions of the final collapse.

In the introduction, we pose three main questions, which we want to answer in this last part.

• i.

  Sudden onset. The sudden onset of collapse sinkhole formation in 2004 points to a change in a boundary condition close to the site, e.g. changes in the mechanical and/or hydrological situation. We are able to explain the sinkhole collapse with a drop in groundwater table, reducing the buoyant force, which stabilises the existing subsurface voids.

• ii.

  Mechanical failure. Our mechanical models of a void in the subsurface, embedded in the typical set-up of the locality (peri-glacial sand, glacial till, and soluble chalk at 20 m depth), indicate stability of the layers above an existing subsurface void when the water table is high.

  A drop in groundwater reduces the buoyant force, which counteracts the gravitational force of the overburden, and the overburden can exceed its mechanical threshold and collapse into the void below, creating a cover-collapse sinkhole.

• iii.

  Void origin. We have speculated about a void origin either shallow in the peri-glacial sand (6–9 m depth) or deeper in the chalk (20–22 m depth). We have shown that the shallow void might be detectable by gravity measurements, but a deeper void is not really detectable with geophysical measurements.

  The high electrical conductivity, together with an elevated hydraulic permeability, which we measured with direct-push-based methods at about 22.7–25 m depth, points to dissolution on top of the chalk layer, thus pointing to void spaces developing at around or below 20 m depth.

Code availability

Simple calculations have been performed with the open-source programming language Python (https://www.python.org, last access: 27 February 2023). The discrete-element model has been run with the open-source YADE package (https://yade-dem.org, last access: 27 February 2023, ). Geoelectric inversions are based on the commercial Res2DInv package (https://www.aarhusgeosoftware.dk/res2dinv, last access: 27 February 2023, , see also ).

Data availability

The direct-push-based data are the property of UFZ Leipzig (please contact Ulrike Werban or Thomas Vienken for details).

Author contributions

GK designed the set-up of the numerical experiments and participated in writing the manuscript. DR designed and ran the distinct-element model runs and participated in writing the manuscript. UW performed the direct-push-based measurements and participated in writing the manuscript. TV performed the direct-push-based measurements and participated in writing the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We acknowledge funding from the BMBF within the SIMULTAN project. Figures were prepared using GMT software and the PARAVIEW software .

Financial support

This research has been supported by the Bundesministerium für Bildung und Forschung (BMBF) within the SIMULTAN project under research grants 03G0843G (Georg Kaufmann, Douchko Romanov) and 03G0843F (Ulrike Werban, Thomas Vienken).We acknowledge support from the Open Access Publication Initiative of Freie Universität Berlin.

Review statement

This paper was edited by Kei Ogata and reviewed by two anonymous referees.