http://www.press.ierek.com

ISSN (Print: 2357-0849, online: 2357-0857)

International Journal on:

Environmental Science and Sustainable Development

DOI: 10.21625/essd.v4i2.558

The Sinkhole Occurrence Risk Mitigation in Urban Areas for the
Historic Salt Mine

Agnieszka Malinowska1, Ryszard Hejmanowski1, Artur Guzy1, Andrzej Kwinta2, Paweł
Ulmaniec3

1Faculty of Mining Surveying and Environmental Engineering, University of Science and Technology, Cracow,

Poland
2Agriculture University, Cracow, Poland
3Wieliczka Salt Mine, Wieliczka, Poland

Abstract
The present research focuses on the definition of a novel methodology for sinkhole risk assessment above shallow
salt mines. The research were carried out on the area above the Wieliczka salt mine, a World Heritage site. The
study of vertical stresses on the basis of a theoretical state of rock mass deformation in the area of test chambers
was performed. Furthermore, the risk of chamber collapse due to ventricular stress exceeding the limit specified
in the zone were calculated based on the arch pressure theory. The final stage of the research consists of spatial
analysis that leading to the identification of chambers potentially influenced by other risk factors. The research
shown in the article strongly suggests that combined spatial analysis with geotechnical analysis may lead to reliable
sinkhole risk assessment methodology.

© 2019 The Authors. Published by IEREK press. This is an open access article under the CC BY license
(https://creativecommons.org/licenses/by/4.0/). Peer-review under responsibility of ESSD’s International Scien-
tific Committee of Reviewers.

Keywords

salt deposit; sinkhole; risk management; convergence; GIS

1. Introduction

1.1. Background

In 1960 a sinkhole of ca. 100 m diameter was formed as a consequence of a goaf in the Schmidt salt chamber at
the Wieliczka salt mine. Within 15 minutes, the terrain had subsided and as a result, two houses collapsed. Similar
events have occurred in other parts of the world, frequently causing serious damage. The problem of sinkhole-
prone chambers is extremely complex and, despite many years of re-search, is not yet fully understood. It is known
that old and shallow workings have to be backfilled. Moreover, some mines such as the one found in Wieliczka,
are UNESCO World Heritage site, and therefore their preservation and maintenance for fu-ture generations is a
priority.

Regardless, safety measures must be implemented to provide assurance for the citizens of Wieliczka, as the oldest
part of the mine rests under the town center. The course of the deformation process in the saline rock mass is

pg. 85

https://creativecommons.org/licenses/by/4.0/


Malinowska / Environmental Science and Sustainable Development, ESSD

influenced by the geomechanical properties of various types of salt present in the Wieliczka salt mine. It is also
caused by a number of mining and geological factors which lead to sink-hole-formation (Çanakcı, Güllü, 2009;
Kleczkowski, 1993; Przybyło, 1980). The studies presented in this paper focus on determining the most important
risk factors responsible for the formation of weaker zones in the rock mass and, consequently, which en-courage the
formation of caving collapse process. In the past, most evaluations on the risk of sinkhole formation concentrated
on empirical solutions or numerical geomechanical calculations (Ryncarz, 1993). However, those methods turned
out be unreliable for a number of reasons.

It is the critical evaluations of previous methods of estimating the sinkhole risk that have led to the concepts
presented in this paper. Widely agreed upon approaches that the sinkhole risk evaluation methods should be
principally based on a spatial analysis en-compassing the present technical state of the chambers, mutual bedding,
applied exploitation systems and interrelations between the chambers.

1.2. Geology of the study area

Wieliczka is located in the south-east of Poland. The saliferous formation is as a narrow strip covering an area
of 7 km to the east-west and 0.9 km to the north-south. The depth of the extracted deposit reaches 370 meters.
Underground mining activity in this region first took place in the 16th century. Salt mining has been carried out
in about 2000 chambers on 9 levels, an area which amounts to 7.5x106 m3. The most intensely developed part of
the town itself, which has a population of 20,000, is located directly above the oldest and shallowest part of the
mine. Because of complex and diverse geological conditions, various salt extraction methods have had to be used
in this region. The salt deposit is divided into two strata. The Lower strata are strongly inclined with a 15◦-70◦

inclination, and are folded. The upper layer consists of loams and salt solids and has been strongly dispersed by
tectonic movements. Five types of salt have been recorded in the analyzed area (Kleczkowski, 1993; Przybylo,
1980):

1. Beds of green salt,

2. Solids of green salt,

3. Shaft-salt,

4. “spiz”-salt,

5. Solids of “zuber”-salt located in the loam bed.

Various exploitation methods were applied with regard to the salt-type, as well as to the geological and hydro-
geological conditions. The salt strata were isolated from water infiltration by a marl and loam bed. Above the
loam and marl bed sandy loams formed, thanks to which water infiltration was able to take place. In those areas,
quicksand beds formed. The northern part of the salt was strongly distributed and cracked, thereby allowing water
infiltration to the chambers in this area was possible.

1.3. Sinkhole risk mitigation in urban areas

At the first stage of the investigation was to design the architecture of functional database. Currently, Geographical
It has been repeatedly shown after many years of research and observation that the main factors which can be
significantly attributed to the formation of sinkholes are: (Mancini et al., 2009; Ryncarz, 1993; Saustowicz, 1955;
Whittaker, Reddish, 1989).

– Height of the chambers

pg. 86



Malinowska / Environmental Science and Sustainable Development, ESSD

Global studies and also events in Poland show that the height of the chambers is one of the most important factors
in causing the development of sinkholes. The height of the chambers in the Wieliczka salt mine ranges between 4
and 50 meters.

– Thickness of the roof ledge

The thickness of the roof ledge is defined as the difference of the terrain surface datum and workings’ roof datum.
The thicknesses of the roof ledge in the study area vary depending on the depth (the roof ledge for the chambers
on the first level ranges between 23 and 300 meters).

– Quaternary sediments and strata overlaying loam beds

Quaternary strata can be a factor which increases the probability in the occurrence of a sinkhole through the
possible infiltration of water to the chamber when the continuity of the strata has been disrupted. The quaternary
strata in the research area are thick with inserts and water-bearing lenses. Water infiltrates through fractures formed
as a result of rock mass relaxation, and do so along vertical workings.

– Type of chambers

The salt exploitation mode is directly connected with the type of salt and the way it is deposited. Salt exploitation in
the Wieliczka area was carried out employing various methods. Generally the chambers can be divided into those
which are made in a block of salt, those which were made in a salt bed, those which are leached, and finally, holed
chambers. Most sinkhole-prone chambers are the shallowest chambers of irregular shape and varying thickness of
salt cover, insulating against water influxes.

– Distribution of chambers with regards their interrelations

The interrelation of the chambers is the most complex issue, requiring advanced computer applications for analyz-
ing spatial interrelations between the chambers. The most hazardous are the complexes of chambers on neighboring
horizons, especially when the size of the chambers is considerable.

The afore mentioned factors will be the criterion of the final selection of sinkhole-prone chambers.3.1. Risk map
development.

2. Methods for assessing the risk of sinkhole occurrence

2.1. Assumptions and general scheme of research

At the beginning of the study was obvious that the complicated geological conditions in the rock mass and his-
torical excavations do not allow the use of any classical numerical calculation methodology to assess the sinkhole
occurrence. The goal of the research was to combine the empirical approach for subsidence prediction with the
geomechanical method for strain calculation and with geographical information system (GIS). Those general as-
sumption allowed to create an own, original method for sinkhole occurrence above the historical salt mine.

The present research result consists of four stages (Fig.1).

The first stage consists of the rock mass deformation analysis in the area of the shallowest salt chambers using
Knothe’s generalized geometrical-integral model (Knothe, 1953)(Fig.1a). The results of those calculations made
it possible to find vertical movement s and strains ε x, ε y, ε z caused by the convergence of chamber workings in
the deeper horizons of the mine.

pg. 87



Malinowska / Environmental Science and Sustainable Development, ESSD

Figure 1. Scheme of research methodology

The second stage of the research facilitated the estimation of deformations and vertical stresses in the area of walls
in shallow chambers (Fig. 1b). With those calculations it was possible to separate chambers for which vertical
stresses were similar or exceed the average tensile stress and compressive stress values, and those which were
regarded as critical. Depending on the risk of the exceeding critical values, categories of differing risk have been
determined.

The third stage of the research evaluates the range of the cracking zone above chambers with pressure arch theory
(Fig.1c). Depending on the depth which could be reached by the cracking zone, four sinkhole-prone zones were
found.

The forth step was the GIS integration of the analytical results with risk factors which could not be analyzed due
to their qualitative characteristics (Fig. 1 d). Finally, the chambers which generate a potential risk of sinkhole
occurrence on the surface have been selected.

2.2. Predicted deformations in the vertical plane

Deformations in the vertical plane were predicted using Knothe’s theory. According to this theory the exploitation
of an elementary volume of deposit, dV and horizontal plane, dP generate an elementary subsidence at point A.
This subsidence can be described as:

dsA = f (x,y)dP (1)

where:

-dsis the elementary subsidence at point A,

-f(x,y) is the influence function,

-dP is the area of an element of the chamber.

Employing the principle of superposition it has been assumed that vertical displacement at point A, as denoted by
sA is the sum of elementary subsidence coming from all elementary volumes in the exploitation field area:

SA =
∫∫

p f (x,y)dP (2)

pg. 88



Malinowska / Environmental Science and Sustainable Development, ESSD

Assuming various forms of the influence function f(x,y), we can define various prediction models for displacements
and deformations. In Knothe’s theory (Knothe 1953), the following generalized influence function for a flat state
of displacements is assumed:

f (x) = smax
h
√

2π
ex p(−h2x2) (3)

where:

-smax= a g is the maximum final subsidence [m],

-a is the exploitation coefficient, depending on the way the post-exploitation void has been filled out,

-g is the thickness of exploitation,

-h is the parameter of the influences dispersion.

Following a parametrization of the influence function, i.e. assuming the limited range of the influence of the
exploitation and introduction of parameter R, the function (3) can assume the form (4):

f (x) = smax
1
r

ex p(−π
x2

R2
) (4)

The parameter R is the radius of the influences dispersion, also called the range of the main influence. This
parameter is connected with the strength characteristic of the rock mass, through the dependence (5):

tgβ =
H
R

(5)

where:

-H is the depth of exploitation [m],

-tgβ is the parameter of rock mass corresponding to strength properties of rocks.

Knothe’s classic theory only enables us to predict deformations in a vertical plane. For the sake of practicality, the
theory had to be supplemented by theoretical calculations in a horizontal plane. Studies by Budryk (Budryk, 1953)
turned out to be crucial for predicting horizontal displacements and strains in Poland.

For the sake of the presented method, it was necessary to determine the state of strains at points of the salt panel
sidewalls. Those strains for all panels were calculated from the relation (6):

εx =
∂ ux
∂ x

, εy =
∂ uy
∂ y

, εz =
∂ uz
∂ z

(6)

where ux, uy, uz=s, are the respective horizontal movements (Budryk, 1953).

Knothe’s principles (see model above) were somewhat generalized for the saline rock mass. Each panel localized
at 9 levels of the Wieliczka salt mine was numerically introduced by its geometrical simplification. In the process
of convergence, specific elements of all panels were generating displacements and deformations at the analyzed
points of the sidewalls in other panels. In this way, it was possible to determine definitely the state of strains for
each sidewall of the salt panel.

2.3. Determining the stresses in vertical plane

In determining the risk of panel stability loss under the influence of deformation of other panels on the basis of
strains determined from Knothe’s general theory, the stresses were defined with the Hook’s general law. On the
basis of deformations and parameters of rocks making up the rock mass, the stresses can be calculated from the
dependence:

σx =
2G

1−2ν
[(1−ν)εx + ν (εy + εz)]

σy =
2G

1−2ν
[(1−ν)εy + ν (εx + εz)]

σz =
2G

1−2ν
[(1−ν)εz + ν (εy + εx)]

(7)

pg. 89



Malinowska / Environmental Science and Sustainable Development, ESSD

where:

-G is the shear modulus,

-v is the Poisson number.

As vertical stresses have the biggest influence on the stability of sidewalls, further analyses were focused only on
those stress values. Due to the fact that all panels stay within the influence of the lower panels, they needed to be
grouped according to the vertical stress. It should be emphasized that these are stresses are exclusively caused by
the deformation of the rock mass coming from the contraction of other panels with respect to those analyzed. This
categorization was introduced for the needs of the saline rock mass analyses presented in this paper. This was the
author’s solution. The following (risk) categories for vertical stresses were proposed:

->Category 1 – deformation stresses within the critical value limits,

->Category 2 – stresses exceed critical values less than twice,

-> Category 3 – stresses exceed critical values twice,

->Category 4 – stresses exceed critical values three times,

->Category 5 – stresses exceed critical values four times.

In total, 283 panels met the criteria assumed for the vertical stress analyses.

The results obtained reveal that a relatively large number of panels could be subjected to the negative impact of
deformations caused by other, deeper-located panels. It should be stressed that the selected panels could have been
back-filled or protected otherwise in view of the expected deformations. This would confirm the assumed stress
analysis methods used in this paper.

2.4. Determining the range of the cracking zone

Determining the feasibility of discontinuous deformation occurrence on the terrain surface over the exploitation
panels is related with defining the range of the direct impact of the cracking zone on the near-surface rock mass
layers. In a given situation, the rock mass in the near-surface zone consists of a few layers of varying geomechan-
ical properties and it is important to indicate which layers will house the cracking zone. Immediately under the
ground surface we have a clayey sand layer with a sliding tendency. Thus, if the cracking zone reaches that layer,
discontinuous deformations in the form of cones will appear on the surface. Then, a clay-gypsum layer follows in
the geological profile. If the cracking zone reaches that layer, its water tightness will be disturbed and a massive
water influx can occur, causing a washing off the sands. In this case, various fractures may be formed.

A further risk analysis of 283 preselected panels lies in using the pressure arch theory. In this way, panels which
are a potential source of discontinuous deformation on the terrain surface could be selected from among the panels
in danger of vertical stress higher than the strength of the rocks.

Taking into account the similar character of data illustrating the geometry of the workings, as well as the build and
geomechanical properties of the rock mass, Sałustowicz’s simplified model was applied (Sałustowicz, 1955). In
this model, also referred to as the pressure arch theory, the assumption was made that the cracking zone takes the
form of an ellipse in the vertical cross-section. The size of the semi-axis of that ellipse depends on the dimension
of the working and the properties of the neighboring rocks. When rocks are made up of a low tensile strength rock
mass, the tensile stresses will have a decisive significance on the range of the cracking zone.

An example scheme presented in Fig. 2 was analyzed.

The depleted panel had the following dimensions:

radius r,

height ge

The semi-axes of the ellipse were denoted as a and b. There were additionally marked with the height coordinates

pg. 90



Malinowska / Environmental Science and Sustainable Development, ESSD

of the terrain surface Zpow , the top of the cracking zone Zsp, the roof of the panel Zst and the floor of the analyzed
panel Z0.

As already mentioned, the prerequisite of cracking zone formation is to exceed the tensile strength of value Rr by
stressing the level of layers making up the panel roof. This condition can be written as follows:

−pz +
(

1 + 2
a
b

)
px ≥ Rr (8)

where pz and px are primary stresses on the top node of the ellipse.

After substituting the dependence of horizontal vs. vertical pressures and determining the boundary value from
(8), the dependence for the semi-axes has been obtained in the form below:

a
b
=

Rr(1−ν)− pz(2ν −1)
2pzν

(9)

The primary vertical pressure at the height of the panel roof was determined from the dependence:

pz =−γ(Zpow −Zst) (10)

Figure 2. Estimating the range of cracking zone over the salt chamber

where:

-g is the average specific weight of the deposited over chamber rock mass.

In line with Fig. 2, the ellipse equation was introduced. In the assumed coordinates system it has the following
form (11): (x

b

)2
+
( z

a

)2
= 1 (10)

After substituting the coordinates of point P (from Fig. 2) to (11) and after transforming it, we get:

a =

√
4
(a

b

)2
r2 + g2e (12)

Taking into account (9) we have:

a =

√
4
(

Rr(1−ν)− pz(2ν −1)
2pzν

)2
r2 + g2e (13)

pg. 91



Malinowska / Environmental Science and Sustainable Development, ESSD

Thus, the vertical coordinate of the cracking zone top node can be ultimately obtained from the dependence (14):

Zsp = Z0 +
ge
2
+ a = Z0 +

ge
2
+

√
4
(

Rr(1−ν)− pz(2ν −1)
2pzν

)2
r2 + g2e (14)

According to the above dependence, the vertical coordinate of the cracking zone for the 283 previously selected
and potentially endangered panels, was determined.

2.5. Determining of the cracking zone

In compliance with the theoretical assumptions presented in the previous sub chapter, the top node of the cracking
zone was calculated. For doing so it was necessary to determine the required geometrical and geomechanical
parameters. The assumption was made that the rock mass was layered; the parameters of each layer were assumed
experimentally and on the basis of literature data (Hwalek, 1971). Based on the theoretical dependencies obtained,
the location of the top node of the cracking zone was determined over each of the analyzed panels. Prior to
determining its degree, the risk of discontinuous deformations occurrence on the surface should first be classified.
The principle, according to which the risk zones were grouped, has been graphically presented in Fig. 3.

Figure 3. The classification of hazard zones due to location of cracks zone apex

The analysis of the obtained results reveals that most of the analyzed panels belong to risk zone one, and only two
cases to category four.

On the basis of the obtained results relating to the formation of a cracking zone and with the results defining the
potential risk of stability loss, it is possible to determine the shallow exploitation panels which can hypothetically
generate discontinuous deformations in the terrain surface.

The below table (Scheme 1) cross-juxtaposes the number of panels grouped in specific risk categories of stabil-
ity loss, (as a result of vertical stresses) with the risk categories of discontinuous deformations (resulting from
the localization of the top node of the cracking arch). In this way fifty-five sinkhole-prone chambers have been
distinguished.

Scheme 1: Analysis of the number of salt panels with respect to risk categories and cracks zones

It can be concluded from the above table that dozens of panels require more detailed analyses. Panels belonging
to hazard zone II can generate discontinuous deformations in the terrain surface mostly in the form of fractures,
whereas the cracking zones over such panels may bring about significant consequences for the safety of the mine’s

pg. 92



Malinowska / Environmental Science and Sustainable Development, ESSD

operation. Breaking the continuity of water-bearing horizons may result in an uncontrollable influx of water to the
mine. Three exploitation panels require special attention. They belong to risk zones three and four as discontinuous
deformations in the form of cones could appear above them.

2.6. GIS analysis

The next step was to match the theoretical results with mining and geological risk factors which were discussed in
chapter three. Global research devoted to the estimation of mining risks, has proven that geographical information
systems are an optimal tool for integration and the analysis of varying geological and mining factors (Malinowska,
Hejmanowski, 2010; Mancini et al., 2009) For this purpose, the selected risk factors were integrated using GIS.
The spatial analyses were mainly based on dependencies stemming from the Boolean logic. As a result of the
analysis performed, 15 chambers were selected as potential sinkhole risks for the terrain surface (Fig. 4).

Figure 4. Map of sinkhole hazard occurrence

3. Results and discussion

The author’s method assess the risk of discontinuous deformations on the terrain surface. This has been done
by modeling the vertical stresses on the basis of the theoretical state of rock mass deformation in the area of the
studied chambers is presented in the paper. For chambers potentially at risk of stresses over the boundary values,
the cracking zone was established using the pressure arch theory. Consequently, the number of chambers was
narrowed down to 55 chambers, which were potentially at risk of discontinuous deformations of the surface. By
supplementing theoretical analyses with results of additional studies of mining and geological factors, 15 sinkhole-
prone chambers were eventually selected. The final stage of the studies was conducted using spatial analyses, from
which chambers, in danger of being influenced by other risk factors were able to be distinguished.

4. Acknowledgments

This research has been financed by the National Sciene Center, grant No.UMO-2014/15/B/ST10/04892.

pg. 93



Malinowska / Environmental Science and Sustainable Development, ESSD

5. References

1. Brudnik K. (2010). The complex hydrogeology of the unique Wieliczka Salt Mine. Polish Geological
Review 58, 787–796 (in Polish).

2. Delle Rose M., Federico A., Parise M. (2004). Sinkhole genesis and evolution in Apulia, and their interrela-
tions with the anthropogenic environment.Nat Hazard Earth Sys,4, 747–755.

3. Malinowska, A., & Hejmanowski, R. (2010). Building damage risk assessment on mining terrains in Poland
with GIS application. International Journal of Rock Mechanics and Mining Sciences,47(2), 238-245.

4. Kleczkowski A. (1993). Groundwater in the vicinity of Cracow – potential and threats. Archives of the
Faculty of Hydrogeology and Geological Engineering of AGH University of Science and Technology, 11,
35–38 (in Polish).

5. Malinowska, A. A., & Dziarek, K. (2014). Modelling of cave-in occurrence using AHP and GIS. Natural
Hazards and Earth System Sciences,14(8), 1945-1951.

6. Malinowska A.A., Matonóg A. (2016). Spatial analysis of mining-geological factors connected with discon-
tinuous deformation occurrence. Acta Geodyn Geomater,14, 159–172.

7. Malinowska A.A., Misa, R., Tajduś, K. (2018). Geomechanical modeling of subsidence related strains
causing earth fissures. Acta Geod Geomater, 15, 197–204.

8. Malinowska, A. (2011). A fuzzy inference-based approach for building damage risk assessment on mining
terrains. Engineering Structures,33(1), 163-170.

9. Mancini, F., Stecchi, F., & Gabbianelli, G. (2009). GIS-based assessment of risk due to salt mining activities
at Tuzla (Bosnia and Herzegovina). Engineering Geology,109(3-4), 170-182.

10. Ryncarz, T. (1978). Analysis of the impact of filling ventricle excavations at levels I-III of the Wieliczka
Salt Mine on surface deformation and stability of lower-lying cavern subject to protection. Cracow: AGH
University of Science and Technology (in Polish).

11. Ryncarz, T. (1993). Outline of rock mass physics. Katowice: Silesian Technical Publisher (in Polish).

12. Saaty, T. L. (1980). The analytic hierarchy process. New York: McGraw-Hill.

13. Wieworka, P., Prunier-Leparmentier, A., Lembezat, C., Vanoudheusden, E., & Vernoux, J. (2009). 3D geo-
logical modelling at urban scale and mapping of ground movement susceptibility from gypsum dissolution:
The Paris example (France). Engineering Geology,105(1-2), 51-64.

14. Wieworka, W., o’Byrn, K. (2012). Selection of backfilling technology works in the Ksawer caverns complex
of the Wieliczka Salt Mine.AGH Journal of Mining and Geoengineering,36, 107–115 (in Polish).

15. Yilmaz, I. (2007). GIS based susceptibility mapping of karst depression in gypsum: A case study from Sivas
basin (Turkey). Engineering Geology,90(1-2), 89-103.

pg. 94


	Introduction
	Background
	Geology of the study area
	Sinkhole risk mitigation in urban areas

	Methods for assessing the risk of sinkhole occurrence
	Assumptions and general scheme of research
	Predicted deformations in the vertical plane
	Determining the stresses in vertical plane
	Determining the range of the cracking zone
	Determining of the cracking zone
	GIS analysis

	Results and discussion
	Acknowledgments
	References