Brief communication: Stalagmite damage by cave ice flow quantitatively assessed by fluid–structure interaction simulations

Abstract

Mechanical damage to stalagmites is commonly observed in mid-latitude caves. Former studies have identified thermoelastic ice expansion as a plausible mechanism for such damage. This study builds on these findings and investigates the role of ice flow along the cave bed as a possible second mechanism in stalagmite damage. Utilizing fluid–structure interaction models based on the finite-element method, forces created by ice flow are simulated for different stalagmite geometries. The resulting effects of such forces on the structural integrity of stalagmites are analyzed and presented. Our results suggest that structural failure of stalagmites caused by ice flow is possible, albeit unlikely.

1 Introduction

Secondary carbonate deposits such as stalactites and stalagmites, collectively referred to as speleothems, are abundant in caves worldwide and grow over extended periods of time, typically many thousands of years, thereby recording valuable proxy information about climate and vegetation outside the cave. Speleothems are therefore highly sought archives of paleoenvironmental change e.g.,. Despite the fact that caves are typically well protected from processes such as weathering or erosion, some caves, particularly in the mid-latitudes, show naturally occurring mechanical damage to speleothems. Several processes have been discussed that cause damage to stalagmites and stalactites, including movement of stalagmites growing on a soft-sediment substrate, gravity-driven breakdown of cave ceilings, strong earthquakes, and the former presence of ice in these subterranean cavities. The latter process has been identified as a main cause for stalagmites being sheared off their bases e.g.,, as well as for sheet-like deposits known as flowstones showing a fractured appearance (e.g.,; see Fig. ). Key to this interpretation is the coexistence of damaged speleothems with so-called cryogenic cave carbonates . The former are indicators of the past presence of perennial ice in these caves and can be reliably dated using radiometric techniques e.g.,. In addition, damaged speleothems are often overgrown by younger calcite of Holocene age, which proves the glacial age of the damage e.g.,. In a recent study, used numerical modeling combined with laboratory measurements to quantitatively assess which processes in former ice-filled caves led to damage of stalagmites. The results show that internal deformation within a cave ice body due to ice flow under gravity cannot fracture stalagmites, even on steep slopes. In contrast, thermoelastic stresses within an ice body due to temperature changes reach values partly exceeding the tensile strength of stalagmites. In this contribution we expand on the analysis of by (a) examining not only the external forces exerted by ice flow but also their internal effect within stalagmites and (b) also including smaller stalagmite dimensions. Structural failure of stalagmites caused by ice flow was investigated using a fluid–structure interaction (FSI) model e.g.,.

Figure 1

(a, b) Examples of stalagmites damaged by ice during the last glacial period as shown by the presence of cryogenically formed calcite deposits in these caves. (a) Two parallel stalagmites detached from their base along subhorizontal cracks (arrows), Obir Caves, Austria. Width of image 0.4 m. (b) Light-brown stalagmite broken near its base and tilted. Younger white calcite (of likely Holocene age) is growing over the lower part of this damaged speleothem and forms a smaller stalagmite to the left. Also note dark-brown flowstone in the foreground, likely fractured due to the former presence of ice. Shatter Cave, Mendip Hills, UK. Width of image 1.1 m. (c) Example of a modern ice cave with floor ice and thick ice deposits on the slope in the background. Feuertal Eishöhle, Austria. Images: Christoph Spötl (a, b) and Harald Zeitlhofer (d).

[Figure omitted. See PDF]

2 Methods

In order to be consistent with the study of , we chose to model the same ice domain with $-\mathrm{7}\le x\le \mathrm{14}$ $\mathrm{m}$, $-\mathrm{3.5}\le y\le \mathrm{3.5}$ $\mathrm{m}$ and the stalagmite location at $x=y=\mathrm{0.0}$ $\mathrm{m}$ (cf. Fig. 3 in ). However, we examined two ice thicknesses ($h_{\mathrm{ice}}=\mathrm{1.0}$ and $\mathrm{2.0}$ $\mathrm{m}$), two slope angle values (${\mathit{\beta }}_{\mathrm{ice}}=\mathrm{10}$ and $\mathrm{40}$°), three stalagmite diameters ($d_{\mathrm{s}}=\mathrm{0.15},\mathrm{0.25}$ and $\mathrm{0.4}$ $\mathrm{m}$), two stalagmite heights ($h_{\mathrm{s}}=\mathrm{0.5}$ and $\mathrm{1.0}$ $\mathrm{m}$), and full-slip and no-slip conditions at the cave base. We did not utilize any glaciological sliding laws, as we wanted to investigate the two end-member states (frictionless full sliding vs. no-slip). In total, 14 different configurations, which are summarized in Table , were examined. The cave walls (sidewalls of the ice domain) were assumed to be no-slip boundary conditions, and the ice surface was a stress-free boundary, allowing for free movement of ice. Stalagmites were modeled as cylinders with a fixed diameter $d_{\mathrm{s}}$, except for configurations 13 and 14, which exhibited a conical shape ($d_{\mathrm{s}}^{\mathrm{base}}=\mathrm{0.25}$ $\mathrm{m}$ and $d_{\mathrm{s}}^{\mathrm{top}}=\mathrm{0.1}$ $\mathrm{m}$). Our numerical mesh typically contains 100 000 tetrahedral elements with mesh element sizes varying between 0.3 $\mathrm{m}$ far away from the stalagmite and 0.02 $\mathrm{m}$ within the stalagmite.

2.1 FSI model

We assumed the deformation of ice to be fluid-like, with a non-Newtonian rheology following Glen's flow law e.g.,. Driven by gravity, the cave ice flows downslope and exerts forces on the stalagmite's outer surface. The structural response of the stalagmite to these external forces was modeled as linearly elastic. In , laboratory experiment results were presented that identified the elastic properties of the stalagmite material. As we are interested in the stress distribution just before structural failure, linear elasticity e.g., is a reasonable choice for our model. Using the chosen FSI model approach, we simulated both the forces exerted by ice deformation (fluid) and the response of the stalagmite (solid structure) to these forces with a finite-element model. A summary of the relevant model parameters is given in Table .

Table 1

Model parameters.

Symbol	Value	Units	Description
$A$	$\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{24}}$	${\mathrm{Pa}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$	Glen's flow law parameter for temperate ice
$n$	$\mathrm{3}$	–	Glen's flow law non-linearity number
$E_{\mathrm{i}}$	$\mathrm{9.5}\times {\mathrm{10}}^{\mathrm{9}}$	$\mathrm{Pa}$	Young's modulus for ice
$E_{\mathrm{s}}$	$\mathrm{6.41}\times {\mathrm{10}}^{\mathrm{10}}$	$\mathrm{Pa}$	Young's modulus for stalagmites
$g$	$\mathrm{9.81}$	$\mathrm{m}{\mathrm{s}}^{-\mathrm{2}}$	Gravitational acceleration
${\mathit{\nu }}_{\mathrm{i}}$	$\mathrm{0.33}$	–	Poisson ratio for ice
${\mathit{\nu }}_{\mathrm{s}}$	$\mathrm{0.27}$	–	Poisson ratio for stalagmites
${\mathit{\rho }}_{\mathrm{i}}$	$\mathrm{910}$	$\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$	Density of ice
${\mathit{\rho }}_{\mathrm{s}}$	$\mathrm{2670}$	$\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$	Density of stalagmite material

We developed this simulation strategy using the open-source software Elmer FEM e.g.,, which includes well-established FSI routines e.g.,. To include the non-linear flow behavior of ice in our model, we used Elmer/Ice e.g.,, which builds upon Elmer FEM.

Elmer FEM applies a partitioned approach to solving FSI problems, thus solving the equations governing non-linear fluid flow and elastic structural response with two different numerical solvers. Forces are coupled at the fluid–structure boundary (stalagmite outer surface), and the model performs iteration until a balance between fluid stress and the elastic response of the solid is found.

Figure 2

Maximum principal stress ${\mathit{\sigma }}_{\max }$ for the conical stalagmite of simulation 13 in ${\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}$ (displacements are scaled by a factor of 1000). Red colors on the upstream face indicate tensile stresses (${\mathit{\sigma }}_{\max }>\mathrm{0}$). The $z$ axis of the coordinate system is aligned vertically upward against gravity, while the $x$ axis aligns with the downstream direction of the ice flow.

[Figure omitted. See PDF]

Figure 3

Maximum principal stress ${\mathit{\sigma }}_{\max }$ in the stalagmite pertaining to simulation 09 in ${\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}$ (displacements are scaled by a factor of 100). This simulation indicates stresses significantly exceeding the tensile strength of the stalagmite on its upstream side. The $x$ axis points in the downstream direction of the ice flow.

[Figure omitted. See PDF]

3 Results

Solutions from the FSI simulations were analyzed for maximum (principal) tensile stresses ${\mathit{\sigma }}_{\max }$ in the stalagmite to identify cases where its tensile strength ($\mathrm{4.3}\times {\mathrm{10}}^{\mathrm{6}}\pm \mathrm{1.5}\times {\mathrm{10}}^{\mathrm{6}}$ $\mathrm{Pa}$; ) is exceeded. According to Rankine's theory , this would result in structural failure of the stalagmite. Typically, the maximum tensile stresses were found to be caused by bending of the stalagmite and occurred at the upstream side of its base in the vertical direction (see Fig. ). The results of the simulations are summarized in Table .

Table 2

Maximum principal stress max ${\mathit{\sigma }}_{\max }$ and maximum shear stress max $\left|{\mathit{\tau }}_{zx}\right|$ for all simulations. ${\mathit{\beta }}_{\mathrm{ice}}$ is the cave slope angle, $h_{\mathrm{ice}}$ the ice thickness, $h_{\mathrm{s}}$ the stalagmite height and $d_{\mathrm{s}}$ the stalagmite diameter. Boundary conditions (BCs) at the ice bed are listed as well. Bold simulation numbers indicate configurations for which the stalagmite would likely break due to the arising tensile stresses.

Sim. no.	${\mathit{\beta }}_{\mathrm{ice}}$	$h_{\mathrm{ice}}$	$h_{\mathrm{s}}$	$d_{\mathrm{s}}$	Base BC	max ${\mathit{\sigma }}_{\max }$	max $\left|{\mathit{\tau }}_{zx}\right|$
	(°)	($\mathrm{m}$)	($\mathrm{m}$)	($\mathrm{m}$)		(${\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}$)	(${\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}$)
01	10	2	1	0.4	no-slip	1.22	0.273
02	10	2	1	0.4	slip	3.53	0.789
03	10	2	1	0.25	no-slip	3.55	0.745
04	10	2	1	0.25	slip	11.1	2.36
05	10	2	0.5	0.15	no-slip	2.73	0.588
06	10	2	0.5	0.15	slip	15.1	3.27
07	10	1	0.5	0.15	no-slip	1.51	0.325
08	10	1	0.5	0.15	slip	13.9	3.01
09	40	2	1	0.25	no-slip	26.1	6.58
10a	40	2	1	0.25	slip	86.9	22.2
11	40	2	0.5	0.15	no-slip	22.4	5.75
12a	40	2	0.5	0.15	slip	105	27.6
13b	10	2	1	0.25	no-slip	3.09	0.712
14b	10	2	1	0.25	slip	10.3	2.4

^a Maximum ${\mathit{\sigma }}_{\max }$ on the downslope side of the stalagmite exceeds the tensile strength of ice (cf. Sect. ). ^b The configurations exhibiting a conical stalagmite shape.

For moderately sloping caves (${\mathit{\beta }}_{\mathrm{ice}}=\mathrm{10}$°), only configurations exhibiting full-slip conditions at the ice–cave base interface (no bed friction) indicate failure (simulations 04, 06, 08 and 14). In the case of steeply sloping caves (${\mathit{\beta }}_{\mathrm{ice}}=\mathrm{40}$°), all investigated configurations (simulations 09 to 12) exhibit maximum tensile stresses significantly exceeding the stalagmite's tensile strength.

We also investigated the effect of ice temperature, as all our simulations assume temperate ice (cf. $A$ in Table ). However, reducing the ice temperature to $-\mathrm{5.0}$ °C only results in a reduction in the maximum tensile stress of less than $\mathrm{0.5}$ ‰ relative to the temperate case.

Ice thickness $h_{\mathrm{ice}}$, and thereby how much ice overflows the stalagmite, significantly affects ${\mathit{\sigma }}_{\max }$. Comparing simulations 05 ($h_{\mathrm{ice}}=\mathrm{2}\mathrm{m}$) and 07 ($h_{\mathrm{ice}}=\mathrm{1}\mathrm{m}$) reveals a reduction in ${\mathit{\sigma }}_{\max }$ by about 45 % in the case of no-slip conditions. In the case of full-slip conditions (simulations 06 and 08), however, only a small reduction in ${\mathit{\sigma }}_{\max }$ by about 8 % for the same $\mathrm{50}\mathit{\%}$ reduction in $h_{\mathrm{ice}}$ is observed.

The stalagmite geometry is found to significantly influence the observed maximum tensile stresses. This is especially true for the stalagmite base diameter $d_{\mathrm{s}}$ and by extension its section modulus $S=\frac{\mathit{\pi }{d}_{\mathrm{s}}^{\mathrm{3}}}{\mathrm{32}}$. Even though a reduction in $d_{\mathrm{s}}$ reduces the surface area exposed to ice flow, ${\mathit{\sigma }}_{\max }$ for simulations 03 and 04 ($d_{\mathrm{s}}=\mathrm{0.25}$ m) is nearly 3 times higher than for simulations 01 and 02 ($d_{\mathrm{s}}=\mathrm{0.4}$ m) due to a disproportionate reduction in $S$. Given the same ice thickness of $h_{\mathrm{ice}}=\mathrm{2}\mathrm{m}$, small stalagmites ($d_{\mathrm{s}}=\mathrm{0.15}$ m and $h_{\mathrm{s}}=\mathrm{0.5}$ m; simulations 05 and 06) are more prone to failure than larger stalagmites ($d_{\mathrm{s}}=\mathrm{0.4}$ m and $h_{\mathrm{s}}=\mathrm{1}$ m; simulations 01 and 02). Conically shaped stalagmites (configurations 13 and 14) exhibit comparable, albeit slightly lower, maximum tensile stresses than cylindrical stalagmites (configurations 03 and 04).

For our simulations we assumed ice frozen to the stalagmite and thus a direct force coupling at the ice–stalagmite interface. Therefore, the modeled ice flow can potentially exert tensile forces on the downslope side of the stalagmite. By comparing the maximum principal tensile stress ${\mathit{\sigma }}_{\max }$ on the downslope ice–stalagmite interface with reported values for ice tensile strength ($\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$–$\mathrm{3.1}\times {\mathrm{10}}^{\mathrm{6}}$ Pa; ), we identified two configurations (configurations 10 and 12) for which failure of the ice itself would likely occur. For these configurations the reported maximum ${\mathit{\sigma }}_{\max }$ values are definitely an overestimation as they assume full stress coupling between the ice and stalagmite on the downslope side. In reality, at such high stresses, the ice itself would fail and cracks would form which would decouple the ice body from the stalagmite at the downstream side. Our simulations do not account for such crack formations, and thus the results for configurations 10 and 12 have to be considered with caution.

High shear stresses within the stalagmite could also lead to structural failure; hence we report maximum absolute values of the base shear stress $\left|{\mathit{\tau }}_{zx}\right|$ in Table . However, maximum $\left|{\mathit{\tau }}_{zx}\right|$ amounts to less than 30 % of maximum ${\mathit{\sigma }}_{\max }$ for all investigated configurations. We therefore consider it unlikely that the stalagmite would fail due to the presence of shear stresses.

4 Discussion and conclusion

Ice caves exhibit complex geometries which comprise galleries, chambers and shafts. In addition, the floor of chambers and galleries is typically covered by coarse debris including variably large and angular blocks. These topographic intricate features render free-sliding (i.e., full-slip) interface conditions between the ice and the cave bed highly unlikely in reality. Most importantly, cave ice bodies are typically frozen to their substrate. They owe their existence to the strong cooling of the rock mass surrounding the cave as a consequence of strong (winter) cooling by airflow or the presence of permafrost e.g.,.

Any zones of full-slip would be compensated for by ice being required to flow around and/or over complex bed topographies. Hence we regard our results with basal full-slip conditions (simulation numbers 02, 04, 06, 08, 10, 12 and 14) as highly unlikely upper-end-member cases for ice cave conditions. Conditions close to the no-slip simulations are regarded as most likely.

We found large tensile stresses on the upslope (i.e., upstream) side of stalagmites to be the main cause for structural failure. In contrast, a large influence of shear stresses near the stalagmite base is ruled out due to their comparatively small magnitude of less than 30 % of the observed tensile stresses.

Our results suggest that only in steeply sloping caves with large ice bodies (e.g., 14 $\mathrm{m}$ of cave ice, 2 $\mathrm{m}$ thick, upslope of a stalagmite) can forces exerted by ice flow become large enough to damage moderately large stalagmites. Only in one case with no-slip conditions (simulation number 09; see Fig. ) did we find a configuration in which a moderately large stalagmite would be damaged.

This study confirms that vertical uplift of stalagmites due to thermoelastic stresses induced by expanding ice, as previously reported , is still the most likely physical explanation for widespread damage to speleothems in mid-latitude caves during the cold periods of the Pleistocene. Even though our results rule out lateral ice flow as a main damage mechanism, which was previously often assumed, in extreme conditions (i.e., steeply sloping, large, ice-filled caves), large tensile stresses caused by ice flow can occur in stalagmites, which ultimately can also cause structural damage.

Code availability

All fluid–structure interaction simulations in this contribution were carried out with the Elmer FEM and Elmer/Ice software packages .

Author contributions

CS and AHJ conceived the study. AHJ performed the numerical simulations and, together with PH, analyzed the results. All authors contributed to writing the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We would like to thank the two anonymous reviewers as well as the editor in charge for their work during the publication process. This research was funded in whole or in part by the Austrian Science Fund (FWF) grant DOIs 10.55776/P31874 and 10.55776/P35877 . For open-access purposes, the authors have applied a CC BY public copyright license to any author-accepted manuscript version arising from this submission.

Financial support

This research has been supported by the Austrian Science Fund (grant nos. P318740 and P358770).

Review statement

This paper was edited by Nanna Bjørnholt Karlsson and reviewed by two anonymous referees.