The initial stages of cave formation: Beyond the one-dimensional paradigm

The solutional origin of limestone caves was recognized over a century ago, but the short penetration length of an undersaturated solution made it seem impossible for long conduits to develop. This is contradicted by field observations, where extended conduits, sometimes several kilometers long, are found in karst environments. However, a sharp drop in the dissolution rate of CaCO_3 near saturation provides a mechanism for much deeper penetration of reactant. The notion of a “kinetic trigger” - a sudden change in rate constant over a narrow concentration range - has become a widely accepted paradigm in speleogenesis modeling. However, it is based on one-dimensional models for the fluid and solute transport inside the fracture, assuming that the dissolution front is planar in the direction perpendicular to the flow. Here we show that this assumption is incorrect; a planar dissolution front in an entirely uniform fracture is unstable to infinitesimal perturbations and inevitably breaks up into highly localized regions of dissolution. This provides an alternative mechanism for cave formation, even in the absence of a kinetic trigger. Our results suggest that there is an inherent wavelength to the erosion pattern in dissolving fractures, which depends on the reaction rate and flow rate, but is independent of the initial roughness. In contrast to one-dimensional models, two-dimensional simulations indicate that there is only a weak dependence of the breakthrough time on kinetic order; localization of the flow tends to keep the undersaturation in the dissolution front above the threshold for non-linear kinetics.

💡 Research Summary

The paper revisits the long‑standing problem of how kilometer‑scale solutional caves can develop in limestone despite the short penetration depth of undersaturated water. Traditional speleogenesis models rely on a “kinetic trigger”: a sharp drop in the CaCO₃ dissolution rate as the solution approaches saturation. These models are fundamentally one‑dimensional (1D), assuming that the dissolution front is planar and that flow and solute transport occur only along a single axis. Under this assumption, the front quickly becomes saturated, limiting further penetration and requiring the kinetic trigger to sustain long conduits.

The authors demonstrate that the planar front assumption is intrinsically unstable even in a perfectly uniform fracture. By performing a linear stability analysis, they show that infinitesimal perturbations to the fracture aperture generate a feedback loop: local increases in aperture accelerate the flow, which in turn raises the local undersaturation and dissolution rate. This positive feedback amplifies the perturbation, causing the front to evolve into a sinusoidal wave rather than remaining flat. The most unstable wavelength λ emerges from a balance between advection, diffusion, and reaction and can be expressed as λ ≈ √(D·q/k), where D is the molecular diffusion coefficient, q the mean volumetric flux, and k the kinetic constant governing dissolution. Crucially, λ depends only on reaction and flow parameters and is independent of the initial surface roughness.

To test the theory, the authors conduct two‑dimensional (2D) numerical simulations that couple Navier‑Stokes flow with a convection‑diffusion‑reaction equation for Ca²⁺ and CO₃²⁻. The fracture is initialized as a uniform aperture with only minute random height variations. Simulations are run for both linear (first‑order) and non‑linear (higher‑order) kinetics. The results confirm the analytical predictions: within a few hundred pore‑volume injections, the initially flat dissolution front breaks up into a series of narrow, high‑flux channels. These channels concentrate flow, keep the fluid undersaturated ahead of the front, and allow the dissolution front to advance deep into the fracture.

A key outcome is that the “breakthrough time” (the moment a continuous conduit spans the fracture) shows only a weak dependence on kinetic order in the 2D simulations. In contrast, 1D models predict a strong sensitivity because the kinetic trigger must keep the solution far from saturation. In the 2D case, flow localization maintains sufficient undersaturation even when the reaction rate slows near saturation, effectively bypassing the need for a kinetic trigger.

The study therefore proposes an alternative, geometry‑driven mechanism for cave initiation: the inherent instability of the dissolution front in a flowing system. This mechanism generates a characteristic spacing of dissolution channels that is set by the interplay of reaction rate and flow, not by pre‑existing fracture roughness. The findings have several implications. First, the existence of long conduits in natural karst does not necessarily require a sharp kinetic transition; the instability alone can produce deep penetration. Second, any realistic modeling of speleogenesis must move beyond 1D approximations and incorporate at least two‑dimensional flow‑reactive coupling. Third, the predicted wavelength provides a potential field observable—channel spacing—that could be used to infer subsurface reaction kinetics and hydraulic conditions. Finally, the concept of a self‑selected erosion pattern may be relevant to other geological and engineering contexts where reactive flow in fractures or porous media occurs, such as CO₂ sequestration, acid stimulation in oil reservoirs, and the weathering of engineered stone.

In summary, the paper overturns the prevailing 1D kinetic‑trigger paradigm by demonstrating that a planar dissolution front is fundamentally unstable. The resulting localized flow and dissolution enable deep, kilometer‑scale cave formation without invoking a special kinetic transition, thereby offering a more robust and physically grounded explanation for the development of extensive karst networks.