## Abstract

Precipitation of calcium carbonate from water films generates fascinating calcite morphologies that have attracted scientific interest over past centuries. Nowadays, speleothems are no longer known only for their beauty but they are also recognized to be precious records of past climatic conditions, and research aims to unveil and understand the mechanisms responsible for their morphological evolution. In this paper, we focus on crenulations, a widely observed ripple-like instability of the the calcite–water interface that develops orthogonally to the film flow. We expand a previous work providing new insights about the chemical and physical mechanisms that drive the formation of crenulations. In particular, we demonstrate the marginal role played by carbon dioxide transport in generating crenulation patterns, which are indeed induced by the hydrodynamic response of the free surface of the water film. Furthermore, we investigate the role of different environmental parameters, such as temperature, concentration of dissolved ions and wall slope. We also assess the convective/absolute nature of the crenulation instability. Finally, the possibility of using crenulation wavelength as a proxy of past flows is briefly discussed from a theoretical point of view.

## 1. Introduction

Water-driven morphological patterns are widespread in geological fluid mechanics, providing a plethora of examples, from glaciology [1–3] to river and sea morphodynamics [4,5] and karst environments [6,7]. One of the more remarkable cases concerns the deposition of sinters from free-surface flows and, in particular, the precipitation of calcite from dripping water [8,9].

Depending on many chemical and physical factors, the process of calcium carbonate deposition generates a wide variety of striking and sometimes weird morphologies, called speleothems. Remarkable examples are stalactites, stalagmites, draperies, helictites and flowstones. The current effort in studying and modelling the evolution of the karst environment is twofold. Firstly, understanding of the fundamental interactions between hydrodynamics and geochemistry that are capable of shaping mineral surfaces over a wide range of spatial scales (10^{−4}–10^{2} m) [10] is interesting. Speleothems are some of the most spectacular and interesting examples of natural beauty, and they are a clear example of the ability of a free surface flow to interact with chemical transport processes in order to carve a dissoluble substrate. Secondly, mineral deposition interacts with several industrial processes [11,12] and, more importantly, can be a precious proof of past climatic conditions. Measurable properties of the speleothems (mineral composition, calcite crystals type, etc.) can be related to the dripping water chemistry and to the surrounding environmental conditions in which the deposition occurred [13]. However, as pointed out by Fairchild *et al.* [14], speleothem records cannot yet be given a unique climatological meaning, and new methods for relating currently measurable properties with past conditions are constantly required.

The focus of this work is crenulations: ripple-like waves (figure 1) that are widely observed over stalactites and flat vertical surfaces [7,8] and characterized by a wavelength ranging from centimetres to a tens of centimetres. Despite their pervasive occurrence, crenulations have received little scientific attention, and only recently two of us [7] have demonstrated that crenulations are the result of an interfacial morphological instability of the dynamical system composed of the calcite surface and the water falling film. The aim of this work is therefore twofold. Firstly, we extend knowledge about the fundamental interactions between hydrodynamics and geochemistry by unveiling the physical and chemical mechanisms that induce the formation of crenulations. Secondly, we show that crenulations have a potential for being used as additional proxies for estimating the chemical composition and rate of past flows, thus playing an important role in paleo-climatic reconstructions.

The past studies that led to the current knowledge about the water-driven morphological evolution of calcite surfaces focused on three fundamental issues: (i) the dependence of the macroscopic calcite precipitation rate on water chemistry, (ii) the role of the hydrodynamic conditions of the precipitating solution, and (iii) the mathematical modelling of speleothem shape evolution.

The first issue was addressed using indirect methods: monitoring of the changes in the solution chemistry during calcite surface growth allowed the growth kinetics to be linked to chemical and physical parameters such as temperature, supersaturation state and pH. There are three major categories of models for crystallization: surface complexation models that take into account the reactions involving surface speciation [15]; summation of the elementary reaction theories that describe growth rate as a function of multiple elementary reactions [16]; chemical affinity models developed in terms of free energy changes of precipitation reactions [17]. In particular, the precipitation model of Plummer [16] was successfully tested against laboratory experiments [18] and is now considered to be the most reliable model for assessing the precipitation and deposition rates of calcite (and thus the corresponding morphological modification of the mineral surface) [19]. For this reason, the Plummer model is incorporated in the present modelling approach.

The second issue (role of the falling film hydrodynamics) was studied by coupling calcite precipitation models with chemical diffusion profiles typical of several flow conditions, and allowed understanding of how calcite precipitation occurs under static, laminar or turbulent films. [19,20]. More recently, numerical and experimental analysis are exploring the role of perturbations applied to laminar, free surface flow in developing instabilities that may ultimately lead to the formation of large-scale morphologies [21].

The third issue (mathematical modelling of speleothems shape evolution) allowed the mechanisms underling the formation of several speleothems to be elucidated. Calcite terraces have been widely investigated over the last decades, both from a theoretical [22] and a numerical [21] point of view, and a robust understanding of the involved processes has been achieved. Stalactites and stalagmites have also attracted the interest of scientists who have developed mathematical laws for the growth and the evolution of such speleothems [6,8,23]. At much lower spatial scales, researchers are also interested in unveiling the mechanisms involved in the formation and growth of dendrites in free surface flows [11]. The study of the interfacial stability of a calcite surface and a water film flowing on it has been addressed through an approximated depth-averaged hydrodynamic model only [24] (i.e. Dressler's equations), with no success in the prediction of crenulation formation.

In this work, we focus on the modelling of the interactions between (i) a thin-film containing a solution of calcium carbonate and flowing over an inclined surface, whose dynamics are described by the Orr–Sommerfeld equations (i.e. the linearized version of the Navier–Stokes equations), written for free surface flows, and (ii) a surface of calcite whose morphological evolution is driven by the geochemistry of the calcium carbonate–water ternary system. In a previous letter [7], it was demonstrated that crenulations are the result of a morphological instability of the film-surface interface. In this work, we pursue the following goals: (i) to describe a number of modelling and mathematical issues that in the preliminary letter were only briefly introduced. The explanation of the mathematical model and the techniques for its solution can be useful for the study of other karst morphologies. We also perform a minor algebraic correction to former results; (ii) to demonstrate the marginal role of carbon dioxide transport in the instability that leads to the formation of crenulations; (iii) to illustrate the chemo-physical mechanisms that induce the interface between the calcite surface and the water film to loose its stability, focusing in particular on the dynamics of the the free surface of the water film; (iv) to investigate the role of the different chemical and physical parameters on the crenulation dynamics and, in particular, how modifications of such parameters change the pattern wavelength; and (v) to study the convective and absolute nature of the crenulation instability.

## 2. Modelling aspects

The film-induced morphological evolution of a calcite surface is the result of the coupling between chemical and physical processes. The evolution of the wall morphology is in fact driven by the calcite deposition. The rate of calcite deposition is given by a number of different chemical reactions. The rate of such reactions depends on the reagent concentrations at the fluid–solid interface. Finally, such concentrations are influenced by the film flow field, which, in turn, is affected by the wall geometry. Therefore, a comprehensive mathematical description of such interactions requires to model (i) the chemical reactions occurring at the fluid–solid interface (rate of precipitation of calcite), (ii) the film flow dynamics, (iii) the transport of chemicals in the fluid film, and (iv) the morphological modifications induced by the deposition of calcite. In particular, we propose a mathematical model in which (i) the rate of precipitation of calcite is evaluated by the celebrated Plummer–Wigley–Parkhurst (PWP) equation [16], (ii) the laminar Navier–Stokes equations are used for modelling the flow dynamics, (iii) the convection–diffusion equations describe the transport of chemicals in the fluid, and (iv) the morphological modifications induced by a calcite flow are evaluated with an Exner–Stefan-like approach [5].

The model will allow us to pursue three fundamental goals. Firstly, we will demonstrate that crenulations are due to a morphological instability of the film-wall system. Previous studies gained important insights about the formation of large-scale speleothems, but were unable to explain crenulations formation due to the approximations in the flow or chemical adopted models [24]. Secondly, we will describe the physical and chemical mechanisms that drive the formation of crenulations. Finally, we will show the marginal role of carbon dioxide transport in the formation of such karst morphologies.

Henceforth, model equations are written in dimensionless form. The longitudinal unperturbed fluid velocity at the free surface, ^{−4}–10^{−1}]. Under the assumption of a quasi-flat bottom, the flow is laminar [25], so that *θ* is the bed mean slope with respect to the horizon (figure 2), *ν* is the fluid kinematic viscosity and *g* is the gravity acceleration. A dimensionless right-handed Cartesian frame {*x*,*z*} indicating the stream-wise and normal-to-the-bed coordinates, respectively, is adopted. Note that in order to address the fundamental dynamics of the crenulation instability, the cylindrical symmetry was assumed. The stream is bounded by the free surface *z*=*H*(*x*,*t*) and by the water–calcite interface *z*=*η*(*x*,*t*), where *t* denotes the time.

The concentration of calcium and bicarbonate will be scaled with the incoming calcium concentration *X*_{i}] denotes dimensional concentration in mol m^{−3} and subscript *i* refers to the *i*th chemical species. Carbon dioxide concentration will be scaled with its equilibrium value at the free surface *p*_{c} is the partial pressure. Accordingly, the chemical concentrations of the main species are *c*_{i} are the corresponding dimensionless concentrations.

### (a) Hydrodynamic model

Water film dynamics are described by the dimensionless Navier–Stokes and continuity equations for uncompressible flows
**u**={*u*,*w*} is the velocity vector, the subscript comma denotes partial derivation in all equations, ∇={∂/∂*x*,∂/∂*z*}, *P* is the pressure and **f**=2/*R*{1,−*cotθ*} is the body force vector.

The boundary conditions are [26]
*H* and *η* refer to the free surface and the bottom, respectively, and **T**=*PR***I**−2**D**, where **I** and **D** are the identity matrix and the rate of strain tensor, respectively. Finally, the Kapitza number [27] reads Ka=*l*^{2}_{c}*g*^{2/3}*ν*^{−4/3}, where *l*_{c} is the capillary length. Equations (2.3*a*–*c*) state (*a*) the kinematic condition, (*b*) the null shear stress and (*c*) the occurrence of surface tension at the free surface (*z*=*H*), while equations (2.3*d*–*e*) state (*d*) the impermeable surface and (*e*) the non-slip fluid velocity at the bottom (*z*=*η*).

### (b) Chemical model

In an open-system with respect to CO_{2}, the chemical reactions responsible for the removal or deposition of calcite are [28,29]
*a*) the dissociation of water into hydrogen- and hydroxyl ions, (*b*) the physical dissolution of carbon dioxide in water, (*c*,*d*) the pH-dependent conversion of carbon dioxide into hydrogen and bicarbonate, (*e*) the dissociation of bicarbonate into hydrogen- and carbonate ions and ( *f*–*i*) the dissolution of calcite, also dependent on the solution pH. Notice that the total amount of available CO_{2}, namely ^{2+} ion dissolved, one CO_{2} molecule is taken from the solution and such molecule is replaced from the carbon dioxide in the atmosphere, as the system is open with respect to CO_{2} by the free surface. Note that (2.5) holds true only for typical karst conditions. In fact, calcite dissolves at low pH, with the production of CO_{2}, and is precipitated at high pH, with the consumption of CO_{2}. Charge equilibrium is achieved if the concentrations of the different chemicals satisfy the following relation
*pH*<8; it follows that [OH^{−}] and ^{2+}] and *z*_{i} is the ionic charge [19]. By considering the ternary system *CaCO*_{3}–*H*_{2}*O*–CO_{2} and the relaxed electroneutrality condition, the ionic strength in natural karst water can be safely assumed as *I*=3[Ca^{2+}]. The activity (measured in mol m^{−3}) of the chemical species, (*X*_{i})=*γ*_{i}[*X*_{i}], is defined by the activity coefficient *γ*_{i}. For karst waters [29] in which *I*<10 mol m^{−3}, *γ*_{i} can be evaluated through the extended Debye–Hückel theory [30] as *I*<10^{−3}, and decrease to values between 0.6 and 0.9 for *I*∼10^{−2}, namely the ionic strength of a saturated solution of calcium carbonate in pure water at 298 K [29].

At the calcite surface, the dissolution or deposition processes are characterized by the three fundamental reactions (2.4*f*–*h*) that macroscopically lead to a flux of calcium ions. Dreybrodt [19] and Reddy *et al.* [31] found that a good quantification of both the dissolution and precipitation processes [32] is given by the Plummer–Wigley–Parkhurst (PWP) equation [16], which relates the calcium ion flux with the local activities of several chemical species evaluated at the calcite–fluid interface. The PWP relation reads
*z*=*η* indicates that activities are evaluated at the calcite–fluid interface, *κ*_{1} to *κ*_{3} are the same rate constants reported in (2.4*f*–*h*) and *κ*_{4} is a function of *κ*_{4}′, *κ*_{4}′′ and *κ*_{4}′′′. To obtain a tractable version of (2.8), one has to relate the solution pH to the incoming calcium concentration ^{2+}]. After some algebra, one obtains
*ρ*_{i} are reported in the electronic supplementary material. We stress that the approximation of relating the solution pH to the incoming calcium concentration *c* will demonstrate the dominant role of the hydrodynamic processes in the crenulation formation. This allows to disregard some chemical details in the calcite deposition modelling.

Consider now the morphological alterations induced on the wall by a depositing calcite flux. The calcite flux modifies the bottom elevation according to the relation *ϱ*∼10^{−4} m^{3} mol^{−1} is the molar volume of calcite. Such relation describes the dynamics of the wall, and plays the same role as the Exner and Stefan equations in other morphodynamic problems [2,5]. It must be noted that this macroscopic relation is valid only at dimensional scales larger than approximately 1 mm. Below this scale, the deposition of calcite cannot be longer considered uniform and does not lead to the formation of a smooth surface. Rather, the calcite surface consists of small, individual crystals protruding out [21]. Anyway, the fluid flows under laminar conditions. It follows that as long as the amplitude of the roughness is much smaller the thickness of the film, then the small protuberances plays no significant role in determining the flow resistance [33].

By introducing the scaling

### (c) Chemical transport model

The transport of the *i*th chemical dissolved in the laminar water film is described by the convection–diffusion–reaction equation
*c*_{i} is the concentration (we recall that subscript 1 and 2 refer to calcium and carbon dioxide, respectively), *r*_{i} is a source or sink term and *b*) reported a wrong sign). If ion-pairing is neglected, calcium does not react with other species and the source–sink term *r*_{1} is zero. At the free surface (*z*=*H*), no calcium is leaving or entering in the solution (no-flux, equation (2.12*a*)). At the wall (*z*=*η*), the flux of calcium is directly proportional to the calcite flux (2.12*b*). In particular, when calcium is leaving (entering) the water film—i.e. deposition, *f*>0 (dissolution, *f*<0)—then the calcium concentration increases (decreases) from the water–calcite interface towards the free surface of the water film.

The dimensionless net production of carbon dioxide can be written as
*k*_{1} and *k*_{2}. At the bottom, *z*=*η*, the carbon dioxide flux is null, see equation (2.12*d*), while at *z*=*H*, the concentration of dissolved carbon dioxide is in equilibrium with the gas partial pressure in the atmosphere (2.12*c*). The assumption of carbon dioxide in Henry's law equilibrium at the free surface is commonly adopted [29], and is acceptable because if we set a fixed concentration of carbon dioxide at the top layer, then the CO_{2} flux that exists from the water to the air is not hampered, as the gradient of the carbon dioxide concentration at the free surface is different from zero. This non-zero flux corresponds to the carbon dioxide that leaves the water solution (degassing). We also stress that the carbon dioxide concentration at the top of the water film is in Henry's equilibrium with the carbon dioxide partial pressure that occurs close to the water film. This value can be different from the CO_{2} concentration far from the free surface. Anyway, we will demonstrate (by a detailed sensitive analysis that spans orders of magnitude and reported in §4*e*) that crenulation formation is practically no sensitive to the external carbon dioxide partial pressure and, therefore, the assumptions introduced to model the carbon dioxide flux through the water air surface and its transport in air have a negligible effect on the global dynamics of the system.

## 3. Stability analysis

In the previous section, the mathematical model describing the film-induced spatio-temporal evolution of a calcite surface has been defined. It consists of (i) the Navier–Stokes equations (2.2) for the modelling of the water film dynamics, flanked by the boundary conditions (2.3); (ii) the equation (2.9)—which quantifies the flux of calcite at the calcite–liquid interface as a function of the calcium ions and carbon dioxide concentrations, evaluated at the calcite–liquid interface—and equation (2.10), which relates the temporal evolution of the bed with a given depositing calcite flux, and (iii) the transport equations (2.11), with the boundary conditions (2.12*a*–*d*), that model the chemical transport in the water film.

In the following, the temporal derivatives of all equations, except in equation (2.10), are disregarded. This hypothesis (so-called quasi-steady approximation) is customary in many morphological problems [5]. It holds true when the morphological temporal evolution of the wall is much slower than the temporal dynamics of the fluid. The rigorous justification for this hypothesis is reported in the electronic supplementary material.

In the real water–calcite karst system, one should also note that the time-averaged calcite deposition rate is not constant along *x*. In fact, the deposition of calcite occurring while the film is flowing causes the average concentration of calcium to be progressively reduced from the initial concentrations at the source point. However, the film degassing and the loss of calcium occur on longitudinal scales (greater than 10 *cm*) much larger than the crenulation scale [7]. Therefore, in order to study the local stability of the water–calcite interface, these non-parallel effects can be disregarded [34].

Finally, we perform the transformation of variables *ζ*=[*z*−*η*(*x*,*t*)]/*D*(*x*,*t*). In this way, the domain [*η*,*η*+*D*=*H*] is mapped in the rectangular domain [0,1]. This allows a correct setting of the boundary conditions at the free surface, that is a key point for the detection of the instability, as it will be shown in the following sections.

In order to study the stability of the flow-calcite interface to an infinitesimal perturbation, the dynamical system is forced with a harmonic disturbance of the bottom elevation
*η*_{0}(*t*) is the uniform growth of the bottom elevation resulting from the uniform calcite deposition *α* is the longitudinal wavenumber. The complex frequency reads *ω*=*ω*_{R}+*iω*_{I}, where (*ω*_{I}) and (*ω*_{R}/*α*) are the perturbation growth rate and celerity, respectively. Hereafter, subscript I (R) refers to the imaginary (real) part.

Equation (3.1) states that the deposition of calcite induces an upward translation of the bottom *η*_{0}(*t*). Such translation acts over very long times, i.e. *dη*_{0}(*t*)/*dt*<1. It follows that the alteration of the average bottom elevation does not induce any significant acceleration to the water film. From a mathematical point of view, the only effect of retaining *η*_{0}(*t*) in (3.1) is to consider a rigid upward translation of the whole system. Therefore, *η*_{0}(*t*) can be set to zero without any loss of generality.

The response of the governing equations to the bottom disturbance is then found by the following normal mode ansatz
*ϵ*<1), and the governing equations can be linearized around the basic state.

The unperturbed basic state can be found by collecting the terms *x* and over the time. In order to evaluate *i*=2) and its corresponding boundary conditions (2.12*c*–*d*), namely
*k*_{1} and *k*_{2} are much smaller than 1, the solution of the previous equation can be approximated to
*μ*=*k*_{1}−*k*_{2}(*a*_{1}+*a*_{2}). We finally recall that from the horizontal (*x*−) and vertical (*ζ*−) components of (2.2*a*) and the corresponding boundary conditions (2.3), the Nusselt's profiles for laminar flows

Let us consider now the first order by collecting the terms

In order to find *a*) has to be considered. It reads
*α*_{i}={*α*^{2},*α*^{2}−*k*_{1}}.

Let us consider first the transport of calcium (i.e. *i*=1). The linearized reactive term in the transport equation reads *a*,*b*) is *d*/*dz*→1/*D*(*x*,*t*) *d*/*dζ*. Implementation of this boundary condition is a fundamental step, but critical, as it involves also the unknown calcite flux perturbation

Let us consider now carbon dioxide (i.e. *i*=2). The reactive term of the transport equation (2.13*a*) is linearized to

It follows that

Finally, we substitute (3.9), (3.13) and (3.14) in the equation (3.6), obtaining the dispersion relation
*D*(*ω*,*α*) relates the growth rate (*ω*_{I}) and the celerity (*ω*_{R}/*α*) of the crenulations forming on the water–calcite interface to the wavenumber *α*. By just observing the dispersion relation, one recognizes that crenulation dynamics are influenced by (i) the water chemistry, through the term *a*–*b*) and depends on the wavenumber to a great extent (the complete procedure for the evaluation of

## 4. Results

### (a) Morphological instability of the water–calcite interface

For a given set of parameters, the linear dynamics of the fluid–calcite interface in the wavenumber space can be summarized by three scalar quantities. The first one is the maximum growth rate of the perturbation, namely

In figure 3, the growth rate is plotted as a function of the wavenumber for *R*=10^{−3} and for a fixed set of benchmark parameters, well representative of the karst environment. We consider the typical conditions **l** is the vector of control parameters and subscript ‘*B*’ denotes benchmark conditions. By setting *θ*_{B}=*π*/2, we firstly focus our analysis on vertical surfaces (e.g. stalactites) where crenulations are very likely to be observed.

The curve *ω*_{I}(*α*) exhibits two maxima, the first one at *α*≃10^{−3} and a second one at *α*≃10^{−1}. Recalling that the dimensional wavelength and the dimensionless wavenumber are related through ^{−2} m). Instead, the second maximum selects a pattern with a wavelength of the order of 0.5 mm, which is the typical size of calcite dendrites (several tens of micrometres [35]). Such second maximum falls in a zone where the validity of the proposed model is questionable, and it is here reported only for completeness. It can be observed, in fact, that the unstable zone on the right arises for wavelength smaller than approximately 2 mm. As anticipated in §2*b*, at such longitudinal scales the deposition of calcite cannot be considered uniform. In fact, the sinter surface consists of a multitude of small individual crystals protruding out. Therefore, for *α*>10^{−2} (*α*<10^{−2}, and we will not consider unstable zones occurring for *α*>10^{−2} to be physically relevant for crenulation dynamics.

In figure 4, the dependence of the crenulation dynamics on the Reynolds number is illustrated. Reynolds number varies from 10^{−4} up to 10^{−1}, a range that comprises most of the flow conditions found in the karst environment. It can be observed (figure 4*a*) that the maximum growth rate is always positive (the calcite-flow interface is always unstable to crenulations) and it decreases when the flow rate (i.e. *R*) increases. In figure 4*b*, the dynamics of the crenulation phase velocity can be observed. In vertical surfaces, crenulations always migrate downstream (*R*. Figure 4*c* demonstrates the very small influence of *R* in defining the most amplified wavenumber: while *R* spans three orders of magnitude (*R*∈[10^{−4},10^{−1}]), ^{−4}. Dimensionally, instead, while the unitary flow *q* spans in the range [10^{−4},10^{−1}] *l* *h*^{−1} m^{−1}, the pattern wavelength (i.e. the most amplified wavelength *d*).

These results confirm the ability of the model to capture the typical wavelengths observed in karst environment. For *R*∼10^{−3}, it predicts the occurrence of centimetric waves, namely crenulations, characterized by downstream migration (in vertical surfaces) and a relatively fast growth rate. For *R*∼10^{−1}, the model predicts the formation of longer waves (approx. 10^{−1} m) that migrate downstream, and with a growth rate much lower than crenulations. Also these dune-like long waves (figure 1*d*) are common morphologies on flowstones (i.e. the inclined side walls of the caves).

### (b) The marginal role of carbon dioxide transport

We now explore the role played by the chemical and physical phenomena responsible of the morphological instability of the water–calcite interface and of the selection of the most unstable wavelength. To this aim, we numerically quantify the different terms appearing in the dispersion relation (3.15).

Let us first compare the numerical value of coefficients *ρ*_{1}, *ρ*_{2} and *ρ*_{3}, which link the rate of calcite deposition to the bottom concentrations *c*_{1} and *c*_{2} by the PWP equation (2.9). By considering the typical numerical values of these coefficients (table 1), it can be observed that the term (*ρ*_{1}+2*ρ*_{3})≃10^{−5} exceeds of three orders of magnitude the term *ρ*_{2}≃10^{−8}. We now compare the magnitudes of terms

We have now the elements to detect the dominant terms of the dispersion relation (3.15). Let us first consider the numerator. The real and the imaginary part of *χ* while its imaginary part is 10^{6} times smaller than the imaginary part of the other term

Summing up, the terms related to the carbon dioxide transport *ρ*_{2} that is very small, and the dispersion relation can be approximated as
*ω*_{I}) and celerity (*ω*_{R}/*α*) are proportional, through the coefficient (*ρ*_{1}+2*ρ*_{3}), to the perturbation of calcium concentration at *ζ*=0,

We finally recall that after using (3.14), equation (4.1) can be rearranged in the final form
*A*=(*ρ*_{1}+2*ρ*_{3})/*χ*, showing that the growth rate is proportional to the film depth perturbation

### (c) The key mechanisms driving crenulation formation

Relation (4.1) shows that the formation of crenulations can be understood by investigating how the film depth perturbation, *a*).

Preliminary, we have to elucidate how the free surface and, in turn, the water depth respond to a bottom perturbation. We, therefore, focus on the phase angle, *ϕ*, between the flow depth perturbation induced by a bottom topography and the topography itself. Such angle is linked to the location of the minimum flow depth (*c*) and three cases can occur (i) the most severe shallowing occurs over the downstream side of a bump (i.e. where *π*/2,*π*[; (ii) *ϕ* falls in the range ]*π*,3*π*/2[; and (iii) *π*/2] or [3*π*/2,2*π*], depending on the sign of *d*) that for the typical values *R*=10^{−3} and *θ*=*π*/2 occurring in karst systems, the phase angle between *π*/2, for all the wavenumber considered. Therefore, the most severe film thinning always occurs over the downstream side of the bumps. The black dotted (black dashed) line in figure 5*d* shows that a reduction (increase) of the Reynolds number induces only a small downward (upward) translation of the phase-wavenumber curves. Also a reduction of the wall slope (dotted grey line) does not affect significantly the behaviour. Instead, in the case of overhanging walls, it can be observed (dashed grey line) that for low wavenumbers, the phase angle switches from values slightly greater than *π*/2 to values slightly smaller than 3*π*/2. Physically, the most severe film shallowing occurs on the upstream side of the protrusions. The switching of the position of min[*θ* exceeds *π*/2 causes several changes in the macroscopic behaviour of the crenulations pattern (e.g. the switch of the migration direction) that will be described in the following sections.

We now investigate how the calcium concentration at the film bottom is altered by a water depth perturbation, *f* reports the spatial behaviour over one wavelength of the perturbed water depth *e*, right part), the corresponding calcium profile (black line) becomes more steep with respect to the unperturbed profile (grey line, see also figure 5*a*). The consequence of a steeper calcium profile (see dashed black and grey lines in figure 5e) is that more calcium ions are pumped through diffusion from the film towards the wall. This high flux of calcium ions, in turn, rises the calcium concentration at the wall. It follows that the calcite flow rate results increased (see linearized PWP equation (3.6) and also the bottom elevation growth rate increases (see the bottom evolution equation (2.9)). The opposite happens if an increase of the water depth occurs (*e*, left part): the calcium profile becomes less steep than the unperturbed calcium profile and the concentration

At the beginning of this section, we have shown that the minimum water depth always occurs over a bump (*concentration mechanism*, in order to underline the concentration of calcium at the bottom induced by the film thinning. We recall that hydrodynamically induced alterations of the chemical gradients were demonstrated to play a key role in the generation of other calcite morphologies [21].

Consider now the feedback of the calcite flux on the calcium concentration. In figure 5*h*, the spatial behaviour over one wavelength of the calcite flux, *g*) and, in particular, near the wall. According to equations (3.6) and (2.9), also the calcite flow rate and the bottom elevation growth rate are reduced. The opposite happens if a lower rate of calcite deposition occurs ( *g*, left): the calcium concentration increases and, in turn, also the calcite flow rate *depletion mechanism*, plays a stabilizing role. In fact, at the locations in which an excess of calcite precipitation occurs, a reduction of the calcium concentration is induced. Such reduction of calcium concentration, in turn, halts the calcite deposition excess.

Examine now the full response of the system to perturbations *i* and 6. The film responds to an initial bottom perturbation (I) with a film depth perturbation (II). Recall (figure 5*d*,*i*) that the maximum film thinning (i.e. min[

The same mechanisms are also found if the convection of chemicals is restored. In such case, we observe only a slight downstream phase lag between the water depth perturbation and the calcite perturbed flow. Anyway, the key processes inducing the instability—concentration mechanism and competition with the depletion mechanisms—remain the same.

### (d) The key role of the wall slope

In the previous section, we have demonstrated the crucial role of the film fluid dynamics in determining the stability of the calcite-flow interface. It follows that, being the flow characteristics (in particular, the response of the film depth to a bottom perturbation) influenced to a great extent by the Reynolds number and by the wall average slope, these quantities impact the dynamics of crenulations. In figure 7*a*, the marginal stability curves (i.e. points where the growth rate *θ*−*R* plane. It can be observed that for relatively high Reynolds numbers (*R*>10^{−2}) crenulation formation is only possible on walls almost vertical, or steeper (*θ*>0.9*π*/2) while for progressively lower values of *R* the instability occurs even for significantly lower slope, so that crenulations can appear when the wall average slope is as low as *π*/4. With a further reduction of *R*, the system becomes stable to crenulations and, for a given *θ*, the marginal stability curves exhibit a threshold value *R*_{c} that separates the zone of stability (*R*<*R*_{c}) from the zone of instability (*R*>*R*_{c}). Note that the critical Reynolds number, *R*_{c}, above which crenulations form is highly dependent on temperature and water chemistry, as shown by the different lines of figure 7*a*–*b*, evaluated by changing some parameters from the benchmark condition.

In figure 7*b*, the locus of the points that have a null phase velocity, namely *θ*>*π*/2), a very common condition found in stalactites. This confirms the result obtained by Camporeale & Ridolfi [7], whereas the downstream migration for *θ*<*π*/2 is instead a novel finding.

The effect of slope is now investigated by the ratios *θ*=*π*/2. Such ratios are evaluated for different Reynolds numbers and only for unstable conditions (i.e.

It can be observed (figure 7*c*) that lowering the mean slope of the calcite wall induces a strong reduction of the growth rate that eventually becomes negative. On the contrary, in overhanging slopes the growth rate is amplified of more than one order of magnitude. This feature is more pronounced for high Reynolds numbers. The wall slope plays a strong influence also in selecting the most unstable wavelength (figure 7*d*): for low Reynolds numbers (*R*<10^{−3}) the pattern wavelength increases up to five times by reducing the slope below *π*/2, while weaker changes are observed if *θ* is increased above *π*/2. For higher Reynolds numbers (*R*>10^{−2}) the pattern wavelength is weakly influenced by a reduction of the slope below *π*/2, while it reduces up to one-fourth by increasing *θ* above *π*/2.

Finally, the wall slope affects the crenulation phase velocity (figure 7*e*). We have previously demonstrated (figure 7*b*) that the crossing of the wall slope threshold *θ*=*π*/2 induces a switch in the direction of migration of crenulation. Figure 7*f* shows that high Reynolds numbers (*R*∼10^{−1}) induce a dramatic increment in the absolute value of the phase velocity.

### (e) Impact of water chemistry on crenulation dynamics

The effect of water chemistry and temperature is here investigated systematically. To this aim, we follow the same approach used in the previous section by focusing on the ratios *l*_{B}, to the new value, *l*. The ratios are evaluated by changing the parameters of the benchmark set one by one. The analysis has been repeated for different Reynolds numbers, but no appreciable differences were observed, and all the curves evaluated for different *R* collapse on the same line.

Figure 8*a* shows that an increase in the concentration of calcium, ^{2} leads to an increment of the maximum growth rate, *p*_{c} increases, the maximum growth rate *b*). From a physical point of view, the higher is the degree of over-saturation of the solution, the faster crenulations form and it can also be observed (line *r*_{α}) that the most amplified wavenumber undergoes significantly changes (the ratio *r*_{α} is in the range [0.5,1.6]). It follows that the chemical composition of the incoming water has an appreciable effect in selecting the dominant wavelength. Finally, the phase velocity of the most unstable wavenumber has a behaviour very similar to the growth rate; in fact,

Let us examine now the role of thermal changes. The temperature modifies a number of physical and chemical parameters of the system (water density, viscosity, surface tension, equilibrium constants of the chemical reactions, solubility products, etc.), but the main effects are a change in the solubility of calcite (that decreases with temperature) and in the solubility of carbon dioxide (decreasing with the temperature, as film degassing is promoted). The final result is therefore an increment of the growth rate of the instability with temperature. From a quantitative point of view, variations of temperature of a few tens of degrees cause changes of the maximum growth rate and of the phase velocity up to two order of magnitudes, while the selection of the most amplified wavenumber is less influenced (*r*_{α} remains in the range [0.9,1.5], figure 8*c*). Such response of the system to changes of the water chemistry and temperature is not surprising. In §4*c*, we have in fact demonstrated that the dilution mechanism is the ultimate responsible of the destabilization of the film–calcite interface. In §4*b*, we have also shown that the strength of dilution mechanism is proportional to *p*_{c} and temperature *T*.

It is now instructive to investigate the role of the water chemistry activity in the crenulation dynamics. In particular, we consider the introduction of other chemical species in the *b*) that the introduction of new ions alters the solution total activity *I*. In particular, the addition of the chemical *i* induces in the solution activity an increment *X*_{i}] of ions that have high ionic charge (e.g. *z*_{i}=3−7) can significantly increase the solution activity and reduce the activity of the single different chemicals. In order to quantify the role of additional ions, we introduce a coefficient *Υ*_{I}≥1 that multiplies the solution activity *I*=3[Ca^{2+}] adopted in the previous analysis. In this way, we model an increment of the solution activity caused by the presence in the water of generic elements other than calcium- and bicarbonate ions and water. It can be observed (figure 8*d*, bold line) that an increment of solution activity (increasing *Υ*_{I}) has a moderate influence on the maximum growth rate *b*). Physically, the increment of the solution activity causes the reduction of the activity coefficients *γ*_{i} and, despite the concentrations [*X*_{i}] are unaffected, the activity (*X*_{i}) of the chemicals involved in calcite deposition is reduced. In turn, the unperturbed calcite deposition rate *I* cause changes of *d*). Also the most amplified wavenumber (figure 8*d*, dashed line) results poorly affected by the solution activity, being *r*_{α} in the range [0.75,1].

Finally, the effect of the alteration of the liquid surface tension is explored. This point is not practically relevant, since a reduction of water surface tension is only possible through the action of a surfactant. Nevertheless, this analysis is theoretically instructive as it further clarifies the significance of the film flow dynamics on the crenulation development. To this aim, we introduce the coefficient 10^{−1}≤*Υ*_{γ}≤10^{1} that multiplies the surface tension *γ* (accordingly, the capillary length *l*_{c} is multiplied by *e*, bold line) as well as a relevant impact on the selection of the most amplified wavenumber (figure 8*e*, dashed line). This is a further proof that film flow dynamics, and in particular, the free surface response to a bottom topography, plays a key role also in selecting both wavenumber of the pattern (a common feature in many free surface flow-induced morphologies [36]) and the pattern growth rate.

## 5. Convective/absolute nature of the crenulation instability

The assessment of the convective/absolute nature of an instability is a key information in unstable open flow hydrodynamic systems [34], as well as in morphodynamic problems [37,38]. A system is convectively unstable if the response to an impulsive perturbation increases in time but migrates and decays to zero at all the spatial locations, and it is absolutely unstable if the response grows exponentially in time at all spatial locations. Convectively unstable systems behave as ‘noise amplifiers’, displaying extrinsic dynamics, as in the absence of continuous forcing the response decays back to zero, whereas absolutely unstable systems are characterized by intrinsic dynamics and behave as ‘oscillators’ [39,40]. From a practical point of view, the knowledge of the type of instability is fundamental for correctly understanding field observations and for an appropriate set up of numerical simulations and experiments.

In order to discriminate between these two kinds of instabilities, one needs to determine if the wave with the zero group velocity is growing (absolute case) or decaying (convective case). This is an easy task only in a limited number of cases, where the dispersion relation is provided analytically. Unfortunately, this is not the case for crenulations, which require a refined modelling of the flow field along the non-homogeneous vertical direction. In order to circumvent this difficulty, we will apply the so-called cusp map method [38,41] to the present problem.

By definition, the complex absolute wavenumber *α*_{0} characterizes the long time behaviour of equation (3.15) at *x* fixed, and it corresponds to the solution of the saddle point condition, *D*(*ω*,*α*)=*ω*_{,α}(*α*)=0, provided that the causality principle is satisfied, i.e. that *α*_{0} is a pinch point. Moreover, *ω*_{0I}=Im[*ω*(*α*_{0})] is the associated absolute growth rate. If only real wavenumbers are considered, equation *ω*_{I,α}=0 is satisfied by the wavenumber *ω*_{0I}<0; on the contrary, the flow is absolutely unstable if *ω*_{0I}>0. Recalling that the *j*th spatial branch of the dispersion relation *ω*_{,α}=0 are pinch points only if at least two spatial branches *α*_{m}(*ω*_{0I}) and *α*_{n}(*ω*_{0I}) pinching in *α*_{0} are well confined within opposite *α*_{I} half-planes when *ω*_{I} is increased [40].

The pinch point criterion usually requires a mapping from the complex frequency plane to the complex wavenumber plane, thus the dispersion relation has to be solved for *α* as a function of *ω*. Unfortunately, in many physically relevant cases (e.g. the water–calcite interface evolution herein investigated), the dispersion relation is a transcendental function of *α* while it is polynomial only in *ω*. In order to circumvent this difficulty, Kupfer *et al.* [41] refined a work by Derfler [42] and developed a technique—called cusp map method—for the assessment of the convective/absolute nature of the instability requiring only a *α*→*ω* mapping.

The cusp map method basically follows two conceptual steps (i) to detect the points that have null group velocity in the complex *ω*-plane and (ii) to determine whether these zero group velocity complex frequencies are actually pinch points in the complex wavenumber plane or not. The former task is accomplished recalling that (in the frequency plane) a point *ω*_{0} that satisfies *D*(*ω*_{0},*α*)=*D*_{,α}(*ω*_{0},*α*)=0 (the saddle point condition) and *D*_{,αα}(*ω*_{0},*α*)≠0, has a local map (*ω*−*ω*_{0})∼(*α*−*α*_{0})^{2}. From a topological point of view, this implies that when a curve lying in the complex *α*-plane and passing trough *α*_{0} is mapped into the complex *ω*-plane, it displays a cusp-like singular point at the branch point *ω*_{0}. Kupfer [41] suggested to map the contour *α**_{I} from zero. The mapping of all the contour points with *Ω*(*α**_{I}). A branch point is obtained when *Ω*(*α**_{I}) displays a cusp-like singularity (occurring at *ω*_{0}), and this occurs exactly when *α**_{I}=*Im*(*α*_{0}). Finally, the check for the pinch point condition is performed by considering the relative position of *ω*_{0} with respect to *α**_{I}=0 into the the *ω*-plane (i.e. the usual temporal branch). It can be demonstrated that branch points *ω*_{0} that are ‘covered’ an odd number of times by

The results obtained by the application of the cusp-map method to the crenulation instability are provided in figure 9 for different values of *R* and *θ*. The curve *Ω*(*α**_{I}) are reported in the *ω*-plane, with *α**_{I} progressively lowering from zero. Left panel refers to the case *θ*=*π*/2 while the central panel refers to (*θ*=1.1⋅*π*/2).

For *θ*=*π*/2, lowering *α**_{I} from 0 causes a progressive lowering of the contours until the whole curve *Ω*(*α**_{I}) lies in the half plane *ω*_{I}<0. Even though *α**_{I} is progressively reduced, no singular points (clue of a zero-group-velocity wave) can be found. It follows that for *θ*=*π*/2 there are no unstable waves (instability would require *Ω*(*α**_{I}) to lie in the half plane *ω*_{I}>0) that exhibit null group velocity. Therefore, the calcite-flow interface is convectively unstable.

For *θ*=1.1⋅*π*/2, lowering *α**_{I} from zero causes a progressive sharpening of the contours *Ω*(*α**_{I}) that ultimately leads to the occurrence of a singular point at *ω*_{0}. In order to show more clearly the selected cusp, an enlarged view of the zone where the cusp is forming (white box in the central part surrounding the branch point *ω*_{0}) is reported in the right part of each panel. Once the branch point is detected, its position is compared with respect to *ω*_{0} is ‘covered’ only once by *ω*-plane therefore correspond to pinching points in the *α*-plane. It follows that, being *ω*_{0I}>0, the calcite-flow interface is absolutely unstable for any *R* here considered.

The dependence of the nature of the instability on *θ* has been thoroughly scrutinized, revealing that *θ*=*π*/2 is a critical value for any value of the other control parameters: for *θ*>*π*/2 crenulations are absolutely unstable while for for *θ*≤*π*/2 the instability is convective.

This is a key result. From the definition of convective instability and from the behaviour of other convectively unstable systems [43], it is in fact known that exists a critical length *L*_{c} of the domain below which the macroscopic detection of the instability is not possible. Moreover, the lower is the perturbation growth rate, the higher becomes this critical length *L*_{c}. We also recall that variations of the chemo-physical control parametrs (e.g. a wall slope reduction, figure 7*c*) can cause dramatic reductions of the growth rate of crenulations. It follows that in some cases (e.g. wall slope below *π*/2, little over-saturation of the water), the critical length *L*_{c} is very large. In particular, *L*_{c} can exceed the speleothem size. As a consequence, the pattern of crenulations cannot be macroscopically detected. This justifies the difficulty of observing crenulations in some caves (little over-saturation of the water) or in some specific speleothems (wall slope below *π*/2).

Consider now the emblematic case of one stalactite (wall slope above *π*/2) hanging over its associated stalagmite (*θ*<*π*/2): the total flow rates of the films dripping on both speleothems are equal. Moreover, stalagmites usually have a radius much larger than the associated stalactite. It follows that stalagmites, in principle, exhibit smaller values of Reynolds number than stalactites. These conditions lead the stalagmites to fall either in the stable region of figure 7*a* or, at most, in the convectively unstable regime but with a very low growth rate (*θ*<*π*/2), both cases being unfavourable to pattern detection. However, provided a sufficient length, the pattern detection remains possible, as testified by the occurrence of regular dune-like patterns on long flowstones.

## 6. Discussion and concluding remarks

In this work, we have presented the results of a detailed physical- and chemical-based mathematical model that predicts the formation of crenulations as the result of a morphological instability of the water film-calcite interface. We have firstly demonstrated that carbon dioxide transport can be neglected, and that only the concentration of calcium at the wall plays a relevant role in determining the stability of the system, namely the occurrence or not of crenulations. This aspect allows the development of future simplified modelling approaches that only focus on calcium transport, even though the complete set of geochemical reactions is fundamental for obtaining the basic state solutions. The possibility of disregarding the carbon dioxide transport allowed us to develop a simplified dispersion relation that shows the physiochemical processes involved in the crenulations instability. In particular, we have distinguished two main processes that affect the calcium concentration at the wall: the *depletion* and the *concentration* mechanisms. The former is related to the consumption of calcium ions near the wall, as a result of calcite deposition, and is showed to be stabilizing. The latter is instead destabilizing, and it is related to the alteration of the calcium concentration depth-profile induced by the free surface dynamics. Similar mechanism where also found in the analysis of the formation of ripples over icicles [2,44], where the heat flux is the key process driving the morphological evolution of the ice wall. We also suspect that the ripple-like patterns that often shape silica sinters [45] are induced by the fluid free surface, but a rigorous analysis that explains how the film dynamics interact with the cooling processes responsible of the sinters formation [46] is still lacking.

We have then investigated the role of many physical (e.g. Reynolds number, average slope) and chemical (e.g. calcium ion concentration, carbon dioxide partial pressure and temperature) parameters involved in the crenulation dynamics, demonstrating that their variations in realistic ranges cause very significant changes in some key control features. A synthesis of the key results is reported in table 3, where the effect of the increment of a given parameter (top row) on the growth rate, phase velocity and most unstable wavenumber is reported. It is a key point to note that the wavelength of the instability (that is easily measurable and stored in speleothems stratigraphy) is influenced especially by the Reynolds number and by the wall slope. Such strong dependence of the pattern wavelength on *R* and *θ* opens the possibility to associate a particular wavelength to past flow conditions that occurred at the time of the calcite deposition. Differently, many parameters affects the growth rate. In particular, we observe high growth rate only for very over-saturated water. It should therefore not be surprising that in several caves where water over-saturation is weak, crenulations are not visible: the film-calcite interface is unstable, but the amplitude of the crenulation pattern grows so slowly that cannot be detected. Anyway, nonlinear analysis are required for gaining further insights about the temporal evolution of the crenulation amplitude.

The convective/absolute nature of the instability has finally been assessed: crenulations are invariably absolutely unstable for *θ*>*π*/2 (typical wall slope in stalactites) while they are convectively unstable for *θ*≤*π*/2. The evaluation of the convective/absolute nature of the instability allows the assessment of the robustness of the results obtained from these palaeo-hydraulic reconstructions. In particular, the use of convectively unstable crenulations’ features requires some caution, because they can result from the amplification of some external disturbance rather than the mirror of an intrinsic dynamics that can be well associated to some flow or environmental characteristics [39]. Differently, the reconstructions of palaeo-flows well fit with the case of absolutely unstable crenulations. Anyway, the actual use of crenulation wavelength for the interpretation and estimation of palaeo-flows requires further experimental efforts to validate the theoretical outcomes here presented.

## Data accessibility

This manuscript does not contain primary data and as a result has no supporting material associated with the results presented.

## Author contributions

All the authors equally contributed in developing the the mathematical model, performing the analysis, interpreting the results and writing the paper. All authors gave final approval for publication.

## Funding statement

This work was supported by the Politecnico di Torino.

## Conflict of interests

The authors have no competing interests.

- Received January 16, 2015.
- Accepted February 19, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.