Hostname: page-component-848d4c4894-wg55d Total loading time: 0 Render date: 2024-05-30T20:24:17.579Z Has data issue: false hasContentIssue false

Boundary-layer approach to the linear karren instability

Published online by Cambridge University Press:  20 December 2023

Simeon Djambov*
Affiliation:
Laboratory of Fluid Mechanics and Instabilities, EPFL, Lausanne CH-1015, Switzerland
François Gallaire
Affiliation:
Laboratory of Fluid Mechanics and Instabilities, EPFL, Lausanne CH-1015, Switzerland
*
Email address for correspondence: simeon.djambov@epfl.ch

Abstract

The present paper concerns the linear fate of transverse perturbations in a gravity-driven, thin-film flow over a soluble substrate. We propose a reduced-order model, based on a boundary-layer treatment of the solute transport and a depth-integration of the Stokes equations, using two extended lubrication methodologies found in the literature. We obtain a closed-form dispersion relation, which we compare to a previous, fully resolved analytical investigation (Bertagni and Camporeale, J. Fluid Mech., vol. 913, 2021, A34). The results allow us to distil the essential physical mechanisms behind the instability.

Type
JFM Rapids
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2023. Published by Cambridge University Press.

1. Introduction

Linear karren patterns are dissolution features adorning the bare, inclined faces of soluble rocks, e.g. limestone (figure 1a). They appear as regularly spaced grooves following the incline and their hydrodynamical origin is undisputed (Ford & Williams Reference Ford and Williams2007; Ginés et al. Reference Ginés, Knez, Slabe and Dreybrodt2009). In contrast to fracture-induced forms, these grooves are caused by flows of rainwater or snow melt over the exposed rock surfaces. Their streamwise length as well as their transverse periodicity – observed in nature to span over length scales from millimetres to metres – are questions of particular interest, as they may shed light on past climatic conditions.

Figure 1. (a) Natural karren in the Bernese Alps (top) and in the Jura (bottom). (b) Sketch of the problem at hand.

Recent laboratory experiments by Guérin et al. (Reference Guérin, Derr, Courrech du Pont and Berhanu2020) showed that a uniform film flow carves an array of parallel rills on its gypsum substrate. Although the film thicknesses they probed were larger than in natural conditions and the flows presented wall-induced turbulence, this evidence highlighted the theoretical lack in the study of karren formation.

Bertagni & Camporeale (Reference Bertagni and Camporeale2021) provided the first theoretical analysis in the framework of linear stability theory. They used a complete description of the laminar hydrodynamical flow and, crucially, addressed the dissolution process by fully resolving the solute transport within the flow. They demonstrated the system's linear instability to transverse perturbations over a substantial range of wavenumbers and revealed the importance of transverse diffusion, both of momentum and solute, for the stabilisation of larger wavenumbers.

Here, we propose a boundary-layer approach to this linear karren instability. We consider a creeping, gravity-driven, streamwise-invariant flow. We model the dissolution process, by the use of Lévêque's (Reference Lévêque1928) boundary-layer solution, which profoundly simplifies the solute transport problem and enables us, together with the solution to the linearised Stokes equations, established by Bertagni & Camporeale (Reference Bertagni and Camporeale2021), to write a dispersion relation in closed form. The latter is remarkably faithful to the complete, tridimensional treatment of the solute transport by Bertagni & Camporeale (Reference Bertagni and Camporeale2021).

Finally, we propose to reduce the hydrodynamical description using two extended lubrication models, found in the literature. The revised lubrication theory, recently proposed by Devauchelle, Popović & Lajeunesse (Reference Devauchelle, Popović and Lajeunesse2022), may lack ‘the mathematical rigour of an order expansion’, notably as it includes the transverse momentum diffusion, while neglecting other effects, a priori of the same asymptotical weight, but is a very intuitive depth-integrated model. The weighted-residual integral boundary-layer methodology (see Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Sheid and Velarde2012) and references therein), on the other hand, has the advantage of being asymptotically consistent, although somewhat tedious to manipulate. Both produce sets of depth-integrated streamwise and transverse momentum conservation equations, which further emphasise the simple underlying physics of the coupling in the problem at hand: the streamwise component of the flow interacts with the dissolution process, while the transverse flow serves to flatten the free surface. Their results are validated against the exact linearised Stokes equations.

2. Problem formulation and modelisation

We investigate the laminar, incompressible, gravity-driven flow of a thin film of Newtonian liquid (water) of density $\rho$, kinematic viscosity $\nu$ and surface tension $\sigma$ over a soluble substrate (e.g. limestone or gypsum), inclined by an angle $\theta$ with respect to the horizontal, and the solute transport inside the film, with diffusion coefficient $D$ (figure 1b). The substrate material's insolubility $\gamma =\rho _s/c_{sat}$ is defined as a dimensionless solid density ratio between its crystallised and dissolved states, where $\rho _s$ is the rock's density and $c_{sat}$ is the equilibrium, i.e. saturated, mass concentration in the liquid.

2.1. Governing equations

The system is governed by the continuity, the Stokes and the advection–diffusion equations. We make the problem dimensionless by rescaling all lengths by the characteristic film thickness $H_N$, the pressure p by a hydrostatic pressure scale $\rho g\sin \theta H_N/2$ and the velocity ${\boldsymbol{u}}$ by a Nusselt velocity scale $g \sin \theta H_N^2/(2\nu )$, where g is the gravitational acceleration. Lastly, we rescale the concentration c by its saturated value $c_{sat}$. Assuming that the concentration field adapts instantly to the flow, i.e. in a quasi-stationary regime, the governing equations read

(2.1)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}
(2.2)\begin{gather}\boldsymbol{\nabla}p=\nabla^2\boldsymbol{u}+\boldsymbol{f}, \end{gather}
(2.3)\begin{gather}{Pe} \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} c=\nabla^2 c, \end{gather}

where $\boldsymbol {f}=2\hat {\boldsymbol {x}}-2\cot \theta \hat {\boldsymbol {y}}$ is the gravity forcing. The Péclet number ${Pe}=g\sin \theta H_N^3/(2 \nu D)$ is the product of the Reynolds and Schmidt numbers ${Sc}=\nu /D=10^3$, and can therefore be large even for creeping flows.

As the considered water films are usually formed by rainwater runoff or snow melt, the water is assumed solute-free at the inlet:

(2.4a)\begin{equation} c|_{x=0}=0. \end{equation}

In general, the dissolution process is governed by a competition between the surface reaction and diffusion rates. For simple dissolution reactions, the comparison is relatively straightforward: the dissolution of very soluble minerals, e.g. salt, are overwhelmingly diffusion-controlled. In contrast, the dissolution of calcite involves a much more sophisticated chemical reaction, the kinetics of which can vary substantially with the conditions. It is usually accepted that calcite's dissolution in very undersaturated water is also transport-controlled (Ford & Williams Reference Ford and Williams2007). This means that the reaction at the solid–liquid interface is quasi-instantaneous; the water is nearly saturated at contact with the rock,

(2.4b)\begin{equation} c|_{y=\eta}=1, \end{equation}

and the dissolution rate is thus regulated by how fast diffusion can evacuate the solute. The solid mass conservation equation then governs the geomorphological evolution:

(2.4c)\begin{equation} \gamma {Pe} \partial_t \eta=\boldsymbol{\nabla}c\boldsymbol{\cdot} \boldsymbol{\nabla}[y-\eta]|_{y=\eta}. \end{equation}

We remind the reader that the dimensionless parameter $\gamma$ measures the insolubility of the substrate material; it is large for weakly soluble materials, e.g. $\gamma =49\ 091$ for calcite (Bertagni & Camporeale Reference Bertagni and Camporeale2021). Moreover, there is no solute flux through the liquid–air interface

(2.4d)\begin{equation} \boldsymbol{\nabla}c\boldsymbol{\cdot} \boldsymbol{\nabla}[y-\eta-h]|_{y=\eta+h}=0. \end{equation}

From (2.4c), we observe that the geomorphological time scale is much larger than the hydrodynamical one. The substrate's evolution – being controlled by diffusion, ${Pe}\gg 1$, and undergoing a dramatic reduction in its solid density, $\gamma \gg 1$ – is vastly slower than the flow velocities, which enables us to use the classical no-slip and no-penetration boundary conditions at the substrate:

(2.5a)\begin{equation} \boldsymbol{u}|_{y=\eta}=\boldsymbol{0}. \end{equation}

The kinematic boundary condition describes the evolution of the film's thickness:

(2.5b)\begin{equation} \partial_t h = \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}[y-\eta-h]|_{y=\eta+h}. \end{equation}

Finally, the liquid–air interface is free of tangential stresses and undergoing a capillary jump in the normal stress:

(2.5c)\begin{equation} {\boldsymbol{\mathsf{T}}}\boldsymbol{\cdot} \boldsymbol{\nabla}[y-\eta-h]|_{y=\eta+h} ={-}{We} {\mathcal{K}} \boldsymbol{\nabla}[y-\eta-h], \end{equation}

with the stress tensor ${\boldsymbol{\mathsf{T}}}=-p{\boldsymbol{\mathsf{I}}}+(\boldsymbol {\nabla } \boldsymbol {u}+\boldsymbol {\nabla } \boldsymbol {u}^T)$, the curvature

(2.6)\begin{equation} {\mathcal{K}}=\boldsymbol{\nabla}\boldsymbol{\cdot} \left[\frac{\boldsymbol{\nabla}[y-\eta-h]}{\|\boldsymbol{\nabla}[y-\eta-h]\|}\right], \end{equation}

and the Weber number ${We}=2\sigma /(\rho g\sin \theta H_N^2)$. The gradient expressions $\boldsymbol {\nabla }[y-\eta -h]$ and $\boldsymbol {\nabla }[y-\eta ]$ are the upward-pointing vectors, normal to the upper and lower boundaries of the liquid domain, respectively. They are not normalised in the kinematic and dissolution conditions ((2.5b), (2.4c)), as $\partial _t h$ and $\partial _t\eta$ are themselves projections of the normal interface velocities onto $\hat {\boldsymbol {y}}$.

2.2. Flat-rock basic flow

The stationary, transverse-invariant basic flow is that of a uniform liquid film falling over an inclined flat rock,

(2.7a,b)\begin{equation} h_b=1,\quad \eta_b=0, \end{equation}

involving a half-Poiseuille streamwise velocity profile and a hydrostatic pressure field

(2.7c,d)\begin{equation} \boldsymbol{u_b}=y(2-y)\hat{\boldsymbol{x}},\quad p_b=2\cot\theta(1-y). \end{equation}

Polyanin et al. (Reference Polyanin, Kutepov, Vyazmin and Kazenin2001, pp. 130–132) provide an exact solution to the advection–diffusion equation (2.3) for this basic flow, ignoring the streamwise diffusion $(\partial _x^2 c_b)$, which is negligible for large Péclet numbers ${Pe}$ (figure 2). Near the inlet, the solute forms a thin boundary layer, outside of which the water remains pure. The thickness of this layer grows due to diffusion and the water's aggressivity towards the substrate decreases. Downstream, the accumulation of solute saturates the film.

Figure 2. Two-dimensional basic concentration field $c_b$ and wall-normal profiles at $x/{Pe}=\{10^{-3}, 10^{-2}, 10^{-1}, 1\}$. The orange solid lines represent the exact solution, while the black dashed lines depict Lévêque's (Reference Lévêque1928) boundary-layer solution.

Since the basic concentration field and, crucially, its gradient at the substrate are a function of the streamwise coordinate $x$, flat rock cannot be an exact solution. Bertagni & Camporeale (Reference Bertagni and Camporeale2021) remedied this by isolating the dissolution rate induced by the uniform film flow $V$ from the evolution of the substrate $\partial _t\eta$. They formulated the dissolution law (2.4c) as

(2.8)\begin{equation} \gamma {Pe}(\partial_t \eta-V)=\boldsymbol{\nabla}c\boldsymbol{\cdot} \boldsymbol{\nabla}[y-\eta]|_{y=\eta},\quad \text{with }\gamma {Pe} V={-}\partial_y c_b|_{y=0}. \end{equation}

2.3. Solutal Lévêque boundary-layer dissolution model

Let us focus on the region near the inlet, where the water is almost pure of solute. The concentration interacts with the velocity field only very close to the wall, in a solutal boundary layer, where the velocity has an almost purely shear profile. Under these conditions, $x\ll {Pe}$, and keeping only the streamwise advection and the wall-normal diffusion, Lévêque (Reference Lévêque1928) reduced the advection–diffusion equation (2.3) to

(2.9)\begin{equation} {Pe} \tau (y-\eta)\partial_x c=\partial_y^2c, \end{equation}

where the advecting streamwise velocity is asymptotically expanded using its shear rate at the substrate $\tau =\partial _y u|_{y=\eta }$, and subject to boundary conditions

(2.10ac)\begin{equation} c|_{x=0}=0, \quad c|_{y=\eta}=1, \quad \partial_y c|_{y\to\infty}=0. \end{equation}

Lévêque (Reference Lévêque1928) obtained a self-similar solution, which reads

(2.11)\begin{equation} c(x,y)=\frac{\varGamma\left(\dfrac{1}{3},\dfrac{{Pe} \tau (y-\eta)^3}{9 x}\right)}{\varGamma(1/3)},\end{equation}

with $\varGamma (a,s)=\int _s^\infty t^{a-1}\,\mathrm {e}^{-t}\,\mathrm {d}t$ and $\varGamma (a)\equiv \varGamma (a, 0)$ – the incomplete and complete gamma functions, respectively. For the basic flow, $\tau _b=2$ and $\eta _b=0$, this solution is compared to the exact one in figure 2 and is found to be fairly valid up to $x/{Pe}\lesssim 0.1$, as already noted by Bertagni & Camporeale (Reference Bertagni and Camporeale2021).

By injecting this solution (2.11) into the dissolution law (2.8), we obtain the following evolution equation:

(2.12)\begin{equation} \partial_t\eta=V-(\gamma{Pe})^{{-}1}\|\boldsymbol{\nabla}[y -\eta]\|^2 \partial_y c|_{y=\eta}=V (1-\|\boldsymbol{\nabla}[y-\eta]\|^2(\tau/\tau_b)^{1/3}), \end{equation}

with $V=(6{Pe}/x)^{1/3}/(\gamma {Pe}\varGamma (1/3))$ and where we have used the saturation condition (2.4b) to rewrite $\boldsymbol {\nabla }c\boldsymbol {\cdot } \boldsymbol {\nabla }[y-\eta ]|_{y=\eta }=\|\boldsymbol {\nabla }[y-\eta ]\|^2\partial _y c|_{y=\eta }$. This model indeed ignores the in-plane concentration diffusion, which is important at small wavelengths, but should capture well the behaviour of the boundary layer's thickness with respect to the streamwise velocity.

2.4. Depth-integrated hydrodynamical models

At this stage, only the streamwise shear rate at the rock surface $\tau$ is necessary for closing the problem with the dissolution model (2.12). As mentioned above, Bertagni & Camporeale (Reference Bertagni and Camporeale2021) achieved an exact solution to the linearised Stokes equations for this problem, which can be directly applied here for the purposes of linear stability analysis, as we will see in § 3.1. However, with the purpose to validate fully depth-averaged models, we set out to generalise two extended lubrication methodologies and obtain three closed-form expressions for the dispersion relation. We believe that depth-averaged models could in future prove very useful when it comes to studying the nonlinear evolution of linear karren patterns. They can suitably be used as surrogates for the full Stokes equations with a free interface, which would require an interface-capturing or interface-tracking multiphase solver. A Wolfram Mathematica script, performing the mathematical derivations from this section, is available as supplementary material at https://doi.org/10.1017/jfm.2023.932.

2.4.1. Revised lubrication model

Firstly, following Devauchelle et al. (Reference Devauchelle, Popović and Lajeunesse2022), and generalising their model to free-surface flows, we make a half-Poiseuille ansatz for the in-plane velocity profiles,

(2.13)\begin{equation} u(y,z)\hat{\boldsymbol{x}}+w(y,z)\hat{\boldsymbol{z}} =\frac{3\boldsymbol{q}(z)}{h(z)}F_0\left[\frac{y-\eta(z)}{h(z)}\right], \end{equation}

with the local flow rate vector $\boldsymbol {q}=q_x\hat {\boldsymbol {x}}+q_z\hat {\boldsymbol {z}}$. The $y$-dependence of the in-plane velocity field is isolated into a semi-parabolic structure

(2.14)\begin{equation} F_0[\zeta ]=\zeta -\zeta ^2/2, \end{equation}

ensuring the no-slip condition (2.5a) on the wall and a leading-order, long-wave version of the stress-free dynamic boundary condition (2.5c) on the water–air interface ($\partial _y u|_{y=\eta +h}=\partial _y w|_{y=\eta +h}=0$). The mass conservation equation reads

(2.15)\begin{equation} \partial_t h+q_z'=0, \end{equation}

where the prime $'$ denotes a differentiation with respect to $z$.

The integration of the wall-normal component of the momentum equations (2.2) produces the hydrostatic pressure field

(2.16)\begin{equation} p(y,z)=2\cot\theta(\eta+h-y)+{We}{\mathcal{K}}, \end{equation}

which includes the capillary contribution of the curved interface and where, in the lubrication spirit, we have ignored the wall-normal momentum's diffusion $\nabla ^2 v$.

Finally, using the ansatz (2.13), we integrate the in-plane components of the Stokes equations (2.2) through the film and obtain the following momentum conservation equations:

(2.17) \begin{align} \left(1+\eta'^2-h'\eta'-h'^2+\frac{h(h''+\eta'')}{2}\right) \boldsymbol{q} & +h(h'+\eta')\boldsymbol{q}'-\frac{h^2\boldsymbol{q}''}{3} \nonumber\\ & =\dfrac{h^3}{3}(2\hat{\boldsymbol{x}}-2\cot \theta(h'+\eta')\hat{\boldsymbol{z}}-{We}{\mathcal{K}}'\hat{\boldsymbol{z}}). \end{align}

There are three sources of depth-integrated momentum: the streamwise projection of the gravitational acceleration and the transverse hydrostatic and capillary pressure gradients, generated by deformations of the free surface.

The revision, proposed by Devauchelle et al. (Reference Devauchelle, Popović and Lajeunesse2022), essentially consists in allowing $\boldsymbol {q}$ to keep an autonomous dependence on $z$, instead of it being adiabatically slaved to $h$ and $\eta$, as in the classical lubrication theory (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Sheid and Velarde2012), where

(2.18)\begin{equation} \boldsymbol{q}=\frac{h^3}{3}(2\hat{\boldsymbol{x}} -2\cot\theta(h'+\eta')\hat{\boldsymbol{z}}-{We}{\mathcal{K}}'\hat{\boldsymbol{z}}). \end{equation}

Evaluating the streamwise shear rate at the substrate from the ansatz (2.13) produces

(2.19)\begin{equation} \tau=\frac{3 q_x}{h^2}. \end{equation}

2.4.2. Weighted-residual integral boundary-layer model

Secondly, following the derivations gathered in the monograph of Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Sheid and Velarde2012), we aim at a weighted-residuals formulation for the free-surface thin-film flow over a variable substrate morphology $\eta (z)$. We rewrite the Stokes equations (2.2) as a set of streamwise-invariant boundary-layer equations, which are asymptotically consistent up to second order (so as to keep the in-plane momentum diffusion) in a thin-film order parameter, separating in-plane and wall-normal length scales,

(2.20a)\begin{gather} 2+\partial_y^2 u+\partial_z^2u=0, \end{gather}
(2.20b)\begin{gather}-2\cot\theta-\partial_y p+\partial_y^2 v=0, \end{gather}
(2.20c)\begin{gather}-\partial_z p+\partial_y^2 w+\partial_z^2 w=0, \end{gather}

with the dynamic boundary conditions now reading

(2.21a)\begin{gather} \partial_y u|_{y=h+\eta}=\partial_z[h+\eta] \partial_z u|_{y=h+\eta}, \end{gather}
(2.21b)\begin{gather}\partial_y w|_{y=h+\eta}=4\partial_z[h+\eta] \partial_z w|_{y=h+\eta}-\partial_z v|_{y=h+\eta}, \end{gather}
(2.21c)\begin{gather}p|_{y=h+\eta}=2\partial_y v|_{y=h+\eta}-{We}\partial_z^2[h+\eta]. \end{gather}

In order to account for higher-order deviations from the half-Poiseuille profile, the in-plane velocity ansatz is corrected with two additional polynomial functions, weighted by their associated coefficients $\boldsymbol {r}=r_x\hat {\boldsymbol {x}}+r_z\hat {\boldsymbol {z}}$ and $\boldsymbol {s}=s_x\hat {\boldsymbol {x}}+s_z\hat {\boldsymbol {z}}$:

(2.22)\begin{equation} u(y,z)\hat{\boldsymbol{x}}+w(y,z)\hat{\boldsymbol{z}}= \frac{3(\boldsymbol{q}-\boldsymbol{r}-\boldsymbol{s})}{h}F_0 \left[\frac{y-\eta}{h}\right]+\frac{45 \boldsymbol{r}}{h}F_1\left[\frac{y-\eta}{h}\right]+\frac{210 \boldsymbol{s}}{h}F_2\left[\frac{y-\eta}{h}\right], \end{equation}

where

(2.23a,b) \begin{align} F_1[\zeta]&=\zeta-\frac{17\zeta^2}{6}+\frac{7\zeta^3}{3}-\frac{7\zeta^4}{12},\nonumber\\ F_2[\zeta]&=\zeta-\frac{13\zeta^2}{2}+\frac{57\zeta^3}{4}- \frac{111\zeta^4}{8}+\frac{99\zeta^5}{16}-\frac{33\zeta^6}{32}, \end{align}

which are chosen orthogonal to one another, and to $F_0[\zeta ]$, with respect to the scalar product $\langle a|b\rangle =\int _\eta ^{h+\eta } a b\,\mathrm {d} y$, and satisfy the leading-order boundary conditions. The local flow rate vector remains $\boldsymbol {q}$ and hence the mass conservation equation (2.15) still holds. The corrections $\boldsymbol {r}$ and $\boldsymbol {s}$ are at least of order 1 in the thin-film parameter, so that their second derivatives with respect to $z$ fall out of the boundary-layer equations (2.20). Since their first derivatives do not appear in the equations either, $\boldsymbol {r}$ and $\boldsymbol {s}$ are adiabatically slaved to $\boldsymbol {q}$, $h$ and $\eta$.

The pressure field is integrated from the wall-normal boundary-layer equation (2.20b):

(2.24)\begin{equation} p=2\cot\theta(h+\eta-y)-\partial_z w -\partial_z w|_{y=h+\eta}-{We}(h''+\eta''), \end{equation}

where we have used the continuity equation (2.1) to replace $\partial _y v$.

Finally, we perform a Galerkin projection of the in-plane boundary-layer momentum equations (2.17) onto the basis functions themselves,

(2.25a)\begin{gather} \left\langle \left.F_i\left[\frac{y-\eta}{h}\right]\right|2 +\partial_z^2u\right\rangle+F_i[1]\partial_y u|_{y=h+\eta} +\left\langle \partial_y^2F_i \left. \left[\dfrac{y-\eta}{h}\right]\right|u\right\rangle=0, \end{gather}
(2.25b)\begin{gather}\left\langle \left.F_i\left[\frac{y-\eta}{h}\right]\right| -\partial_z p+\partial_z^2 w\right\rangle+F_i[1] \partial_y w|_{y=h+\eta}+\left\langle \partial_y^2F_i \left.\left[\frac{y-\eta}{h}\right]\right|w\right\rangle=0 \end{gather}

with $i=\{0,1,2\}$, and where we have integrated by parts the wall-normal diffusion terms in order to introduce explicitly the second-order dynamic boundary conditions (2.21a,b). The projections onto $F_0$ readily produce the following momentum conservation equations:

(2.26a) \begin{equation} \left(1+\eta'^2+\frac{h'\eta'}{2}-\frac{3h'^2}{10}+ \frac{h(23h''+15\eta'')}{40}\right)q_x+\frac{2h h' q_x'}{5}- \frac{2h^2q_x''}{5}=\frac{2h^3}{3}, \end{equation}
(2.26b)\begin{align} \left(1+2\eta'^2+h'\eta'-\frac{8h'^2}{5}+\frac{h(24h'' +15\eta'')}{10}\right)q_z&+\frac{9h h' q_z'}{5}- \frac{9h^2q_z''}{5} \nonumber\\ & =\dfrac{h^3}{3}({-}2\cot\theta(h'+\eta')+{We}(h'''+\eta''')), \end{align}

while the corrections $\boldsymbol {r}$ and $\boldsymbol {s}$ can be computed from the other residual equations $i=\{1,2\}$.

The streamwise shear rate is then obtained from the corrected ansatz (2.22):

(2.27)\begin{equation} \tau=\left(1+\frac{51h'\eta'}{64}-\frac{h'^2}{320}+ \frac{h(53h''-75\eta'')}{640}\right)\frac{3q_x}{h^2}+ \frac{69h'-315\eta'}{320}\frac{q_x'}{h}-\frac{q_x''}{5}. \end{equation}

3. Linear stability analysis

3.1. Linearised Stokes equations

As mentioned in § 2.4, Bertagni & Camporeale (Reference Bertagni and Camporeale2021) provide an exact solution in closed form to the hydrodynamical problem ((2.1), (2.2)), linearised around the basic flow (2.7c,d), with boundary conditions (2.6c,d), linearised and flattened around the flat-rock solution (2.7a,b). Thanks to the streamwise invariance, the continuity equation enables them to write the cross-stream components of the perturbation velocity field in terms of a scalar streamfunction $\varphi$, separate from the streamwise perturbation velocity. After a normal mode decomposition

(3.1)\begin{equation} (\boldsymbol{u};h,\eta)=(\boldsymbol{u}_b;1,0) +\epsilon(\tilde{u}(y)\hat{\boldsymbol{x}}-\mathrm{i} k\tilde{\varphi}(y)\hat{\boldsymbol{y}}+\partial_y \tilde{\varphi}(y)\hat{\boldsymbol{z}};\tilde{h}, \tilde{\eta})\exp[\mathrm{i}(k z-\omega t)]+\textrm{c.c.}, \end{equation}

with real transverse wavenumber $k$ and angular frequency $\omega$, the imaginary part of which is the temporal growth rate, their solution reads in the present notation

(3.2a,b) \begin{equation} \left.\begin{array}{@{}l@{}} \tilde{u}=2\operatorname{sech} k \left((\tilde{h}+\tilde{\eta})\dfrac{\sinh[k y]}{k}-\tilde{\eta}\cosh[k-k y]\right),\\ \tilde{\varphi}=\mathrm{i}(2\cot\theta+{We} k^2) \\ \qquad \times\dfrac{(k y\coth[k y]-1)(k\sinh k+\cosh k)-k^2 y\cosh k}{k^2(2k^2+\cosh[2k]+1)} \sinh[k y](\tilde{h}+\tilde{\eta}). \end{array}\right\} \end{equation}

The linearised and flattened kinematic boundary condition and dissolution model ((2.5b), (2.12)) can then form the system of solvability equations

(3.3a)\begin{gather} \omega \tilde{h}=k\tilde{\varphi}|_{y=1}, \end{gather}
(3.3b)\begin{gather}\mathrm{i}\omega \tilde{\eta}=\dfrac{V}{6} (\partial_y \tilde{u}|_{y=0}+\tilde{\eta}\partial_y^2 u_b|_{y=0}), \end{gather}

which is singular and hence allows for non-trivial solutions when

(3.4)\begin{equation} \omega_\pm^{lS}(k)=\dfrac{\omega_0^{lS}}{2}\pm \sqrt{\left(\dfrac{\omega_0^{lS}}{2}\right)^2-\dfrac{V(\cosh k\sinh k-k)(1-k\tanh k)}{3k(2 k^2+\cosh[2k]+1)}(2\cot\theta+{We} k^2)},\end{equation}

with

(3.5)\begin{equation} \omega_0^{lS}={-}\mathrm{i}\left\{\dfrac{\cosh k\sinh k-k}{k(2 k^2+\cosh[2k]+1)}(2\cot\theta+{We} k^2)+\dfrac{V}{3}(\operatorname{sech} k+k\tanh k-1)\right\}. \end{equation}

3.2. Linearised depth-integrated models

The two depth-integrated models, developed in § 2.4 to replace the Stokes equations, instead reduce the problem to four scalar variables $(q_x,q_z;h,\eta )$, governed by four equations ((2.17) or (2.26), (2.15), and (2.12) with either (2.19) or (2.27)). The perturbed state is expressed in the form

(3.6)\begin{equation} (q_x,q_z;h,\eta)=(2/3,0;1,0)+\epsilon (\tilde{q}_x,\tilde{q}_z;\tilde{h},\tilde{\eta}) \exp[\mathrm{i}(k z-\omega t)]+ \textrm{c.c.} \end{equation}

For the revised lubrication model, presented in § 2.4.1, the linearised perturbation evolution equations then read

(3.7a)\begin{gather} (3+k^2)\tilde{q}_x-(6+k^2)\tilde{h}-k^2\tilde{\eta}=0, \end{gather}
(3.7b)\begin{gather}(3+k^2)\tilde{q}_z+\mathrm{i}k(2\cot\theta+{We} k^2) (\tilde{h}+\tilde{\eta})=0, \end{gather}
(3.7c)\begin{gather}k \tilde{q}_z-\omega\tilde{h}=0, \end{gather}
(3.7d)\begin{gather}V(3\tilde{q}_x/2-2\tilde{h})-3\mathrm{i}\omega\tilde{\eta}=0. \end{gather}

We obtain the following dispersion relation:

(3.8)\begin{equation} \omega_\pm^{rl}(k)=\dfrac{\omega_0^{rl}}{2}\pm \sqrt{\left(\dfrac{\omega_0^{rl}}{2}\right)^2- \dfrac{V k^2 (1-2 k^2/3)}{(3+k^2)^2}(2\cot\theta+{We} k^2)}, \end{equation}

with $\omega _0^{rl}=-\mathrm {i}k^2(2\cot \theta +{We} k^2+V/2)/(3+k^2)$.

For the weighted-residuals model, presented in § 2.4.2, the linearised perturbation evolution equations instead read

(3.9a)\begin{gather} (3+6k^2/5)\tilde{q}_x-(6+23k^2/20)\tilde{h}-3k^2\tilde{\eta}/4=0, \end{gather}
(3.9b)\begin{gather}(3+27k^2/5)\tilde{q}_z+\mathrm{i}k(2\cot\theta+{We} k^2) (\tilde{h}+\tilde{\eta})=0, \end{gather}
(3.9c)\begin{gather}k \tilde{q}_z-\omega\tilde{h}=0, \end{gather}
(3.9d)\begin{gather}V(5\tilde{q}_x/4-(3/2-5k^2/384)\tilde{h}+23 k^2\tilde{\eta}/128) -3\mathrm{i}\omega\tilde{\eta}=0. \end{gather}

The dispersion relation becomes

(3.10)\begin{equation} \omega_\pm^{wr}(k)=\dfrac{\omega_0^{wr}}{2}\pm \sqrt{\left(\dfrac{\omega_0^{wr}}{2}\right)^2-\dfrac{V k^2 (1-k^2(9+k^2)/15)}{(3+6k^2/5)(3+27k^2/5)}\left(2\cot\theta+{We} k^2\right)}, \end{equation}

with

(3.11)\begin{equation} \omega_0^{wr}={-}\mathrm{i}k^2\left\{\dfrac{2\cot\theta+{We} k^2}{3+27k^2/5} +\dfrac{21V(1+46k^2/315)}{128+256k^2/5}\right\}. \end{equation}

3.3. Results

Figure 3 presents these dispersion relations ((3.4), (3.8), (3.10)) for different inclination angles $\theta$. In the presented range, the weighted-residuals dispersion relation is visually indistinguishable from the linear Stokes one. Both gravity and surface tension are destabilising at small $k$. The instability is neutralised linearly when $k\to 0$ for $\theta \ne {\rm \pi}/2$ and quadratically for a vertical substrate. Larger wavenumbers, $k\tanh k>1\Leftrightarrow k^2>1.439$ for the linear Stokes model, $k^2>3/2$ for the revised lubrication model, or $k^2>(\sqrt {141}-9)/2\approx 1.437$ for the weighted-residuals model, are stabilised by the transverse transfer of momentum, emulated strikingly well by the diffusion terms in the two extended lubrication models ((2.17), (2.26)).

Figure 3. Dispersion relations ((3.4), (3.8), (3.10)) compared to Bertagni & Camporeale (Reference Bertagni and Camporeale2021). Here ${Pe}=10^3$, $\gamma =49\ 091$ (calcite), $x=10$, $\theta =\{{\rm \pi} /8, {\rm \pi}/4, 3{\rm \pi} /8, {\rm \pi}/2\}$, ${We}\approx 4210\times \sin ^{-1/3}\theta$.

These simplified models appear to capture in a distilled form most of the necessary physics to compare rather favourably with the much more involved and complete analytical treatment of Bertagni & Camporeale (Reference Bertagni and Camporeale2021). Deviations are observed only for the cutoff wavenumber, which can be accounted for by the lack of transverse solute diffusion in Lévêque's (Reference Lévêque1928) boundary-layer dissolution model (see their figure 5c).

In the spirit of Bertagni & Camporeale (Reference Bertagni and Camporeale2021), figure 4(a) illustrates the influence of the Péclet number ${Pe}$ on the dispersion relation (3.4). Our analysis retrieves a behaviour rather faithful to what they report. We remind the reader that ${Pe}={Sc}{Re}$ with the Schmidt number ${Sc}=\nu /D$. Figure 4(b) shows the dispersion relation for different materials, to be compared again with the results of Bertagni & Camporeale (Reference Bertagni and Camporeale2021). More soluble minerals, e.g. salt, exhibit a shorter plateau of larger growth rate.

Figure 4. (a) Contour plot of the growth rate ${\textrm {Im}}[\omega _+]$ (3.4) in the {$k$, ${Pe}$} space. The most unstable wavenumber $k_m$ is displayed by a dashed line. The striped area indicates negative growth rates. Here $\gamma =49\ 091$ (calcite), $x/{Pe}=0.05$, $\theta ={\rm \pi} /4$, ${We}\approx 4.7\times 10^5\times {Pe}^{-2/3}$. (b) Dispersion relation (3.4) for several materials. Here ${Pe}=10^3$, $x=10$, $\theta ={\rm \pi} /4$, ${We}\approx 4726$.

3.4. Instability mechanism

We briefly elucidate here the instability mechanism, explained by Guérin et al. (Reference Guérin, Derr, Courrech du Pont and Berhanu2020) and depicted in figure 7 of Bertagni & Camporeale (Reference Bertagni and Camporeale2021). Let the substrate be perturbed by a sinusoidal wave in the transverse direction $z$. Then, over the troughs (i) the film gets thicker – on the one hand because of a hydrostatic pressure gradient, which tends to bring water back down in the troughs, and on the other hand because of a capillary pressure gradient, which tends to keep the free surface flat; (ii) the flow accelerates; (iii) the boundary layer becomes thinner; and (iv) the dissolution is faster than on the crests.

For shorter wavelengths, transverse transfer of momentum becomes important and the flow acceleration in step (ii) is hindered. Inside the narrower troughs the flow is confined and in fact slowed down, thus stabilising the system.

For very long wavelengths, the flattening of the free surface, step (i), becomes a growth-limiting factor and the instability is neutralised. This explains the qualitatively different asymptotic behaviour of the growth rate for $\theta ={\rm \pi} /2$: without a wall-normal component of the gravitational acceleration, surface tension is the only effect restoring the flat free surface. This is also why in the rigid-lid approximation, where $\tilde {h}+\tilde {\eta }=0$ is imposed, bypassing step (i) altogether, the unstable plateau reaches $k=0$, as highlighted by Bertagni & Camporeale (Reference Bertagni and Camporeale2021).

4. Conclusion

The dispersion relation displays a large plateau of almost equally unstable wavenumbers in the interval $10^{-3}\lesssim k\lesssim 10^{-1}$ (figure 3). This impedes to some degree the model's predictive power, since a large range of wavelengths – spanning from ten to approximately a thousand times the film thickness $H_N$ – can grow with very similar rates. On the other hand, this feature may be consistent with the wide-ranging length scales spanned by natural linear karren forms. The competition for the emerging pattern's ultimate wavelength would be settled by more subtle mechanisms. Finite amplitude and stochastic effects, such as the roughness and/or grain distribution of the initial rock surface, or raindrop impacts on the falling film, might play a role in the wavelength selection.

It is clear that the classical lubrication theory could not predict a cutoff wavenumber – while the film is thicker over the troughs, the flow would always be faster, and the troughs would dissolve more, independent of $k$. In this sense, both extended lubrication models, because they account for the transverse transfer of momentum, have proven to be suitable. Moreover, their striking adjacency underlines the place of physical insight in simplifying mathematically cumbersome systems. Finally, their validation against the exact Stokes equations, at least in the linear regime, holds a promise for their future use in nonlinear analyses.

Supplementary material

Supplementary material are available at https://doi.org/10.1017/jfm.2023.932.

Acknowledgements

The authors are grateful to T. Salamon and M. Chaigne for fruitful discussions.

Declaration of interests

The authors report no conflict of interest.

References

Bertagni, M.B. & Camporeale, C. 2021 The hydrodynamic genesis of linear karren patterns. J. Fluid Mech. 913, A34.Google Scholar
Devauchelle, O., Popović, P. & Lajeunesse, E. 2022 Viscous transfer of momentum across a shallow laminar flow. J. Fluid Mech. 932, A32.Google Scholar
Ford, D. & Williams, P. 2007 Karst Hydrogeology and Geomorphology. John Wiley and Sons.Google Scholar
Ginés, A., Knez, M., Slabe, T. & Dreybrodt, W. 2009 Karst Rock Features. Karren Sculpturing. Carsologica, vol. 9. Založba ZRC.Google Scholar
Guérin, A., Derr, J., Courrech du Pont, S. & Berhanu, M. 2020 Streamwise dissolution patterns created by a flowing water film. Phys. Rev. Lett. 125, 194502.Google Scholar
Kalliadasis, S., Ruyer-Quil, C., Sheid, B. & Velarde, M.G. 2012 Falling Liquid Films. Applied Mathematical Sciences, vol. 176. Springer.Google Scholar
Lévêque, A. 1928 Les lois de la transmission de chaleur par convection. PhD thesis, Faculté des Sciences de Paris.Google Scholar
Polyanin, A.D., Kutepov, A.M., Vyazmin, A.V. & Kazenin, D.A. 2001 Hydrodynamics, Mass and Heat Transfer in Chemical Engineering. Topics in Chemical Engineering, vol. 14. CRC.Google Scholar
Figure 0

Figure 1. (a) Natural karren in the Bernese Alps (top) and in the Jura (bottom). (b) Sketch of the problem at hand.

Figure 1

Figure 2. Two-dimensional basic concentration field $c_b$ and wall-normal profiles at $x/{Pe}=\{10^{-3}, 10^{-2}, 10^{-1}, 1\}$. The orange solid lines represent the exact solution, while the black dashed lines depict Lévêque's (1928) boundary-layer solution.

Figure 2

Figure 3. Dispersion relations ((3.4), (3.8), (3.10)) compared to Bertagni & Camporeale (2021). Here ${Pe}=10^3$, $\gamma =49\ 091$ (calcite), $x=10$, $\theta =\{{\rm \pi} /8, {\rm \pi}/4, 3{\rm \pi} /8, {\rm \pi}/2\}$, ${We}\approx 4210\times \sin ^{-1/3}\theta$.

Figure 3

Figure 4. (a) Contour plot of the growth rate ${\textrm {Im}}[\omega _+]$ (3.4) in the {$k$, ${Pe}$} space. The most unstable wavenumber $k_m$ is displayed by a dashed line. The striped area indicates negative growth rates. Here $\gamma =49\ 091$ (calcite), $x/{Pe}=0.05$, $\theta ={\rm \pi} /4$, ${We}\approx 4.7\times 10^5\times {Pe}^{-2/3}$. (b) Dispersion relation (3.4) for several materials. Here ${Pe}=10^3$, $x=10$, $\theta ={\rm \pi} /4$, ${We}\approx 4726$.

Supplementary material: File

Djambov and Gallaire supplementary material 1
Download undefined(File)
File 71.6 KB
Supplementary material: File

Djambov and Gallaire supplementary material 2
Download undefined(File)
File 102.2 KB