2. Formulation

Let us consider a pendent film which flows along the underside of a wall. For the sake of simplicity, the domain is separated into two parts: a curved and a flat region. The upstream curved region is assumed to have a cylindrical shape, with centre of curvature $O$ and constant radius of curvature $\tilde{\mathscr{R}}$ , and it is joined to the flat region downstream (see figure 2). In the curved region, we will adopt a cylindrical reference frame $\{{\it\theta},\tilde{r},{\tilde{y}}\}$ , corresponding to the azimuthal, radial and spanwise coordinates, respectively. Note that the curvilinear coordinate is $\tilde{s}=\tilde{\mathscr{R}}{\it\theta}$ , so the alternative reference frame $\{\tilde{s},\tilde{r},{\tilde{y}}\}$ will also be adopted in the following. In the flat region, we will use the standard Cartesian reference frame $\{\tilde{x},{\tilde{y}},\tilde{z}\}$ , corresponding to the streamwise, spanwise and wall-normal coordinates, respectively. In this latter region, the slope of the wall with respect to the vertical axis reduces to a constant value, i.e.  ${\it\theta}\rightarrow {\it\vartheta}$ (see figure 2) and, for unperturbed conditions, the depth and the surface velocity are respectively provided by Nusselt’s solutions

(2.1a,b ) $$\begin{eqnarray}\tilde{D}_{0}=\left(\frac{2{\it\nu}^{2}\,\mathit{Re}}{g\cos {\it\vartheta}}\right)^{1/3},\quad \tilde{U} _{0}=\left(\frac{g\cos {\it\vartheta}\,\mathit{Re}^{2}\,{\it\nu}}{2}\right)^{1/3},\end{eqnarray}$$

where $g$ is the gravitational acceleration, ${\it\nu}$ is the viscosity and $\mathit{Re}=\tilde{U} _{0}\tilde{D}_{0}/{\it\nu}$ is the Reynolds number. Henceforth, tilde refers to the dimensional notation.

In order to tackle the morphological instability, a wall evolution equation is necessary. With this aim, it is convenient to focus on the flat region and to decompose the film thickness into the sum of the Nusselt depth, $\tilde{D}_{0}$ , and a (small) local pattern-induced perturbation, ${\it\epsilon}\tilde{d}(\tilde{x},{\tilde{y}},\tilde{t})$ . The quantity ${\it\epsilon}$ is a small perturbative parameter that represents the amplitude of the perturbation, so $\tilde{d}$ is of order one. At length scales much larger than a single crystal, it can be shown that the growth rate of both calcite–water and ice–water interfaces is linearly related to the depth perturbation ${\it\epsilon}\tilde{d}$ .

Let us first consider the karst environments. In the work by Short et al. (Reference Short, Baygents, Beck, Stone, Toomey and Goldstein2005a ), the stalactite shape was derived by demonstrating that the calcite deposition rate can be written as

(2.2) $$\begin{eqnarray}{\it\varrho}\frac{\text{d}\tilde{{\it\eta}}}{\text{d}\tilde{t}}=a\{\tilde{D}_{0}+{\it\epsilon}\tilde{d}(\tilde{x},{\tilde{y}},t)\}\{k_{-}^{\prime }[\text{Ca}^{2+}]-k_{+}H[\text{CO}_{2}]_{\infty }+k_{0}[\text{H}^{+}]\},\end{eqnarray}$$

where ${\it\varrho}$ is the ratio of molar mass to density for calcite, $\tilde{{\it\eta}}$ is the wall perturbation, square brackets refers to ion concentration, $H$ is Henry’s constant, $k_{-}^{\prime }$ , $k_{+}$ and $k_{0}$ are pH-dependent reaction constants, and the subscript of the carbon dioxide concentration refers to far atmospheric values. The coefficient $a$ accounts for geometrical effects but is very close to unity for thin films, so it can be omitted in the present context (as also remarked by Short et al. (Reference Short, Baygents and Goldstein2005b )). The term related to calcium concentration is formally space-dependent; however, as already pointed by Short et al. (Reference Short, Baygents, Beck, Stone, Toomey and Goldstein2005a ), the calcium concentration profile is nearly flat, because of the very low deposition rate, so there is no significant reduction of calcium along the length of the speleothem. On the other hand, if one were interested in modelling sub-centimetre transverse ripples (i.e. crenulations), it would be important to consider calcium concentration as harmonically perturbed in the streamwise direction, as pointed out by Camporeale & Ridolfi (Reference Camporeale and Ridolfi2012a ). For our purposes, the whole content of the second curly brackets – which we denote ${\it\varrho}\mathscr{F}_{K}$ – can be considered instead independent of hydrodynamics. In this way, (2.2) can be rewritten in the generalized form of an evolution equation,

(2.3) $$\begin{eqnarray}\frac{\partial \tilde{{\it\eta}}}{\partial \tilde{t}}=A+{\it\epsilon}B\tilde{d},\end{eqnarray}$$

where $A=\mathscr{F}_{K}\tilde{D}_{0}$ and $B=\mathscr{F}_{K}$ .

In the ice context, the master equation for wall evolution is instead the celebrated Stefan equation, which states the heat balance at a freezing–melting interface and reads

(2.4) $$\begin{eqnarray}{\it\rho}_{i}\mathscr{L}\frac{\text{d}\tilde{{\it\eta}}}{\text{d}t}=[({\it\kappa}_{i}\boldsymbol{{\rm\nabla}}\tilde{T}_{i})_{I^{-}}-({\it\kappa}_{w}\boldsymbol{{\rm\nabla}}\tilde{T}_{w})_{I^{+}}]\boldsymbol{\cdot }\boldsymbol{n}_{I},\end{eqnarray}$$

where ${\it\kappa}_{w,i}$ are the thermal conductivity of water and ice, respectively, ${\it\rho}_{i}$ is ice density, $\mathscr{L}$ is the latent heat, $\boldsymbol{n}_{I}$ is the normal versor to the liquid–solid interface, and $\tilde{T}_{w,i}$ are respectively water and ice temperature. Similarly to the depth, these latter can be decomposed into the sum of a base state plus a pattern-induced perturbation, namely $\tilde{T}_{w,i}=\tilde{T}_{0_{w,i}}(\tilde{z})+{\it\epsilon}\tilde{{\it\Theta}}_{w,i}(\tilde{x},{\tilde{y}},\tilde{z})$ . Additionally, we can introduce the normalized vertical coordinates ${\it\zeta}=(\tilde{z}-\tilde{{\it\eta}})/\tilde{D}$ and ${\it\zeta}_{i}=(\tilde{z}-\tilde{{\it\eta}})/\tilde{S}$ , where $\tilde{S}$ is an assigned depth in the ice. In this way the liquid and the solid domains are rectangularized so that ${\it\zeta}\in [0,1]$ and ${\it\zeta}_{i}\in [-1,0]$ , respectively. One therefore obtains that $\tilde{T}_{0_{w,i}}$ are linearly distributed over the film depth, i.e. 

(2.5a,b ) $$\begin{eqnarray}\tilde{T}_{0_{w}}=\tilde{T}_{S}{\it\zeta},\quad \tilde{T}_{0_{i}}=\tilde{T}_{b}{\it\zeta}_{i},\end{eqnarray}$$

where $\tilde{T}_{S}$ is a prescribed free surface temperature (in equilibrium with air), and $\tilde{T}_{b}$ is the ice temperature at $\tilde{z}=-\tilde{S}$ (i.e.  ${\it\zeta}_{i}=-1$ ). Since the Prandtl number of water at melting conditions is equal to $\mathit{Pr}=13$ and $\mathit{Re}$ can be assumed to be of order 10⁻² or less, the Péclet number, $\mathit{Pe}=\mathit{Pr}\,\mathit{Re}$ , is $O(10^{-1})$ or less. This justifies the neglect of thermal convection in the heat equation of the liquid phase. The perturbations $\tilde{{\it\Theta}}_{w}$ and $\tilde{{\it\Theta}}_{i}$ therefore satisfy the Laplace equations ${\rm\nabla}^{2}\tilde{{\it\Theta}}_{w,i}=0$ , with homogeneous boundary conditions that ensure the melting temperature at the liquid–solid interface and the constancy of the overall heat flux throughout the domain (i.e. the heat flux induced by perturbations is zero). With the use of the normalized coordinates, after linearizing with respect to ${\it\epsilon}$ and Fourier-transforming from $\{\tilde{x},{\tilde{y}}\}$ to the horizontal wavevector $\tilde{\boldsymbol{k}}\equiv \{\tilde{{\it\alpha}},\tilde{{\it\beta}}\}$ , where $\tilde{{\it\alpha}}$ is the streamwise wavenumber and $\tilde{{\it\beta}}$ is the spanwise wavenumber, the diffusive thermal problem becomes

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\frac{\text{d}^{2}}{\text{d}{\it\zeta}^{2}}-\tilde{D}_{0}^{2}\tilde{k}^{2}\right]\tilde{{\it\Theta}}_{w}^{\ast }+\tilde{D}_{0}\tilde{k}^{2}\tilde{T}_{s}(\tilde{D}_{0}+\tilde{d}^{\ast }{\it\zeta})=0, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\frac{\text{d}^{2}}{\text{d}{\it\zeta}_{s}^{2}}-\tilde{S}_{0}^{2}\tilde{k}^{2}\right]\tilde{{\it\Theta}}_{i}^{\ast }+\tilde{S}_{0}\tilde{D}_{0}\tilde{k}^{2}\tilde{T}_{b}(1+{\it\zeta}_{s})=0, & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{{\it\Theta}}_{w}^{\ast }(0)=\tilde{{\it\Theta}}_{i}^{\ast }(0)=\left.\frac{\text{d}\tilde{{\it\Theta}}_{w}^{\ast }}{\text{d}{\it\zeta}}\right|_{{\it\zeta}=1}=\left.\frac{\text{d}\tilde{{\it\Theta}}_{i}^{\ast }}{\text{d}{\it\zeta}}\right|_{{\it\zeta}=-1}=0. & \displaystyle\end{eqnarray}$$

$\tilde{S}_{0}=\tilde{S}-{\it\epsilon}$

$\tilde{k}=|\tilde{\boldsymbol{k}}|$

${\it\epsilon}$

(2.9) $$\begin{eqnarray}{\it\rho}_{i}\mathscr{L}\frac{\text{d}\tilde{{\it\eta}}^{\ast }}{\text{d}t}=-\frac{\mathscr{F}_{I}}{\tilde{D}_{0}}+{\it\epsilon}\left\{\frac{\mathscr{F}_{I}}{\tilde{D}_{0}^{2}}\tilde{d}^{\ast }-\frac{\tilde{k}^{2}}{2}[(2\tilde{D}_{0}+\tilde{d}^{\ast }){\it\kappa}_{w}\tilde{T}_{s}-\tilde{S}_{0}{\it\kappa}_{i}\tilde{T}_{b}]\right\}+O({\it\epsilon}^{2}),\end{eqnarray}$$

where $\mathscr{F}_{I}={\it\kappa}_{w}\tilde{T}_{s}+{\it\kappa}_{i}\tilde{T}_{b}$ . At this point, it is instructive to perform an order-of-magnitude analysis of (2.9). The flow conditions considered herein imply that $\tilde{D}_{0}$ varies in the range $[10{-}50]~{\rm\mu}\text{m}$ (see also table 1), while the morphological patterns we are considering have centimetre length scale, namely $\tilde{k}^{2}\sim 10^{-7}~{\rm\mu}\text{m}^{-2}$ . Since by definition $\tilde{d}=O(1)$ and $\tilde{S}_{0}$ has the same order of magnitude as $\tilde{D}_{0}$ , it is evident that the second term in the curly brackets is much smaller than the first one. The back-transform of the resulting approximate equation falls again into the form of (2.3), provided the following hold:

(2.10a,b ) $$\begin{eqnarray}A=-\frac{\mathscr{F}_{I}}{{\it\rho}_{i}\mathscr{L}\tilde{D}_{0}},\quad B=\frac{\mathscr{F}_{I}}{{\it\rho}_{i}\mathscr{L}\tilde{D}_{0}^{2}}.\end{eqnarray}$$

We note that, depending on the signs of the terms $\mathscr{F}_{K}$ and $\mathscr{F}_{I}$ , both karst and ice environments can be modelled as either depositional or dissolutive (remember that $\tilde{T}_{0_{i}}$ is in fact negative by construction). However, considering the formative conditions of the morphological patterns that superimpose over the net evolution of the liquid–solid substrate, it reasonable that karst flutings occur in the same conditions that have generated the supporting speleothem, i.e. a depositional environment. On the contrary, the presence of a melt flow during ice fluting morphogenesis suggests that, in the absence of other water sources, the icefall is melting. Based on such considerations and without any loss of generality, we will assume henceforth that $\mathscr{F}_{K}$ and $\mathscr{F}_{I}$ are both positive quantities.

Table 1. Typical values of dimensional and dimensionless quantities.

At this point, all lengths and velocities can be made dimensionless by using the Nusselt solutions (2.1a,b ), respectively, and the dimensionless time follows accordingly as $t=\tilde{t}\tilde{U} _{0}/\tilde{D_{0}}$ . Quantities without a tilde are therefore intended to be dimensionless. The evolution equation (2.3) recast in general dimensionless form is

(2.11) $$\begin{eqnarray}{\it\gamma}\frac{\text{d}{\it\eta}}{\text{d}t}=\pm 1+{\it\epsilon}d(\tilde{x},{\tilde{y}},t),\end{eqnarray}$$

where the positive (negative) sign refers to the karst (ice) case and ${\it\gamma}=\tilde{U} _{0}/B\tilde{D}_{0}$ is the ratio of the morphological to hydrodynamic time scales. In this unified approach, ${\it\gamma}$ is the only element that discriminates between the karst and ice problem, henceforth referred to with subscripts $K$ and $I$ , respectively. Assuming $\tilde{T}_{s}=0.1\,^{\circ }\text{C}$ , $\tilde{T}_{b}=-0.01\,^{\circ }\text{C}$ and standard conditions for speleothem formation, i.e. $\text{d}\tilde{{\it\eta}}_{K}/\text{d}t=\tilde{c}\sim O(\text{cm}~\text{century}^{-1})$ (e.g. Short et al. Reference Short, Baygents and Goldstein2005b ), one obtains

(2.12a,b ) $$\begin{eqnarray}{\it\gamma}_{I}=10^{4}\,\mathit{Re},\quad {\it\gamma}_{K}=\frac{\tilde{U} _{0}}{\tilde{c}},\end{eqnarray}$$

which furnish ${\it\gamma}_{I}$ of order 10² and ${\it\gamma}_{K}$ of order 10⁷ or larger (see table 1 for a summary of typical values of dimensional and dimensionless quantities). Without any loss of generality, the validity of (2.11) can also be extended to the curved region, provided the function $d$ is set to depend on $({\it\theta},y,t)$ .

The high values of the ratio between the morphological and hydrodynamic time scales, ${\it\gamma}$ , allows one to adopt the so-called quasi-steady approximation, where the flow field $\boldsymbol{u}$ is assumed to be slaved to the domain deformation induced by changes in either ${\it\eta}$ or $D$ . This assumption is customary in morphodynamic problems (Camporeale & Ridolfi Reference Camporeale and Ridolfi2012b ) and allows us to neglect the time derivative in momentum equations. The main modelling effect of this approximation is to disregard shear waves (Tollmien–Schlichting waves), which are always stable for the conditions herein considered ( ${\it\theta}>0$ and $\mathit{Re}\,\ll 1$ ).

Bearing in mind the above remarks, the governing equations are the dimensionless form of Navier–Stokes equations for incompressible fluids (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012)

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{{\rm\nabla}}p=\frac{1}{\mathit{Re}}{\rm\nabla}^{2}\boldsymbol{u}+\boldsymbol{F}, & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

$\boldsymbol{F}=(g\tilde{D}_{0}/\tilde{U} _{0}^{2})\{\sin {\it\theta},0,\cos {\it\theta}\}$

$\mathscr{I}(x,y,t)\equiv z-{\it\eta}(x,y,t)=0$

$\mathscr{S}(x,y,t)\equiv z-{\it\eta}(x,y,t)-D(x,y,t)=0$

(2.15) $$\begin{eqnarray}\frac{\text{D}\mathscr{S}}{\text{D}t}=0,\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\frac{\text{D}\mathscr{I}}{\text{D}t}=\left\{\begin{array}{@{}ll@{}}0\quad & (\text{dissolution{-}precipitation}),\\ \displaystyle {\it\rho}_{iw}\frac{\partial {\it\eta}}{\partial t}\quad & (\text{melting{-}freezing}),\end{array}\right.\end{eqnarray}$$

(2.17a−c ) $$\begin{eqnarray}(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}\boldsymbol{\cdot }\boldsymbol{t})_{z={\it\eta}+D}=0,\quad (\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}\boldsymbol{\cdot }\boldsymbol{n})_{z=D}+\mathit{Ka}\left(\frac{2}{\mathit{Re}^{5}\,\cos ^{2}{\it\theta}}\right)^{1/3}\mathscr{K}_{\mathscr{ S}}=0,\quad (\boldsymbol{t}\boldsymbol{\cdot }\boldsymbol{u})_{z={\it\eta}}=0,\end{eqnarray}$$

where $\text{D}/\text{D}t=\partial /\partial t+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}$ is the material derivative, $p$ is pressure, ${\it\rho}_{iw}$ is the density ratio of ice to water, $\boldsymbol{n}$ and $\boldsymbol{t}$ are the unit normal and tangent vector to a generic surface, respectively, $\unicode[STIX]{x1D64F}=p\unicode[STIX]{x1D644}-2\,\mathit{Re}^{-1}\,\unicode[STIX]{x1D63F}$ is the dimensionless Newtonian stress tensor (with $\unicode[STIX]{x1D644}$ and $\unicode[STIX]{x1D63F}$ the identity matrix and the rate-of-strain tensor, respectively) and $\mathscr{K}_{\mathscr{S}}=\{[(\unicode[STIX]{x1D644}-\boldsymbol{n}_{S}\boldsymbol{n}_{S})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}]\boldsymbol{\cdot }\boldsymbol{n}_{S}\}/2$ is the local mean curvature of the free surface (while $\mathscr{K}_{\mathscr{I}}$ is the corresponding value at the interface). The parameter $\mathit{Ka}=l_{c}^{2}g^{2/3}{\it\nu}^{-4/3}$ is the so-called Kapitza number (Chang & Demekhin Reference Chang and Demekhin2002), representing the ratio between surface tension and viscous forces, which depends on fluid properties ( $l_{c}\sim 2{-}3~\text{mm}$ is the water capillary length). Henceforth, if not further specified, we will assume the melting conditions for the ice and standard conditions, $T=20\,^{\circ }\text{C}$ , for karsts. It is also worth noting that thermocapillary effects can be neglected, since the Marangoni number $\mathit{Ma}=\mathit{Ka}{\it\delta}_{{\it\sigma}}$ is very small, where ${\it\delta}_{{\it\sigma}}\sim 1.8\times 10^{-4}$ is the temperature-induced (relative) variation in the surface tension (see table 1 and also refer to Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012, par. 9.5.3)). Since calcium concentration does not affect surface tension, solutal capillarity is also negligible.

Equation (2.15) is the kinematic condition for the free surface (responsible for surface waves), while equations (2.16) are the kinematic condition for the liquid–solid interface, where the source term on the right-hand side that appears in the melting–freezing case accounts for the meltwater flux entering the film (Hutter Reference Hutter1983, p. 136). Equations (2.17a,b ) state the tangential and normal components of momentum conservation at the free surface, i.e. the so-called dynamic conditions. We emphasize that the last term in (2.17b ) takes into account the effect of surface tension. Equation (2.17c ) states the no-slip condition. In summary, the complete dimensionless mathematical problem is provided by (2.11) and (2.13)–(2.17).

3. Parallel stability theory

As a first step, we neglect the curvature of the wall and perform a stability analysis in the planar region of the wall depicted in figure 2. As already done in the thermal problem, Prandtl’s mapping $z={\it\eta}+{\it\zeta}(1+d)$ is also adopted in the fluid dynamics problem, and a normal mode perturbation is assumed by following the decomposition of film thickness already introduced in the previous section, with the ansatz

(3.1) $$\begin{eqnarray}\{{\it\eta},D\}=\left\{\pm \frac{t}{{\it\gamma}},1\right\}+{\it\epsilon}\hat{d}\left\{\frac{1}{{\it\gamma}},1\right\}\text{e}^{{\it\lambda}t+\text{i}({\it\alpha}x+{\it\beta}y)}+\text{c.c.},\end{eqnarray}$$

where $\text{i}=\sqrt{-1}$ , ${\it\lambda}$ is the complex growth factor and the hat refers to the perturbation mode. As will be evident in the sequel, the time dependence of the base state of ${\it\eta}$ has no effects on the equations governing the perturbations, since it is linear in time and can be virtually eliminated by a rigid wall-normal translation of the reference frame ( ${\it\zeta}\rightarrow {\it\zeta}\pm t/{\it\gamma}$ ).

3.1. The morphodynamic Orr–Sommerfeld model

From standard arguments (Chang & Demekhin Reference Chang and Demekhin2002), it is well known that the linearized form of the problem (2.13)–(2.17) is equivalent to the Orr–Sommerfeld problem extended to open boundary domains. Accordingly, let us set $k^{2}=|\boldsymbol{k}|^{2}={\it\alpha}^{2}+{\it\beta}^{2}$ , where ${\it\alpha}$ and ${\it\beta}$ are the streamwise and spanwise components of the wavevector, respectively, and define a suitable streamfunction ${\it\phi}={\it\phi}_{0}({\it\zeta})+\hat{{\it\phi}}({\it\zeta})\exp [{\it\lambda}t+\text{i}({\it\alpha}x+{\it\beta}y)]$ . The constraint of the continuity equation allows the velocity perturbation modes $\{\hat{u} ,\hat{v},{\hat{w}}\}$ to follow the relations (Camporeale & Ridolfi Reference Camporeale and Ridolfi2012b )

(3.2a,b ) $$\begin{eqnarray}{\it\alpha}\hat{u} +{\it\beta}\hat{v}=k\frac{\text{d}\hat{{\it\phi}}}{\text{d}{\it\zeta}}+{\it\alpha}(1+\hat{d})u_{0}({\it\zeta}),\quad {\hat{w}}=-\text{i}k\hat{{\it\phi}},\end{eqnarray}$$

where $u_{0}=\text{d}{\it\phi}_{0}/\text{d}{\it\zeta}=2{\it\zeta}-{\it\zeta}^{2}$ is the base-state Nusselt velocity profile. The final linearized set of equations reads

(3.3) $$\begin{eqnarray}\left(\frac{\text{d}^{4}}{\text{d}{\it\zeta}^{4}}-2k^{2}\frac{\text{d}^{2}}{\text{d}{\it\zeta}^{2}}+k^{4}\right)\hat{{\it\phi}}=\text{i}{\it\alpha}\,\mathit{Re}\left[u_{0}\left(\frac{\text{d}^{2}}{\text{d}{\it\zeta}^{2}}-k^{2}\right)+2\right]\hat{{\it\phi}},\end{eqnarray}$$

(3.4a,b ) $$\begin{eqnarray}k\hat{{\it\phi}}-\text{i}{\it\chi}{\it\gamma}^{-1}\hat{d}=0,\quad k\frac{\text{d}\hat{{\it\phi}}}{\text{d}{\it\zeta}}+2{\it\alpha}=0\quad ({\it\zeta}=0),\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \left(k\frac{\text{d}^{2}}{\text{d}{\it\zeta}^{2}}+k^{3}\right)\hat{{\it\phi}}-2{\it\alpha}(1+\hat{d})=0\quad ({\it\zeta}=1), & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}{\it\zeta}}\left(\frac{\text{d}^{2}}{\text{d}{\it\zeta}^{2}}-3k^{2}-\text{i}{\it\alpha}\,\mathit{Re}\right)\hat{{\it\phi}}-{\it\Omega}(1+\hat{d})={\it\lambda}\,\mathit{Re}\frac{\text{d}\hat{{\it\phi}}}{\text{d}{\it\zeta}}\quad ({\it\zeta}=1), & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & (\text{i}{\it\lambda}-{\it\alpha})(1+\hat{d})=k\hat{{\it\phi}}\quad ({\it\zeta}=1), & \displaystyle\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle & {\it\gamma}{\it\lambda}=\hat{d}, & \displaystyle\end{eqnarray}$$

where ${\it\Omega}=\text{i}k[\mathit{Ka}(2/\mathit{Re}^{2}\,\cos {\it\vartheta})^{1/3}k^{2}-2\tan {\it\vartheta}]$ and ${\it\chi}=1-{\it\rho}_{i}/{\it\rho}_{w}$ in the I case and ${\it\chi}=1$ in the K case. Note that the quantity ${\it\Omega}$ is the sum of a term proportional to the Kapitza number, accounting for capillarity effects, and one accounting for gravitational effects, which are competing mechanisms for the stability problem, as will be shown in the Discussion. Additionally, in order to obtain (3.4), we have substituted (3.8) in the kinematic equation at the wall. In the following, the set of equations (3.3)–(3.8) will be referred as the ‘morphodynamic Orr–Sommerfeld’ (MOS) model.

3.2. Numerical solution

The MOS model (3.3)–(3.8) is solved numerically through a spectral Galerkin technique by adapting the procedure developed in Camporeale, Canuto & Ridolfi (Reference Camporeale, Canuto and Ridolfi2012) and Camporeale & Ridolfi (Reference Camporeale and Ridolfi2012b ), to which the reader is referred for additional details. Although analytical results will be provided in the following, a consistent numerical solution is in fact decisive in order to distinguish the exact boundaries for pattern formation parametrically.

Essentially, the differential system can be recast into a generalized algebraic eigenvalue problem of the form $\mathscr{S}\{{\it\phi}_{i}\}={\it\lambda}_{i}\mathscr{M}$ with $\mathscr{M}$ and $\mathscr{S}$ as the mass and stiffness matrices, respectively. With this in mind, a modal representation of the solution is adopted, where the function $\hat{{\it\phi}}$ is expanded in the (truncated) spectral form

(3.9) $$\begin{eqnarray}\hat{{\it\phi}}=\mathop{\sum }_{i=-3}^{N}{\it\phi}_{i}{\it\Phi}_{i}({\it\xi}),\end{eqnarray}$$

with ${\it\Phi}_{i}$ being a set of trial functions, whereas ${\it\phi}_{i}$ and $\hat{d}$ represent the new unknowns. For numerical convenience, the vertical coordinate is further mapped into the range ${\it\xi}\in [-1,1]$ , according to ${\it\xi}=2{\it\zeta}-1$ . The original ordinary differential equation problem is multiplied by a set of (smooth) test functions with zero first derivative at ${\it\xi}=\pm 1$ and zero value at ${\it\xi}=-1$ and integrated over the domain $[-1,1]$ . The repeated use of integration by parts allows the fourth derivative in (3.3) to be reduced to a second-order derivative; in addition, the dynamic boundary condition (3.6) is incorporated in the mass and stiffness operators in the so-called weak form, through the boundary term that arises from the integration by parts. The correct choice of the trial and test functions enables the use of the remaining boundary conditions to be set in a strong form, that is, as additional rows in the ultimate algebraic system.

We consider the set of polynomials

(3.10) $$\begin{eqnarray}{\it\Phi}_{i}=\sqrt{\text{i}+\frac{3}{2}}\left(\frac{L_{i+3}-L_{i+1}}{(2\text{i}+3)(2\text{i}+5)}-\frac{L_{i+1}-L_{i-1}}{(2\text{i}+1)(2\text{i}+3)}\right),\quad i\in [1,N-3],\end{eqnarray}$$

where $L_{i}({\it\xi})$ denotes the $i$ th Legendre polynomial. These functions are obtained by integrating twice each Legendre polynomial, while enforcing zero boundary conditions at ${\it\xi}=\pm 1$ for the function and its first derivative. The set of functions in (3.10) completes the set of solutions of the classical Orr–Sommerfeld problem, where the boundary conditions are all homogeneous (e.g. Shen Reference Shen1994). In order to accommodate non-vanishing boundary conditions, we need to add further low-degree polynomials, which are

(3.11a,b ) $$\begin{eqnarray}{\it\Phi}_{-3}={\textstyle \frac{1}{4}}({\it\xi}+1)^{2}({\it\xi}-1),\quad {\it\Phi}_{-2}={\textstyle \frac{1}{4}}(3-{\it\xi}^{2}-2{\it\xi}),\end{eqnarray}$$

(3.12a,b ) $$\begin{eqnarray}{\it\Phi}_{-1}={\textstyle \frac{1}{2}}(1-{\it\xi}^{2}),\quad {\it\Phi}_{0}={\textstyle \frac{1}{4}}(3+2{\it\xi}-{\it\xi}^{2}).\end{eqnarray}$$

It is straightforward to verify that trial and test spaces for the present problem are arranged in the following way, respectively,

(3.13a,b ) $$\begin{eqnarray}{\it\Phi}_{i}=\{{\it\Phi}_{-3},{\it\Phi}_{-2},{\it\Phi}_{-1},{\it\Phi}_{0},{\it\Phi}_{i}\}^{\text{T}},\quad \quad V=\{{\it\Phi}_{0},{\it\Phi}_{i}\}^{\text{T}}.\end{eqnarray}$$

The trial functions ${\it\Phi}_{-3}$ and ${\it\Phi}_{-2}$ become unnecessary in the particular cases ${\it\gamma}\rightarrow \infty$ and ${\it\alpha}=0$ , respectively, being the homogeneous boundary conditions at the wall.

The above specifications allow one to obtain the mass and stiffness matrices (not reported for the sake of the space). The adoption of a Galerkin discretization with the aforementioned Legendre modal basis and the use of integration by parts reduces the numerical complexity due to the fourth- and third-order derivatives of the Orr–Sommerfeld–Squire problem. This trick eliminates pointwise values of the higher derivative in the weak form of the equations, thus increasing the accuracy of the results, and avoids the occurrence of spurious eigenvalues. Unlike the classical Orr–Sommerfeld problem, the quasi-steady approximation yields a void mass matrix, except for three lines corresponding to non-homogeneous equations (3.6)–(3.8). The eigenvalue problem is therefore reduced to the solution of a $3\times 3$ linear system.

3.3. The Stokes model: analytical solution

In the context of the Stokes regime ( $\mathit{Re}\,\ll 1$ ), the right-hand sides of (3.3) and (3.6) can be neglected. This choice is numerically justified by the fact that the range of unstable wavenumbers and the considered Reynolds numbers yield the term ${\it\alpha}\,\mathit{Re}$ to be less than 10⁻³. Furthermore, from (3.8), it follows that ${\it\lambda}\,\mathit{Re}=O({\it\gamma}^{-1}\,\mathit{Re})$ , which is ${<}10^{-3}$ for the ice case and ${<}10^{-8}$ for the karst case. We also recall that dropping the right-hand side of (3.3) provides a biharmonic equation that is exact on the ${\it\beta}$ -axis, while the asymptotic case ${\it\gamma}\rightarrow \infty$ provides the classical Orr–Sommerfeld problem with no wall evolution. In summary, an analytical approach encompassing two fundamental approximations will be considered: neglecting convective terms in the momentum equation and the right-hand side of (3.6).

The thus derived differential Stokes problem allows a straightforward but cumbersome explicit solution for $\hat{{\it\phi}}(1)$ , as a function of the unknown $\hat{d}$ and parameters ${\it\alpha}$ , $k$ , $\mathit{Re}$ and ${\it\vartheta}$ . By substituting the expression for $\hat{{\it\phi}}(1)$ in the linearized kinematic condition (3.7), and combining with the wall evolution equation (3.8), we obtain a quadratic characteristic polynomial for ${\it\lambda}$ , with two analytical solutions:

(3.14) $$\begin{eqnarray}{\it\lambda}_{\pm }=-{\it\sigma}_{1}\pm ({\it\sigma}_{1}^{2}-2{\it\sigma}_{2})^{1/2},\end{eqnarray}$$

where

(3.15) $$\begin{eqnarray}\displaystyle {\it\sigma}_{1} & = & \displaystyle \frac{\text{i}{\it\alpha}}{2}+\frac{1}{2{\it\gamma}}-\frac{k\cosh k+\sinh k\left(2\text{i}k^{2}-R{\it\alpha}\right)}{k{\it\gamma}\left(2k^{2}+\cosh (2k)+\text{i}{\it\alpha}\mathit{Re}+1\right)}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{2k^{3}{\it\Omega}-k^{2}{\it\Omega}\sinh (2k)+\text{i}{\it\alpha}\left(4k^{4}+\text{i}\left(2k^{2}+1\right){\it\alpha}\mathit{Re}\right)-\text{i}{\it\alpha}^{2}\mathit{Re}\cosh (2k)}{4k^{4}\left(2k^{2}+\cosh (2k)+\text{i}{\it\alpha}\mathit{Re}+1\right)},\qquad\end{eqnarray}$$

(3.16) $$\begin{eqnarray}\displaystyle {\it\sigma}_{2} & = & \displaystyle \frac{\text{i}{\it\alpha}}{2{\it\gamma}}-\frac{2k^{3}{\it\Omega}-2k^{2}\cosh k\left(4k^{2}{\it\alpha}+{\it\Omega}\sinh k+2\text{i}{\it\alpha}^{2}\mathit{Re}\right)}{4\text{i}{\it\gamma}k^{4}\left(2k^{2}+\cosh (2k)+\text{i}{\it\alpha}\mathit{Re}+1\right)}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{{\it\alpha}\left(4k^{4}+\text{i}{\it\alpha}\mathit{Re}+2\text{i}k^{2}\right)-\text{i}{\it\alpha}^{2}\mathit{Re}(\cosh (2k)-4k\sinh k)}{4\text{i}{\it\gamma}k^{4}\left(2k^{2}+\cosh (2k)+\text{i}R{\it\alpha}+1\right)}.\end{eqnarray}$$

The values of ${\it\lambda}_{\pm }^{(r)}\equiv \text{Re}[{\it\lambda}_{\pm }]$ provide two modes of the growth rate of the perturbations that are related to the wall and free surface dynamics. The feature of the Stokes solution will be compared to the numerical solution of the MOS problem in the following section.

The critical condition on the ${\it\beta}$ -axis is obtained by setting ${\it\alpha}={\it\sigma}_{2}=0$ , leading to

(3.17) $$\begin{eqnarray}{\it\beta}_{c}=(2\,\mathit{Re}\,\sec {\it\vartheta})^{1/3}\sqrt{\frac{\sin {\it\vartheta}}{\mathit{Ka}}}.\end{eqnarray}$$

The square root in the above formula shows that fluted perturbations are stable when ${\it\vartheta}<0$ and this condition does not depend on the parameter ${\it\gamma}$ . Thus, the unstable domain on the ${\it\beta}$ -axis is not affected by the wall evolution dynamics, which is a first clue of HLM, as will be explained in the sequel. A synthetic result is provided by writing (3.17) in terms of the critical wavelength $L_{c}=2{\rm\pi}\tilde{D}_{0}/{\it\beta}_{c}$ . Quite curiously, this quantity does not depend on $\mathit{Re}$ and can be written in the compact form as $L_{c}/l_{c}=2{\rm\pi}\sqrt{\csc {\it\vartheta}}$ .

3.4. Results

The above mathematical model involves four dimensionless numbers ( $\mathit{Re}$ , ${\it\alpha}$ , ${\it\beta}$ ,  ${\it\vartheta}$ ) plus the calcite depositional rate, $c$ , which is given in $\text{cm}~\text{century}^{-1}$ for the sake of clarity. An extensive exploration of the results in the parameter subspace $\{\mathit{Re},{\it\alpha},{\it\beta},|{\it\vartheta}|,c\}<\{1,1,1,{\rm\pi}/2,40\}$ demonstrates that the real part of the least stable eigenvalue (henceforth referred to as ${\it\lambda}^{(r)}$ ) is positive under the condition ${\it\vartheta}>0$ only. The theory therefore provides instability only in inverted conditions, in agreement with observations (e.g. figure 1). This result has been obtained through both the numerical and analytical solution and for the karst and ice cases.

However, some important differences have emerged between the analytical and numerical solutions, which are worth remarking upon. Firstly, the (approximate) analytical solution localizes the maximum growth rate on the ${\it\beta}$ -axis systematically, which means that Squire’s theorem is always violated at any $\mathit{Re}$ . Conversely, the numerical solution is able to distinguish between two patterns: fluting and ripple. In the former case, the maximum growth rate is on the ${\it\beta}$ -axis; in the latter case, it is on the ${\it\alpha}$ -axis. If we define the function $\mathit{Re}_{t}({\it\vartheta};{\it\gamma})$ as the slope-dependent value of the Reynolds number for the transition between these two modes, then fluting is dominant at $\mathit{Re}<\mathit{Re}_{t}$ , while ripples dominate at $\mathit{Re}>\mathit{Re}_{t}$ . The shift between these two modes is sharp, that is, there is no hybrid bulb-like patterns with ${\it\alpha}$ and ${\it\beta}$ both not null. As stated earlier, the analytical solution is only valid in the fluting regime, which is also the focus of the present work.

The above remarks are well summarized in figure 3(a) where, for the K case, we observe that the fluting regime extends to higher values of $\mathit{Re}$ with increase in the calcite deposition rate, $c$ , and the transition curve has a maximum at ${\it\vartheta}=[1{-}10]^{\circ }$ and $\mathit{Re}=O(10^{-2})$ . In the I case (figure 3 b), the transition occurs at $\mathit{Re}=O(1)$ . This remarkable difference between the ice and karst cases is confirmed by the substantial absence of ripple-like rings on ice stalactites of icefalls when the Reynolds number is low.

Figure 3. Boundaries between fluting ( $\mathit{Re}<\mathit{Re}_{t}$ ) and ripple ( $\mathit{Re}>\mathit{Re}_{t}$ ) regimes as a function of ${\it\vartheta}$ : (a) karst case at different calcite deposition rates; (b) ice case. Points A, B and C refer to some particular cases that are commented upon in the text.

It is also instructive to discuss a few particular cases. Figure 4(a) shows the numerically computed contour map of ${\it\lambda}^{(r)}$ in the ${\it\alpha}{-}{\it\beta}$ plane for the particular case $\text{A}_{K}$ , marked in figure 3(a) (K case, $\mathit{Re}=0.0115$ , ${\it\vartheta}=2^{\circ }$ , ${\it\gamma}=5.4\times 10^{7}$ ). The outer contour line refers to the critical condition, namely ${\it\lambda}^{(r)}=0$ (dark blue curve). According to (3.17) the theoretical value of ${\it\beta}_{c}$ is $9\times 10^{-4}$ , which is in perfect agreement with the intercept of the outer contour with the ${\it\beta}$ -axis. Since this case belongs to the fluting regime, the maximum growth rate occurs at ${\it\alpha}=0$ . There is in fact a very narrow internal boundary layer with high numerical values reported on the ${\it\beta}$ -axis. A close-up of this thin layer is reported in figure 4(b) by using the analytical solution (3.14)–(3.16). The numerical solver described in § 3.2 provides an almost identical plot (not reported), supporting the fact that, very close to the ${\it\beta}$ -axis ( ${\it\alpha}$ being of order 10⁻⁸ or less in this case), the Stokes approximation theory provides reliable results. A posteriori, it is thus evident that the occurrence of such a strange narrow boundary layer – which is crucial for the existence of the fluting regime – is not a numerical artifact, since it is confirmed by the Stokes approximation theory. This theory is in fact exact on the ${\it\beta}$ -axis (see § 3.3). Incidentally, we note that the width of this layer scales as  ${\it\gamma}^{-1}$ .

Figure 4. (a) Contour plot of ${\it\lambda}^{(r)}$ from blue to red, for the case $\text{A}_{K}$ ( $\text{maximum}=1.2\times 10^{-8}$ ). (b) Close-up of a thin stripe around the ${\it\beta}$ -axis (numerical and analytical computations provide identical plots).

Furthermore, far from the ${\it\beta}$ -axis, but on the ${\it\alpha}$ -axis, two other relative maxima are also located. In the case of the fluting (ripple) regime, these maxima are lower (higher) than the peaks on the ${\it\beta}$ -axis. Differently from the ${\it\beta}$ -axis, these additional structures are not correctly provided by the the Stokes approximation theory, regardless of the Reynolds number, and therefore the ripple regime is not detectable analytically. This latter finding is much more evident if the behaviours on the ${\it\alpha}$ - and ${\it\beta}$ -axes are directly compared against the growth rate (see figure 5). In the fluting regime (case $\text{A}_{K}$ , panel a), the maximum along the ${\it\beta}$ -axis (thick solid line) is larger than that along the ${\it\alpha}$ -axis, when this latter is computed both numerically (thin line) and analytically (dashed line). On the contrary, in the ripple regime (case $\text{B}_{K}$ , panel b), the thin solid line is higher than the thick one, but the dashed line still remains lower. Therefore, the analytical prediction of the curve ${\it\lambda}^{(r)}({\it\alpha})$ always underestimates the true behaviour and mirrors the behaviour on the ${\it\beta}$ -axis, i.e. the outer contour line is nearly circular. The failure of the analytical model on the ${\it\alpha}$ -axis is a consequence of the two approximations contained in the Stokes problem, which both contribute to reducing the value of ${\it\lambda}^{(r)}$ . The two effects can be separated by solving the problem while neglecting the convective terms but not the right-hand side of (3.6). The result of this latter choice is reported by the dot-dashed line in figure 5, which still remains remarkably different from the dashed line where both approximations have been made.

Figure 5. The growth rate plotted against the ${\it\alpha}$ -axis (thin lines) and ${\it\beta}$ -axis (thick line). Plots against the ${\it\alpha}$ -axis: numerical computation (solid thin line), analytical computation (dashed thin line), numerical computation without convective terms (dash-dotted thin line). (a) Case $\text{A}_{K}$ and $c=5~\text{cm}~\text{century}^{-1}$ . (b) Case $\text{B}_{K}$ and $c=5~\text{cm}~\text{century}^{-1}$ .

The above picture is completed by reporting the results for the ice case (figure 6). Figure 6(a) shows the numerical result concerning point $\text{A}_{I}$ of figure 3(b) ( ${\it\gamma}=10^{3}$ ), where one can observe a large boundary layer centred on the ${\it\beta}$ -axis (the width of which is again of order ${\sim}{\it\gamma}^{-1}$ ) and the absence of secondary peaks. Owing to the circularity of the outer contour line, once again the estimation given by the analytical formula (3.17) corresponding to ${\it\beta}_{c}\sim 0.005$ provides a good approximation for the critical longitudinal wavenumber, ${\it\alpha}_{c}$ . Figure 6(b) (point $\text{B}_{I}$ , ${\it\gamma}=10^{4}$ ) refers to a condition that is close to the fluting–ripple transition. In this case the circularity of the outer contour line is completely lost, as is evident by the distorted scales of the picture, but maxima are still on the ${\it\beta}$ -axis and the ${\it\gamma}^{-1}$ scaling of the internal boundary layer is still valid. Incidentally, we note that, if the behaviour along the ${\it\alpha}$ -axis is plotted (not shown), then a very weak relative maximum appears, but it is several orders of magnitude smaller than the dominant one on the ${\it\beta}$ -axis. In figure 6(c) (point $\text{C}_{I}$ , ${\it\gamma}=2\times 10^{4}$ ), the scenario is completely different, since it occurs in the ripple regime. The higher contour lines in fact appear centred on the ${\it\alpha}$ -axis at ${\it\alpha}\sim 0.03$ , corresponding to a dimensional wavelength of nearly 5 cm. We observe that the plot reported in figure 6(c) extends to larger values of ${\it\alpha}$ and ${\it\beta}$ than the previous cases, which means that the emerging patterns involve shorter wavelengths. In this case, the theoretical critical wavelength (3.17) does not predict correctly the border of the instability on the ${\it\alpha}$ -axis.

Figure 6. Contour plot from blue to red of ${\it\lambda}^{(r)}$ for the cases $\text{A}_{I}$ (a),  $\text{B}_{I}$ (b) and $\text{C}_{I}$  (c). The maximum value is also reported: (a)  $1.1\times 10^{-6}$ ; (b)  $3.2\times 10^{-6}$ ; (c)  $3.5\times 10^{-4}$ .

It is worth noting that the maximum growth rate for the I case is numerically higher than for the K case; this is quite reasonable, since the time scales of the wall evolution are different in the two cases, as accounted for by the parameter ${\it\gamma}$ , which has the same order of magnitude as ${\it\lambda}^{-1}$ (see (3.8)).

At this point, we sum up the main results of the above analysis. (i) At the Reynolds numbers under investigation, the mathematical problem formulated in § 3.1 predicts pattern formation at only inverted conditions ( ${\it\vartheta}>0$ ). (ii) Provided the latter condition is satisfied, the results depend on the slope angle ${\it\vartheta}$ , as testified by the shape of the curve $\mathit{Re}_{t}({\it\vartheta})$ . (iii) At $\mathit{Re}<\mathit{Re}_{t}$ , the selected patterns are longitudinally oriented, i.e. flutes. (iv) Since the Stokes approximation is exact on the ${\it\beta}$ -axis, the developed theoretical formulas can describe the fluting regime correctly. On the contrary, in the ripple regime ( $\mathit{Re}>\mathit{Re}_{t}$ ) the numerical solution developed in § 3.2 is necessary. (v) These remarks are general and valid for both karst and ice problems, the only difference being quantitative and concerning the value of $\mathit{Re}_{t}$ for the fluting–ripple transition and the width of the internal boundary region located around the ${\it\beta}$ -axis. (vi) There is no bulb-like pattern formation corresponding to maxima out of the ${\it\alpha}$ -axis or ${\it\beta}$ -axis.

The fundamental result is that the present model turns out to be a suitable tool for explaining fluting formation, provided $\mathit{Re}<\mathit{Re}_{t}$ . The selected transverse wavenumbers are of order 10⁻³, which corresponds to spanwise wavelengths of centimetres, whereas the dimensional time scale of the instability is $\tilde{{\it\tau}}\sim \tilde{D}_{0}/\tilde{U} _{0}{\it\lambda}_{+}^{(r)}({\it\beta}_{m})$ , which yields $\tilde{{\it\tau}}_{K}=14$  days and $\tilde{{\it\tau}}_{I}=27$  min (for points $\text{A}_{K}$ and $\text{B}_{I}$ , respectively). The typical (spatial and temporal) scales of pattern inception are therefore captured for both problems (see also Martin-Perez, Martin-Garcia & Alonso-Zarza Reference Martin-Perez, Martin-Garcia and Alonso-Zarza2012; Chen Reference Chen2014). Additionally, non-negligible (albeit non-dominant) values of ${\it\lambda}^{(r)}$ at ${\it\alpha}\neq 0$ are a clue to the longitudinal component of the instability, and the well-known presence of draperies in the K case and their absence in the I case support these results. We infer that this feature is dependent on the value of parameter  ${\it\gamma}$ .

The dots in figure 7(a) show the selected transverse wavelength $\tilde{L}_{y}$ , as numerically computed, versus wall slope. Note that, because of the $\mathit{Re}$ dependence of $\tilde{D}$ , the dimensional outcome of $\tilde{L}_{y}$ is unaffected by Reynolds number (which varies over two orders of magnitude in the figure). This feature is shared with speleothem dynamics (Camporeale & Ridolfi Reference Camporeale and Ridolfi2012a ) and could justify the pattern regularity and the lack of noise in the shapes observed, since any seasonal or inter-annual variation in the flow does not force significant changes on the selected length scale. However, we will see in § 4 that a weak dependence on the Reynolds number arises when wall curvature is accounted for. Furthermore, the results are also ${\it\gamma}$ -independent (both cases are in fact reported in figure 7) which is a first symptom of HLM, a feature that will be definitively proved in the next section.

Figure 7. (a) Dimensional transverse wavelength versus ${\it\vartheta}$ ( $\mathit{Re}=[10^{-3}{-}10^{-1}]$ , $c=[1{-}5]~\text{cm}~\text{century}^{-1}$ ): MOS model (dots, either ice and karst case); fixed-bed model ( $\tilde{L}_{y}^{(H)}$ , continuous line); and non-parallel theory ( $\tilde{L}_{y}^{(G)}$ ), for $\mathit{Re}=0.03$ (full symbols) or 0.3 (open symbols) with $\tilde{\mathscr{R}}=2~\text{m}$ (squares) or 10 m (triangles), and $T=273~\text{K}$ (left row) or $T=288~\text{K}$ (right row). (b) Normalized growth factor versus ${\it\beta}$ for the case $\text{A}_{K}$ , where arrows indicate the value given by (3.26).

Although the focus of the present work is on the fluting regime, some remarks on ripple formation are due. In the karst problem, the emergence of ripples, also called crenulations, confirms the findings by Camporeale & Ridolfi (Reference Camporeale and Ridolfi2012a ), where a novel model successfully predicted the formation of ripple-like undulations on the surface of stalactites, by considering hydrodynamics, calcium and carbon dioxide geochemistry in detail. Despite the fact that the present model does not consider at all the dynamics of solutes (apart from the phenomenological evolution equation (2.11)), its ability to predict crenulations is confirmation of the recent remarks by Vesipa, Camporeale & Ridolfi (Reference Vesipa, Camporeale and Ridolfi2015), where a thorough sensitivity analysis has shown that a correct modelling of free surface dynamics is sufficient to provide successful predictions.

As far as the ice problem is concerned, the ripple regime can be assimilated to the so-called icicles. Past theoretical works have explained the formation of icicles by pure water falling films (e.g. Ueno et al. Reference Ueno, Farzaneh, Yamaguchi and Tsuji2010), but the recent experimental findings by Chen & Morris (Reference Chen and Morris2013) have questioned whether pure water is able to trigger icicle formation. The experiments in fact required the presence of dissolved salts in order to see appreciable instabilities. Likewise, the reliability of the ripple regime that is herein theoretically predicted could be questioned too. Yet, it is worthwhile to compute the value of the Reynolds number corresponding to the experiments conducted by Chen & Morris (Reference Chen and Morris2013). They reported a water flow rate $q=3~\text{g}~\text{min}^{-1}$ or less. Assuming a centimetre stalactite diameter, ${\it\phi}$ , this is equivalent to $\mathit{Re}=q/2{\rm\pi}{\it\phi}{\it\nu}=O(10^{-2})$ . When reported in figure 3(b), this value falls in the fluting regime, which could explain why the experiments did not produce icicles. If confirmed by further experiments made at higher values of Reynolds number, this argument might resolve the debate opened by the experimental results by Chen & Morris (Reference Chen and Morris2013). Nevertheless, it remains to be explained why dissolved salt is a destabilizing factor.

3.5. The fixed-bed model and the HLM conjecture

We now demonstrate that the geometrical structure of the fluting instability is completely driven by film hydrodynamics of inverted flows. To this end, we neglect the contribution of the wall evolution equation (2.3) and consider standard fixed-bed conditions (i.e. no slip and impermeability), in place of (3.4). This is formally equivalent to setting ${\it\gamma}^{-1}\rightarrow 0$ . Since fluting formation is addressed here, henceforth we will focus only on spanwise wavelength selection ( ${\it\alpha}=0$ ). The model (3.3)–(3.8) therefore reduces to

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\text{d}^{4}}{\text{d}{\it\zeta}^{4}}-2{\it\beta}^{2}\frac{\text{d}^{2}}{\text{d}{\it\zeta}^{2}}+{\it\beta}^{4}\right)\hat{{\it\phi}}=0\quad (0\leqslant {\it\zeta}\leqslant 1), & \displaystyle\end{eqnarray}$$

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{{\it\phi}}=\frac{\text{d}\hat{{\it\phi}}}{\text{d}{\it\zeta}}=0\quad ({\it\zeta}=0), & \displaystyle\end{eqnarray}$$

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\text{d}^{2}}{\text{d}{\it\zeta}^{2}}+{\it\beta}^{2}\right)\hat{{\it\phi}}=0\quad ({\it\zeta}=1), & \displaystyle\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}{\it\zeta}}\left(\frac{\text{d}^{2}}{\text{d}{\it\zeta}^{2}}-3{\it\beta}^{2}\right)\hat{{\it\phi}}-{\it\Omega}=0\quad ({\it\zeta}=1), & \displaystyle\end{eqnarray}$$

(3.22) $$\begin{eqnarray}\displaystyle & \text{i}{\it\lambda}-{\it\beta}\hat{{\it\phi}}=0\quad ({\it\zeta}=1). & \displaystyle\end{eqnarray}$$

The above mathematical problem refers to a purely hydrodynamic fixed-bed modelling, which has fully lost the morphological mode of the MOS and Stokes models, as well as the ability to catch the temporal dynamics of pattern formation. Notwithstanding this, our aim is to show that the geometrical features of the selected patterns are preserved through the action of the hydrodynamic mode, which is the essence of the HLM.

The solution of (3.18)–(3.21) is

(3.23) $$\begin{eqnarray}\displaystyle \hat{{\it\phi}}={\it\Omega}\frac{2{\it\beta}^{2}{\it\zeta}(\text{e}^{2{\it\beta}}-\text{e}^{2{\it\beta}{\it\zeta}})+{\it\beta}\text{e}^{2{\it\beta}}[({\it\zeta}-1)\text{e}^{2{\it\beta}{\it\zeta}}+({\it\zeta}+1)(\text{e}^{2{\it\beta}{\it\zeta}}+1)](\text{e}^{2{\it\beta}}+1)(\text{e}^{2{\it\beta}{\it\zeta}}-1)+{\it\beta}{\it\zeta}-1}{4{\it\beta}^{3}\text{e}^{{\it\beta}({\it\zeta}+1)}(2{\it\beta}^{2}+\cosh (2{\it\beta})+1)}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The use of (3.22) yields the hydrodynamic growth factor (hereafter denoted by  $H$ )

(3.24) $$\begin{eqnarray}{\it\lambda}_{r}^{(H)}=\frac{({\it\beta}-\sinh {\it\beta}\cosh {\it\beta})[2^{1/3}\mathit{Ka}\,\mathit{Re}^{-2/3}\,{\it\beta}^{2}\sec {\it\vartheta}^{1/3}-2\tan {\it\vartheta}]}{{\it\beta}[1+2{\it\beta}^{2}+\cosh (2{\it\beta})]}.\end{eqnarray}$$

In figure 7(b), the result provided by (3.24) is compared to the growth rates obtained by the MOS problem for the K case, under the conditions of case $\text{A}_{K}$ . Since the focus here is on the dependence of the growth factor on ${\it\beta}$ , the ordinate has been rescaled by the maximum value of each curve, ${\it\lambda}_{m}$ . It is remarkable that, despite some small differences in the shape, the values of the largest unstable wavenumber, ${\it\beta}_{c}$ , and the wavenumber that maximizes the instability, ${\it\beta}_{m}$ , coincide among the three curves. The same behaviour can be obtained if comparisons are made with the ice case. It follows that the hydrodynamic theory evaluates the same wavelengths selected as the morphodynamic theories, regardless of the value of  ${\it\gamma}$ .

Accordingly, we take advantage of (3.24) to analytically derive formulas for: (i) the critical wavenumber ${\it\beta}_{c}$ , by setting ${\it\lambda}_{r}^{(H)}({\it\beta}_{c})=0$ , and (ii) the most unstable wavenumber ${\it\beta}_{m}$ , by setting $\text{d}{\it\lambda}_{r}^{(H)}/\text{d}{\it\beta}=0$ . After Taylor-expanding in ${\it\beta}$ (owing to the reasonable assumption ${\it\beta}\ll 1$ ), one obtains

(3.25) $$\begin{eqnarray}\frac{\partial {\it\lambda}_{r}^{(H)}}{\partial {\it\beta}}=\frac{4}{3}{\it\beta}\left[-{\it\beta}^{2}\left(\frac{2^{1/3}\mathit{Ka}\sec {\it\vartheta}^{1/3}}{\mathit{Re}^{2/3}}+\frac{18}{5}\tan {\it\vartheta}\right)+\tan {\it\vartheta}\right]+O({\it\beta}^{5}).\end{eqnarray}$$

From (3.24) and (3.25), we arrive at

(3.26a,b ) $$\begin{eqnarray}{\it\beta}_{c}^{(H)}=\frac{(2\,\mathit{Re}\,\sec {\it\vartheta})^{1/3}}{\sqrt{\mathit{Ka}}}\sqrt{\sin {\it\vartheta}},\quad {\it\beta}_{m}^{(H)}\sim \left(2^{4/3}\mathit{Ka}\sin {\it\vartheta}(\,\mathit{Re}\,\sec {\it\vartheta})^{2/3}+\frac{18}{5}\right)^{-1/2}.\end{eqnarray}$$

Recalling that $\mathit{Re}\,\ll 1$ , (3.26a ) gives

(3.27) $$\begin{eqnarray}{\it\beta}_{m}=\frac{(\mathit{Re}\,\sec {\it\vartheta})^{1/3}}{2^{1/6}\sqrt{\mathit{Ka }\csc {\it\vartheta}}}\quad \Longrightarrow \quad \tilde{L}_{y}^{(H)}=2{\rm\pi}l_{c}\sqrt{2\csc {\it\vartheta}}.\end{eqnarray}$$

For the particular case ${\it\vartheta}={\rm\pi}/2$ , this latter result corresponds to the Rayleigh–Taylor instability, namely $\tilde{L}_{y}=2{\rm\pi}l_{c}\sqrt{2}\sim 8.8l_{c}$ (Charru Reference Charru2001). The behaviour of $\tilde{L}_{y}^{(H)}$ is reported in figure 7(a) with a continuous line. Remarkably, despite its simplicity, equation (3.27) completely agrees with the MOS theory (dots). In order to better appreciate this comparison, the relative difference between $\tilde{L}_{y}^{(H)}$ and the wavelength selected by the MOS theory has been computed for $\mathit{Re}=[10^{-3}{-}10^{-1}]$ and $c=[1{-}5]~\text{cm}~\text{century}^{-1}$ , obtaining a relative error always smaller than 10⁻⁴. This latter test furnishes the ultimate proof of the HLM conjecture.

4. Non-parallel stability theory

According to the local stability analysis developed so far, $\tilde{L}_{y}$ diverges for vertical flows (i.e. when ${\it\vartheta}$ approaches zero), so no wavelength is selected close to the critical condition and the formulas reported in the previous section are only valid for ${\it\vartheta}>0$ . Another criticism is that the flow-induced morphological channelization is triggered at the upper boundary of the unstable zone, namely where the azimuthal coordinate ${\it\theta}$ is zero. These shortcomings can be definitively remedied through a refined non-parallel stability analysis, which also encloses the radius of curvature $\tilde{\mathscr{R}}$ of the wall and non-parallelism of the flow field, in a neighbourhood of the vertical tangent point (figure 2). In order to simplify the analysis, we will axiomatically assume that, in the limit of small curvature, the HLM is also valid for the non-parallel case. Therefore, the analysis will be focused on the fixed-bed case, in the fluting regime, with no distinction between the karst and ice cases.

If all quantities are scaled with film depth and surface velocity at ${\it\theta}=0$ , a new parameter emerges, namely the dimensionless curvature ${\it\mu}=(\tilde{D}/\tilde{\mathscr{R}})_{{\it\theta}=0}$ . Since ${\it\mu}\ll 1$ , we are allowed to simplify several relationships that will be derived in the following. It is convenient to adopt cylindrical coordinates $\{{\it\theta},r,y\}$ , as depicted in figure 2. Accordingly, the velocity vector becomes $\boldsymbol{u}=\{u_{{\it\theta}},u_{r},u_{y}\}$ and the mathematical problem is recast in the form

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{r}\frac{\partial p}{\partial {\it\theta}}-\frac{1}{\mathit{Re}}\left({\rm\nabla}^{2}u_{{\it\theta}}+\frac{2}{r^{2}}\frac{\partial u_{r}}{\partial {\it\theta}}-\frac{u_{{\it\theta}}}{r^{2}}\right)-2\frac{\sin {\it\theta}}{\mathit{Re}}=0, & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial p}{\partial r}-\frac{1}{\mathit{Re}}\left({\rm\nabla}^{2}u_{r}-\frac{2}{r^{2}}\frac{\partial u_{{\it\theta}}}{\partial {\it\theta}}-\frac{u_{r}}{r^{2}}\right)-2\frac{\cos {\it\theta}}{\mathit{Re}}=0, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial p}{\partial y}-\frac{1}{\mathit{Re}}\left(\frac{\partial ^{2}u_{y}}{\partial {\it\theta}^{2}}+\frac{\partial ^{2}u_{y}}{\partial y^{2}}+\frac{\partial ^{2}u_{y}}{\partial r^{2}}\right)=0, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial r}(ru_{r})+r\frac{\partial u_{y}}{\partial y}+\frac{\partial u_{{\it\theta}}}{\partial {\it\theta}}=0, & \displaystyle\end{eqnarray}$$

${\rm\nabla}^{2}\equiv r^{-2}\partial _{{\it\theta}{\it\theta}}+r^{-1}\partial _{r}(r\partial _{r})+\partial _{yy}$

(4.5) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}=0\quad (r=\mathscr{R}), & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle D_{t}+\frac{u_{{\it\theta}}}{r}D_{{\it\theta}}+u_{y}D_{y}-u_{r}=0\quad (r=\mathscr{R}+D), & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}\boldsymbol{\cdot }[\boldsymbol{t}_{x}\,\boldsymbol{t}_{y}]=0\quad (r=\mathscr{R}+D), & \displaystyle\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle p-2\,\mathit{Re}^{-1}(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{n})+2^{1/3}\mathit{Ka}\,\mathit{Re}^{-5/3}\mathscr{K}_{S}=0\quad (r=\mathscr{R}+D), & \displaystyle\end{eqnarray}$$

$\boldsymbol{n}$

$\boldsymbol{t}_{x}$

$\{r^{-1}D_{{\it\theta}},-1,D_{y}\}$

$\{1,r^{-1}D_{{\it\theta}},0\}$

$\boldsymbol{t}_{y}=\boldsymbol{t}_{x}\times \boldsymbol{n}$

$D$

(4.9) $$\begin{eqnarray}\unicode[STIX]{x1D63F}=\left(\begin{array}{@{}ccc@{}}\displaystyle \frac{1}{r}\frac{\partial u_{{\it\theta}}}{\partial {\it\theta}}+\frac{u_{r}}{r} & \displaystyle \frac{1}{2}\left(-\frac{\partial u_{{\it\theta}}}{\partial y}-\frac{1}{r}\frac{\partial u_{y}}{\partial {\it\theta}}\right) & \displaystyle \frac{1}{2}\left(\frac{u_{{\it\theta}}}{r}-\frac{\partial u_{{\it\theta}}}{\partial r}-\frac{1}{r}\frac{\partial u_{r}}{\partial {\it\theta}}\right)\\ \displaystyle -\frac{\partial u_{{\it\theta}}}{\partial y}-\frac{1}{r}\frac{\partial u_{y}}{\partial {\it\theta}} & \displaystyle \frac{\partial u_{y}}{\partial y} & \displaystyle \frac{1}{2}\left(-\frac{\partial u_{y}}{\partial z}-\frac{\partial u_{r}}{\partial y}\right)\\ \displaystyle \frac{1}{2}\left(\frac{u_{{\it\theta}}}{r}-\frac{\partial u_{{\it\theta}}}{\partial z}-\frac{1}{r}\frac{\partial u_{r}}{\partial {\it\theta}}\right) & \displaystyle \frac{1}{2}\left(-\frac{\partial u_{y}}{\partial z}-\frac{\partial u_{r}}{\partial y}\right) & \displaystyle \frac{\partial u_{r}}{\partial r}\end{array}\right).\end{eqnarray}$$

Prandtl’s mapping is posed in the form $r={\it\zeta}d+{\it\mu}^{-1}$ , thus the solution ansatz of a generic quantity $f$ (e.g. depth, velocity and pressure) is recast in a global-mode fashion,

(4.10) $$\begin{eqnarray}f({\it\theta},{\it\zeta},y,t;{\it\mu})=f_{0}({\it\theta},{\it\zeta};{\it\mu})+{\it\epsilon}\hat{f}({\it\theta},{\it\zeta};{\it\mu})\exp [{\it\lambda}^{(G)}t+\text{i}{\it\beta}y],\end{eqnarray}$$

where the explicit dependence on ${\it\theta}$ accounts for non-parallel effects, whereas the basic state $O({\it\epsilon}^{0})$ can be further expanded in powers of ${\it\mu}$ up to the first order, namely $f_{0}=f_{00}({\it\theta},{\it\zeta})+{\it\mu}f_{01}({\it\theta},{\it\zeta})$ . In the present non-parallel theory, (4.10) is a generalization of (3.1). Note that in a standard global stability analysis the azimuthal coordinate ${\it\theta}$ also appears in the exponential, through a multiple scaling approach, and the problem is then tackled through the WKBJ technique (Huerre & Rossi Reference Huerre, Rossi, Godreche and Manneville1998). Since we are only interested to the longitudinally parallel perturbations (flutes), this trick can be avoided. Moreover, a detailed order-of-magnitude analysis of continuity equation (4.4) suggests that, at the basic state, the radial component of velocity is at most of order ${\it\mu}$ , and its correct ansatz at the order $O({\it\epsilon}^{0})$ is therefore $u_{r0}={\it\mu}u_{r01}({\it\theta},{\it\zeta})+{\it\mu}^{2}u_{r02}({\it\theta},{\it\zeta})$ .

4.1. Problem $O({\it\epsilon}^{0})$

The equations at leading order pose ${\it\epsilon}={\it\mu}=u_{y00}=0$ and therefore read

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial ^{2}u_{{\it\theta}00}}{\partial {\it\zeta}^{2}}+2(1+D_{00})^{2}\cos {\it\theta}=0, & \displaystyle\end{eqnarray}$$

(4.12) $$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Re}\frac{\partial p_{00}}{\partial {\it\zeta}}-2(1+D_{00})\sin {\it\theta}=0, & \displaystyle\end{eqnarray}$$

(4.13) $$\begin{eqnarray}\displaystyle & \displaystyle (1+D_{00})\frac{\partial u_{{\it\theta}00}}{\partial {\it\theta}}-{\it\zeta}D_{00}^{\prime }\frac{\partial u_{{\it\theta}00}}{\partial {\it\zeta}}+\frac{\partial u_{r01}}{\partial {\it\zeta}}=0, & \displaystyle\end{eqnarray}$$

(4.14) $$\begin{eqnarray}\displaystyle & u_{{\it\theta}00}=u_{r01}=0\quad ({\it\zeta}=0), & \displaystyle\end{eqnarray}$$

(4.15) $$\begin{eqnarray}\displaystyle & D_{00}^{\prime }u_{{\it\theta}00}-u_{r01}=0\quad ({\it\zeta}=1), & \displaystyle\end{eqnarray}$$

(4.16) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{{\it\theta}00}}{\partial {\it\zeta}}=p_{00}=0\quad ({\it\zeta}=1), & \displaystyle\end{eqnarray}$$

where the prime refers to the ${\it\theta}$ derivative. The above set of equations is readily solved, yielding

(4.17a−d ) $$\begin{eqnarray}D_{00}=(\sec {\it\theta})^{1/3},\quad u_{{\it\theta}00}=\frac{u_{0}}{D_{00}},\quad u_{r01}=\frac{u_{0}{\it\zeta}}{3}\tan {\it\theta},\quad p_{00}=-\frac{u_{0}^{\prime }D_{00}}{\mathit{Re}}\sin {\it\theta},\end{eqnarray}$$

where $u_{0}=2{\it\zeta}-{\it\zeta}^{2}$ is the Nusselt parabolic profile, which was already derived in the parallel problem.

After substituting (4.17) and collecting the terms of order $O({\it\epsilon}^{0}{\it\mu}^{1})$ , one obtains the problem

(4.18) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial ^{2}u_{{\it\theta}01}}{\partial {\it\zeta}^{2}}+4(\cos {\it\theta})^{2/3}D_{01}+\frac{2}{3}[5-6{\it\zeta}+(\sec {\it\theta})^{2}]=0, & \displaystyle\end{eqnarray}$$

(4.19) $$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Re}\frac{\partial p_{01}}{\partial {\it\zeta}}-2D_{01}\sin {\it\theta}+2(3{\it\zeta}-2)\frac{\sin {\it\theta}}{3\cos {\it\theta}^{2/3}}=0, & \displaystyle\end{eqnarray}$$

(4.20) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{{\it\zeta}}{3}[D_{01}(9{\it\zeta}-14){\it\Upsilon}_{2}+6D_{01}^{\prime }(\cos {\it\theta})^{1/3}+{\it\Upsilon}_{1}]+\frac{\partial u_{r02}}{\partial {\it\zeta}}=0, & \displaystyle\end{eqnarray}$$

(4.21) $$\begin{eqnarray}\displaystyle & u_{{\it\theta}01}=u_{r02}=0\quad ({\it\zeta}=0), & \displaystyle\end{eqnarray}$$

(4.22) $$\begin{eqnarray}\displaystyle & \displaystyle D_{01}^{\prime }-\frac{2\tan {\it\theta}}{7}D_{01}+{\it\Upsilon}_{3}=0\quad ({\it\zeta}=1), & \displaystyle\end{eqnarray}$$

(4.23a,b ) $$\begin{eqnarray}\frac{\partial u_{{\it\theta}01}}{\partial {\it\zeta}}-(\cos {\it\theta})^{1/3}=1,\quad \mathit{Re}\,p_{01}-\frac{2}{3}{\it\Upsilon}_{2}-\mathit{Ka}\left(\frac{4}{\mathit{Re}}\right)^{2/3}=0\quad ({\it\zeta}=1),\end{eqnarray}$$

where

(4.24a ) $$\begin{eqnarray}\displaystyle {\it\Upsilon}_{1}=\frac{[2(5-2{\it\zeta})\sec {\it\theta}^{2}-(18{\it\zeta}^{2}-28{\it\zeta}+7)]\sin {\it\theta}}{3(\cos {\it\theta})^{4/3}}, & & \displaystyle\end{eqnarray}$$

(4.24b,c ) $$\begin{eqnarray}\displaystyle {\it\Upsilon}_{2}=\frac{\sin {\it\theta}}{(\cos {\it\theta})^{2/3}},\quad {\it\Upsilon}_{3}=\frac{\cos (3{\it\theta})+9\sin {\it\theta}}{36(\cos {\it\theta})^{11/3}}. & & \displaystyle\end{eqnarray}$$

The complete solution of (4.18)–(4.23) is

(4.25) $$\begin{eqnarray}\displaystyle & \displaystyle D_{01}=\frac{7[33(\cos {\it\theta})^{8/21}-8\sec ^{2}{\it\theta}-25]}{600(\cos {\it\theta})^{2/3}}, & \displaystyle\end{eqnarray}$$

(4.26) $$\begin{eqnarray}\displaystyle & \displaystyle u_{{\it\theta}01}=-\frac{1}{3}{\it\zeta}[18D_{01}({\it\zeta}-2)(\cos {\it\theta})^{2/3}+({\it\zeta}-2)(\sec {\it\theta})^{2}-{\it\zeta}(2{\it\zeta}-5)-7], & \displaystyle\end{eqnarray}$$

(4.27) $$\begin{eqnarray}\displaystyle & \displaystyle p_{01}=\frac{10D_{01}}{\mathit{Re}}({\it\zeta}-1)\sin {\it\theta}+\frac{12\mathit{Ka}(\cos {\it\theta})^{2/3}+2^{2/3}((4-3{\it\zeta}){\it\zeta}+1)\,\mathit{Re}^{2/3}\sin {\it\theta}}{3\,\mathit{Re}(\mathit{Re}\sin 2{\it\theta}\csc {\it\theta})^{2/3}}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(4.28) $$\begin{eqnarray}\displaystyle u_{r02} & = & \displaystyle \frac{{\it\zeta}^{2}[72({\it\zeta}-3)\cos ^{11/3}({\it\theta}){D_{01}}^{\prime }]}{108(\cos {\it\theta})^{10/3}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{{\it\zeta}^{2}\sin {\it\theta}[(36(7-3{\it\zeta})D_{01}\cos ^{8/3}{\it\theta}+({\it\zeta}(27{\it\zeta}-56)+21)\cos 2{\it\theta}+{\it\zeta}(27{\it\zeta}-40)-39)]}{108(\cos {\it\theta})^{10/3}}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The sum $D_{00}+{\it\mu}D_{01}$ provides an approximate solution of the depth, up to order ${\it\mu}$ , for the non-parallel base state problem and, to the author’s knowledge, is a novel result for the present geometry. It should be noted that the first term corresponds to the uniform Nusselt’s solution computed at a varying slope ${\it\theta}(s)\neq {\it\vartheta}$ . The flow is actually well approximated by this solution (also called normal solution) when ${\it\mu}$ is small, namely when ${\it\theta}$ changes slowly in the streamwise direction due to the small curvature of the wall. The normal solution has a minimum $D_{00}=1$ at the vertical tangent point, ${\it\theta}=0$ , and it results from a balance between the streamwise component of the gravity force, the pressure gradient and the basal shear stress. At the order ${\it\mu}$ , normal viscous stresses acting on the vertical planes play a role in smoothing the free surface curvature, and they slightly reduce its deviation from the value obtained at the vertical tangent point. This means that $D_{01}$ always gives a negative contribution. By using the relations in (4.17) and (4.25) it is straightforward to compute the relative deviation between the normal $D_{00}$ solution and its correction ${\it\mu}D_{01}$ , which is at most ${\sim}0.15{\it\mu}$ for $|{\it\theta}|<{\rm\pi}/4$ . Since in the present analysis ${\it\mu}$ will be considered of order 10⁻³ or less, the correction to the normal solution is not relevant, despite its theoretical value and applicability to other problems. One could therefore neglect the effects of this term in the following development of the perturbed problem. However, we prefer to keep it in the following derivation, in order to pose the analysis in a more general way.

4.2. Problem $O({\it\epsilon})$ : the global eigenvalue

The problem at order $O({\it\epsilon})$ is much more involved, but it can be addressed quite easily by recognizing that the term related to the mass forces is wholly absorbed at the leading order and the Stokes flow problem, written at order ${\it\epsilon}$ , reduces to the search for a solution of the equation ${\rm\nabla}^{2}\hat{\boldsymbol{u}}=\mathit{Re}^{-1}\,\boldsymbol{{\rm\nabla}}\hat{p}$ posed in a cylindrical reference frame. This allows us to use the Papkovich–Neuber solutions (Tran-Cong & Blake Reference Tran-Cong and Blake1982) to the perturbed quantities,

(4.29a,b ) $$\begin{eqnarray}\hat{\boldsymbol{u}}=\boldsymbol{{\rm\nabla}}(\boldsymbol{r}\boldsymbol{\cdot }{\it\Psi}+{\it\Psi}_{0})-2{\it\Psi},\quad \hat{p}=2\,\mathit{Re}^{-1}\,\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\Psi},\end{eqnarray}$$

where $\boldsymbol{r}=r\{\cos {\it\theta},\sin {\it\theta}\}$ , and ${\it\Psi}=\{{\it\Psi}_{{\it\theta}},{\it\Psi}_{r}\}$ and ${\it\Psi}_{0}$ are harmonic potentials ( ${\it\Psi}_{y}$ is set to zero because of the transverse invariance), which satisfy the following equations:

(4.30) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial ^{2}{\it\Psi}_{{\it\theta}}}{\partial {\it\theta}^{2}}-{\it\Psi}_{{\it\theta}}+r\left[\frac{\partial }{\partial r}\left(r\frac{\partial {\it\Psi}_{{\it\theta}}}{\partial r}\right)+r\frac{\partial ^{2}{\it\Psi}_{{\it\theta}}}{\partial y^{2}}\right]+2\frac{\partial {\it\Psi}_{r}}{\partial {\it\theta}}=0, & \displaystyle\end{eqnarray}$$

(4.31) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial ^{2}{\it\Psi}_{r}}{\partial {\it\theta}^{2}}-{\it\Psi}_{r}+r\left[\frac{\partial }{\partial r}\left(r\frac{\partial {\it\Psi}_{r}}{\partial r}\right)+r\frac{\partial ^{2}{\it\Psi}_{r}}{\partial y^{2}}\right]-2\frac{\partial {\it\Psi}_{{\it\theta}}}{\partial {\it\theta}}=0, & \displaystyle\end{eqnarray}$$

(4.32) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial ^{2}{\it\Psi}_{0}}{\partial {\it\theta}^{2}}+r\left[\frac{\partial }{\partial r}\left(r\frac{\partial {\it\Psi}_{0}}{\partial r}\right)+r\frac{\partial ^{2}{\it\Psi}_{0}}{\partial y^{2}}\right]=0. & \displaystyle\end{eqnarray}$$

The above equations are solved by invoking once again the global-mode ansatz

(4.33) $$\begin{eqnarray}{\it\Psi}_{j}={\it\epsilon}\hat{{\it\Psi}}_{j}({\it\theta},r,y)\exp [{\it\lambda}^{(G)}t+\text{i}{\it\beta}y]\quad (j=0,{\it\theta},r).\end{eqnarray}$$

After substituting (4.33) in (4.30)–(4.32), using the Prandtl mapping ( $r:\rightarrow {\it\zeta}$ ) and expanding up to the first order in ${\it\epsilon}$ and ${\it\mu}$ , the three harmonic equations for the perturbation mode of potentials reduce to

(4.34) $$\begin{eqnarray}\frac{\partial ^{2}\hat{{\it\Psi}}_{j}}{\partial {\it\zeta}^{2}}+{\it\mu}D_{00}\frac{\partial \hat{{\it\Psi}}_{j}}{\partial {\it\zeta}}-{\it\beta}^{2}D_{00}^{2}\hat{{\it\Psi}}_{j}=0\quad (j={\it\theta},r,0).\end{eqnarray}$$

By recalling that $D_{00}=(\sec {\it\theta})^{1/3}$ and ${\it\theta}={\it\mu}s$ , the approximate solutions of (4.34) are

(4.35) $$\begin{eqnarray}\hat{{\it\Psi}}_{j}(s,{\it\zeta})=C_{j}^{\pm }(s)\exp [\pm ({\it\beta}-{\it\mu}/2){\it\zeta}+O({\it\mu}^{2})]\quad (j=0,{\it\theta},r).\end{eqnarray}$$

The coefficients $C_{j}^{\pm }$ are obtained by imposing no slip and impermeability at ${\it\zeta}=0$ and the dynamic conditions at ${\it\zeta}=1$ (i.e. the O( ${\it\epsilon}$ ) order version of (4.5) and (4.7), written in the new coordinate system), which read, respectively,

at ${\it\zeta}=0$

(4.36) $$\begin{eqnarray}\left[\begin{array}{@{}ccc@{}}{\it\mu} & \displaystyle -\frac{1}{2}\frac{\partial }{\partial s} & \displaystyle -{\it\mu}\frac{\partial }{\partial {\it\zeta}}\\ 0 & \displaystyle -{\it\mu}+\frac{\partial }{\partial {\it\zeta}} & 0\\ 0 & 1 & {\it\mu}\end{array}\right]\left(\begin{array}{@{}c@{}}\hat{{\it\Psi}}_{{\it\theta}}\\ \hat{{\it\Psi}}_{r}\\ \hat{{\it\Psi}}_{0}\end{array}\right)=0,\end{eqnarray}$$

at ${\it\zeta}=1$

(4.37) $$\begin{eqnarray}\displaystyle & & \displaystyle \left[\begin{array}{@{}ccc@{}}\displaystyle {\it\mu}-\frac{\partial }{\partial {\it\zeta}} & \displaystyle \left[{\it\mu}-1+\left(\frac{1}{{\it\mu}}-2{\it\mu}\right)\frac{\partial }{\partial {\it\zeta}}\right]\frac{\partial }{\partial s} & \displaystyle \left[\left(1-\frac{{\it\mu}}{2}\right)\frac{\partial }{\partial s}-{\it\mu}\right]\frac{\partial }{\partial {\it\zeta}}\\ 0 & \displaystyle (5{\it\mu}^{2}-{\it\mu}-1)\frac{\partial }{\partial {\it\zeta}}+\frac{{\it\mu}s^{2}}{6}\frac{\partial }{\partial {\it\zeta}} & \displaystyle {\it\mu}\frac{\partial }{\partial {\it\zeta}}\\ \displaystyle \frac{\partial }{\partial s} & \displaystyle {\it\mu}+\left(1+{\rm\mu}s\frac{\partial }{\partial s}-\frac{1+{\it\mu}}{{\it\mu}}\frac{\partial }{\partial {\it\zeta}}\right)\frac{\partial }{\partial {\it\zeta}} & \displaystyle -\frac{\partial ^{2}}{\partial {\it\zeta}^{2}}\end{array}\right]\left(\begin{array}{@{}c@{}}\hat{{\it\Psi}}_{{\it\theta}}\\ \hat{{\it\Psi}}_{r}\\ \hat{{\it\Psi}}_{0}\end{array}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\left(\begin{array}{@{}c@{}}{\it\mu}\\ 0\\ \displaystyle \frac{\mathit{Ka}{\it\beta}^{2}}{(2\,\mathit{Re})^{2/3}}\end{array}\right)\hat{D}.\end{eqnarray}$$

After substituting (4.35) in (4.36) and (4.37) and solving, it follows that $C_{j}^{\pm }=c_{j}^{\pm }\hat{D}(s)$ , with $j=\{{\it\theta},r,0\}$ . If a truncated series expansion in ${\it\mu}$ and ${\it\beta}$ is performed (recall in fact that ${\it\beta}\ll 1$ ), the coefficients $c_{j}^{\pm }$ read

(4.38a,b ) $$\begin{eqnarray}c_{{\it\theta}}^{\pm }=\pm \frac{{\it\mu}}{24{\it\beta}}\left[{\it\beta}^{2}\left(\frac{\sqrt[3]{2}\mathit{Ka}s}{\mathit{Re}^{2/3}}-3\right)+6\right],\quad c_{r}^{\pm }=\mp \frac{{\it\beta}\mathit{Ka}[{\it\beta}(2\mp 8{\it\mu})+{\it\mu}+2]}{2^{5/3}\,\mathit{Re}^{2/3}},\end{eqnarray}$$

(4.39) $$\begin{eqnarray}c_{0}^{\pm }=\pm \frac{{\it\beta}\mathit{Ka}\{{\it\beta}[\pm 8-{\it\mu}(109{\it\mu}+32)]\pm 3(4-29{\it\mu}){\it\mu}+8\}}{16~2^{2/3}{\it\mu}\,\mathit{Re}^{2/3}}.\end{eqnarray}$$

It is worth recalling that $\hat{D}(s)$ is the global mode of depth perturbation and it plays the role of the unknown eigenfunction of the problem. After plugging all this in the kinematic condition (4.6) and applying the aforementioned transformations, we finally get a first-order differential eigenvalue problem relating $\hat{D}(s)$ and ${\it\lambda}^{(G)}$ in the form

(4.40) $$\begin{eqnarray}\left[\frac{\partial }{\partial s}+\frac{\mathscr{F}(s)+{\it\lambda}^{(G)}}{1+{\it\varrho}s}\right]\hat{D}=0,\end{eqnarray}$$

where the coefficient ${\it\varrho}$ and the function $\mathscr{F}(s)$ depend in a complicated manner upon ${\it\mu}$ and ${\it\beta}$ . Since (4.40) is formally solvable by a separation of variables, its solution is proportional to the boundary condition $\hat{D}(0)$ , which can therefore be absorbed by the prefactor ${\it\epsilon}$ , by putting $\hat{D}(0)=1$ . It is also evident that the eigenvalue must satisfy the relation ${\it\lambda}^{(G)}=-\mathscr{F}(0)-\hat{D}^{\prime }(0)$ , whence a continuous spectrum of eigenvalues parametrized by the quantity $\hat{D}^{\prime }(0)$ emerges. A reasonable choice is to put $\hat{D}^{\prime }(0)=0$ since, in this manner, it can be shown that the solution (4.40) is the most regular one (i.e. the average value of $\hat{D}^{\prime }(s)$ is minimized). We therefore arrive at our main goal, namely the solution of the global eigenvalue

(4.41) $$\begin{eqnarray}{\it\lambda}^{(G)}=-\mathscr{F}(0)=\frac{{\it\beta}(63{\it\mu}+8)\cosh {\it\beta}-[{\it\beta}^{2}(77{\it\mu}+8)+4({\it\mu}+2)]\sinh {\it\beta}}{8(2\,\mathit{Re})^{2/3}}{\it\beta}\mathit{Ka}\text{e}^{-{\it\mu}/2}.\end{eqnarray}$$

Inspection of the dependence of (4.41) on ${\it\beta}$ , in the range ${\it\beta}\ll 1$ , shows that there is a maximum (the selected wavenumber) at ${\it\beta}_{m}^{(G)}\sim \sqrt{177{\it\mu}}/4\sqrt{2}$ and a zero value (the cutoff) at ${\it\beta}_{c}^{(G)}=\sqrt{2}{\it\beta}_{m}^{(G)}$ . This amazing result translates into a dimensional selected wavelength that reads

(4.42) $$\begin{eqnarray}\tilde{L}_{y}^{(G)}\sim 3{\it\nu}^{1/3}\left(\frac{\mathit{Re}}{g}\right)^{1/6}\tilde{\mathscr{R}}^{1/2}.\end{eqnarray}$$

The latter formula provides a finite wavelength selection at finite values of dome radius of curvature. A novel weak $\mathit{Re}$ dependence has appeared in (4.42), through a one-sixth power law. The effect of the capillary length, which is crucial in planar conditions, is therefore replaced by $\tilde{\mathscr{R}}$ when curvature is considered. Physically speaking, the effect of surface tension is overwhelmed by streamwise momentum diffusion, because of non-parallel effects. The symbols in figure 7(a) provide some computations of (4.42); for example, $\tilde{L}_{y}^{(G)}=3~\text{cm}$ for $\tilde{\mathscr{R}}=2~\text{m}$ and $\mathit{Re}=0.3$ .

A possible scenario is that the patterns are triggered at wavelength $\tilde{L}_{y}^{(G)}$ around $s=0$ , where ${\it\theta}=0$ and the curvature is maximum, then fluting spacing gradually converges downstream to the value provided by $\tilde{L}_{y}^{(H)}$ , where the wall becomes almost planar. The non-parallel theory presented herein has therefore resolved the inconsistency observed with the parallel theory: finite wavelengths can be selected at ${\it\theta}=0$ .

6. Concluding remarks

Elongated organ-pipe morphologies are widespread in ice and karst environments. The former case is very common over icefalls, the latter case bedecks the cave walls, wherein speleothems represent an important source of palaeo-climatic data. A novel approach relating fluid mechanics to the morphogenesis of such structures could therefore be of paramount importance. Through a unified modelling approach, we have solved the morphodynamic problem in terms of a stability analysis, and the typical features of karst and ice flutings have been shown to be hydrodynamically locked. The results obtained reject Squire’s theorem when the Reynolds number is below a critical value that depends on the slope. The correct identification of the transition between flutings and ripples has been possible thanks to the numerical solution of the complete Orr–Sommerfeld-like problem, through an advanced spectral Galerkin technique.

Although the occurrence of the HLM is not essential to the instability discussed here, it has allowed us to explain that the morphogenesis of such patterns is completely driven by hydrodynamics. Note that the HLM would be a trivial finding if (2.3) was replaced by an algebraic relationship. The result is instead counterintuitive because the presence of a phenomenological evolution equation is expected to introduce new dynamics into the system, as in the case of ice ripples and crenulations (Camporeale & Ridolfi Reference Camporeale and Ridolfi2012a ,Reference Camporeale and Ridolfi b ). Differently, in the present case, it appears that the wall dynamics of fluted patterns (which occur at inverted conditions) are fully slaved to hydrodynamics in a very wide range of the morphological parameter ${\it\gamma}$ , thus embedding completely different environmental conditions, such as icefall and caves.

Another advantage in demonstrating the HLM comes from adopting the framework of the hydrodynamic stability analysis for Stokes flow. Indeed, parallel and non-parallel stability analyses of inverted creeping flows have provided several analytical results about the selected wavelength, which match well with observations. The parallel theory readily provides the wavelengths at various slopes, but fails close to verticality. The non-parallel theory gives a local result only in a neighbourhood of verticality.

It is worth noting that in the parallel case there is no dependence of flute spacing on Reynolds number; instead it appears as a result of the non-parallel theory, namely at the head of flutes, through a one-sixth power law. Such a dependence of pattern wavelength on flow characteristics could be of palaeo-climatic interest. Nowadays, karst environments are in fact acquiring relevance in the field of palaeo-climate reconstructions. Several pieces of research have demonstrated that mineral deposits constituting speleothems can record past climatic conditions: measurable properties of the speleothems (e.g. mineral composition, calcite crystals) can in fact be related to dripping water and to the surrounding environmental conditions in which the deposition occurred (McDermott, Mattey & Hawkesworth Reference McDermott, Mattey and Hawkesworth2001). However, as pointed out by Fairchild et al. (Reference Fairchild, Smith, Baker, Fuller, Spotl, Mattey and McDermott2006), speleothem records cannot yet be given a unique climatological meaning.

The analytical theories provided in this work are limited to Stokes flows, namely when the Reynolds number is much smaller than unity, which is the norm for the phenomena herein investigated. Although an improvement to larger values of $\mathit{Re}$ could be desirable, it is instructive to show that formula (4.42) is quite robust nevertheless. A suggestive example is provided by figure 1(g): a close observation of the fountain’s features allows one to assume $\mathscr{R}\sim 0.2~\text{m}$ , while the mean annual temperature is $T=15\,^{\circ }\text{C}$ . The evaluation of a precise value of the Reynolds number is prohibitive, so it is reasonable to assume a conservative range of potential values such as $\mathit{Re}=[20{-}500]$ . With this setting, formula (4.42) returns $\tilde{L}_{y}^{(G)}=[5{-}8.5]~\text{cm}$ , which is fairly close to the actual value, despite the fact that the Stokes approximation would be no longer valid at such conditions.

More generally, other boundary value problems involving interface dynamics that follow (2.3) should support the same results (e.g. photosynthetic biofilms), as long as the linear theory is valid. Nonlinearities can contribute to modify the final picture, possibly introducing pattern coarse-graining. We have neglected front dynamics, since the related fingering instability develops with time scales too short to have a morphological significance in the present context, as is evident from numerical simulations (Lin, Kondic & Filippov Reference Lin, Kondic and Filippov2012).