1 Introduction

The Apulian region (south-eastern Italy) is broadly characterized by the presence of underground cavities in soft rock formations. These caves generally originated as a result of the extraction processes of soft rock, used as construction material, according to the so-called “room and pillar” technique. In those cases where cave systems are located rather close to the ground surface, the stability is marginal and the eventual collapse would imply interaction with overlying structures and infrastructure (Parise and Lollino 2011; Castellanza et al. 2018). In fact, in the past, urban development and infrastructure building operations have rarely been planned considering the possibility of interaction with abandoned underground cavities, therefore taking into account the assessment of the corresponding sinkhole hazard (Delle Rose et al. 2004; Parise 2010; Waltham and Lu 2007; Zhou and Beck 2011). Over time, the risk of sinkholes affecting built-up areas has become a real problem due to an accelerating urban expansion as well as the evolution of rock weathering and related degradation processes over time (Lollino et al. 2013; Fazio et al. 2017; Castellanza et al. 2018). Due to a combination of geological, hydrological, and seismic characteristics, the anthropogenic underground cavity systems are characterized by a high degree of failure susceptibility and represent a serious threat to the built heritage, due to the sudden and unexpected occurrence of sinkholes, often triggered by rainfall events and even moderate earthquakes (di Santolo et al. 2015). Generally, in the literature, the risk due to the presence of underground cavities is associated with weathering effects induced by climatic factors. These act on porous rocks generating water-induced weakening processes, with the formation of micro-fractures and consequent reduction in strength over time (Auvray et al. 2004; Castellanza and Nova 2004; Grgic et al. 2006; Andriani and Walsh 2007; Ciantia and Hueckel 2013; Ciantia et al. 2015a, b). The evaluation of the stability of underground cavities becomes highly intricate in the context of “room and pillar” cave systems. These systems present a challenging geometry problem, intricate boundary and loading conditions, and a dynamic evolution of both boundary conditions and rock material properties over time. This complexity required more advanced methods to address these multifaceted challenges. Numerical modelling has emerged as a potent tool for investigating the stability of the rock masses surrounding the cave, with the study of the stress–strain conditions, as well as the resultant displacement field caused by specific loading conditions or alterations in boundary conditions. Bekendam (1998) investigated the stability of calcarenite and limestone mine pillars in the Netherlands using two-dimensional elasto-plastic finite element (FE) models, incorporating time-dependent creep processes. Similarly, Parise and Lollino (2011) employed the same 2D FE approach to emphasize the role of degradation processes in limestone and calcarenite rocks surrounding caves in Southern Italy, particularly in the formation of sinkholes. Ghabezloo and Pouya (2008) conducted finite element (FE) analyses to examine the stability of the roof in limestone caves in France, focusing on the degradation of tensile strength caused by karst processes. Whereas, Diederichs (2003) explored rock fracture mechanisms and the overall collapse of caves using the distinct element method.

However, in several areas, such as the Apulia region, the seismic action could even represent a possible triggering factor for cave failure, as discussed in detail by de Lucia et al. (2022), which highlights the role of the dynamic signal and the cave geometry in the dynamic stability of single cavities or cave systems. Zhao et al. (2021) have also recently proposed a numerical investigation of the seismic performance of underground loess caves with traditional dwellings in China and the results have been compared with the observations derived from shake table tests. Nevertheless, the role of dynamic stresses on man-made underground cave stability has been frequently underestimated or assessed only with qualitative approaches (Sharma and Judd 1991). The vast majority of studies were focused on the modification of the seismic motion at the ground surface as a consequence of the presence of the cavities, as discussed by Lee and Karl (1992), Smerzini et al. (2009), Conte and Dente (1993) to mention a few. Genis and Gercek (2003) first explored the effect of cavity stability using numerical approaches, considering the effect of wave direction and peak acceleration on the enlargement of the yield zone around the cave with respect to the static condition. Genis and Aydan (2007) have instead detected a threshold peak acceleration value capable of triggering yield zone enlargement with respect to the static case. However, as specifically concerns the “room and pillar” underground systems, the assessment of the dynamic behaviour of single pillars or internal septa is of crucial importance for sinkhole hazard analysis. In these systems, the septum is a long wall that divides two adjacent cavities, whereas the pillar is a single column in the middle of the cavity. Indeed, these cave structures are rock remains left to support the weight of overburdened material between adjacent underground cavities. Due to the high levels of stress acting at large depths, pillar failure represents a widespread hazard within underground environments (Cook 1976). A limited number of scientific works have been specifically dedicated to the numerical investigation of underground pillar behaviour in static conditions. Mortazavi et al. (2009) first presented a numerical analysis of the rock pillar failure mechanisms in underground openings using both perfectly plastic and brittle constitutive models, comparing the results against field data recorded in Canadian mines. Elmo and Stead (2010) have proposed a hybrid finite-discrete model to describe the behavior of naturally fractured pillars, along with the effect of the implementation of a fracture network approach. Wang et al. (2011) have investigated the mechanism of failure of both serial and parallel pillar models, while Wang and Cai (2021) have shown from a numerical point of view the effect of time-dependent spalling on the progressive failure of pillars. However, few literature studies have investigated the effect of interaction between adjacent cavities. In terms of dynamic behavior, Genis and Aydan (2008) demonstrated with a 3D numerical model that, at low acceleration peak values, tensile failures develop along the cave roof, whereas with larger acceleration values shear failure generates within pillars. Gercek (2005), Conte and Dente (1993), and Landolfi (2013) have highlighted that the presence of cavities at short distances induces larger risk conditions under dynamic loading, with amplification or deamplification of the seismic input in the pillar zone, depending on the wave characteristics and cave geometry. de Lucia et al. (2022) have instead investigated the interaction distances between twin and parallel cavities under dynamic loading conditions, also specifying the minimum distance beyond which interaction is practically negligible. However, the response of single underground pillars under dynamic loading has not been specifically studied so far.

This paper aims to provide a detailed numerical insight into the dynamic behaviour of internal elements of underground caves, such as septa or pillars, in accordance with the regional seismicity recorded in the study area and the typical mechanical properties of the soft rock forming the caves. The study particularly addresses the behavior of septa/pillars, focusing on the effects of the corresponding geometrical features.

In this paper, an overview of the geotechnical properties of the soft rocks involved in these processes is first described, followed by a discussion on the assumptions made in the numerical finite element models.

For these septum/pillar geometries, the evolving stress–strain field is observed under dynamic inputs, both under 2-D and 3-D configurations. To quantitatively assess the septum stability, a septum stability index SSI is also proposed according to the results of 2D analyses. The SSI has been compared with the factor of safety FoS proposed by de Lucia et al. (2022), to highlight potentialities and differences in the evaluation of the stability of such structures under dynamic loading.

The numerical results are pointed out and discussed in detail to derive some indications of the dynamic behaviour of the internal rock elements under study.

2 Geotechnical properties of soft calcarenite rock

A huge number of underground cavities, of interest in the present work, have been excavated within soft calcarenites, which outcrop extensively in the Apulian region. In dry conditions, the uniaxial compressive strength (UCS) of calcarenites typically ranges from 1.2 to 6 MPa, while in saturated conditions, it is generally half of dry values (Andriani and Walsh 2010). The rocks with an UCS in the range of 0.25–25 MPa can be classified as soft to extremely soft rocks (ISRM 1979; Dobereiner and De Freitas 1986). This factor, along with strength values sufficient to support the excavation, has favoured the ease of extraction and workability of such rock. The latter was extensively mined for use as a building material, resulting in the creation of underground cavities that were abandoned over time. The calcarenite units cropping out in Apulia can be divided into three main geological groups:

  • Upper Miocene calcarenites,

  • Upper Pliocene-Lower Pliocene calcarenites,

  • Middle Pleistocene-Upper Pleistocene calcarenites (named as Gravina calcarenites, Andriani e Walsh 2010).

Most of the Apulian underground cavities are excavated in the central part of the region, where the Gravina calcarenite crops out. Thus, the mechanical characterization of this soft rock is here described. The Gravina calcarenites have a low and variable dry uniaxial compression strength, generally between 1 and 9 MPa (Coviello et al. 2005; Andriani and Walsh 2010; Ciantia et al. 2015a, b). The tensile strength generally varies in the range of 0.2–2.0 MPa (Coviello et al. 2005; Ciantia et al. 2015a, b), so the ratio between uniaxial compression and tensile strength is approximately equal to 3–12. The mechanical strength is generally observed to depend on porosity, which varies between 26 and 55% (Andriani and Walsh 2010; Andriani et al. 2019) The degree of saturation has a significant influence on the rock strength: under saturation conditions, the uniaxial compression strength strongly decreases down to 50%, or in some cases even more, in the range of 0.3–6 MPa (Coviello et al. 2005; Andriani and Walsh 2010; Ciantia et al. 2015a, b).

The petrographic and structural features of this rock formation are highly variable. Such large rock heterogeneity, taking place at all scales of analysis, significantly influences the corresponding mechanical response. Also, the weathering processes due to chemical and physical agents and the drying-wetting cycles induce a reduction in the mechanical strength of the calcarenite (Andriani and Walsh 2010, Lollino 2021).

Based on the well-known technique of linearization of the Hoek–Brown failure criterion, considering a range of the cavity depth between 2 and 20 m and the aforementioned parameter for the calcarenite (with intact rock condition and mi = 4), the parameters of the Mohr–Coulomb failure criterion are set for the present study. In particular, rock cohesion, c′, can be assumed in a range of 100–200 kPa, according to the different degrees of rock cementation. Friction angle, ϕ′, can be considered to vary between 30 and 35°, whereas the dilation angle has been imposed equal to ψ’ = 0° to reduce the risk of strength overestimation. From the results of the uniaxial compression test, Young’s modulus, E′, ranges between 100 and 300 MPa at large strains. From typical values of shear wave velocity measured in calcarenite deposits through down-hole tests, instead, E′ values larger than 1000 MPa have been indirectly calculated at small and very small strain levels (Gallipoli and Lupo 2012). A Poisson ratio equal to 0.3 can be set for the calcarenite rock.

3 Finite-element modelling

3.1 Numerical modelling assumptions

To evaluate the influence of a seismic input on the shear behaviour of a septum within an underground cave system, a large number of parametric analyses were conducted on ideal schemes of rectangular twin cavities. The analyses were carried out using finite-element simulations performed with PLAXIS 2D software (Brinkgreve et al. 2022a). In most real cases, the site conditions are characterized by a homogeneous half-space in any direction with a horizontal ground surface. In general, plane-strain conditions, taking into account simplified cavity geometric features, can be assumed to be representative of the most frequent in situ conditions of the cavity systems existing in the Apulian region. Thus, the length of the caves, i.e. the size in the 3rd direction—normal to the domain plane, is frequently larger than those in the transverse plane. Under such conditions, the use of two-dimensional simulations is justified. However, to investigate also the behaviour of a single pillar, which is supposed to be more unstable than the septum configuration, three-dimensional FE analyses were carried out and the corresponding results have been compared with those obtained from 2D analyses. The 3D numerical simulations were performed with PLAXIS 3D (Brinkgreve et al. 2022b).

All the numerical domains have a width of 200 m and a height of 86.5 m. These dimensions are large enough to reduce any boundary effects during the static phase based on sensitivity analyses and previous studies conducted by Perrotti et al. (2018). In particular, the domain height is such that the first natural frequency of the site domain (\({V}_{S}/4H\)) is similar to the predominant frequency of the seismic input, described thereafter. So that the most unfavourable condition for stability is selected.

The discretization mesh of the finite-element model is performed using 15-node triangular elements. To reach a sufficient calculation accuracy in the interested area, a finer mesh was chosen in the area around the cave and the septum, while, a coarser mesh was set in the area far from the septum. In detail, a coarseness coefficient of 0.2 was set in the septum area, compared to a value of 1 for distant areas. The latter indicates the triangular element size for a finer mesh with respect to the target element size, which is a function of the outer geometry dimensions of the model. However, to respect the propagation condition of the highest frequencies of the seismic input, the maximum size of the mesh triangular element was set to comply with the well-known Kuhlemeyer and Lysmer relationship (1973):

$${h}_{{\text{max}}}=\frac{{V}_{s}}{\begin{array}{c}\left(8\div 10\right){f}_{{\text{max}}}\\ \end{array}}$$
(1)

where Vs is the average value of the velocity of the shear wave measured in the calcarenite rock (i.e. 900 m/s) and \({f}_{{\text{max}}}\) is the maximum frequency of the input signal (i.e. 10 Hz). Thus, \({h}_{{\text{max}}}\) could vary from 9 to 11.25 m, accordingly.

The numerical analyses were performed according to different calculation stages:

  • First stage: the initialization of the geostatic stress state, within the rock mass domain, was done by assigning a coefficient of earth pressure at rest, \({K}_{0}\), equal to 1. Therefore, the vertical stresses generated were calculated through gravitational forces, while the horizontal stresses were calculated using the specified value of \({K}_{0}\). Higher values of \({K}_{0}\) have generally been observed to increase the stability of cavity roofs, whereas, for shallow rock masses, lower values of \({K}_{0}\) are unlikely (Brown and Hoek 1978; de Vallejo and Hijazo 2008). Furthermore, the studied area is characterised by not very high tectonic thrusts, which does not justify higher \({K}_{0}\) values.

  • Intermediate stages: the excavation process was performed by deactivating the internal clusters of the cavity through a static plastic phase, according to the chosen constitutive model.

  • Final stage: the seismic shaking was assigned with the application of an acceleration time history at the bottom of the model.

The boundary conditions applied to the domain for each stage were as follows:

  • For the static stage, the vertical boundaries had horizontal zero displacements (\({u}_{x}=0\)) and were left free in the vertical direction, while the mesh bottom was fixed in x–y-directions (\({u}_{x}={u}_{y}=0\)) and the ground surface had free displacement.

  • For the dynamic stage, free-field boundary conditions were applied to the lateral sides of the model to simulate the waves distancing in the far domain. These conditions consist of two dashpots, in the normal and shear directions, which allow reducing the reflection of the incident waves on the boundary. Instead, a compliant base consisting of viscous boundary and x-direction prescribed displacements was implemented at the model base to simulate wave propagation into deep soil with minimal reflection at the bottom boundary, for cases with non-high impedance contrast. This condition combines viscous boundaries and line-prescribed displacements. While the viscous boundaries allow the downward propagating compressive and shear waves to be absorbed. the input signal is transferred to the main domain by applying equivalent normal and shear forces. (Brinkgreve et al. 2022a, b).

For the dynamic phase, in PLAXIS 2D, an implicit integration scheme is used according to the Newmark formulation (1959). The displacement and the velocity at the time \(t+\Delta t\) depend on the coefficients \({\alpha }_{N}\) and \({\beta }_{N}\), which controls the accuracy of the numerical time integration. The coefficients must respect the modified Newmark condition proposed by Hilber et al. (1977) to obtain a stable solution:

$${\alpha }_{N}=\frac{{\left(1+\gamma \right)}^{2}}{4} ; \quad{\beta }_{N}=\frac{1}{2}+ \gamma$$
(2)

where γ represents a numerical dissipation parameter, which varies between 0 and 1/3. In this study, no numerical damping is considered (\(\gamma = 0\)), since a dissipation related to the soil damping has been assigned (i.e. Rayleigh formulation). When \(\gamma = 0\) the coefficients \({\alpha }_{N}\) and \({\beta }_{N}\) coincide with those of the original Newmark method (\({\alpha }_{N}\)= 0.25 and \({\beta }_{N}\)= 0.5).

As regards the constitutive model of the calcarenite rock, an elastic-perfectly plastic model with the Mohr–Coulomb failure envelope was chosen. This choice is motivated by the lack of experimental data in the literature concerning the dynamic characterization of these rocky soils under study, with an accurate description of the non-linear dynamic and cyclic behavior. However, with the aim of studying the failure mechanism involving such structures, an elastic-perfectly plastic model can accurately represent the ultimate state behaviour of a rock mass exposed to relatively strong motions, similar to those assumed in this work. In addition, particular attention must be given when using such a constitutive model in dynamic analyses, selecting stiffness values that accurately predict the shear wave velocity of the soils. This leads to an overestimation of shear stiffness values that remain constant until failure. For this purpose, in the analyses presented here, a Rayleigh damping formulation has been employed to account for dissipation in elastic behavior and wave propagation, simulating hysteric energy dissipation in the soil. To avoid overestimation of the rock shear strength, a non-associated flow rule (\(\Psi {\prime}\hspace{0.17em}\)= 0°) was adopted (Vermeer and De Borst 1984).

In the Rayleigh formulation, the damping matrix is defined as follows:

$$\left[C\right]={\alpha }_{R}\left[M\right]+{\beta }_{R}\left[K\right]$$
(3)

where the damping matrix \([C]\) is a linear combination of the mass matrix \([M]\) and the stiffness matrix \([K]\), with frequency-dependent coefficients \({\alpha }_{R}\) and \({\beta }_{R}\). These coefficients have been calculated following the relationship proposed by Chopra (1995) and Clough and Penzien (1995):

$$\left\{\begin{array}{c}{\alpha }_{R}\\ {\beta }_{R}\end{array}\right\}=\frac{2D}{{\omega }_{m}+{\omega }_{n}}\left\{\begin{array}{c}{\omega }_{m}{\omega }_{n}\\ 1\end{array}\right\}$$
(4)

where \(D\) is the damping associated with the two angular frequencies (\({\omega }_{m} , {\omega }_{n}\)) that are related to the target frequency interval [\({f}_{m}- {f}_{n}\)]. The selection of the target frequency interval is a delicate stage of the model calibration: different references can be found in literature (Park and Hashash 2004; Amorosi et al. 2010; Tsai et al. 2014). However, a small damping rate was adopted in the presented study, for which the frequency interval is selected in correspondence with the greater energy content of the Fourier spectra at different depths. For the calcarenite rock, the damping was assumed constant and equal to 0.5%, as generally chosen in the literature (Baranello et al. 2003; Fierro et al. 2020).

All the parameters used in the numerical analyses are summarized in Table 1.

Table 1 Parameters used in the numerical analyses
Full size table

A real accelerogram, consistent with the Apulian historical seismicity, was selected to assess the dynamic effects on the septum/pillar between two adjacent cavities. The earthquake is characterized by a peak ground acceleration PGA = 0.2 g, a predominant frequency fp = 2.6 Hz, and a low-pass filter with a cut-off frequency of 10 Hz, as shown in Fig. 1. For the dynamic analysis, the maximum time step should be set by following Eq. (5) to obtain reliable numerical results (Brinkgreve et al. 2022a, b):

$$\Delta {t}_{{\text{max}}}=\frac{{l}_{{\text{min}}}}{{V}_{S}}$$
(5)

where \({l}_{{\text{min}}}\) is the minimum distance between the apexes of the triangular elements in all the defined mesh. To fulfill the constraint, the analyses were performed by setting a time step equal to 0.01 s.

Fig. 1
figure 1

a Time history of accelerations applied at the model base b Fourier spectra

Full size image

To validate the 2D model geometry, the adopted boundary conditions and the viscous formulation, a local seismic response analysis of the 2D elasto-plastic model without cavity was compared with those obtained from an equivalent linear analysis with STRATA software (Kottke et al. 2019). The two analyses highlighted the same response in terms of accelerogram at the ground surface, as shown in Fig. 2.

Fig. 2
figure 2

Validation procedure in terms of accelerogram at the ground surface between 2D elasto-plastic analysis without cavity (PLAXIS 2D) and 1D equivalent linear analysis (STRATA)

Full size image

3D finite-element analyses were also carried out on the most relevant cases treated in the 2D numerical analyses. The 3D model was defined by a parallelepiped volume characterized by a width and height equal to that of 2D analyses (i.e. 200 and 86.5 m respectively), while the length was set equal to 200 m. Also, in this case, the model size is much larger than the cavity dimension to exclude boundary effects. The 3D analyses were carried out using the same boundary conditions and the same three calculation steps set for the 2D analyses. The dynamic phase in the 3D model fulfills the mesh condition expressed by Eq. 1. An elastic perfectly plastic constitutive model with the Mohr–Coulomb failure criterion was chosen to simulate the calcarenite mechanical behaviour and the parameters assigned are reported in Table 1, while the viscous damping follows the Rayleigh formulation as previously described. In the following, 2D analyses will explicitly refer to septa (i.e. long wall dividing adjacent cavities), in accordance with plane-strain conditions, whereas 3D analyses will investigate the behaviour of both pillars and septa.

3.2 2D finite-element modelling results

Two-dimensional parametric analyses were carried out to explore the influence of the dynamic input on the shear behaviour of a septum, assuming different septum slenderness conditions. The 2D finite-element study was proposed for ideal schemes of the septum, as shown in Fig. 3.

Fig. 3
figure 3

Geometric scheme of the septum in 2D configuration

Full size image

In particular, keeping fixed the cavity width \(L\) and the roof thickness t, the width \(i\), the height \(h\), and the slenderness ratio \(i/h\) of the septum were singularly varied. The width of the septum is supposed to be smaller than the interaction distance between adjacent cavities under dynamic loading, as determined by de Lucia et al. (2022). In fact, the cited study identified the maximum thickness of the septum for the two cavities to interact under dynamic loading. This width is such that for thicker septa, no interaction in the septum is observed and the two cavities behave independently under dynamic conditions. For septum widths i less than approximately 1.5L, where L is the length of the cavity, dynamic interaction is present for both shallow and deep cavities. Therefore, the septum widths were considered to be well below this threshold. Based on typical values observed in situ and reported in the literature, the slenderness ratio \(i/h\) was varied in the range of 0.3–1.7.

In detail, the dynamic influence was investigated for the following septum configurations:

  1. a)

    Shallow septum (\(t=6 \,\text{m}\)) dividing wide cavities (\(L=15-20 \,\text{m}\))

  2. b)

    Deep septum (\(t=20 \,\text{m}\)) dividing wide cavities (\(L=15- 20 \,\text{m}\))

  3. c)

    Deep septum (\(t=20 \,\text{m}\)) dividing narrow cavities (\(L=5 \,\text{m}\))

In the hypothesis of shallow septum dividing wide cavities (configuration a), assuming the mechanical properties reported in Table 1 (c′=200 kPa, \(\Phi \hspace{0.17em}\)= 35°, Vs=900 m/s), the main results are presented in Fig. 4 in terms of plastic shear points.

Fig. 4
figure 4

Plastic shear points (2D analysis) for the configuration of shallow septum (t = 6 m) dividing wide cavities (L = 15–20 m)—configuration a

Full size image

It should be pointed out that to avoid shear failure in the septum during the static phase, it is necessary to have septa with a width \(i\) larger than or equal to 3 m, regardless of the \(i/h\) ratio. For septa with \(i\hspace{0.17em}\)≥ 3 m, the dynamic stability of the septa is different, depending on the slenderness ratio \(i/h\). In detail:

  • For \(\frac{i}{h}<1.5\), the shear behavior depends on both the slenderness ratio and the absolute value associated with \(i\) and \(h\):

  • For \(\frac{i}{h}\le 0.5\), a threshold value of the septum width \(i\) equal to 4 m is identified: for \(i\) < 4 m, the septum fails already in the static phase; for \(i\) ≥ 4 m the rupture affects the apexes or the inner zones of the septum in the static phase, and then, in the dynamic phase, it propagates in the septum—less diffusely as the septum thickness \(i\) increases—and also in the roof zone.

  • For \(0.5 <\frac{i}{h} < 1. 5\) again a threshold value is identified for the septum width of i = 4 m. It distinguishes two different types of septum shear behaviour: for \(i\) < 4 m, in the static phase, the shear failure affects the area of the septum at the apexes (\(i\) values close to 4 m), or along the vertical walls of the septum (\(i\) values close to 3 m) and in the area above the septum starting from the vertices. While, in the dynamic phase, shear failure propagates within both the septum and the roof area. For \(i\) ≥ 4 m, in the static phase, failure involves the area of the roof above the septum starting from the apexes and, in the dynamic phase, it propagates only in the roof vertices and the roof

  • For \(\frac{i}{h}\ge 1.5\), the shear behavior is practically the same, regardless of the absolute value assigned to the geometric features \(i\) and \(h\): in the static phase, the shear failure affects the area of the roof immediately above the septum, starting from the vertices of the septum; in the dynamic phase, the failure propagates mainly in the roof area.

Under the aforementioned assumptions in terms of septum geometrical configuration and mechanical properties of the soil, the septa with a slenderness ratio of \(i/h \le 0.5\) (slender septa) are stable in static conditions, but the slender septa are the most affected by the influence of seismic loading, with a significant generation of plastic points within the septa, according to the typical \(x\)-shaped mechanism.

As regards the case of deep septa sided by wide cavities (configuration b), the main results are presented in Fig. 5.

Fig. 5
figure 5

Plastic shear points (2D analysis) for the configuration of the deep septum (t = 20 m) dividing wide cavities (L = 15–20 m)—configuration b

Full size image

It is possible to notice two different shear behaviours of septa in accordance with the slenderness ratio \(i/h\). In particular:

  • For \(\frac{i}{h} < 1.5\), regardless of the absolute values assigned to the width and the height of the septum, it always exhibits shear failure in the static phase.

  • For \(\frac{i}{h} \ge 1.5\) the shear behavior is always the same: in the static phase, the shear failure points develop mainly along the septum walls, according to spalling mechanisms, whereas, in the dynamic phase, the shear failure propagates significantly within the septum. The dynamic failure also develops within the roof, starting from the outer vertical-cavity walls, according to a less diffuse mechanism as the septum width \(i\) increases.

It comes out that, with the aforementioned assumptions of deep septa and wide cavities, the dynamic instability generated by the chosen seismic motion is observed only for septa with a slenderness ratio \(i/h \ge 1.5\) (squat septa). For such configurations, shear instability starts already in the static phase and becomes significantly more pronounced in the dynamic phase, along with the generation of plastic shear points affecting mainly the septum. For \(i/h < 1.5\), the septum is unstable already in the static phase.

Finally, concerning configuration c, keeping the depth of the roof constant, the cavity width \(L\) is also seen to influence the behaviour of the septum. In fact, for narrow cavities (L = 5 m) the static and dynamic behavior in the septum is different from the case of wider cavities (L = 20 m). The main results for this configuration are presented in Fig. 6.

Fig. 6
figure 6

Plastic shear points (2D analysis) for the configuration of the deep septum (t = 20 m) dividing narrow cavities (L = 5 m)—configuration c

Full size image

The figure shows that:

  • For \(\frac{i}{h} < 1\) and for absolute values of i ≤ 3 m, the shear failure points affect the septum in the static phase and increase in the dynamic phase; whereas, for absolute values of i > 3 m, the septum is stable in the static phase, while, in the dynamic phase, the typical \(x\)-shaped shear failure develops.

  • For \(\frac{i}{h} \ge\) 1, in the static phase, the septum exhibits shear failure points at the vertices, whereas shear failure points develop in the dynamic phase as a function of threshold values of \(i\), which varies according to \(i/h\) ratios. For \(i/h=1\), the threshold value for achieving dynamic shear failure in the septum is i ≤ 4 m, while, for \(\frac{i}{h}=1.5,\) the threshold value is around 3 m.

Similarly, for deep septa dividing narrower cavities, the most unstable geometries for shear failure are associated with slenderness ratios \(i/h < 1\).

The results described above are shown in Fig. 7 in terms of the contours of the deviatoric strain. These plots confirm the results in terms of shear failure points, but also allow a comparison between shallow and deep septa, respectively dividing wide or narrow cavities.

Fig. 7
figure 7

Deviatoric strain contours for geometric configurations (a), (b), and (c) with slenderness ratios of 0.5 and 1.5

Full size image

In particular, with the same values assigned to the septum width and height, i.e. slenderness ratio \(i/h=0.5\), the most severe configuration, with the concentration of large deviatoric strains in the dynamic phase, is observed to be that for the shallow septum dividing wide cavities (configuration a). For the deep septum dividing wide cavities (configuration b) a shear collapse occurs already in the static phase, with deviatoric strain levels that are lower than those achieved for configuration a under dynamic conditions. For the deeper septum with narrower cavities (configuration c), the level of deviatoric strains within the septum under dynamic conditions is comparable with that of the static phase in configuration b.

For slenderness ratio \(i/h=1.5\), the most severe situation is confirmed to be that for deep septum dividing wide cavities (configuration b). For shallow septa (configuration a), a concentration of deviatoric strains involves the roof, whereas, for deep septa flanked by narrow cavities (configuration c), a negligible difference is observed between static and dynamic stages.

An interesting perspective of the strain evolution of the septum under dynamic loading is provided by the plot of vertical settlement at the top of the septum (see point A in Fig. 3) during the seismic loading history as a function of time, as shown in Fig. 8. The figure indicates the different amounts of settlement calculated in the different septum configurations, in accordance with the evolution of strains above described. For all the considered geometrical configurations, a clear increment of the vertical displacement calculated at the end of the loading history is observed for slender septa with respect to thicker ones. Moreover, for shallow septa dividing wide cavities (configuration a), large vertical displacements have been calculated (3–10 mm) for the case of slenderness ratio equal to 0.5 Lower displacement values, approximately equal to 2–4 mm, have been simulated for deep septa dividing wide cavities (configuration b), and even lower than 0.5 mm for deep septa dividing narrow cavities (configuration c).

Fig. 8
figure 8

Vertical displacement at the top of septum against dynamic time, for geometric configurations a, b, and c with slenderness ratios of 0.5 and 1.5

Full size image

In order to assess the stability of the septum from a quantitative point of view, a septum stability index (SSI) was proposed. In detail, the SSI is calculated by considering only the shear stress state of the septum, at the end of both the static and dynamic phases. The index represents the ratio between the average shear-resistant force and the shear-mobilized force of the septum and is obtained analytically according to the following formulas:

$$SSI= \frac{\sum_{e=1}^{n}{{\tau }_{{\text{max}}}}_{e}{A}_{e}}{\sum_{e=1}^{n}{{\tau }_{{\text{mob}}}}_{e}{A}_{e}}$$
(6)
$${{\tau }_{{\text{max}}}}_{e}={\text{avg}}\left({{\tau }_{{\text{max}}}}_{i}\right)$$
(7)
$${{\tau }_{{\text{mob}}}}_{e}={\text{avg}}\left({{\tau }_{{\text{mob}}}}_{i}\right)$$
(8)

with \(i\in e\) and \(e\in {\text{septum}}\)

where \({\tau }_{{\text{max}}}\) is the maximum shear strength obtained when the Mohr’s circle is expanded to contact the Mohr–Coulomb shear failure envelope, while keeping fixed the centre of the Mohr’s circle. Following Eq. (7), \({{\tau }_{{\text{max}}}}_{e}\) is the average value of all 12 values recorded by the stress points i inside the triangular finite element e belonging to the septum, while \({\tau }_{{\text{mob}}}\) is the radius of the Mohr’s circle and represents the maximum value of shear stress. Following Eq. (8), \({{\tau }_{{\text{mob}}}}_{e}\) is the average value of all 12 values recorded by the stress points i inside the triangular finite element e belonging to the septum. \({A}_{e}\) is the area of the triangular element e belonging to the septum, for which \({{\tau }_{{\text{max}}}}_{e}\) and \({{\tau }_{{\text{mob}}}}_{e}\) were computed. The index is calculated by considering only the shear stress state, i.e. assuming that such structures mainly collapse by shear failure.

Figure 9 shows the trends of SSI as a function of septum width i, calculated under both static and dynamic conditions; the figure indicates that the effect of the earthquake on septum stability tends to fade as the septum width increases, with an asymptotic trend. This is true except for configuration b, where the dynamic effect becomes more evident as the width of the septum increases. In fact, this is the case where the septum experiences the largest stress level (thick roof and adjacent wide cavities), so it is close to failure already in the static phase for smaller septum width i. Thus, for configuration b the static stability increases by increasing i, and the effect of the earthquake is more strongly resented.

Fig. 9
figure 9

Static and dynamic SSI as a function of the septum width i, for geometric configurations a, b, and c, with slenderness ratios of 0.5, 1, and 1.5

Full size image

The results shown in Fig. 9, tend to be in good agreement with those reported in Figs. 4, 5, 6, and 7 in terms of shear failure points and deviatoric strains. Depending on the width-to-height ratio and the width of the septum, different trends can be distinguished. In detail:

  • for configuration a, it is observed that the slender septa (\(i/h\le 0.5)\), although statically stable, are most affected by the earthquake loading, decreasingly with increasing i, as already observed from the numerical results in terms of shear plastic points. For squat septa (\(i/h\ge 1.5)\) the dynamic instabilizing effect is instead not captured by the SSI, since the rupture mainly affects the roof. For the intermediate case (\(i/h=1)\) there is a slight dynamic effect on the septa with i < 4 m, whereas for i ≥ 4 m the earthquake does not influence the septum stability since the shear failure again affects only the roof;

  • for configuration b, the SSI confirms the qualitative results previously observed and captures a more instabilizing effect caused by the earthquake, with the septum being stable in the static phase;

  • for configuration c, the dynamic instabilizing effect can be detected depending on the ratio i/h and the septum width i. The SSI trend indicates that for \(i/h=1.5\) the dynamic effect is most evident for i ≤ 3 m, while for \(i/h=1\) the dynamic effect is clearly observed for i ≤ 4 m. In this configuration, for slender septa, no dynamic effect is observed unlike what was evident from the plastic points. In fact, for these geometries, the septum is already very unstable under static conditions.

These results highlight that the simple analysis of the numerical results in terms of failure points does not indicate the proximity of the stress point Mohr’s circle to the failure envelope, since it just represents the points that have reached the failure; the calculation of SSI becomes critical to study the stress state of the septum and the corresponding safety margin against shear failure.

A comparison of the proposed index SSI with the shear factor of safety (FoS) calculated in the pillar, proposed by de Lucia et. al (2022), is presented. The FoS is a ratio of the shear resisting force to the shear mobilized force and is calculated by dividing the thickness of the septum into vertical slices. Figure 10 shows the FoS in the cavity septum area for each slice, as a function of the distance from the central septum axis, and the value of the index SSI, calculated for the static and dynamic phases. The SSI and FoS are valued for a septum with a width i of 5 m, and a slenderness ratio of 0.5, in the three different configurations (a, b, c §3.1).

Fig. 10
figure 10

Static and dynamic SSI and FoS (de Lucia et al. 2022) comparison as a function of the distance from the central septum axis

Full size image

The results highlight that the SSI is almost an average value of the FoS calculated for each septum slice, and it is always lower than the maximum value of FoS. Therefore, if the FoS provides an accurate and local indication of the area most prone to failure, the SSI provides a rough indication of the global safety margin of the septum against shear failure. Accordingly, SSI and FoS are complementary tools in the evaluation of septum shear stability.

3.3 Sensitivity analyses

A sensitivity analysis was conducted to evaluate the parameters of greatest influence on the dynamic shear behaviour of the septum. In detail, several parametric analyses were performed on a reference geometric configuration, varying the calcarenite shear wave velocity Vs, the calcarenite cohesion, and the earthquake applied to the model base in terms of PGA and frequency content.

Unchanged from the other properties of the analysis, the velocity Vs was considered to be 500 m/s, 700 m/s and 900 m/s. From the results in Fig. 11a, in terms of SSI as Vs changes, it is observed that Vs does not influence static SSI, while SSI decreases slightly with Vs increases. However, the Vs influence on the septum dynamic behaviour can be considered negligible.

Fig. 11
figure 11

a Static and dynamic SSI as a function of calcarenite versus b Static and dynamic SSI as a function of calcarenite cohesion

Full size image

Figure 11b shows the results of the parametric analysis conducted by varying the cohesion (i.e. c′ = 170 kPa; 185 kPa; 200 kPa). It can be seen that as the cohesion increases, the septum becomes more statically stable while the dynamic SSI does not vary, resulting in a greater instabilizing effect due to the dynamic input.

The influence of different seismic inputs was evaluated by subjecting the same cavity geometry to three different accelerograms, shown in Fig. 12a–b. The accelerograms vary in energy content and predominant frequency but have the same PGA. From the results in Fig. 12c, it can be observed that accelerograms with similar frequency content, same PGA, and varying in predominant frequency, have negligible influence on the dynamic stability of the septum geometry considered.

Fig. 12
figure 12

Sensitivity analysis for accelerograms a accelerations time history b Fourier spectrum c Static and dynamic SSI as a function of accelerograms

Full size image

The SSI calculation was also performed for an earthquake out of scale for the typical Apulian seismicity. Such an earthquake was obtained by scaling the input accelerogram so that only the intensity of the earthquake varied, increasing the PGA to 0.3 g and leaving the other characteristics unchanged. The results obtained are shown in Fig. 13 in terms of SSI values against the septum width. It can be seen that the SSI well captures the effect of a more intense seismic loading.

Fig. 13
figure 13

Static and dynamic SSI as a function of the septum width i, for geometric configurations a, and c, varying the earthquake PGA

Full size image

3.4 Three-dimensional effects

Three-dimensional finite element analyses were also carried out to investigate the role of three-dimensionality in the pillar response under dynamic actions. In particular, three different models were investigated (Fig. 14): 3D model with a long septum dividing twin cavities (a), 3D model with a single pillar in the middle of a shallow cave (b), and 3D model with a single pillar in the middle of a deep cave (c). The geometric characteristics of the different models are summarized below:

  1. a)

    L = 20 m, t = 6 m, i = 4 m, h = 8 m, i/h = 0.5, for the septum configuration;

  2. b)

    L = 20 m, t = 6 m, i = 4 m, h = 8 m, i/h = 0.5, for the pillar configuration;

  3. c)

    L = 20 m, t = 20 m, i = 7.5 m, h = 5 m, i/h = 1.5, for the pillar configuration.

Fig. 14
figure 14

Three-dimensional models of septum and pillars assumed in the study: a long septum dividing twin cavities; b single pillar in the middle of a shallow cave; c single pillar in the middle of a deep cave

Full size image

All the assumptions regarding the seismic loading, the choice of the constitutive model, the problem geometry, the discretization mesh, and the boundary conditions made in the 3D analysis have been kept consistent with those used in the 2D analyses. For all the 3D analyses, the seismic loading has been applied as a horizontal wave acting in the x-direction.

In particular, the numerical results obtained from the 3D model corresponding to configuration a in Fig. 14a should be compared with the results derived from the 2D analysis with the same geometrical features, as reported in Fig. 4. In this case, the 3D model results in terms of plastic points indicate a shear failure mechanism very similar to that of the corresponding 2D analysis, with x-shaped failures developing already in the static phase, which later on propagates within the roof during the dynamic phase from the apexes of the two caves (Fig. 15). However, dealing specifically with the dynamic phase, the 3D numerical results seem to show plastic failures within the cave roof less diffuse than for the case of 2D analysis. Such a difference should be presumably the effect of different stress redistribution in the 3D model, which obviously should be considered more consistent with reality than 2D analysis. This effect is confirmed by observing the contours of deviatoric strains calculated in the 3D model under dynamic conditions reported in Fig. 16, which indicates a cross-section of the 3D model with a strain level lower than that simulated in the 2D analysis (Fig. 7a).

Fig. 15
figure 15

3-D model results in terms of plastic points calculated for configuration a, with long septum dividing twin cavities (t = 6 m, L = 20 m; i/h = 0.5)

Full size image
Fig. 16
figure 16

Contours of deviatoric strains under dynamic conditions calculated in a cross-section of the 3D model (configuration a)

Full size image

As regards configuration b, which is characterized by a single pillar in the middle of a shallow cavity (Fig. 14b), the results of the 3D model reported in Fig. 17 show, in the static phase, diffusion of plastic points both within the pillar but also in the roof in comparison with the 2D model (Fig. 4). This is always due to the stress redistribution which is quite higher concentrated in the roof compared to the plane strain condition. Deviatoric deformations in the static 3D model are also slightly larger than in the static 2D model (Figs. 7a, 18). Regarding dynamic stability, the 3D model shows an increase in plastic shear points only in the roof and not in the pillar, as the 2D septum model shows. Indeed, the deviatoric strains of the 3D pillar show an increase in strain levels compared to the static analysis results but not as significant as those shown by the 2D model (Figs. 6a, 18). This again confirms some stabilizing effect of the third component perpendicular to the seismic input.

Fig. 17
figure 17

3-D model results in terms of plastic points calculated for configuration b, with a single pillar in the middle of a shallow cave (t = 6 m, L = 20 m; i/h = 0.5)

Full size image
Fig. 18
figure 18

Contours of deviatoric strains under static and dynamic conditions, as calculated in a cross-section of the 3D model (configuration b)

Full size image

The same results are confirmed for configuration c, which is characterized by a single pillar in the middle of a deep cavity (Fig. 13c). The 3D model results indicate in the static phase a strong concentration of shear failures within the pillar and the area above due to the larger concentration of vertical stress with respect to the 2D plane-strain condition (Figs. 5, 18 and 19). Then, in the 3D dynamic phase, shear failure highly propagates within the area above the pillar and not in the pillar. In fact, deviatoric strain levels in the pillar are also slightly higher in the 3D dynamic model than in the 3D static model, but this increase is much less pronounced than in the 2D dynamic configuration, as illustrated by comparing Fig. 20 with Fig. 7b.

Fig. 19
figure 19

3-D model results in terms of plastic points calculated for configuration c, with a single pillar in the middle of a deep cave (t = 20 m, L = 20 m; i/h = 1.5)

Full size image
Fig. 20
figure 20

Contours of deviatoric strains under static and dynamic conditions, as calculated in a cross-section of the 3D model (configuration c)

Full size image

4 Results discussion

The first objective of the present work is to study the evolution of the stress–strain field under dynamic conditions for ideal septum geometries in 2D configurations by means of parametric finite element analyses, varying the geometrical feature of the septum. For this purpose, the variation of plastic shear points, deviatoric strains and vertical displacements at the top of the septum were observed at the end of the static and dynamic phases. The numerical results obtained under 2D conditions have shown that for shallow septa, with i/h < 1.5, the shear behaviour depends on both the slenderness ratio and the absolute value of \(i\) and \(h\), whereas, for i/h > 1.5, the shear behaviour is practically the same, regardless of the absolute value assigned to the geometric parameters \(i\) and \(h\). The slender septa have resulted to be the most affected by the influence of seismic loading, with a significant generation of plastic points within the septa, according to the typical \(x\)-shaped mechanism. For deep septa, the dynamic instability generated by the seismic motion is observed only for septa with a slenderness ratio i/h > 1.5. At large depths, the cavity width is also seen to influence the behaviour of the septum, due to a different stress concentration transferred in the septum. For all the geometrical configurations examined, a clear increment of the vertical displacement calculated at the apex of the septa is observed for slender septa, with respect to thicker ones.

A stability index SSI was also suggested in 2D models to assess from a quantitative point of view the stability of the septum under both static and dynamic conditions, so that the effect of the dynamic loading can be clearly detected. The results expressed in terms of stability index are in good agreement with the failure behaviour observed in terms of shear failure points and deviatoric strains.

The parameters, that most influence the behaviour of the septum, are the cohesion of the calcarenite and PGA of the dynamic input. The comparison of the proposed index SSI with the shear-based Factor of Safety (FoS), as introduced by de Lucia et al. (2022), indicates that the SSI and FoS are complementary tools in the evaluation of the septum shear stability.

Three-dimensional FE analyses were also carried out on the most relevant cases treated in 2D analyses. The results have highlighted rather similar results under static conditions, while a stabilizing effect in the pillar of the third component perpendicular to the earthquake is evident in the 3D dynamic phases compared with the 2D case, mainly in terms of deviatoric strains.

5 Concluding remarks

A numerical investigation of the effects of dynamic loading on pillars and septa within underground caves in soft rocks has been proposed in this work. To this purpose, both two-dimensional analyses representing plane-strain configurations of septa dividing adjacent cavities and three-dimensional analyses simulating single pillars in the centre of underground caves have been carried out. Typical values of the static and dynamic properties of soft carbonate rocks outcropping in the Apulia region (South-eastern Italy) have been chosen. Geometrical features largely used to excavate underground environments in the past have been considered for the parametric numerical study. An elastic-perfectly plastic model with the Mohr–Coulomb failure envelope, with a soil damping dissipation according to the Rayleigh formulation, was chosen. It should be also emphasised that the simple constitutive model used in this study allows catching reasonably the ultimate state behaviour of the rock masses subjected to strong motion, despite the limitation of not taking into account the dynamic non-linearity of the materials involved.

For shallow cavities, slender septa were identified as the most susceptible to the impact of seismic loading. In the case of deep cavities, dynamic instability is only observed for rather squat septa, and the cavity width also plays a role in influencing the dynamic behavior. A quantitative assessment of septum stability was conducted through the introduction of a stability index in 2D models. This proposed index facilitated the quantification of the parameters that have the most significant influence on the dynamic stability of these structures. Furthermore, three-dimensional analyses revealed a stabilizing effect in the pillar due to the stress component perpendicular to the earthquake.

Ideal and simplified geometries were taken into consideration in the current work, although different cavity shapes could result in variable stress releases on the septum/pillar, which would require further investigation. Moreover, in order to understand in a more comprehensive way how septa and pillars behave when subjected to seismic wave motion, different accelerometric input data should be taken into account, this being the objective of future scientific efforts.