The Impact of Earthquakes on Dropout Doline (Cover Collapse Sinkhole) Development: A Case Study from the Environs of Mečenčani and Borojevići (Croatia)

Abstract

In this study, the effect of earthquakes, beginning on 28 December 2020, on dropout doline development in the environs of Mečenčani and Borojevići was investigated. For that purpose, the shape of the doline, the inclination of the bearing surface and the rate of doline development were determined. A further analysis involved the characteristics of groundwater in the environs of the dolines and a functional relationship between the depth and the diameter of the dolines was sought. A model is proposed for the failure of the ceiling of cover cavities without support. The intensity of doline development is explained by favourable environment (dual cavity system, low inclination of the bearing surface, the presence and fluctuation of groundwater, etc.), the direct effect of earthquakes (material failure induced by earthquakes) and by their indirect effect (the partial solifluction of the ceiling material, lowered groundwater level).

1. Introduction

In this study, the development of dropout dolines (cover collapse sinkholes) in the environs of Mečenčani and Borojevići (Croatia) to the effect of the earthquakes of Petrinja between 28 December 2020 and 03 March 2021 is interpreted. Croatian researchers compiled a comprehensive documentation on the effect of the Petrinja earthquakes and also described the dolines that were formed during this time [1, 2]. All other investigations are of geomorphological nature [3], those concerned with remediation measures and priorities for immediate action [4], ground displacement using data from orbits of the Sentinel-1 mission action [5, 6]; fault geometry and the coseismic slip distribution [7], earthquake and deformations [8, 9]. To date, not a single investigation considered the effect of earthquakes on dropout doline development in the environs of Mečenčani and Borojevići such as the shape of the doline, the inclination of the bearing surface and the rate of doline formation. It is very important to emphasize that doline development damaged houses in the villages. The novelty of this investigation includes a model for the failure of the ceiling of cover cavities without support, analysis of the characteristics of groundwater in the environs of the dolines and a functional relationship between the depth and the diameter of the dolines.

Subsidence dolines (subsidence sinkholes) develop on unconsolidated, permeable or partly permeable rock (on concealed karst). Their varieties are dropout dolines (cover collapse sinkholes), suffosion dolines and compaction dolines [10, 11]. Dropout dolines are formed by collapse, while suffosion dolines develop by suffosion [10, 12].

The development of dropout dolines is also contributed by earthquakes [13, 14]. Earthquake waves may trigger rock failure and thus, the collapse of rocks that became looser.

Irreversible and reversible (fluctuating) water level changes occur to the effect of earthquakes [15], and they also affect collapses.

2. Description of the area

The Sunja River Valley is situated between two larger geographical/landscape units of Sisak Posavina and Banovina (Banija) (Figure 1) and belongs to the Petrinja-Sunja hilly terrain. Permanent surface watercourses constitute a highly developed parallel drainage network with significant amounts of surface water. The Sunja River is a right-bank tributary of the Sava River and is a part of the Sava River Basin (Danube River Basin).

Figure 1.

Overview map of the area with geographical coordinates 45°16′59″N and 16°25′30″E.

The Sunja River Valley is a covered karst where subsidence dolines occur, but they are also widespread between Hrvatska Kostajnica and Petrinja. Following the earthquakes of 2020–2021, further subsidence (dropout) dolines developed within a relatively short time in relatively great numbers [16, 17]. The altitude of the bearing surface is 175–190 m. The established level of the groundwater is 177–183 m and generally flows in S-E and N-E direction towards the springs Pašino vrelo, Bojanića vrelo and Davidovića vrelo (Figure 2). These springs create lakes in inactive dolines. The spring water of Pašino vrelo originates from Quaternary sediments in 70%, from Badenian limestone in 30%, while at Bojanića vrelo this proportion is 66% and 34%, at Davidovića vrelo 90% and 10% and at PBV-3 well 56% and 44% [18, 19].

Figure 2.

Hydrological map of the environs of the dolines ([18], modified). 1. Gravel, sandy clay, clay, coarse-grained stream load and boulder, 2. Badenian limestone, 3. Sandstone, conglomerate, marl, calcareous and clayey marl, sand, clay, 4. Detected geological boundary, 5. Assumed geological boundary, 6. Erosional-discordance boundary, 7. Detected fracture, 8. Assumed fracture, 9. Fracture series, 10. Doline, 11. Direction of groundwater flow, 12. Water level according to the state of 11 November 2005, 13. Fluvial measurement site, 14. Geological cross section (A-B) and the cross section of its groundwater level (C-D), and geophysical cross section (E-F) 15. Occupied permanent spring with permanent water, 16. Occupied spring with permanent water, for local use, 17. Spring of permanent water with low discharge, 18. Intermittent spring with higher discharge, 19. Intermittent spring, 20. Drilling, 21. Planned drilling, 22. Drilled well, 23. Dug well, 24. Subsurface water level.

The appearance of the springs is the consequence of a pronounced diagonal fault in the Sunja River Valley, which brought into direct contact permeable Upper Badenian predominantly carbonate rocks especially lithothamnium limestones (M₄) and poorly permeable to impermeable Pannonian marls (M₆) [18, 19]. Also, the area consists of Plio-Pleistocene age (Pl-Q) constituted by clay, clayey sand and pebbles and the alluvium (al- Q₂) constituted by pebbles, clay, clayey pebbles, older rocks, deluvial and proluvial sediments (Figures 2 and 3) [18, 20].

Figure 3.

Geological cross section (A-B) [18]. 1. Alluvium of the Sunja River: gravel, clay, sand, boulder, 2. Badenian limestone, 3. Eocene, Ottnangian, Sarmatian, Pannonian clastic sediments, 4. Boundary of beds, 5. Erosional discordance boundary, 6. Fracture, 7. Structural-piezometric drilling.

Geological structures in the investigated area extend mostly in NNW-SSE or NW-SE directions and follow the so-called ‘Dinaric’ strike (NW-SE), with predominantly dip-slip movements which are tectonically disturbed by the intersection of longitudinal NW-SE right-lateral and transverse NE-SW left-lateral faults of different size and importance along the transitional contact zone of the Dinarides and the Pannonian Basin [8]. The most important fault zone in the area is an active fault zone Pokupsko-Banja Luka in the Dinaridic ophiolite zone, Sava zone. According to [21], the strongest earthquake in the Kupa Valley M = 5.8 was recorded in the year 1909 (Mohorovičić discontinuity). The first strong earthquake in the wider Petrinja area was recorded on 28 December 2020. The day after, the second and third ones had magnitudes of M = 5 and M = 6.2 (h = 10 km). These earthquakes occurred on the Hrastovički fault, a segment of the Pokupsko Fault Zone stretching from Jastrebarsko through Pokupsko towards Banja Luka. These earthquakes caused the opening of approximately 100 new sinkholes (dolines) in the wider Mečenčani and Borojevići area.

Water inflow features are dolines, which are covered karst features and occur on the floor of the Sunja Valley. Among the dolines, there are distinguished [20] old dolines (they developed preceding the earthquake), more recent dolines (which developed during the earthquakes) and even newer dolines (which developed after the earthquakes) and primitive dolines and buried dolines (Figures 4 and 5). Old features (42 dolines) have gentle slopes, are covered with vegetation and without water. During the earthquakes, 82 new dolines developed in an area of 4 km² (but doline development also continued after the earthquake activity, and by the beginning of May, their number was 91 and by the beginning of December 2021, they numbered more than 100), which are collapse features with steep slopes (Figure 5). Several boreholes reached greater depths in order to get information on the composition of the cover. In case of a dropout doline, the sediments of the cover are organic soil, sandy lean clay and lean clay in 8-metre thickness [22].

Figure 4.

Distribution of dolines (modified from [20]). 1. Buried doline, 2. Primitive doline, 3. New doline, 4. Old doline, 5. Development date (year that is missing from the map: in 2020, the month is 12 and in 2021, the months are January and February), 6. Boundary of the area marked for the calculation of the development rate.

Figure 5.

New (dropout) subsidence dolines: A. non-narrowing dropout doline with lake (photo taken by Željko Grgić in February 2021), B. dropout doline with lake that is narrowing towards its floor from close-up (photo taken in February 2021 by Ronald Goršić): 1. Traces of liquefaction, 2. Limestone block, 3. Thrown out material, 4. Ragged lawn with traces of primary collapse, 5. Secondary collapse, C. narrowing dropout doline with lake from far view (a) and dropout doline without water (b) (photo taken in January 2021 by Marijan Car, Mario Bačić, Josip Terzić). D. Twin-like doline: Partial depression at the right side of the photo which developed either by the newer collapse of the cover cavity or during the secondary collapse of the side slope (photo taken by Sonja Zlatovič) 1, 2. Partial depressions, 3. Secondary collapse. All dolines were survey measured and monitored by Jeronim Moharić from January to September 2021 for the purpose of this study.

Water outflow sites are springs and spring lakes and the depressions containing water (paleodolines) such as the spring Pašino vrelo.

3. Methods

• Earthquake and precipitation data and the distribution of dolines were analysed.

• The water-level elevation of a lake of a doline with permanent water was compared with the water-level elevation of one of the wells.

• A hydrological profile was made, and the characteristics of water outflow as well as its relation to flow systems was analysed.

• The morphological and hydrological characteristics of the dolines were studied.

• The inclination of the bearing slope was calculated, the average width and depth of the dolines with various ages were compared, and the doline development rate was calculated (for the calculation of rate, density was calculated and for density, the occurrence area of new dolines was confined with tangent lines).

• The shape of dolines was calculated by the quotient of depth and diameter. Shape values were described according to frequency by putting them into classes, a functional relationship was determined between the shape of old and new dolines and the number of dolines belonging to different classes. A conclusion was drawn from the shape of dolines to the shape of their former cavities.

• Putting the new dolines into two groups (at one group, the depth was greater, and at the other group, the width was larger), a function relation was looked for between the width and the depth of the dolines.

• A theoretical study of the failure of cavity ceilings was performed in the following way.

The calculation of the ceiling thickness that can collapse without earthquakes for a given cavity diameter was done and led to the formation of a dropout doline at the surface. Using this relatively simple theoretical model, it is possible to determine the critical ceiling thickness at which collapse will necessarily occur.

The model under consideration is shown in Figure 6. It depicts a cylindrical cavity of diameter D at depth H below the surface. The initial assumption of model is that the ceiling of the cavity is not subject to any supporting forces from below, so the equilibrium of the rock mass above the cavity is ensured by the cohesion force and internal friction between the particles of the overlying aggregate. If these forces cannot counterbalance the weight of the rock mass, collapse will occur.

Figure 6.

Sketch of a theoretical model of dropout dolines. Explanation: 1. Uncollapsed rock mass, 2. Collapsed rock mass, 3. Cavity, 4. Interface between collapsed and intact rock masses, 5. Vertical stress in the rock, 6. Horizontal stress in the rock, 7. Compressive stress along the discontinuity surface, 8. Shear stress along the discontinuity surface, 9. The gravity of the destroyed rock aggregate, 10. The diameter of the cavity, 11. The thickness of the ceiling, 12. The slope of the discontinuity surface, 13. The infinitely thin layer of rock (imaginary) cut from the destroyed aggregate, 14. The depth below the surface of the selected layer, 15. The radius of the cut rock layer (disc), 16. The thickness of the cut rock layer.

Experience from mining and soil mechanics shows that in low-strength rocks such as the alluvium of the Mečenčani and Borojevići areas, collapse occurs along the upward-expanding fracture surfaces that start at the edge of the cavity and bound a conoidal frustum. The rock beds within the conoidal frustum are subjected vertically downwards to the gravitational forces of the Earth’s gravity and upwards to the cohesive and frictional forces of the rock particles outside the fracture surface along the conoidal frustum mantle. These forces are defined below. The rock beds within the stump cone are vertically downwardly influenced by the gravitational forces due to the Earth’s gravity and upwardly by the cohesive and frictional forces of rock particles outside the fracture surface along the conoidal frustum mantle. The formulae defining these forces are derived below.

The weight of the falling rock (G)

$$G=\mathit{\rho gV},$$E1

where ρ is the density of the rock [kg/m³], g is the acceleration due to gravity (g = 9,81 m/s²), V is the volume of the rock mass [m³], which in the case of a conoidal frustum is

$$V=\frac{\pi }{12}H\cdot \left(3D^2+3H\cdot D\mathrm{ctg}\alpha +H^2\mathrm{ct}{\mathrm{g}}^2\alpha \right).$$E2

α is the angle of the fracture surface with respect to the horizontal, which, according to the relationship known from soil mechanics [23]

$$\alpha =45°+\frac{\mathit{\varphi }}{2},$$E3

where ϕ is the angle of internal friction of the rock.

There are two types of forces acting along the conoidal frustum sheath (i.e. the discontinuity surface). One is proportional to the compressive stress (σ) perpendicular to the cone surface, the other proportional to the shear stress (τ) parallel to the cone’s surface. To determine σ and τ, let us imagine a prismatic body in the rock, infinitesimally small, bounded by a horizontal surface, a vertical surface and a surface parallel to the discontinuity surface. This body is loaded from above (in the vertical direction) by the pressure (σ_v):

$${\sigma }_v=\mathit{\rho gz},$$E4

and sideward (on the vertical side of the imaginary prism), a horizontal stress (σ_h) is applied:

$${\sigma }_h=\frac{\nu }{1-\nu }\mathit{\rho gz},$$E5

where z is the depth below-ground surface [m], and ν is the Poisson’s ratio of the rock. (Poisson’s ratio is the relationship between the transverse and longitudinal deformation of rock under uniaxial pressure. Its value varies between 0 and 0.5.) The equilibrium of forces acting on the prism requires that stresses (σ, τ) also occur on the inclined side. These give the stresses along the discontinuity surface. The compressive stress, which is perpendicular to the surface (according to the elementary relationship known from the strength theory), is

$$\sigma ={sin}^2\alpha \cdot {\sigma }_h+{cos}^2\alpha \cdot {\sigma }_v.$$E6

The shear stress (τ) acting parallel to the discontinuity surface at the moment of initiation of collapse consists of the internal friction between the rock particles sliding side by side, which is directly proportional to the sum of the compressive stress on either side of the fracture surface and the specific cohesive force “pulling” the particles together:

$$\tau =\mathrm{tg}\mathit{\varphi }\cdot \sigma +c,$$E7

where c is cohesion [N/m²].

The stresses σ and τ are specific, that is they give the force per unit area of the envelope of cone. Let f_v denote the vertical component of the resultant of these forces, which, according to elementary calculations, is

$$f_v=cos\alpha \cdot \sigma +sin\alpha \cdot \tau .$$E8

According to terms (3)–(7) (after appropriate aggregations)

$$f_v=\left(\nu \frac{cos\alpha {sin}^2\alpha +\mathrm{tg}\mathit{\varphi }{sin}^3\alpha }{1-\nu }+{cos}^3\alpha +\mathrm{tg}\mathit{\varphi }sin\alpha {cos}^2\alpha \right)\mathit{\rho gz}+sin\alpha \cdot c.$$E9

Imagine cutting out a ‘ring’ of height dz from the cone’s mantle at depth z. Its radius (r), is obviously

$$r=\frac{D}{2}+\left(H-z\right)\mathrm{ctg}\alpha .$$E10

Mark the area of the cloak of this ring dA:

$$\mathit{dA}=\frac{2\mathit{r\pi }}{sin\alpha }\mathit{dz}.$$E11

Since f_v is constant along this ring, the resultant of the vertical component of the compressive and frictional forces acting on the rock bed is

$$\mathit{dF}=f_v\cdot \mathit{dA},$$E12

If we add up the forces acting on all the elementary rings covering the envelope of cone of the rock-face before the collapse, we obtain the resultant F of the forces acting upwards (opposite to gravity) on the body, which, taking into account (9), (10) and (11)

$$F=\underset{0}{\overset{h}{\int }}\pi \left[D+2\left(H-z\right)\mathrm{ctg}\alpha \right]\cdot [\left(\nu \frac{cos\alpha sin\alpha +\mathrm{tg}\mathit{\varphi }{sin}^2\alpha }{1-\nu }+\frac{{cos}^3\alpha }{sin\alpha }+\mathrm{tg}\mathit{\varphi }{cos}^2\alpha \right)\mathit{\rho gz}+c]\mathit{dz}$$E13

For ease of writing, let us introduce the symbol

$$\kappa =\frac{1}{6}\left(\nu \frac{cos\alpha sin\alpha +\mathrm{tg}\mathit{\varphi }{sin}^2\alpha }{1-\nu }+\frac{{cos}^3\alpha }{sin\alpha }+\mathrm{tg}\mathit{\varphi }{cos}^2\alpha \right).$$E14

Performing the integration formulated in (13)

$$F=\mathit{\pi \kappa }\left(2H^3\mathrm{ctg}\alpha +3{\mathit{DH}}^2\right)\mathit{\rho g}+\pi \left(H^2\mathrm{ctg}\alpha +\mathit{DH}\right)c.$$E15

Eq. (15) gives the force that ensures the equilibrium of the cavity ceiling. It can be seen that this has a linear (first-order) relationship with cohesion, that is the more cohesive, more frictional rocks have a greater ‘holding power’ at the same depth than rocks with less cohesion or less internal friction. There is also a linear relationship between the retention force and the density of the rock.

The formation of a doline starts when the balance of forces is upset for some reason, that is when

$$F\le G,$$E16

The critical (H) cavity depth at which collapse can occur is given by F = G. Expressions (1) and (14) yield a second-degree equation for H:

$$\left(24\kappa \mathrm{ctg}\alpha -\mathrm{ct}{\mathrm{g}}^2\alpha \right)\rho {\mathit{gH}}^2+\left[12{\sigma }_c\mathrm{tg}\mathit{\varphi }\mathrm{ctg}\alpha +\left(36\mathit{\kappa D}-3D\mathrm{ctg}\alpha \right)\mathit{\rho g}\right]H+\left(12\mathit{Dc}-3D^2\mathit{\rho g}\right)=0$$E17

A numerical analysis of (17) (with realistic data for the Mečenčani and Borojevići area) shows that the first term is several orders of magnitude smaller than the second and third terms and is therefore negligible. The numerical analysis also shows that the first member of the second bracket term is negligible next to the second and third. The solution to the simplified equation for H (taking into account the expression (14) for κ) is given by:

$$H=\frac{1-\nu }{2\nu }\frac{D-\frac{4c}{\mathit{\rho g}}}{sin\alpha \left(cos\alpha +\mathrm{tg}\mathit{\varphi }sin\alpha \right)}.$$E18

H is the critical ceiling thickness (expression 18) at which the roof covering is ‘just’ not yet collapsing. However, at this ceiling thickness, the cavity is already unstable, that is the slightest disturbance (earthquake, rock scour, soil moisture variation, etc.) will cause immediate collapse.

4. Results

The pattern of doline distribution shows their tectonic determination. The fault zone detected in the bedrock [3, 17, 19] directed the karst water flow, the transportation of dissolved material and thus, the direction of cavity formation there. However, the cavities that formed in this way determined the site of the material loss of the cover and the zone in which it took place. The geoelectrical resistance values of the profiles indicate the extent of fracturing and cavity formation of the bedrock below the zone of the dolines.

Based on the earthquake and precipitation data of the area, the following can be established:

• Between 28 December 2020 and 3 March 2021, two foreshocks and 14 aftershocks took place with magnitudes of 3.3–5.2. Their hypocentres were at depths of 8.7–10.9 km and their epicentre was in Strašnik 3 km away from Petrinja in WSW direction [1]. The development of dolines, some of them exactly dated, was mainly related to the earthquakes in January (Figure 4).

• According to the data of the Kostajnica Meteorological station, January was an especially rainy month, with a total of 106.4 mm rainfall, but February was also wet (72.7 mm). In January, there was a daily 6 mm or more rainfall eight times per month (once its value was 21 mm and on another occasion 28 mm).

Taking into consideration the observed dates of doline development (Figure 4, available only for some dolines), the earthquakes may have triggered immediately (on the day of the earthquake) and delayed doline development (days after the earthquakes). An example for immediate doline formation was what happened on 29 December 2020 when two dolines were formed and 30 December 2020 and 31 December 2020 when one doline developed each day. A similar situation was observable on the days of earthquakes of 2 January 2021, 4 January 2021, 5 January 2021 and on 6 January 2021 when one doline was formed each day. On the contrary, a delayed doline formation took place on 23 January 2021 when two dolines developed; however, no earthquake occurred on that day, but they happened only on the days of 7 January 2021, 9 January 2021, 10 January 2021 and 15 January 2021. No doline development coincides with the day of significant amount of rainfall. However, on 29 December 2020 when two dolines developed, the quantity of precipitation was significant. The relation between rainfall (infiltration) and doline development is less obvious. On 29 December 2020, two dolines developed at rainfall of 6 mm, while one doline was formed at precipitation of 12 mm on 30 December 2020. It is more probable that precipitation has a delayed effect. On 26 December 2020 when precipitation was 31 mm, it must have ensured favourable conditions for the development of the above-mentioned dolines on the 29th and 30th of December. The role of precipitation in doline development is indicated by the fact that dolines continued to take shape at the end of 2021. Thus, in December 2021, four new dolines developed probably to the effect of autumn rainfalls and also after a stronger M 3.2 earthquake.

The distribution of depressions is of banded (at some places, of linear) pattern. The band has a NW-SE direction and a length of 3100 m (outside that band only six dolines occur which developed during the earthquakes). Perpendicular to this direction, the width of the band is maximum 800–850 m (Figure 4). The band of dolines intersects the Sunja River in NW, and only some new dolines occur on its left bank. The width of the band covered by old dolines is larger; thus, the appearance of new dolines is more concentrated than that of the old dolines. The occurrence of new dolines at the north-western and south-eastern ends of the band (mainly at the latter) is slightly dispersed. However, older dolines frequently group along this arcuate band. The occurrence of dolines is combined in the band: there are sections made up of only new and only old depressions, but mixed sections are also present. Depressions outside this arcuate band are in groups of an irregular pattern. These may be constituted by homogeneous depressions (made up of only new or only old dolines), but by depressions that developed at various time. The direction of the band and the tectonic structure as well as the concordance of the directions of some fractures (faults) refers to the fact that doline development is controlled by tectonics. The area bearing the dolines (this is the floor of the Sunja River Valley) is of a very low inclination. It is 1.29° along the profile perpendicular to the Sunja River (Figure 2).

The development rate of new dolines is extremely high. If we consider the area with dolines and presume a period from 28 December to the end of March (doline development is probably shorter than this: in March no doline development was likely to the direct effect of the earthquakes, but the last earthquake, as we observed, took place on 3 March) calculating with 93 days, in the bearing area of 0.0113 km², the average rate of doline development is 27.33 dolines per month. Calculating this to 1 km² and 1 month, the rate is 2418.70 doline/month or calculating with the actual area to 1 year, it is 327.96 doline/1 km², calculating to 1 km² and 1 year, 201.56 doline develops, and to 1000 m² and 1 month the development rate is 29.02 doline/month. Calculating with the actual area, the development rate is 0.91 doline per day. Thus, comparing with the development rates of literary data, in the area of Mečenčani and Borojevići, the same number or more dolines per day were formed during a day than in other karst areas during a year. This high rate of doline development can only be explained by the impact of the series of earthquakes.

The morphological and hydrological characteristics of the dolines are the following:

• Dolines that developed preceding the earthquakes have gentle slopes with vegetation.

• Two varieties of new dolines that were formed during the earthquakes can be distinguished morphologically. Dolines belonging to one variety have gentler slopes, narrowing towards their floor, but their margins are made up of arcuate sections (Figure 5B). The other variety involves dolines with steep slopes (sub-vertical) and a wide, expanded floor (Figure 5A). The cover is exposed on the dolines of both varieties, their margin is sharp and the lawn is ragged. The morphology of new dolines refers to the collapse origin. Dolines collapse for two reasons: primarily resulting in their development and secondarily which collapse on their side slopes (Figure 5). The result of secondary collapses is that dolines are widening and their slopes are becoming gentle. At some sites, patchy and zone-like sediment accumulations can be seen beyond the margins (Figure 5B and C). On their side slopes, flow traces (Figure 5B) and the scars and slides of subsequent, secondary collapses occur (Figure 5D). Sometimes, the side slope is constituted by parts of various steepness (parts of low inclination replace steeper slope sections or the floor of variable elevation (Figure 5C and D). They have a circular ground plan, but elongated dolines and twin-like complex dolines also occur (Figure 5D). The surface is flat in the environs of the dolines, and no drainage feature is connected to it.

The depth of new depressions does not exceed 5 m predominantly; thus, they were formed in the unconsolidated superficial deposit (al-Q₂). Therefore, they are subsidence dolines. Only one doline has a depth of 12 m. However, it refers to the fact that the collapse triggering its development also spreads over the Badenian limestone. The cavity of the Badenian limestone must have been without water (or partly without water) at least during the time of the earthquake which triggered its collapse since the collapse of the cavity is only possible if the ceiling of the cavity is not supported. However, the development of some depressions (those with a depth close to 5 m) may have started at the Badenian limestone since the depth of the doline was reduced by the collapsed material.

Among the new ones, there are depressions that have been dry since their development, there may be constantly wet depressions (there is a lake with permanent water on their floor) and there are depressions that are intermittently wet (a lake appears on their floor intermittently). The altitude difference between the water level of the doline with permanent water and the groundwater table is 4.5 m, and the distance between the distance between the doline with lake and the well is 80 m. The similarity in water levels indicates that the water of dolines with lake partly originates from groundwater or the dolines developed from cavities that had at least partly been below the groundwater table. (The reason for the difference of the above water levels may be that the collapse of the cavity increased the water level in the doline and the water level of the lake follows the subsidence of the groundwater level with delay.) The water-level fluctuates in both the permanent and intermittent lakes.

If the altitude of the floor of new dolines is almost the same as the altitude of the cavities from which they developed, then the extent to which the dolines are filled with water or its lack indicates the position of former cavities as compared to the groundwater level. As compared to the groundwater level, the dolines and their cavities may be and may have been of the following types:

• Dolines which are constantly dry (a depression that has been dry since its development can be seen at the right side of Figure 5C), the floor of their cavity has always been above the groundwater level. These are dolines of small depth (shallower than 1–2 m).

• Dolines with intermittent lakes, and their cavity was situated at least partly between the high and low groundwater levels (according to drying traces visible on its floor, its depression lost its water by the observation date of 5 June 2021.) The water level of the lakes changes in the dolines because their floor is in the groundwater fluctuation zone. These dolines are deeper than the former. Their depth ranges from 1 to 2 to 5 m.

• Dolines with permanent lake, the floor of their cavity was below the current groundwater level. An example is the 12-m-deep doline, but probably at least three other dolines whose depth is about 5 m belong to this group. This doline still had a lake at an observation of November 2021.

• Dolines that developed after the earthquakes do not have water intermittently either; thus, they also developed from cavities above high groundwater level.

Dolines of water outflow sites have steep slopes and permanent lakes. These features are not active (paleodolines) such as the spring Pašino vrelo. They are situated at the emergence sites of groundwater and karstwater. Based on the composition of their water, they are drainage sites for them. According to [19], groundwater from limestones and water from Quaternary deposits are mixed. The groundwater is the water of the al-Q₂ beds; in the NE, it follows the Sunja River and is situated at a depth of 1–2 m. Its position as compared to the surface SW from Mečenčani is unknown. However, it can be observed that the altitude difference between the surface and the water level increases farther from the Sunja River (Figure 7). The different proportion of the two waters at various drainage sites refers to the existence of two kinds of flow systems, but also to their partial or widespread relation. Their discharge fluctuation was between 24.4 and 12.4 l/s at the spring of Pašino vrelo. The fluctuation of the groundwater discharge and also the fluctuation of its water level [18, 19] refer to the presence of high and low water levels.

Figure 7.

Position of groundwater. Legend: 1. surface, 2. groundwater level.

The size of dolines that developed at water inflow sites at various times is different. Older dolines are of a greater diameter and a smaller depth (the average diameter is 5.1 m, the average depth is 0.47 m), newer dolines have a smaller diameter, but they are deeper (their average diameter is 3.7 m, and their average depth is 1.5 m). These values are more striking if we only consider the sizes of those with a shape index below 0.3 in case of old dolines (Figure 8): here, the average diameter is 5.73 m, and the average depth is 0.34 m. The latter are the smallest, with a diameter of only 1–2 m, and having steep slopes, they are passage-like and aligned with a collapse yard of small depth. The shape of 34 old dolines is in the 0.0–0.1, 0.1–0.2, 0.2–0.3 classes and only eight dolines have a shape larger than class 0.3 (Figure 8). The shape of 63 new dolines is in classes larger than 0.3. Thus, at old dolines, a larger diameter belongs to a given depth, while in case of young dolines, a smaller diameter belongs to a given depth. The class distribution of the two doline groups is also different. The distribution of old dolines can be described with a decreasing polynomial function, while that of new dolines can be described with a polynomial function with two maximum values (Figure 8).

Figure 8.

Shape distribution of old and new dolines and the functions of shape distributions.

New dolines can be referred into two groups. At one of them, the depth is dominant (the depth is larger than the width), at the other, the width is dominant (the width is larger than the depth) since there is a functional relationship between the width and depth at the members (at every depression) of the group (Figure 9). If the depressions developed from cover cavities, the former cavities can also be put into two types: at one of them, the vertical expansion is dominant, while at the other, the horizontal expansion is dominant. Therefore, the cavities that were situated at a similar depth as compared to the surface were mainly either vertical or horizontal. During the collapse, the cavities preserved their horizontal and vertical dimensions and they were only modified to a small degree. The small number of vertical depressions (there are four) and the large number of depressions with large width prove that great width favours the denudation of the cavity ceiling.

Figure 9.

Functional relationship between the width and depth of new (dropout) dolines. 1. Dolines whose depth exceeds their width, 2. Dolines whose width is larger than their depth, 3. Function fitted to deeper dolines, 4. Function fitted to wider dolines.

In case of depressions whose depth exceeds their diameter, there is a strong relation between the width and the depth at the width and depth function (Figure 9). Depressions becoming wider and wider are more and more expanded vertically. Thus, former wider cavities were probably more and more expanded downwards. At depressions where their width exceeds their depth, in case of those with a small width, depths are less diverse. At depressions whose width is larger, the diversity of their depth increases. Therefore, in the latter case, the cavities from which the depressions developed had an extremely diverse vertical cavity size during their widening. This is possible since the floors of some cavities from which the depressions developed had a deeper and deeper position in case of similar thicknesses of the cavity ceilings.

According to soil mechanics studies in the Mečenčani and Borojevići area (Tomac et al. 2021c), the alluvium has an angle of internal friction of ϕ ≈ 24–28°, a specific cohesive strength c ≈  20–50 kPa, a density of ρ ≈ 1800–2000 kg/m³, a Poisson’s ratio ν ≈ 0.25–0.33 and a gravity acceleration g = 9.81 m/s² (Tomac et al. 2021c). Based on the mean value of these (ϕ = 25°, c = 20 kPa, ρ = 1800 kg/m³, ν = 0.33), the ceiling may stop above a cavity of diameter D = 10 m at a thickness of H = 6.7 m.

To illustrate the relationship (18), the evolution of the ceiling thicknesses at the stability limit as a function of the size of the cavity before collapse for different realistic values of cohesion (c) and internal friction angle (ϕ) is presented in Table 1. (ρ = 2000 kg/m³, ν = 0.33) It can be seen that H is very sensitive to both ϕ and c. During earthquakes, these factors can undergo significant changes, which can result in previously quiescent cavity ceilings losing their stability and collapsing.

D [m]	H [m]
	c = 10 [kPa]			c = 20 [kPa]			c = 30 [kPa]			c = 40 [kPa]
	Φ = 20°	Φ = 25°	Φ = 30°	Φ = 20°	Φ = 25°	Φ = 30°	Φ = 20°	Φ = 25°	Φ = 30°	Φ = 20°	Φ = 25°	Φ = 30°
2	—	—	—	—	—	—	—	—	—	—	—	—
4	2.8	2.5	2.3	—	—	—	—	—	—	—	—	—
6	5.6	5.1	4.6	2.7	2.5	2.3	—	—	—	—	—	—
8	8.5	7.7	7.0	5.6	5.1	4.6	2.7	2.4	2.2	—	—	—
10	11.3	10.3	9.3	8.4	7.7	6.9	5.5	5.0	4.6	2.6	2.4	2.2
12	14.2	12.9	11.7	11.3	10.2	9.3	8.4	7.6	6.9	5.5	5.0	4.5
14	17.0	15.5	14.0	14.1	12.8	11.6	11.2	10.2	9.2	8.3	7.6	6.9
16	19.8	18.1	16.4	16.9	15.4	14.0	14.1	12.8	11.6	11.2	10.1	9.2
18	22.7	20.6	18.7	19.8	18.0	16.3	16.9	15.4	13.9	14.0	12.7	11.5
20	25.5	23.2	21.1	22.6	20.6	18.7	19.7	18.0	16.3	16.8	15.3	13.9

Table 1.

Critical (on a limit of stability) ceiling thickness (H) as a function of cavity size (D) for different cohesion (c) and internal friction (ϕ).

5. Discussion

The morphology and shape of the dolines indicate that the dolines of the area between Mečenčani and Borojevići originated by suffosion, collapse and probably compaction.

Old dolines with a shape smaller than 0.3 have a small depth because of limited material transport (suffosion and/or compaction), a large diameter (widespread material transport can take place) and gentle slopes (as a result of slow subsidence and slope denudation). New dolines with a shape larger than 0.3 are deeper because of collapses. Therefore, old dolines with a shape index smaller than 0.3 are suffosion dolines, while new dolines with a shape larger than this are dropout dolines.

Among old dolines occur dolines with a shape larger than 0.3 (Figure 8). Therefore, old dolines may also have developed by collapse, but their side slopes became gentle subsequently (pseudo-suffosion doline). Among new dolines, there occur dolines with a shape smaller than 0.3. This may refer to the fact that in case of new dolines of this type, collapse already stopped at an early phase or was followed by subsidence.

In the case of some new (dropout) dolines, the reason for shape increase is not width decrease, but ever greater depth. Presuming a cavity ceiling of small thickness (which is a precondition of collapse), this is possible if the former cavities were increasingly expanded downwards. Thus, cavities develop in the zone of groundwater fluctuation; in addition, in case of those with a shape larger than 1, there was a cavity which was originally in the bedrock; then, its development spreads onto the cover.

The relatively large number of old (suffosion) doline refers to the fact that in this area the tendency to doline development might have been significant earlier as well. This can probably be explained by the low inclination of the bearing terrain in addition to cavity formation that is discussed below. In our case, surface inclination on the valley floor, as already mentioned, is low. This can partly be traced back to the fact that on such terrains, the infiltration proportion of precipitation is large; on the other hand, as a result of its low inclination, the groundwater level is in wide expansion close to the surface (and the karstwater level is also close to the bedrock surface). The water level being close to the surface favours widespread subsurface cavity formation, which enables widespread material transport from the cover. However, collapse tendency is high because the cover cavities are close to the surface. The cavity formation in the bedrock favours collapses and the reception of the material transported from the cover.

Regarding the water outflow sites of the area, it can be established that the water level oscillates at Pašino vrelo whose degree alternates between 0.46 and 0.52 m since the discharge of its spring fluctuates between 22.8 and 12.8 l/min [18, 19]. This indicates the fluctuation of groundwater level, that is high groundwater level and low groundwater level. Their range between dry and wet seasons is 2 m [16].

Since the groundwater of the cover and the karstwater of the Badenian limestone are present to various proportions at springs Pašino vrelo, at Davidovića vrelo [18, 19], the two kinds of water are separated from each other at least temporarily, and they are primarily mixed at the water outflow sites. The water of the Badenian limestone functions as artesian water in wet seasons [16]. Therefore, the groundwater and the karstwater flow circulate individually.

A two-storey cavity system developed as a result of the two water bodies and their circulation and the fluctuation of groundwater level. The lower level is in the Badenian limestone, while the upper is in the cover.

The shape and size of cover cavities were different. This is proved by the fact that the depth and shape of dolines differ for the dolines (it is presumed that during their development, the dolines preserve the shape and size of former cavities). The shape of the former cavity can be determined from the shape of the depression. The shape, but also the size and position of the cavity point to a collapse tendency. The wider the cavity (the smaller shape it has) and the closer to the surface, the greater its collapse tendency.

Cover cavities are of various positions as compared to the changing groundwater level, since among the developed dropout dolines there occur dry dolines, dolines with permanent lakes and dolines with intermittent lakes. (According to the observation of 05 June 2021, some dolines lost their water.) Among former cavities there were some which were situated in their whole expansion above the high groundwater level and there were others which were above the low groundwater level in their whole expansion. Cavities, whose lower part was below the groundwater table developed by the partial collapse of the bedrock and the cover above it since the permanent presence of groundwater, hindered suffosion.

The altitude of the groundwater level affects the size and development of the cavities and the way of material transport out of the cavities. The difference between the altitude of the groundwater level and the altitude of the surface increases in SW direction from the Sunja River (Figure 7). However, the extent of the fluctuation of high and low groundwater levels also increases in the same direction since the water supply from the Sunja River into the groundwater and thus, its affect raising the water level, diminishes moving farther from the river. In case of low water level, its water level reducing effect is limited since the water level cannot reach below the water level of the river.

Consequently, the cavities of the cover may have developed in the following manner:

• By suffosion, cavities in the environs of which the cover is at least temporarily without groundwater and it is coarse-grained. Suffosion is possible because a passage develops above the cavity of the bedrock that reaches above the karstwater level and this transports the material of the cover into the karstic cavity.

• By the further collapse of the suffosional embryonic cavities of the cover.

• By the collapse of the bedrock. The cavity of the bedrock is inherited onto the cover, and then to the surface. This rarely observed process must have taken place, for example in the case of the doline with a depth of 12 m. If inheritance happens only onto the cover, a cavity is formed in the cover by collapse and with its later collapse a doline develops. If the cover follows the bedrock collapse in its total expansion, no cavity develops in the cover, but a depression is formed immediately at the surface.

There are a large number of cavities in the upper level of the Badenian limestone according to measurements [18]. The cavities are situated both above and below the water level. Evidence for the existence of cavities above the water level is the already mentioned 12-m-deep-doline. Since it could have developed in a way that it continues in the bedrock, thus, the bedrock cavity had to collapse as well. It is only possible if the cavity reached above the water level.

The material flow of the two-level cavity system may be the following:

• Material transport solely happens by suffosion. This is possible at those cavities that are above the high groundwater level (and at cavity parts which are situated above this level). Suffosion always takes place during the rainy season. In the cover, vertical passages develop between the surface and the cavity and at the portion of the cavity above the high groundwater level. Passage development is associated with suffosion doline formation and with the partial filling up of the cavity. If suffosion material transport takes place from the cavity, the cavity increases and its sediments are transported away.

• Material transport is partly of suffosion and partly of collapse origin. The cavity is between the high and the low groundwater level in its partial or complete expansion. Suffosion only takes place and only for the period until the groundwater level is below the floor of the cavity. In this case, passages may form too. If the groundwater is above the cavity, suffosion is horizontal, which does not generate passages, but it induces material transport in the groundwater. If the amount of precipitation increases, vertical suffosion stops because of the rise of the groundwater level.

• The transport of the cover material is solely of collapse origin if the cavities of the Badenian limestone above the karstwater level collapse. Suffosion may be present accessorily, temporarily which may be expressed by the reworking of the collapse material. The cavity is either in the bedrock or both in the bedrock and in the cover. The cavity of the cover may continue its development by collapse (dropout doline) or by suffosion (suffosion doline).

There is a higher probability of old (suffosion) doline development at sites where:

• the cavity is situated above the high groundwater level in its total expansion,

• the host rock of the cover cavity is less cohesive,

• there is no earthquake.

New (dropout) dolines may develop in the following ways.

• By cover collapse. The cover cavity may be formed by bedrock collapse, which is inherited onto the cover or by suffosion. The cover cavity collapses if the pore water pressure increases in its environs, to the effect of earthquakes if the water level decreases in the cavity, but it may also take place without any external effect if the width of the cavity reaches a critical value and the ceiling is thin.

• By the collapse of the bedrock which is inherited onto the cover in a way that no cavity is formed there, the material of the bedrock ceiling and the cover collapse together. This may happen if the cavity is above the karstwater level, close to the bedrock surface and horizontally developed.

The collapse of cover cavities may take place only during the failure of the ceiling in case of cavities above the high groundwater level (direct earthquake effect) and when it was contributed by other factors too (indirect earthquake effect). The collapses that triggered the development of dolines in the environs of Mečenčani and Borojevići may have been affected by several factors too if the ceiling did not directly collapse because of the failure such as the partial liquefaction of the ceiling material may have happened at cover cavities situated below the high groundwater level. Traces of liquefaction can be seen in Figure 5B. Concomitant phenomena, material shots can be recognized in the environs of some dolines (Figure 5B and C).

In the area of Mečenčani and Borojevići, after the series of earthquakes, the groundwater level became (and stayed) somewhere 30–50 cm lower than earlier in the wells (conversations with locals in 2021). This was probably caused by the irreversible displacement of karstic blocks and the resulting karstwater subsidence spreads onto the groundwater too and/or because of the decrease of cover compaction the groundwater level sank in an irreversible way. The groundwater subsidence reduced the support of cavity ceilings (doline development is without support), which launched a process that is well known in karst literature [24]: the collapse of the ceiling weakened by liquefaction. The thinning out of the ceilings was also enabled by the fact that parts became separated from there, partly because the infiltration of precipitation increased pore water pressure and partly because of liquefaction. However, the infiltrating precipitation also caused an increase in burden, which resulted in the warping of the ceilings. Local people reported on surface subsidence preceding doline formation [16]. This may also have enhanced the separation of parts from the ceilings. To the effect of repeating shocks, the chance of the collapse of thicker ceilings may also increase. Cavities with thin ceilings probably collapsed to the immediate effect of earthquakes, while the cavities with thicker ceilings collapsed due to the delayed effect from the cumulative impact of the shocks.

Dolines, the ceilings of whose cavities had been above the high groundwater level, may have developed during the failure of the cover due to direct earthquake effect by collapse. A favourable condition for this, as already mentioned, is that the cover thickness is less than 1–2 m, and it is not greater than some tens of centimetres. Similarly, the collapse of bedrock cavities was caused by the direct effect of earthquakes.

The relationship (18) highlights several factors that are important for the formation of dropout dolines and even allows a quantitative analysis of these factors.

It can be seen that there is a linear relationship between the thickness (H) of the still stable cavity ceiling and the diameter (D) of the cavity: the main thickness of a cavity with a stable ceiling (i.e. not yet collapsing) is proportionally larger. It is noteworthy that according to this model, there is a minimum cavity diameter below which the ceiling of smaller cavities (due to their own weight) does not collapse. If

$$D<\frac{4c}{\mathit{\rho g}},$$E19

then the numerator of (19) becomes negative, which means that no matter how thin the cavity ceiling is, it is held together by cohesion. (Of course, in this extreme case, other external factors neglected in the present model may still trigger collapse.)

From (18), it can be seen that if the cohesion (c) of the rock is larger, then H is smaller. This is consistent with the experience that more consolidated rocks do not collapse even with thinner ceilings. It also turns out that if the Poisson’s ratio (ν) of the rock is smaller, a larger cover thickness is required to avoid collapse. The reason for this is because, with a smaller Poisson’s ratio, the primary vertical rock pressure (σ_v) results in a smaller horizontal pressure (σ_h) (see relation (4)), which results in a smaller frictional force balancing the cover assembly.

The relationship (18) gives an account of all the effects of earthquakes that lead to the formation of dropout sinkholes. The most important of these (which also play a central role in building damage) are the accelerations induced by the quake waves. These can be calculated from seismology (using the distance from the epicentre or hypocentre of the quake and the magnitude of the quake) and are added to the gravity acceleration. For example, for a 5-magnitude quake, this can increase the value of g by up to 30%. Since in expression (18)g is contained in the denominator of a subtraction term, and an increase in g reduces the subtraction term and therefore increases the value of H. In practice, this means that a cavity whose ceiling thickness still ensures stability under ‘normal’ conditions will no longer comply and collapse due to the increased critical principal strain during an earthquake.

Similarly, the reduction in internal friction and cohesion due to the shear stress of the earthquake increases the value of H: a reduction in both c and ϕ leads to an increase in the value of H according to (18); that is, the effective cover thickness becomes less than the minimum value required to avoid collapse.

6. Conclusions

The karst features of Mečenčani area were classified, the morphology of dolines was described and dolines were put into groups based on their position as compared to the groundwater level. The coincidence of the development of dolines and the time of earthquakes was described. The shape distribution of old and new dolines was also studied.

Statements related to doline development are as follows:

• The variability of doline development was due to the diversity of hosting sediments. Therefore, dolines of various types and different formation date may occur next to each other.

• The dropout dolines may have developed by the collapse of the cover cavity of suffosion origin, by the inheritance of the bedrock cavity onto the cover by collapse, from which a dropout doline develops by collapse too. They may develop by the inheritance of the bedrock cavity onto the surface through the cover.

• Based on the model proposed, the effect of rock mechanics parameters that influence the stability of cavities can be given in numbers.

• Earthquakes may result in collapse directly on the ceiling of the cover or of the bedrock. Collapse can take place because accelerations due to earthquakes are added to acceleration of gravity increasing the weight of the rock mass and decreasing friction and cohesion. Collapse may also take place indirectly by earthquakes when the ceiling loses its support and it becomes thinner by liquefaction.

• During periods without earthquakes primarily suffosion dolines develop, while the periods of earthquakes result in dropout dolines. As a result of the presence of groundwater, suffosion dolines can only develop at cavities situated above lower groundwater level.

• The high development rate during earthquakes was accompanied by cavity formation due to dual water body and groundwater-level fluctuation and attributed to the thin ceiling of the cavities and low slope angle.

• Doline development would have taken place in any case because of the cavities of the bedrock and the cover, but not within such short time, but during several hundreds or thousands of years and probably partly by suffosion.

• Wide cavities situated in the fluctuating groundwater level have the highest probability of dropout doline formation due to earthquakes.

• In contrast to literary data, the development of suffosion dolines in the areas of Mečenčani and Borojevići did not take place by the transport of the cover material into the bedrock, but it also happened by transport into the cover cavity.

The study of the dolines of the Mečenčani area draws attention to the fact that doline development is influenced by several factors. Some factors are only specific of this area (the closeness of the karstwater level to the bedrock), while others are generally valid (groundwater, earthquake). The significance of our study is that the role of various factors in doline development was given in quantities. Thus, a better prognosis for doline development can be given on various covered karsts of the Earth and human made structures can be planned more safely.

Acknowledgments

Many thanks for overall support to Mr. Gordan Hanžek (State Secretary of the Central State Office for Reconstruction and Housing of the Republic of Croatia) and the Ministry of Physical Planning, Construction and State Property. We are very grateful for the help with collection data on the field to Mr. Jeronim Moharić with colleagues (Geo Gauss d.o.o. from Čakovec), Borna Gašpert and Mirko Stanković. The English language was edited by Tatjana Jauk (linguist at Hrvatske vode) and by Dénes Lóczy (University of Pécs), and we thank them a lot for her kind help.

Conflict of interest

There is no conflict of interest.