1 Introduction

Due to urban growth, renewal, and densification, geotechnical working sites tend to come closer and closer to natural or human preexisting structures. The use of explosives and heavy machinery for civil engineering works such as quarrying, tunneling, or soil compaction generates potentially harmful vibrations at short to medium distances (tens to hundreds of meters, (Crispim 1999; Lacave et al. 2003). This is becoming a major issue for projects taking place in dense areas and underlines the need to reduce vibrations as much as possible. Blast-induced ground vibration level may be contained by several techniques. Trenches and holes, precutting blasts may be used to reduce vibration propagation. The explosive charge per blast hole may be reduced by dividing the bench, and delay time may be used between holes to avoid seismic wave superimposition (Mogi et al. 2000). This requires controlling the firing time of the explosive charge placed in each drilled hole, with the use of precise electrical detonators fired successively with sufficient delay (Hoshino et al. 2000; Rustan 1998). For tunneling or quarrying, numerical softwares are more and more frequently used for blasting excavation optimization (Xie et al. 2018).

In many sedimentary regions, concretionary caves represent a major natural heritage to be preserved (Lacave et al. 2003). It yet also appears to be a vibration-sensitive target (Gauchon et al. 2007). The irreversibility of damages caused to speleothems and the uniqueness of each cave site underlines the need for dedicated vibration studies. Until now, the effects of vibrations on speleothems have been mostly addressed in speleoseismological surveys. Speleoseismology consists in detecting and characterizing speleothem anomalies in caves to study earthquake records (Becker et al. 2006). Growth anomalies, deformation, and/or speleothem failure are notably used to derive some earthquake parameters (e.g. magnitude, azimuth, location, date, peak ground acceleration—PGA or peak ground velocity—PGV) (Bokelmann and Gribovszki 2015; Forti 2001; Forti and Postpischl 1980; Gilli et al. 1999; Méjean et al. 2015; Postpischl et al. 1991; Ferranti et al. 2019; Kagan et al. 2017). The accuracy and robustness of the method yet remains debated as most speleothems—except long and slender ones—may withstand earthquakes (Gilli 2005; Lacave et al. 2004). In addition, many other processes can break speleothems, such as erosion and soil creep, underground glaciers and ice creep, floods, mud and debris flows, frost action, decompression, load and slope movements, and anthropogenic causes (Becker et al. 2006).

Precise characterization of speleothem mechanical properties is of primary importance in speleoseismological studies (Konečnỳ et al. 2015). In the literature, static parameters (geometry, density ρ, ultimate strength) as well as elastic and dynamic properties (Young modulus E, natural frequencies fn, damping ξ) have been provided by both in situ measurements and laboratory experiments (Table 1). Few studies only measured the speleothem density, generally through broken samples weighing and volume estimation. Ultimate compressional (σu,c) and tensile (σu,t) strengths were estimated via uniaxial compression or traction tests, Brazilian tests, or bending tests. Dynamic properties have been measured in situ via geophones or laser interferometers. Significant wandering in mechanical parameters has been observed, both between speleothems of the same site as well as between sites (Bednárik 2009). Such variability may originate from differences in crystallization, porosity, boundary conditions, and/or previous mechanical damage (Bednárik 2009; Szeidovitz et al. 2008). A probabilistic approach has been used by Lacave et al. (2004) to describe variations in tensile rupture stress (σu,t) between samples.

Table 1 Literature review of speleothem mechanical properties

Due to their slenderness, speleothems may be subject to significant dynamic amplification phenomena. The dynamic behavior of speleothems has generally been addressed in the literature using an equivalent static approach (ESA) (Cadorin et al. 2001; Gribovszki et al. 2008, 2013a, 2013b; Konečnỳ et al. 2015, see Table 1). The speleothem is considered as a Euler-Bernoulli beam loaded at mid-span with a force F. This force corresponds to the inertial solicitation, taken equal to the speleothem mass multiplied by the peak ground acceleration (PGA). In ESA, the maximal tensile stress is located at the embedded extremity. The ground horizontal acceleration causing the rupture is then only function of the speleothem geometry (length L, internal and external diameters Dint and Dext), density ρ, and ultimate strength σu. For such loading, the speleothem bending rupture depends on the calcite tensile strength σu,t rather than its compressional resistance σu,c.

A generalization to irregular beams was proposed by Lacave et al. (2004) and Cadorin et al. (2001). However, the dynamic resonance that occurs when the ground motion frequency content matches the speleothem natural frequencies cannot be simulated in the ESA. This technique is hence restricted to the “rigid-body” domain, i.e. when the ground motion occurs at lower frequencies than speleothem natural frequencies (the latter generally ranging between tens of Hz to hundreds of Hz). Such conditions are not systematically fulfilled: slender speleothems may indeed experience significant amplification phenomena under earthquakes, allowing rupture at lower acceleration values (Gribovszki et al. 2017; Lacave et al. 2004). Long, slender, tubular speleothems such as soda straws are thus particularly vulnerable to ground motion, making them suitable targets for speleoseismological studies (Becker et al. 2006; Lacave et al. 2000).

In their speleoseismological study, Lacave et al. (2004) adapted Eurocode 8 elastic response spectra. For speleothems with fundamental natural frequency f0 lying in the seismic bandwidth (< 25 Hz), a 4.5 flat amplification factor was applied. This plateau amplification may yet be conservative or unconservative, depending on speleothem mechanical properties as well as ground motion spectrum. In addition, previous studies used only the peak ground acceleration (PGA) to compute the maximal bending stress undergone by speleothems (Cadorin et al. 2001; Gribovszki et al. 2017; Lacave et al. 2004). In contrast, the acceleration full waveform has never been used, although it contains the whole ground motion information. To our knowledge, the study of vibrations induced by geotechnical engineering works located close to a concretionary cave has been addressed only by Lacave et al. (2003). Their approach yet does not simulate entirely the dynamic amplification that speleothems may experience.

Our study focuses on blasting excavation works conducted to restore the pedestrian access to the Choranche cave (Vercors, France). The rock blasting takes place only tens of meters away to very delicate, hollow tubular soda straw speleothems. We first characterize soda straw dynamic parameters and ultimate strength thanks to laboratory and in situ experiments. We then carry out blast tests to record ground motion along the cave wall and characterize the wave attenuation with distance. We use a dynamic finite element modeling (FEM) to estimate the state of stress in the soda straws subjected to the recorded ground motion. Finally, the admissible peak particle velocity (PPV) along the cave wall was estimated, as well as the maximal instantaneous charge of explosive (MIC) as a function of distance.

2 Site description

The study site is located in the Isère department (southeast France), in the Vercors Pre-Alps Mesozoic massif (Fig. 1a). The Choranche cave consists two main rooms (lake and cathedral, Fig. 1 b and c) and several halls and galleries. The cave is included in a > 50 km long karstic network called Gournier-Coufin-Chevaline (Delannoy 1986) and located about 350 m below the surface. This network develops at the interface between massive Urgonian limestone (lower Cretaceous) and underlying impervious marls and calcareous marls from Hauterivian (Delannoy et al. 2009). The cave hosts many speleothems, including very long (up to 3 m), thin, translucent soda straws (Fig. 1b). These hollow cylinders get longer via successive crystallization of calcite rings around water drops hanging at the speleothem’s end (Fairchild and Baker 2012). The mean annual growth rate ranges between 0.09 and 0.92 mm/year (Perrette and Jaillet 2010) so that the longest Choranche soda straws may date from 3200 to 33,000 BP. Most of the soda straws concentrate in the lake room and the soda straw hall close to the cave exit (Fig. 1 b and c).

Fig. 1
figure 1

a Location map. Choranche site is shown with a black star. b Picture taken inside the cave (lake room). Soda straw speleothems range mostly between 1 to 3 m in the picture. c Map of the site. The cave wall is delineated with the thick black line, the rounded side pointing toward the rock mass. The rock mass is shown in white. The lake is shown in blue. Pedestrian pathways are shown in light gray. The soda straws (red dotted area) are located about tens of meters away from the blasting excavation work site (black hatched area). Seismic sensors (GPa, GPb, GPc) installed on the internal cave wall are shown with green squares. Tests shot locations (A, B, C) are shown with orange pentagons. The seismic profile (PS) is shown with the orange segment

This concretional cave represent an exceptional geological inheritance (Delannoy et al. 2009; Gauchon et al. 2007; Reynard et al. 2011). The cave site is opened to the public throughout the year and welcomes about 100,000 visitors per year (Garnier, pers. comm. 2018). A narrow, hanging pedestrian pathway located at the cliff toe gives access to the cave (Fig. 1c). Rehabilitation works were required, consisting in moving the pedestrian pathway closer to the rock wall along about 100 m. As the rock site is not accessible for heavy machinery, explosives were planned for rock slope cutting.

Two drill cores (15 and 50 m) and a geological penetrating radar (GPR) campaign were carried out during fall 2017 and winter 2017–2018 (not shown here). They revealed a major subvertical fracture (N175/subvertical) striking roughly parallel to both the outward cliff and the cave wall (Fig. 1c). The rock mass at depth appeared massive, with a rock quality designation (RQD) > 90% and no significant reflector. In contrast, the external rock skin (3 to 6 m in depth) showed intense fracturing. Such rock mass decoupling and fracturing was expected to significantly reduce the ground motion level into the cave (Kumar et al. 2016). The distance between the blasts and the soda straw ranges about tens of meters (Fig. 1c).

3 Soda straw characterization

3.1 Dynamic properties

Fifteen samples of soda straws were collected into the cave, mostly found lying onto the lake bottom (dataset D1). We measured their length, diameter, and wall thickness using a Vernier caliper. Most of the soda straws show a perfectly regular, cylindrical shape, although a few of them contain enlarged cross-sections. The length of samples ranges between a few centimeters up to about 40 cm. Longer soda straws cannot be found intact in the cave as they are probably smashed into small bits when impacting the ground or the lake. The mean external diameter Dext for regular-shaped soda straws is 4.5 × 10−3 m and average calcite wall thickness is 5 × 10−4 m.

The speleothem density was characterized for a larger speleothem sample coming from the same Choranche cave. The speleothem density ρ = 2500 kg m−3 was computed from sample weighing with a scale and volume estimation from water displacement method. This value falls in the upper range compared with previous studies, probably because of reduced low porosity in soda straws. It agrees very well with values taken for soda straws by Lacave et al. (2000).

We carried out dynamic testing on 4 soda straw samples in the laboratory (dataset CHO_D1) and on 1 soda straw in situ (dataset CHO_D3, Table 2). The length of CHO_D3.01 could not be measured on the field because of many surrounding, fragile soda straw preventing access to the cave ceiling. Due to instrument availability issues, we used an Ometron VH300+ Laser Doppler Vibrometer Type 8329 (www.bksv.com) for the lab experiment and a Polytec PDV-100 (www.polytec.com) for the in situ measurements. Figure 2 a shows the laboratory experimental setup. The laser beam points at the speleothem top, which was painted in metallic silver for better reflection coefficient. For the in situ experiment, the water dripping on the soda straw caused proper laser reflection and a good signal-to-noise ratio. The bottom of the speleothem was embedded into a concrete block with a chemical seal. The speleothem was excited either with compressed air or light hammer impacts on the concrete block. The velocity of the speleothem extremity was recorded during 10 s at 5000 Hz. For the in situ instrumentation, the laser was setup on a pole to reach the height of the soda straw free extremity. The pole position was strongly constrained by the location of the pedestrian pathway in the cave, the lake preventing the instrument setup in the area (Fig. 1c). In addition, it had to be installed without causing any harm to the soda straws, which restricted again the eligible areas. The pole natural frequencies were measured by pointing the laser beam against a wall and were removed from the records. For each shot, the signal was sampled at 2000 Hz.

Table 2 Geometrical characterization and dynamic parameters of Choranche soda straw speleothems
Fig. 2
figure 2

a Laboratory experimental setup for soda straw dynamic characterization. The laser vibrometer is visible at the forefront, pointing at the speleothem’s top (see red laser dot in the zoom inset, delineated with the red circle).Soda straws CHO_D1.01, CHO_D1.02, CHO_D1.06, and CHO_D1.07 are sealed vertically with chemical grout into a concrete block (background). b Soda straw resonant frequencies (f0–f13) as a function of their length L. Soda straws from datasets CHO_D1 (D1.01, D1.02, D1.06, and D1.07) and CHO_D3 are shown with black triangles and black squares, respectively (this study). Measurements from Lacave et al. (2000) are shown with white stars (CHO_LA). Theoretical resonant frequencies (f0–f13) from Eq. 1 are shown as dashed black lines and labeled for f0–f6. The background red color scale shows the spectral content of shot A (see Sect. 4), with a clear spectral peak close to 175 Hz.

For both laboratory and field experiment, the natural frequencies (fn) were pointed on the laser vibrometer record spectra (Table 2). Previous works conducted by Lacave et al. (2000) in the same cave have also been included (dataset CHO_LA, Table 2).

The theoretical natural frequencies of a Euler-Bernoulli beam are given by the analytical formula in Eq. 1 (from Chopra 2001) and drawn with black dashed lines in Fig. 2b.

$$ {f}_{\mathrm{n}}=\frac{X_{\mathrm{n}}^2}{2\pi}\sqrt{\frac{EI}{\ \mu {L}^4}} $$
(1)

With:

fn:

natural frequency at nth mode (Hz);

\( {X}_{\mathrm{n}}^2 \):

mode coefficient (without unit, see Table 3);

E:

Young’s modulus (N m−2);

I:

moment of inertia of the beam (m4), equal to π (Dext4-Dint4)/64 for a hollow cylinder;

μ:

mass per unit length (kg m−1), equal to ρ (π/4) (Dext2-Dint2) for a hollow cylinder;

L:

beam length (m)

Table 3 Mode coefficients Xn2 for Euler-Bernoulli beam natural frequencies

Figure 2b presents the natural frequencies fres as a function of speleothem length L. Experimental tests for both laboratory experiment (CHO_D1, black triangles) and in situ measurements (CHO_LA, white stars) show consistent results for soda straw ranging from a few centimeters up to about 30 cm in length (Fig. 2b). The best fit between theoretical natural frequencies (dashed black lines) and experimental values was obtained for a 60 GPa Young’s modulus E. This value agrees very well with previous studies conducted on soda straws (Lacave et al. 2000). It is about 2 to 3 times greater than for other kind of speleothems (see Table 1), probably because of complex shape, porosity effects, and/or crystallization defaults for the latter.

Natural frequencies of the CHO_D3.01 soda straw have been measured up to the 13th mode (black squares, Fig. 2b). Experimental resonance frequencies and the ratios between successive modes are both perfectly described by Euler-Bernouilli beam theoretical formulas (Eq. 1 and Table 3). We assumed that embedded-free conditions with chemical seal in the laboratory experiment are comparable with calcite bond in caves. We hence estimated the CHO_D3.01 soda straw length at 1.78 m, which is in good agreement with visual estimate in the field. Our results suggest that Choranche soda straws can be accurately modeled as regular, tubular beams for dynamic characterization.

The damping of soda straw was measured via logarithmic decrement technique (Inman 2008, Table 2). The damping ξ is expressed in percent relative to critical damping. The damping ξ0 measured at the fundamental frequency is generally very low (< 2%), with a mean of about 0.5%. This value is compatible with measurements from Lacave et al. (2000) carried out on the same Choranche soda straws (same table). Higher modes show damping ranging between about 0.1 and 1%. In the rest of this study, we used a ξd = 0.1% design damping value for all modes our computations.

3.2 Rupture parameters

Soda straw failure parameter (i.e. ultimate tensile strength σu,t) was measured on another set of samples (dataset D2, Table 4).

Table 4 Results of failure tests conducted on dataset CHO_D2. B.C.: boundary conditions

We used both single-end loading on embedded-free speleothems (EF) and three-point flexural test (3PFT, Table 1). The load P was progressively increased by pouring sand in a recipient attached to the speleothem extremity (EF) or mid-span (3PFT). The tensile stress at failure was computed using the formula given in Eq. 2 with K = 32 (EF) and K = 8 (3PFT), respectively.

$$ {\sigma}_{\mathrm{u},\mathrm{t}}=\frac{K\ P\ L\ {D}_{\mathrm{ext}}}{\pi\ \left({D}_{\mathrm{ext}}^4-{D}_{\mathrm{int}}^4\right)} $$
(2)

The speleothem own-weight was no measured during computation of σu,t as it is negligible compared with the applied load P.

Our results (Table 4) show that σu,t wanders greatly between 2.5 and 16.2 MPa, most of values concentrating between 2.5 and 7.5 MPa. The lowest values might be caused by preexisting microfractures because of the collected samples have experienced failure: they were found lying on the ground or on the bottom of the lake. Tensile strength σu,t measured agree well with previous studies conducted on soda straws (Table 1). Other types of speleothems generally show lower tensile strength because of porosity effects, crystallization defaults, and/or microdamage.

In the rest of this article, we selected the characteristic value of failure strength as the weakest value of our tests, i.e. about 2.5 MPa. We then defined a conservative design value σu,t,d = 1.25 MPa so that soda straws would not endure more than ~ 50% of their ultimate failure. This value corresponds to a ~ 5% rupture probability based on the tests conducted in Milandre cave (Lacave et al. 2003).

4 Seismic measurements

We measured the ground motion velocity in three locations (GPa, GPb, GPc, Fig. 1c) spread out along the cave wall. Three orthogonal 1C, 4.5 Hz geophones were set at each location. Vertical (Z) and transverse (T) sensors were mounted on metallic L-brackets. Geophone coupling with the rock was achieved by pouring plaster into the drill holes. Radial component (R) was facing the outward normal to the cave wall. We used a Geometrics® Geode® acquisition system operating at 4000 Hz during 4 s. The records started as the same time as source triggering.

Blast tests were conducted using Titanobel® Emulstar 8000 UG and Daveydet® electrical detonators. The maximal instantaneous charge (MIC) ranged from 0.013 to 0.2 kg. The shots were triggered individually at one go. Explosive was fired in holes drilled into the rock at three different spots (A, B, C, Fig. 1c). Each firing hole was used only one time.

The ground motion recorded at sensor GPa during blast test A (MIC = 0.2 kg, d ≈ 55 m) is shown in Fig. 3a. Due to instrument malfunction, only the radial component was usable. The signal amplitude was multiplied by \( \sqrt{3} \) as a conservative proxy of dominant motion (Chapot 1981a). This hypothesis was later confirmed during excavation blasts. The signal shown in Fig. 3a is referred as GP3C in the rest of this study. The velocity signal shows a sudden burst in ground velocity, reaching a peak particle velocity (PPV) of about 1 mm/s. The vibration then decays rapidly after about 0.1 s in length. GP3C spectral content (Fig. 3c, black curve) shows a wide, regular spectrum, spreading between 10 and about 500 Hz. The spectrum culminates between 150 and 200 Hz and shows negligible amplitude below 10 Hz and over 500 Hz. The signal shape and spectrum are both consistent with an explosive source fired at short-medium distance (tens of meters) and traveling through a rock medium.

Fig. 3
figure 3

a, b Ground particle velocity for signals GP3C and GP3C10, respectively. Signal GP3C derives from test blast A (MIC = 0.2 kg) recorded at sensor GPa inside the cave (D ≈ 55 m). Synthetic signal GP3C10 corresponds to 10 micro delays randomly distributed between 20 to 30 ms. c Fourier spectra of GP3C (black) and GP3C10 (gray)

5 Speleothem vibration modeling

The Choranche soda straws were modeled as bending beams, without considering the shear deformation and rotational inertia. A homemade elastodynamic 2D code with Euler-Bernoulli finite elements is used to simulate either a modal analysis (computation of mode shapes and natural frequencies) or a transient response analysis (temporal evolution of displacements, velocities, and stresses at each node) of a soda straw (Roux et al. 2014). Linear dynamic differential equations are solved using Newmark’s numerical integration method (Newmark 1959). The geometrical and mechanical parameters of the beam are given in Sect. 3.1. We neglected water drops hanging at the soda straw extremity—when present—as they weigh generally less than 1–2% of the speleothem own weight. During the transient response analysis, the solicitation consists of the ground velocity signal GP3C derived from Shot A recorded at sensor GPa (see Sect. 4 and Fig. 3a). This ground motion is applied transversally at the speleothem embedded extremity. Preliminary tests were performed to ensure the numerical convergence of the solver so that a 10−5 s time step and 20 finite elements per soda straw were selected. At each time step, the motion and stress state at each node of the soda straw was recorded.

We first simulated the CHO_D3.01 (L = 1.78 m) soda straw response during GP3C solicitation. The soda straw free extremity motion is shown in Fig. 4 for both experimental laser measurements (black) and FEM transient analysis (gray). Time signals are shown on the top and corresponding spectra at the bottom. Both time signals show a sharp increase with peak amplitude followed by a smooth decay. The signal duration is consistent between both signals and lasts a few seconds. The amplitude of the numerical tip motion (~ 2 mm/s) yet appears about 4 times greater than the experimental motion (about 0.5 mm/s). This overestimation could originate from multiple factors. Firstly, the soda straw tip motion may not be aligned along the laser line of sight as the laser position was constrained by the lake and the speleothem tip motion direction was not known a priori. We estimate that the angular bias could cause up to 25% motion underestimation in the experimental laser records. Secondly, the solicitation GP3C applied to the numerical soda straw was estimated from records made along the cave wall (sensor GPa, Fig. 1), which is located about 25 m closer to the shot than the soda straw CHO_D3.01. Based on the propagation law later shown in Sect. 6, this distance effect may represent a factor of ~ 2 in amplitude. In the rest of this study, we hence retained this factor 4 as a safety factor.

Fig. 4
figure 4

Motion of CHO_D3.01 soda straw tip (length L = 1.78 m). a In situ laser vibrometer measurements during explosive shot A. b Numerical FEM simulation results with signal GP3C as input. c Spectrogram signal shown in a. d Fourier spectra (black: laser measurements; gray: FEM modeling)

The spectral content for both experimental (black) and numerical (gray) tip motions are shown in Fig. 4b. We observe a very good agreement between experimental and numerical natural frequencies (f0f13, see also Table 2), especially in the range 10–300 Hz, which is most excited by the blast (Fig. 3). These resonance frequencies values also fit almost exactly the Euler-Bernoulli beam equations. The narrowness of the measured spectral peaks reflects a low damping ratio at the natural frequencies, which is consistent with the laboratory characterization (Sect. 3.1) and satisfactorily simulated by FEM. These results suggest the appropriateness of the FEM dynamic model to simulate the response of CHO_D3.01 soda straw, as its spectral content, damping, and transient response are correctly reproduced.

We then modeled the Choranche soda straw population in the lake room with speleothem length L varying between 0.025 and 3.0 m with a step of 0.025 m. For each simulation, each computation time step and each soda straw length, the maximal tensile stress σu,t,max in the speleothem was saved. σu,t,max value takes into account both the static stress due to speleothem own weight and the dynamic stress due to shaking. This study is the first to simulate resonance amplification when applicable, unlike previous work (see Table 1). Figure 5 shows the maximal tensile stress σu,t,max as a function of speleothem length L. The soda straw own weight is shown as dashed black line. We observe that the stress level remains limited during test shot GP3C (black continuous line), peaking at about 0.2 MPa for a 0.15 m long speleothem. For soda straws in the range 0.025–1 m, we observe successive stress peaks at specific lengths (black arrows and gray vertical lines). These spectral peaks correspond to soda straws natural frequencies (f0–f4), which are very close to the maximal spectral content of signal GP3C (~ 150–200 Hz, with a peak at about 175 Hz, see Fig. 3c). For longer soda straws, the maximal stress appears flatter, showing only little increases with length due to the speleothem own weight. In other words, most of short soda straws (L < 1 m) experience greater stress than the longer ones (L > 1 m) under this test blast A because of dynamic amplification.

Fig. 5
figure 5

Maximal tensile stress σu,t,max as a function of speleothem length computed in FEM numerical simulations. Black: signal GP3C as input. Colored areas: signal GP3C10 as input, for several multiplication factors α. Mean value: continuous line, min/max values: filled area. Vertical gray lines and black arrows: natural frequencies associated with stress peaks. Tensile resistance design value σu,t,d = 1.25 MPa derived from laboratory measurements is shown as horizontal dashed red line. Tensile stress due to speleothem own weight is shown as black dashed line (bottom right corner)

As the peak in stress may correspond to a higher mode (i.e. not the fundamental mode f0), the maximal stress is not necessarily located at the speleothem embedded extremity.

This observation is in contradiction with the equivalent static approach ESA (e.g. Cadorin et al. 2001), which states that the maximal stress is due to flexural bending and is located at the embedded extremity of the speleothem. In the ESA, the stress increases with the squared speleothem length (L2) for a given ground acceleration (the speleothem mass that controls the inertia force amplitude scales with L and the momentum this inertia force squares also with L).

We then carried out simulations for delayed blast shots, which are almost always used for blasting excavation works. It consists in firing several explosive charges separated by a few tens of milliseconds, the delayed charges lying inside the same borehole or in separated boreholes. The assumption is that delaying the successive wave will prevent constructive interaction between the wavefields (i.e. equal chances of constructive and destructive interactions). For large delays between explosive charges and stiff media (i.e. high frequency content, short signal duration), a complete decoupling between charges may even be achieved (Chapot 1981a). The quantification of the optimum delay yet remains delicate (see Hoshino et al. 2000) and standard 25 ms delays are commonly used for civil engineering excavation works.

We created a synthetic signal GP3C10, which simulates 10 explosive charges of 0.2 kg MIC (maximal instantaneous charge of explosive). Each charge was separated by a delay randomly ranging between 20 and 30 ms, to take into account variations in detonator accuracy and propagation times. Ten simulations were carried out to estimate the effect of random delay wandering. An example of GP3C10 is shown in Fig. 3b (time signal) and in Fig. 3c (spectrum, gray curve). We observe that GP3C10 signal culminates at about 1.3 mm/s, which is comparable with GP3C peak ground velocity. GP3C10 yet shows a spectrum richer in high frequencies compared with single blast GP3C.

The maximal tensile stress σu,t,max during GP3C10 (Fig. 5b, α = 1) wanders across the green area during the 10 random simulations. Stress average value is shown by the continuous green line. We observe that the peak stress increases of a factor about 3.5 for GP3C10 (green line) compared with GP3C (black line, Fig. 5b), although the signal peak velocity increased only slightly. This effect is probably due to the successive seismic wavefronts, which tend to enrich the signal content in high frequencies (Fig. 3c), hence increasing the ground acceleration level. Such negative impact of delayed charges on speleothem integrity was also pointed out by Lacave et al. (2003), who advised delays as large as 500 ms for underground blasts to avoid superposition effects. This delay is probably greater for soda straws, which have very low internal damping and slender shape. Unfortunately, electrical detonators used for this study were available only with 25-ms time steps.

Varying the input signal amplitude by a coefficient α, we observed that the stress level appeared linearly related to the input signal amplitude, for a given signal (i.e. for a specific frequency content). Figure 5b shows the results for GP3C10 input signal and amplitude factors α = 1, α = 2, and α = 5. The tensile failure design value σu,t,d = 1.25 MPa (see Sect. 3.2) is reached for α = 2, which means a GP3C10 PPVmax of about 2.4 mm/s. The peak in tensile stress occurs for ~ 0.15 m long soda straws, which experience dynamic amplification at f0.

In order to test the influence of different input signals, we conducted these simulations for other test shots B and C, which have different MIC and location (see Fig. 1c). The same processing steps as shot A were applied to derive corresponding GP3C and GP3C10 signals. Maximal tensile stress σu,t,max as a function of ground peak particle velocity (PPV) is shown in Fig. 6, for both single charge case (black dots input signals GP3C) and delayed charge case (blue dots, input signals GP3C10). Results for blast test A (see Fig. 5b) are shown with squares in Fig. 6. A least square linear fit for all blast tests is drawn as a dashed line. Whatever the blast shot and its frequency content, the maximal tensile stress σu,t,max increases with ground velocity for both single shots (Fig. 6a) and delayed charges blasts (Fig. 6c). We yet note that our study lacks from too few test shots and especially from strong motion blast tests, for cost and safety reasons. We estimated approximately the maximal admissible peak ground velocity PPVmax as the abscissa where the regression (dashed lines) intercepts the design resistance σu,t,d = 1.25 MPa (dashed red line). PPVmax is hence estimated at about 8.6 mm/s for single charge blasts and about 2.4 mm/s for delayed shots. These values are considered along the cave wall, i.e. at the soda straw embedded extremity.

Fig. 6
figure 6

a Maximal tensile stress σu,t,max computed by FEM as function of ground peak particle velocity (PPV). Single blast shots are shown in black and delayed blasts are shown in blue. GP3C and GP3C10 for shot A presented in Fig. 3 are shown with squares. Least-square regressions are shown with dashed lines. b Ground peak particle velocity (PPV) as a function of scaled distance (i.e. distance D to the blast divided by square root of the maximal instantaneous explosive charge MIC). c Maximal tensile stress σu,t,max computed by FEM as function of scaled distance. Single blast shots are shown in black and delayed blasts are shown in blue. GP3C and GP3C10 for shot A presented in Fig. 3 are shown with squares. Least-square regressions are shown with dashed lines. For a and c, the maximal tensile stress admissible in the soda straws σu,t,d is shown as dashed red line

6 Estimation of maximal charge

In order to estimate the maximal explosive charge that can be fired for the excavation, we fitted a propagation law on the vibration peak particle velocities (PPV). Measurements from three component geophones, GPa, GPb, and GPc, for all test shots were used in combination with an additional seismic profile (PS) located outside the cave (Fig. 1c). PS profile was composed of ten 4.5 Hz vertical geophones cemented on the rock with plaster.

The measured PPV as a function of scaled distance are shown in Fig. 6b. GP measurements (inside the cave) are shown with white squares, PS profile with black crosses. Peak particle velocity shows a steady decrease with the scaled distance (in log-log axes). A greater spreading is observed for GP measurements, which recorded test shots fired at different locations (Fig. 1c) than for profile PS.

$$ \mathrm{PPV}=K\ {\left(\frac{D}{\sqrt{\mathrm{MIC}}}\right)}^{-b} $$
(3)

with:

PPV:

peak particle velocity (in mm s−1)

MIC:

maximal instantaneous charge (in kg)

D:

distance (in m)

K:

site coefficient (without unit)

b:

site attenuation coefficient (without unit)

The power-law shown in Eq. 3 was then fitted on the measurements (Chapot 1981b; Duvall and Petkof 1958 in Kumar et al. 2016). This formula scales the distance by the square root of explosive maximal instantaneous charge (i.e. \( D/\sqrt{\mathrm{MIC}} \)) to account for the mix of body and surface waves emitted. The exponent b controls the attenuation and ranges generally between 1 and 2 (Kumar et al. 2016). We chose b = 1.8 as best fit on the experimental data. The coefficient K depends both on the site and shot parameters; it tunes the PPV level for unit scaled distance. It ranges commonly between 300 and 6000, with a median value of 2500 (Chapot 1981b). In our case, K ranges between 800 and about 4000 (respectively dotted and continuous lines in Fig. 6b). The PS profile located on the external, fractured rock skin shows the lower K, in agreement with the results from the GPR campaign (Sect. 2). Higher K is observed for most GP measurements located inside the cave, reflecting the greater rock quality at depth. Therefore, we selected K = 4000 as a reasonably safe design value for the estimation of the maximal MIC admissible during the works.

Maximal tensile stress σu,t,max as a function of scaled distance is shown in Fig. 6c for both single shots (black dots) and delayed charges blasts (blue dots). Both series show a clear and consistent decrease in σu,t,max experienced by the speleothems when increasing D and/or reducing MIC. We modified Eq. 3 into Eq. 4 to express the maximal tensile stress σu,t,max as a function of scaled distance (\( D/\sqrt{\mathrm{MIC}} \)). The site coefficient K1 in Eq. 4 is equal to regression slope coefficients from Fig. 6a multiplied by K = 4000. K1 is equal to 580 for single shots and about 2096 for delayed shots.

$$ {\upsigma}_{\mathrm{u},\mathrm{t},\max }={K}_1\ {\left(\frac{D}{\sqrt{\mathrm{MIC}}}\right)}^{-b} $$
(4)

with:

σu,t,max:

maximal tensile stress (in MPa)

MIC:

maximal instantaneous charge (in kg)

D:

distance (in m)

K1:

site coefficient (without unit)

b:

site attenuation coefficient (without unit)

The critical scaled distance (\( D/\sqrt{\mathrm{MIC}} \)) is reached when the power law equals the design resistance of the speleothem σu,t,d = 1.25 MPa (dashed red line). Critical scaled distances range between 23.7 (single shot) and 48.5 m/√kg (blast with delayed charges). Lacave et al. (2003) advised similar minimal scaled distances (i.e. between 50 and 70 m/√kg) when using delayed explosive charges close to speleothems.

7 Discussion

Our study dealt only with perfectly tubular soda straws, which is the case for most of the speleothems found in the lake room. However irregular sections or local wall thickening can be found on some samples in the Choranche cave. In addition, other speleothems are present in the cave such as massive stalactites, stalagmites, columns, flowstone, and drapery. For such speleothems, the dynamic parameters may differ significantly from our study. In addition, the FEM should be performed individually and consider the actual speleothem geometry. This represents a time-consuming task, applicable only for a limited number of speleothems.

During our vibration study, we fired only underloaded shots for safety reasons so that no rock was fractured. Such blasts tend to generate higher vibration amplitudes than excavation blasts for the same explosive charge. This assumption provides an additional safety margin. In addition, test blasts generally consist in isolated shots, while delayed charges are fired for the production phase. Synthetic signals created by repeating pulses e.g. every 25 ms tend to underestimate the acceleration level applied at the speleothem embedded extremity (i.e. soda straw top in our case). We hence applied a random delay in the range 20–30 ms to consider the accuracy of electric detonators (about a few milliseconds, Silva et al. 2018) and variations in source-speleothem distance due to blast hole configuration. The nonperiodicity of the input signal created stress peaks in the soda straws. This effect appears clearly in low-damped structures such as speleothems (and more especially soda straws) for which only little energy decay occurs between successive wave train arrivals (Lacave et al. 2003).

In our study, the finite element model allowed us to simulate the dynamic amplification experienced by soda straws, whatever the excited mode. Signal GP3C and corresponding spectrogram recorded during single blast test A are recalled in Fig. 7a. The maximal tensile stress in the speleothems is shown in Fig. 7c (FEM, black curve; same curve as in Fig. 5). Our study showed that relatively short soda straws (0.1–1 m in length) were the most exposed to nearby rock blasting because they experience higher dynamic amplification. Indeed, the conjunction between peaks in ground motion frequency content (170–180 Hz for blast A, Fig. 3b) and soda straws resonance frequencies gives rise to significant amplification. Higher modes (up to f4 for soda straw length L = 0.9 m) experienced significant amplification under nearby blasting.

Fig. 7
figure 7

a Ground acceleration (top) and spectrogram (bottom) of signal GP3C applied at the speleothem embedded extremity. b Ground acceleration (top) and spectrogram (bottom) of earthquake signal applied at the same location. Spectrogram are normalized between 0 (white color) and 1 (red). Subplots c and d show the maximal tensile stress σu,t,max as a function of speleothem length, for blast test GP3C a and earthquake b as input signals, respectively. Results are drawn for the dynamic FEM (continuous black line), one degree of freedom 1DOF (continuous gray line), and equivalent static approach (ESA) (thick gray dotted line). Tensile resistance design value σu,t,d = 1.25 MPa derived from laboratory measurements is shown as dashed red line. Tensile stress due to speleothem own weight is shown as dashed black line. The top x-axis shows the theoretical fundamental frequency f0 for a speleothem of length L (see Eq. 1 and Fig. 2b)

For comparison purposes, we computed the maximal tensile stress σu,t,max yielded by two other approaches during the same GP3C solicitation. In the first approach, we used a simple 1 degree of freedom (1DOF) spring-mass system to compute the maximal relative mass displacement. The maximal tensile stress at the embedded extremity was estimated using static Euler-Bernoulli beam (EBB) equations for embedded-free flexural test. For speleothems shorter than ~ 0.3 m in length, the maximal tensile stress predicted by 1DOF (gray curve, Fig. 7c) is slightly underestimated but shows consistent shape compared to FEM (black curve). Indeed, the fundamental natural frequency f0 coincides with the predominant ground motion excitation (~ 175 Hz for shot A) for short speleothems (< 0.3 m). In this case, the maximal stress can be approximately estimated using 1DOF dynamic model and EBB flexural equations. In contrast, the stress level yielded by 1DOF is significantly underestimated (up to a factor 5) for longer speleothems. They indeed exhibit significant dynamic motion at higher modes (f1f4, see FEM results in Fig. 7c). Those modes are not simulated by the 1DOF spring-mass model, which yields negligible mass motion (“inertial” behavior) for speleothems longer than ~ 0.6 m (i.e. for f0 < 10 Hz, Fig. 7c).

In the second approach, we used the equivalent static approach (ESA, gray dotted line in Fig. 7c). The maximal tensile stress σu,t,max is correctly estimated for very short soda straws (L < 0.08 m). The rigid-body assumption is fulfilled for such speleothems as their resonance frequencies are considerably higher than the blasting spectral content (see f0 in Fig. 7c, top axis). The tensile stress is yet considerably underestimated for 0.08–0.2 m long soda straws as ESA does not take into account the dynamic amplification that occurs mainly at the fundamental mode for those speleothems (f0 = 175 Hz for L ~ 0.15 m). In contrast, σu,t,max is considerably overestimated by the ESA when L > 0.2 m. Indeed, σu,t,max scales with L2 in ESA as both the inertia force and the bending moment scale individually with L. In contrast, FEM shows lower σu,t,max because the rigid body assumption is not valid anymore and tends to overestimate the stress level.

We also compared the FEM, 1DOF, and ESA methods for a regional earthquake solicitation (Fig. 7 b and d). The historical reference earthquake in the Vercors massif is the 25/04/1962 Corrençon ML = 5.3–5.5, located at about 11 km in distance from the Choranche cave. Unfortunately, its seismic signal was not recorded at that time. We hence selected the 07/04/2014 Barcelonette earthquake as a proxy (ML = 5.2, MW = 4.8, depth about 13 km/surface). The ground motion was recorded at station SURF located at ~ 10.0 km in distance from the epicenter (RAP-RESIF network, Sira et al. 2014). The acceleration and velocity amplitudes reached a peak values PGA = 0.5 m/s2 and PPV = 0.044 m/s. Such event is quite rare for mainland France (once or twice per decade) but remains weak compared with possible earthquakes in the Choranche cave area (e.g. the construction design PGA for new buildings equals 1.6 m/s2 in the Vercors, for comparison).

As earthquake-induced motion is strongly attenuated in underground cavities (Becker et al. 2006; Gribovszki et al. 2017), the earthquake signal amplitude was corrected following Hu and Xie (2005) formula. Site conditions for seismic station SURF belong to the soil over rock class (Hollender et al. 2015), with a mean shear wave velocity over the first 30 m of subsurface Vs30 = 420 m/s. In this case, the expected PGA in a 350-m-deep cave scales at about 1/10th the surface PGA, for both vertical and horizontal components (Hu and Xie 2005). The modified acceleration signal (top) and spectrogram (bottom) are shown in Fig. 7b. The PGA culminates at about 0.05 m/s2, with a signal frequency content ranging between 0.1 and ~ 10 Hz. Results are shown in Fig. 7d for FEM (continuous black line), 1DOF (continuous gray line), and ESA (dotted gray line). Maximal tensile stress curves yielded by FEM and 1DOF show a similar shape for most of the speleothems. 1DOF yet tends to underestimate the stress level as it takes only into account the contribution of the fundamental mode f0 for dynamic amplification. The error is greater for longer speleothems because their upper modes f1–f2 may also lie in the seismic frequency range (i.e. below 10 Hz, see Fig. 2b) contrary to Lacave et al. (2004) assumption. The ESA performs approximately well for soda straws shorter than 0.6 m, which resonance frequencies are greater than the seismic bandwidth (f0 > 10 Hz, Fig. 2b). For such soda straws, the rigid body assumption is valid. The ESA yet underestimates σu,t,max when the dynamic resonance becomes non negligible (0.6 < L < 2.5 m) and overestimates the tensile stress for speleothems longer than 2.5 m.

For a 0.05 m/s2 PGA at soda straw embedded extremity, the maximal tensile stress σu,t,max peaks at about 1.0 MPa (Fig. 7d). Highest tensile stresses (0.5 < σu,t,max < 1 MPa) are observed in speleothems measuring between 1 and 2 m in length, which are poorly simulated by 1DOF or ESA. The maximal induced stress σu,t,max during the earthquake remains smaller than the design value σu,t,d = 1.25 MPa in the FEM simulations so that soda straw rupture is unlikely for the selected earthquake event. Our result suggest that calcite strength design limit σu,t,d is reached when PGA at soda straw embedded extremity exceeds about 0.06 m/s2. Such PGA value may be reached for closer earthquakes, higher source magnitude, and/or lower PGA attenuation with depth, depending on the site conditions.

Our study underlines the need for a proper dynamic amplification modeling when speleothem natural frequencies and solicitation spectrum coincide, even partially. The PGA value itself does not contain enough information as the solicitation frequency content plays a major role. Our results yielded admissible PGA of about 0.06 m/s2 under earthquake and about 2.4 m/s2 for a blast with ten delayed explosive charges (see signal GP3C10; α = 1 in Fig. 5). Such relatively low values result from our adequate dynamic amplification simulation, from the intrinsic vulnerability of the soda straws, and also from the safe assumptions made in terms of calcite tensile strength (σu,t,d = 1.25 MPa) and admissible rupture probability. Under seismic solicitation, Lacave et al. (2004) found slightly higher admissible PGA (0.1–0.2 m/s2) for a σu,t,d distribution peaking at about 3 MPa. The admissible PGA was computed as 5% rupture probability during earthquake for very vulnerable, slender, speleothems. For a slender (3 m long, 0.03 m in diameter)-filled stalagmite, an admissible PGA of 0.5 m/s2 was found by Gribovszki et al. (2008) assuming σu,t,d = 2.61 MPa. For the unique previous case study of speleothems exposed to blasting works, Lacave et al. (2003) defined an admissible PGA of about 1.8 m/s2, at 5% rupture probability for soda straws.

In speleoseismology, the ESA should be restricted to low frequency solicitation (e.g. earthquakes) applied to rigid, massive speleothems with high resonance frequencies. ESA application to soda straws—whose natural frequencies may be as low as a few tenth of Hz or a few Hz—may significantly underestimate or overestimate the induced dynamic stress. The 1DOF approach performs slightly better as it simulates the dynamic amplification occurring at the fundamental mode f0. It yet tends to underestimate the induced tensile stress, especially when higher modes are excited. Application of 1DOF or ESA without caution may induce severe biases in speleoseismological studies such as for earthquake hazard studies.

8 Conclusions

This research work was conducted to preserve remarkable soda straws exposed to nearby rock blasting in Choranche cave. We first characterized the geometric, dynamic, and ultimate properties of the speleothems through laboratory and field tests. Blast tests in situ were used as input for a dynamic 2D finite element code. We showed that the soda straw dynamic response varies considerably depending on the match between ground motion frequency content and speleothem natural frequencies. Short soda straws (L < 1 m) experienced higher stress under blast tests than the longer ones due to dynamic resonance. Opposite conclusions were made for speleothems exposed to a regional earthquake because of lower ground motion frequency content. This underlines the good performance of our dynamic finite element model, whatever the solicitation and excited modes. We show that the omission of dynamic resonance or its simplification as made in previous studies may significantly underestimate the induced load in speleothems.

We defined an admissible vibration threshold along the cave wall and we monitored its observance during the excavation blasts. The works were carried out without causing any harm to the soda straws. Our research work may constitute a milestone for future speleoseismological studies. As urban growth, renewal, and densification processing are increasing in many countries, we expect more and more case studies of concretionary caves exposed to anthropogenic vibrations.