2D numerical analysis of the seismic response of a karst rock mass: importance of underground caves and geostructural details

The presence of discontinuities and faults inside a rock mass has a significant influence on the way of propagation of seismic waves and on their amplification effects on the ground level [1] [2]. The situation will further complicate if one or more underground caves (anthropogenic or natural) are present in the rock mass [3], [4], [5], [6], [7]. The site effect evaluation appears to be complex because of these conditions, so that it cannot be conducted with one-dimensional numerical code. On the contrary, two-dimensional or three-dimensional software, such as UDEC [8], are more appropriate to conduct this kind of analysis because it is possible to model both the presence of an underground cave and the discontinuous and anisotropic characteristics of the geotechnical model. Previous studies mainly referred to models with homogeneous and isotropic rocks, without considering rock mass with faults and underground karstic caves. The local seismic response (LSR) analysis of these models (especially in the professional practice) is usually conducted on 1D models, in which only the stratigraphic variation is considered. In the past few years, studies on LSR have been increasing significantly in literature. They aim to evaluate the influence of reliefs, valleys, faults and underground caves [3], [4], [5], [6], [7], [9], as well as the soil–structure interaction [10], [11]. However, even more recent studies are based on homogeneous and isotropic subsoil models and are mainly conducted with finite element method (FEM), finite difference method (FDM), and boundary element method (BEM) calculation software which are not suitable in the case of discretely fractured and/or stratified rock masses, unlike the discontinuous deformation analysis (DDA) [12], [13] and distinct element method (DEM), the latter used in the present case study [14], [15]. Studies about LSR on real and/or parametric cases [6], [7] were conducted on homogeneous and isotropic subsoil, varying the diameter and the depth of the cave. These studies highlighted the strong influence of these two variables over the amplification factor (FA), mostly in the frequency domain (Fig. 1). Overall, a decrease in the FA factor is observed as the cave depth (H) increases and, considering the same depth, as the cave diameter (D) decreases. Moreover, previous studies underline that diameter’s influence on LSR decreases dramatically for D ≤ 5.6 m, that is, for diameters comparable to case study. Furthermore, the effects of amplification/deamplification of the ground motion above a cave are significantly different, depending on whether the V_p or V_s volume waves hit the cave. In the case of V_s waves, the amplification effects on the surface are concentrated on the cave sides, while a deamplification occurs above the cavity. On the other hand, in the case of V_p waves, the amplification of the ground motion is concentrated at the centre of the cavern and decrease sharply, moving to its boundaries.

Figure 1

The UDEC numerical program was used to simulate the LSR of an anisotropic and discontinuous geotechnical model of a real case study, in order to analyse the effects of the simultaneous presence of underground caves, discontinuities and lithostratigraphic variations.

The case study area is part of the extensive Apulian karst area of the Murge. In particular, it is located in the northwest neighbourhood of the town of Turi (Bari, Italy). From the geological and stratigraphic point of view, the rocks that constitute the substratum can be traced back to the ‘Calcari di Bari’ Mesozoic formation. This is a succession of whitish micritic limestones, dolomite limestones and greyish dolomite with local small inclusions and/or grey leading to black coloured laminations (Fig. 2). These limestones have a thickness of several hundred metres. In the study area, they are stratified with layers of 1-m thickness (standard deviation, 0.5 m), position D = 330°, i = 5° (standard deviation, 1°) and I = 60°. The tectonic structures found in the area are attributable to two phases, characterised by different tension fields. The first phase had a compressive nature and generated a smooth anticlinal fold, whose axis are oriented approximately northwest–southeast. The fold lies a few hundred metres south of the study area which is, therefore, placed on the northeast side of this anticline. On the other hand, the second phase generated two normal faults defining a small graben, with northeast–southwest oriented axis, in which the case study area falls.

Figure 2

However, it must be said that the two normal faults are not directly detectable on the surface, as they are inter-formational and with minimal vertical component of dip-slip.

Their location and orientation stemmed from the correlations between the stratigraphic reconstruction through different borings performed and the electrical tomographies acquired from previous studies. An ephemeral watercourse has set itself, with a south south east - north north west course, along the tectono-karstic depression of the graben, which is locally named as ‘lama’. This depression plays a role in drainage and vehiculation of rainwater. However, the water flow occurs underground (and not on the surface), because of the lithological nature of the rocks in the area and their high permeability (because of fracturing and karst), especially in the most superficial layers of the rock mass. As a result of the carbonatic nature of the rock mass, the landscape morphology of the area is essentially related to the karst phenomenon. It manifests itself both with epigeal morphologies, such as the dolina located mostly south of the study area, and with hypogean morphologies, such as the sinkhole ‘capovento’ located at the west side of the ‘Don Bosco’ nursery school. Other hypogean karst forms have been found a few hundred metres to north north east from the study area, such as the cave of ‘S. Oronzo’ and an adjacent small cave located a few decimetres underneath the surface.

The widespread presence of underground karst caves in the Murge area, the proximity of the S. Oronzo’s cave and the structural arrangement of the area suggest that the ‘capovento’ sinkhole is not isolated, but connected to an underground karstic cavern. This assumption is supported by resistivity values found in previous indirect investigations, provided by the technical department of the municipality of Turi (BA), such as electric tomographies with Wenner–Schlumberger arrays (Fig. 3).

Figure 3

The analysis of the n.1 tomography (the most representative one) showed a low resistivity anomaly, flanked by high resistivity areas, beneath the sinkhole ‘capovento’. The trend of resistivity is confirmed by the graben structural layout, in the centre of which the karstic sinkhole ‘capovento’ develops. The normal faults, positioned on the sides of the graben, are identified by the sharp contrast between high resistivity values on the sides and low values in the centre. Finally, low resistivity values can be explained by two evidences. The first one is the underground water circulation, which is channelled along the tectono-karstic depression ‘lama’, at a depth of about 10–14 m from the surface, along directives from southeast to northwest. The second one is the possible presence of underground karst caves filled with residual clay, known in the area as ‘red soils’. Mechanical surveys, with continuous core drilling and core destruction, were performed with the purpose of clearly identifying the presence of a karstic cave connected to the sinkhole and to potentially reconstruct its geometry (Fig. 4).

Figure 4

In particular, 25 core destruction borings (diameter Φ = 60 mm and depth of 6 m from the ground level) and 2 rotation drills with continuous coring (diameter Φ = 101 mm and depth of 15 m from the ground level) were performed. Both the core destruction borings and the rotation drills with continuous coring had the purpose of finding the hypogeal caves and determining their geometric extension (Fig. 5). Moreover, those with continuous coring have allowed to collect six samples for geotechnical tests (two in the S1 core drill and four in the S2 one) and to evaluate the rock quality designation (RQD) index and its variation along boring verticals (Fig. 6). Furthermore, the planimetric and altimetric characteristics of the nearby S. Oronzo’s cave were considered for the reconstruction of the geometry of the hypogeal caves, which have been identified with the geological surveys.

Figure 5

Figure 6

The S. Oronzo’s cave, indeed, is a karstic one of interlayer type with a planimetric development (91.00 m), which is predominant compared to the altimetric one (20.50 m). Besides, it is articulated on two overlapping karst levels. The deepest level, northwest–southeast oriented, presents cave bottom altitudes between 215.50 and 209.50 m a.s.l., while the more superficial one, oriented north–south and northwest–southeast, presents cave bottom altitudes between 218.00 and 214.00 m a.s.l.

Also, the cave connected to the vertical sinkhole ‘capovento’, adjacent to the ‘Don Bosco’ nursery school, is of the interlayer type, with a prevalently sub-horizontal development on two overlapping karst levels. The development direction of the deepest level is from SE to NW, with a total height of about L = 1.50 m and cave bottom altitudes between 226.90 and 222.30 m a.s.l. The most superficial karst level consists of two ramifications that start from the sinkhole ‘capovento’ and develop one in the northwest direction and the other in the north north east direction. The cave has a height of about L = 4.00 m in correspondence with the sinkhole ‘capovento’ and barely more than L = 1.00 m in correspondence of the two branches that start from it; moreover, the altitudes of the cave bottom are between 225.20 and 234.90 m a.s.l. In the geomechanical model section, used for the LSR evaluation, the compact and not altered limestone was distinguished from the altered and karstic one, depending on the alteration level and the karst increase (Fig. 7).

Figure 7

Considering that the numerical code used for the LSR analysis is a DEM one [8], [14], [15], [17], [18], it was necessary to parameterise the geomechanical model. Therefore, laboratory tests on the six specimens taken from the mechanical borings (conducted at the ‘M. Salvati’ laboratory–DICAR Department, Politecnico di Bari) and geomechanical measurement in situ were conducted. The collected data allowed to classify the rock mass with the Bieniawski RMR’76 indexes [19] and to estimate the friction angle and the cohesion values. Moreover, the GSI values were obtained (Table 1) through the correlations between the geomechanical classification indexes [20], [21]: GSI=RMR'76 (valid when RMR'76≥18){\rm{GSI}} = {\rm{RMR}}'76\,({\rm{valid\ when\ RMR}}'76 \ge 18)

In addition, the Hoek–Brown curvilinear strength criterion [20] was implemented using the values of the mechanical parameters obtained from the laboratory tests and the GSI classification indexes. This criterion considers the strong non-linearity of the strength envelope, which is typical of the rock masses, and estimates the representative parameters, even by means of the Mohr–Coulomb linearisation. The physical–mechanical parameters used for the rocks and the discontinuities of the geomechanical model are summarised in Table 2.

In order to assess the LSR of the case study, the seismic input consists of two accelerograms, respectively, for the X-component and Z-component of seismic input motion. Therefore, according to the seismic hazard of the case study area (Table 3) and the Technical Standards for Italian Constructions (Norme Tecniche per le Costruzioni Italiane [23]), the accelerograms were extracted from the European Strong-motion Database [24] [25] present in the REXEL 3.5 software [26].

These accelerograms are natural and spectrum compatible; moreover, the focal mechanism (normal fault) and the bedrock lithology (carbonatic rocks) were compatible with those of the case study area (Table 4).

Subsequently, the accelerograms previously extracted, which are the most spectrum-compatible and also those that have a scale factor closer to the average value (SF_mean) of the respective combinations of seven earthquakes, were processed with SeismoSignal software [27] to obtain the X-velocity and Z-velocity which were used to define the seismic input for the UDEC numerical simulation (Fig. 8). The seismic waves velocities V_s = 2821 m/s and V_p = 4607 m/s in the rock mass were inferred from previous seismic refraction traverse, provided by the technical department of the municipality of Turi (BA).

Figure 8

To evaluate the LSR of the case study, two discontinuous models with deformable blocks were designed using the UDEC 2D distinct element numerical code. The constitutive model used for rock blocks, discontinuities and faults was the Mohr–Coulomb plasticity one. This model was chosen both because it is one of the most efficient computationally plasticity models and the analysis concerns general rock mechanics issues (e.g. slope stability and underground excavation). Before starting the dynamic analysis, aimed at the LSR evaluation, it was crucial to properly define the boundary conditions and the calculation options: –

presence of a ‘graben’ adjacent to the school building and within which the two overlapping underground caves are located;

The static stress-deformation analysis of the rock mass showed that the model has tensile plasticisation points that tend to cluster on the roof of the upper karst cavity (Fig. 9). These yielding conditions, accompanied by downward displacements of the rock blocks of the karst cavity roof, are even more evident in the results of the dynamic analysis.

Figure 9

According to previous studies [9], [10], and general considerations described above, the LSR simulation of the case study was finally assessed considering two hypotheses: the absence and presence of the two overlapping karst caves. For the simulation, the model was extended for a length of 25–30 m on both the sides of karst caves, the boundary conditions listed above were applied and, in addition, 18 control points were placed, 9 at the base of the model, and the others 9 (aligned with the first ones) on the ground level (Fig. 10). The amplification/deamplification (FA) and the most amplified frequencies on the ground level of the model were analysed to estimate the effects of LSR on the X-accelerations component. More specifically, in each control point on the ground level of the model, the amplification/deamplification (FA) has been calculated using two methods (Table 5): a)

ratio between the values of the maximum X-acceleration (A_Xmax) in the points at the top of the model and in the corresponding ones at the base [9], [10];

Figure 10

The maximum values of FA, for both the calculation hypotheses (with or without the caves) and for the two calculation methods implemented, are found on the right edge of the caves at no. 8 control point, close to the edge of the school building foundations (Fig. 10).

Furthermore, the FA values obtained with both calculation methods showed that, in general, in the control points at the ground level of the model, there is a deamplification of the seismic motion (Table 5). The only exception is represented by the point no. 8 in the hypothesis without caves, where the FA = 1.672 (ratio between the max X-accelerations) would highlight a local effect of seismic amplification. Furthermore, between the three frequency ranges analysed in the simulations, those that are most amplified are 2.0–3.8 Hz (without caves) and 3.8–5.6 Hz (with caves).

Finally, it should be observed that the UDEC dynamic analysis is conducted with a linearly elastic behaviour of the rocks. The values of the seismic FA on the ground level of the model are, therefore, lower or at the most equal to those that would be obtained with a non-linear analysis, because the degradation of the stiffness and shear modules of the rocks, especially of the vegetal and/or landfill soil, increases the seismic amplification on the ground level [2], [15].

Preliminarily to the case study LSR, the stress–strain analysis of the rock mass (Fig. 9) showed that the model presents, already in static conditions, but even more in dynamic ones, tensile plasticisation points that tend to cluster on the roof of the upper karst cavity. Moreover, in dynamic conditions, the downward displacements of the rock blocks of the karst cave roof are compatible with a local collapse of the same. The results of the modelling are in perfect agreement with the real occurrence of a sinkhole that, in the recent past, brought to the light the underlying karst cave in the case study area.

The LSR numerical simulation evaluated by the values of FA obtained with two distinct calculation methods shows a generalised deamplification of seismic motion on the ground level of the model, both in the absence and in the presence of hypogean karst caves. The only exception to this trend is represented by point no. 8 in the hypothesis without caves, where an FA = 1.672 (ratio between the max X-accelerations) should highlight a local effect of seismic amplification. The presence of the lower karst cave, filled by ‘red soils’ clay, and the thin vegetal soil that fills the paleochannel ‘lama’ do not seem to influence the model’s LSR.

The generalised deamplification for the horizontal component of the ground level seismic motion is in agreement with what was reported in the previous studies [1], [15], [28] about the relevant effect of seismic energy absorption operated by the layer discontinuities. Moreover, the effect of seismic deamplification at the top and in axis with the underground caves, already highlighted in the literature for cavities with diameter D > 5.6 m in homogeneous and isotropic substrate models [7], is confirmed by the minimum FA values reached at the point no. 5, with both calculation methods, even for a cave with D < 5.6 m.