2.1. Governing equations

A steady water film flowing on a flat surface performs a semi-parabolic velocity profile (figure 2). This is commonly referred to as the Nusselt solution and it can be readily derived from momentum conservation (Nusselt Reference Nusselt1916). Accordingly, the dimensional film thickness and the surface velocity read

(2.1a,b)\begin{equation} \hat{h}_0 =\left( \frac{3\nu\hat{q}}{g\sin\theta}\right)^{1/3},\quad \hat{u}_s = \left(\frac{9g\hat{q}^2\sin\theta}{8\nu}\right)^{1/3}, \end{equation}

where $\nu =10^{-6}$ m$^2$ s$^{-1}$ is the water kinematic viscosity, $g$ is the gravitational acceleration, $\theta$ is the angle with the horizontal and $\hat {q}=2 \hat {u}_s \hat {h}_0/3$ is the flow rate per unit span. The hat refers to dimensional variables. We use the quantities in (2.1a,b) to scale the governing equations of the problems, which are the Navier–Stokes equations and the advection–diffusion equation for the solute concentration $c$ (scaled with its equilibrium value $\hat {c}_{{eq}}$, see table 1)

(2.2)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

(2.3)\begin{gather}{\textit{Re}} (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}+ \boldsymbol{\nabla} p) = \nabla^2 \boldsymbol{u}+\boldsymbol{f}, \end{gather}

(2.4)\begin{gather}{\textit{Pe}}\,\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} c = \nabla^2 c, \end{gather}

where $\boldsymbol {\nabla }=(\partial _x,\partial _y,\partial _z)$, $\boldsymbol {u}=(u,v,w)$ is the velocity field (figure 2), $p$ is pressure – scaled with $\rho \hat {u}_s^2$, with $\rho$ the water density – and $\boldsymbol {f}=(2,-\delta ,0)$ is the gravity term, with $\delta =2\cot {\theta }$. Because of the slow morphological evolution of the rock surface, the quasi-steady approximation ($\partial _t=0$) is used in (2.2)–(2.4). The dimensionless numbers appearing in (2.2)–(2.4) are the Reynolds and Péclet numbers, which read

(2.5a,b)\begin{equation} {\textit{Re}}=\frac{\hat{u}_s\hat{h}_0}{\nu},\quad {\textit{Pe}}=\frac{\hat{u}_s\hat{h}_0}{D}={\textit{Re}}\,{\textit{Sc}}, \end{equation}

where $D$ is the molecular diffusivity coefficient. For the solute species considered (table 1), an average value of $D\simeq 10^{-9}$ m$^2$ s$^{-1}$ can be assumed (Ford & Williams Reference Ford and Williams2013). Notice that ${\textit {Pe}}$ and ${\textit {Re}}$ are linked through the Schmidt number ${\textit {Sc}}=\nu /D=10^3$ (${\textit {Pe}}=10^3 {\textit {Re}}$).

Table 1.

Dissolution reaction, density and solubility of some representative karst minerals. Solubility is given as the amount of mineral that can be dissolved in water at 25$^{\circ }$ and $p_{\text {CO}_2}=4\cdot 10^{-4}$ atm (Dreybrodt Reference Dreybrodt2012; Ford & Williams Reference Ford and Williams2013).

2.2. Hydrodynamic boundary conditions

The boundary conditions are defined on the rock–water interface $\eta (x,z,t)$ and water–air interface $\eta (x,z,t)+h(x,z,t)$ (figure 2). It is thus convenient to introduce the vertical coordinate

(2.6)\begin{equation} \zeta=\frac{y-\eta(x,z,t)}{h(x,z,t)}, \end{equation}

so that the water domain is always defined between $\zeta =0$ (rock–water interface) and $\zeta =1$ (water–air interface). Notice that $y$ and $\eta$ are measured relative to the position of the initially flat rock surface and they are scaled with $\hat {h}_0$.

At the rock–water interface ($\zeta =0$), the velocity field satisfies the no-slip and impermeability conditions, which are respectively

(2.7a,b)\begin{equation} u=w=0,\quad v=0, \end{equation}

where we neglect the influence of the slow dissolution process in the vertical hydrodynamic velocity. At the water–air interface ($\zeta =1$), the continuity of the stress is ensured by the so-called dynamic conditions

(2.8a,b)\begin{equation} \boldsymbol{n}\cdot {\boldsymbol{\mathsf{T}}} \cdot \boldsymbol{n} ={-}{\textit{We}}\,\mathcal{K},\quad \boldsymbol{n}\cdot {\boldsymbol{\mathsf{T}}} \cdot \boldsymbol{t}=0, \end{equation}

where $\boldsymbol {n}=[\boldsymbol {\nabla }(y-\eta -h)]/n$ and $\boldsymbol {t}$ are the versors normal and tangent to the free surface ($n=\lbrace {1+[\partial _x(h+\eta )]^2+[\partial _z(h+\eta )]^2\rbrace }^{1/2}$ is the versor normalization). Here, ${\boldsymbol{\mathsf{T}}}$ is the stress tensor (${\boldsymbol{\mathsf{T}}}_{ij}=p \delta _{ij} - (\partial _{x_i} u_j +\partial _{x_j} u_i)/{\textit {Re}}$). The term ${\textit {We}}\,\mathcal {K}$ accounts for the normal stress induced by the surface tension $\sigma$, where ${\textit {We}}$ is the Weber number and $\mathcal {K}$ is the mean curvature of the free surface (e.g. Chang & Demekhin Reference Chang and Demekhin2002)

(2.9)\begin{gather} {\textit{We}}=\frac{\sigma}{\rho \hat{h}_0 \hat{u}_s^2}=\frac{{\textit{Ka}}}{{\textit{Re}}^{5/3}\sin{(\theta)}}, \end{gather}

(2.10)\begin{gather}\mathcal{K}=\partial_x \left[\frac{\partial_x(h+\eta)}{n}\right]+\partial_z \left[\frac{\partial_z(h+\eta)}{n}\right]. \end{gather}

From (2.9) we notice that ${\textit {We}}$ is a function of ${\textit {Re}}$ and the angle $\theta$ through the Kapitza number ${\textit {Ka}}=2^{1/3} \sigma /(g^{1/3} \nu ^{4/3} \rho )$; ${\textit {Ka}}$ only depends on the water properties and is thus constant for our purposes.

The temporal evolution of the free surface is described by the kinematic condition in $\zeta =1$

(2.11)\begin{equation} \partial_t h = \boldsymbol{u}\cdot\boldsymbol{n}=v-u (h+\eta)_x - w (h+\eta)_z, \end{equation}

where, as in (2.7a,b), we neglected the influence of $\eta _t$ in the hydrodynamics. Equation (2.11) ensures the water mass conservation and will be later used as the first solvability equation.

2.3. Concentration boundary conditions

Differently from the hydrodynamic quantities, the concentration $c$ has also a longitudinal dependence due to the downstream solute accumulation within the water film (see figure 2). Basically, this makes the concentration field a three-dimensional problem. In fact, $c$ depends on $x$ because of the downstream solute accumulation, on $y$ due to the normal-to-wall effect of the dissolution occurring at the rock–water interface and on $z$ because of the transverse nature of linear karren forms.

Upstream, as the water film is generally produced by rain or snow melt, it might be considered free of solute

(2.12)\begin{equation} c\vert_{x=0}=0. \end{equation}

With the dissolution of the rock surface, solute molecules are dissociated from the solid through chemical reactions and then transported into the film bulk by diffusion. In general, the slowest between these two processes (surface reaction and transport by diffusion) controls the dissolution kinetics of a karst rock. Very soluble minerals, such as salt, have a rapid dissolution reaction and thus the kinetics is always diffusion controlled (Ford & Williams Reference Ford and Williams2013). For relatively less soluble minerals, such as gypsum or calcite, numerous studies have demonstrated that in very undersaturated waters, such as those commonly found in exposed karst environments, a diffusion-controlled dissolution kinetics prevails (see the review on dissolution kinetics by Morse & Arvidson (Reference Morse and Arvidson2002), and references therein). As equilibrium is approached, there is a transition region to a surface-controlled dissolution kinetics. Moreover, for calcite, which dissolves in water enriched with CO$_2$, the slow uptake of CO$_2$ from the atmosphere may also become a limiting factor as equilibrium is approached (Buhmann & Dreybrodt Reference Buhmann and Dreybrodt1985; Kaufmann & Dreybrodt Reference Kaufmann and Dreybrodt2007). Notice that CO$_2$-conversion does not influence salt and gypsum dissolution as they directly dissociate in pure water (table 1).

Since we wish to model incipient karren formation and linear karren forms have been observed to develop upstream first, i.e. where the water is very undersaturated (Glew & Ford Reference Glew and Ford1980; Ginés et al. Reference Ginés, Knez, Slabe and Dreybrodt2009; Slabe, Hada & Knez Reference Slabe, Hada and Knez2016), we consider a dissolution kinetics that is diffusion controlled. This implies that surface reactions are fast compared to diffusion and that the solute accumulates in a very thin diffusion boundary layer at the solid–liquid interface that is nearly saturated (e.g. Ford & Williams Reference Ford and Williams2013). By also imposing no flux of solute between water and air, the vertical boundary conditions eventually read

(2.13a,b)\begin{equation} \partial_\zeta c\vert_{\zeta=1} = 0,\quad c\vert_{\zeta=0}=1. \end{equation}

The morphological evolution of the rock surface is thus regulated by the flux of solute that moves into the bulk of the water film through diffusion, i.e. $\rho _s (\hat {V}+\partial _{\hat {t}}\hat {\eta })= D \partial _{\hat {y}} \hat {c}$, where $\rho _s$ is the rock density and $\hat {V}$ is the dissolution rate induced by the uniform hydrodynamic flow ($\hat {V}<0$). In dimensionless form, the same equation reads

(2.14)\begin{equation} \gamma\,{\textit{Pe}}\,h(V+\partial_t \eta) = \partial_\zeta c, \end{equation}

where $\gamma =\rho _s/\hat {c}_{eq}\gg 1$ is related to the mineral considered (table 1). Equation (2.14) will be used as the second solvability equation.

3.1. Hydrodynamics

The Navier–Stokes equations (2.2) and (2.3) reduce to

(3.1a,b)\begin{equation} u_0''={-}2,\quad p_0'={-}\delta/{\textit{Re}}, \end{equation}

with boundary conditions

(3.2a–c)\begin{equation} u_0\vert_{\zeta=0}=0,\quad u_0'\vert_{\zeta=1}=0,\quad p_0\vert_{\zeta=1}=0, \end{equation}

where $'$ denotes differentiation in $\zeta$. The solutions are the Nusselt semi-parabolic profile for the longitudinal velocity and the hydrostatic pressure profile

(3.3a–c)\begin{equation} u_0=(2-\zeta)\zeta,\quad v_0=w_0=0,\quad p_0=\frac{\delta}{{\textit{Re}}} (1-\zeta). \end{equation}

3.2. Concentration

The advection–diffusion equation (2.4) reduces to

(3.4)\begin{equation} {\textit{Pe}}\,u_0(\zeta)\partial_x c_0 = c_0'', \end{equation}

where the longitudinal diffusion has been neglected since its effect is smaller than that of longitudinal advection. The initial and boundary conditions (2.12) and (2.13a,b) read

(3.5a–c)\begin{equation} c_0\vert_{x=0}=0,\quad c_0\vert_{\zeta=0}=1,\quad c_0'\vert_{\zeta=1}=0. \end{equation}

Equation (3.4) highlights the coupling between the concentration field $c_0$ and the velocity profile $u_0(\zeta )$, as well as the $x$-dependency of the concentration distribution due to the downstream solute accumulation. The solution to this problem is resumed in Polyanin et al. (Reference Polyanin, Kutepov, Kazenin and Vyazmin2001, pp. 130–132) and it is here more extensively reported. For $x\ge 0$, we seek the solution in the form of the series (Davis Reference Davis1973)

(3.6)\begin{equation} c_0=1-\sum_{m=0}^{\infty} A_m F_m(x)G_m(\zeta), \end{equation}

where $A_m, F_m(x),G_m(\zeta )$ are unknown functions that must be determined. Substituting the expansion (3.6) into (3.4) and then separating the variables, we obtain

(3.7)\begin{gather} \partial_x F_m(x)+\frac{\alpha_m^2}{{\textit{Pe}}}F_m(x)=0, \end{gather}

(3.8)\begin{gather}G_m''(\zeta)+\alpha^2_m u_0(\zeta) G_m(\zeta)=0. \end{gather}

The solution to the $x$-problem (3.7) is readily given by

(3.9)\begin{equation} F_m(x)=\exp{\left(-\alpha^2_m\frac{x}{{\textit{Pe}}}\right)}. \end{equation}

Instead, (3.8) defines a Sturm–Liouville problem, whose boundary conditions, that are obtained after replacement of (3.6) into the second and third equations in (3.5a–c), are

(3.10a,b)\begin{equation} G_m\vert_{\zeta=0}=0,\quad G_m'\vert_{\zeta=1}=0. \end{equation}

The solution of the Sturm–Liouville problem for the eigenfunctions $G_m$ (up to a constant factor) is

(3.11)\begin{equation} G_m(\zeta)=\exp{\left[-\frac{\alpha_m}{2}(1-\zeta)^2\right]} {\Phi}\left(\frac{1}{4}-\frac{\alpha_m}{4},\frac{1}{2};\alpha_m(1-\zeta)^2\right), \end{equation}

where ${\Phi }(\cdot ,\cdot;\cdot )$ is the degenerate hypergeometric function (Abramowitz, Stegun & Romer Reference Abramowitz, Stegun and Romer1988). Substituting (3.11) into the first of the boundary conditions (3.10a,b) provides the transcendental equation for the eigenvalues $\alpha _m>0$

(3.12)\begin{equation} {\Phi}\left(\frac{1}{4}-\frac{\alpha_m}{4},\frac{1}{2};\alpha_m\right)=0. \end{equation}

For the definition of the coefficients $A_m$, the first step is to substitute the series (3.6) into the first equation in (3.5a–c), which yields

(3.13)\begin{equation} \sum_{m=0}^{\infty} A_m G_m(\zeta)=1. \end{equation}

Multiplying (3.13) by $u_0G_n$, with $n\neq m$, then integrating between $\zeta =0$ and $\zeta =1$, and using the orthogonality condition $\int _0^1 u_0 G_m G_n d\zeta =0$, one obtains

(3.14)\begin{equation} A_m=\frac{\int_0^1 u_0 G_m(\zeta)\,\textrm{d}\zeta}{\int_0^1 u_0 G_m(\zeta)^2\,\textrm{d}\zeta}. \end{equation}

Notably, instead of solving numerically (3.12) and (3.14), it is possible to evaluate the constants $A_m$ and $\alpha _m$ with the approximated relationships (Polyanin & Nazaikinskii Reference Polyanin and Nazaikinskii2015)

(3.15)\begin{gather} \alpha_m=4m+1.68 \quad (m=0,1,2\cdots), \end{gather}

(3.16a,b)\begin{gather}A_0=1.2,\quad A_m=({-}1)^m 2.27\alpha_m^{{-}7/6} \quad (m=1,2,3\cdots), \end{gather}

whose maximum error is less than 0.2 %. It is also worth mentioning that, for $x/{\textit {Pe}} \rightarrow 0$, an easier asymptotic solution for $c_0$ can be obtained (Polyanin et al. Reference Polyanin, Kutepov, Kazenin and Vyazmin2001)

(3.17)\begin{equation} c_0={\Gamma}\left(\frac{1}{3},\frac{2 \zeta ^3}{9 X}\right)/{\Gamma} \left(\frac{1}{3}\right) \quad \mbox{for}\ x/{\textit{Pe}}=X\rightarrow 0, \end{equation}

where $\Gamma (\cdot )$ and ${\varGamma }(\cdot ,\cdot )$ are the complete and incomplete gamma functions (Abramowitz et al. Reference Abramowitz, Stegun and Romer1988); $X=x/{\textit {Pe}}$ is a convenient scaling that maps the dissolution domain of the problem between $X=0$, where the water is free of solute, and $X=1$, where the water film is basically saturated. In fact, some easy algebra shows that the time scale $\hat {\tau }$ needed by water to go from $X=0$ to $X=1$ is $\hat {\tau }\sim \hat {h}_0^2/D$, i.e. equivalent to the time scale required by diffusion to involve the full water depth.

The spatial trends of the exact solution (3.6) and the asymptotic solution (3.17) for $c_0$ are shown in figure 3.

Figure 3.

Vertical profiles of $c_0$ at different $X$ (log scale). The solid red lines refer to the exact solution (3.6). The dotted-black lines refer to the asymptotic solution (3.17), which is shown to be valid up to $X\sim 0.1$. The coloured areas evidence the saturation ratio of water in the solute.

3.3. Dissolution

The dissolution rate $V$ can be evaluated by combining the evolution equation for the rock surface (2.14) and the solution for $c_0$

(3.18)\begin{equation} V=\frac{1}{\gamma\,{\textit{Pe}}}c_0' \vert_{\zeta=0}. \end{equation}

It is also convenient to define the saturation ratio $s$ as

(3.19)\begin{equation} s=\int_0^1 c_0\,\textrm{d}\zeta, \end{equation}

where $s=0$ indicates pure water and $s=1$ corresponds to saturated water (ion pairing is neglected). Because the concentration has been scaled with its saturation value, $s$ is also equivalent to the dimensionless depth-averaged concentration. A graphical evidence of the saturation ratio of the water film is given by the light red areas in the panels of figure 3.

The trend $V=V(X)$, obtained after replacing (3.6) in (3.18), is reported in figure 4(a). Notice that the quantity $\gamma {\textit {Pe}}\,\lvert V \rvert$ is independent of the rock type and the film hydrodynamics. The dissolution rate decreases with $X$ as the water loses chemical aggressivity by increasing its saturation ratio $s$ (right $y$-axis). Because clear water enhances dissolution, half-saturation ($s=0.5$) is reached in the upstream part of the flow ($X\sim 0.2$) and the rate of solute accumulation progressively decreases downstream. Although at $X=1$ the flow is not completely saturated (coloured area in figure 3d) at those levels of saturation – and even before ($s>0.9$) – dissolution is known to drop drastically because of the inhibitory action of other chemical species in the rock texture (Kaufmann & Dreybrodt Reference Kaufmann and Dreybrodt2007; Ford & Williams Reference Ford and Williams2013).

Figure 4.

Dissolution of an initially flat rock. (a) Longitudinal trends of the scaled dissolution rate $\gamma {\textit {Pe}} \lvert V\rvert$ (blue) and the saturation ratio $s$ (red); (b) $\gamma {\textit {Pe}} \lvert V\rvert$ versus $s$, from the previous panel, highlighting a nonlinear relationship with a change of dissolution rate around $s\sim 0.3$.

Between the saturation ratio $s$ and the dissolution rate, there is a nonlinear relationship (figure 4b). This nonlinearity, which here arises thanks to the spatially dependent solution for the concentration, highlights the limits of the common assumption of a diffusive flux linearly related to the saturation ratio, i.e. $V\propto (s-1)$. Indeed, the dissolution of the rock surface is better specified by the local gradient of concentration at the rock surface ($c_0'\vert _{\zeta =0}$), rather than by the average concentration gradient along the film thickness ($s-1$ dimensionally scales with $(\hat {c}-\hat {c}_{{eq}})/\hat {h}_0$). Furthermore, it is remarkable to notice that the trend in figure 4(b) is qualitatively similar to experimental evidence for calcite dissolution (Buhmann & Dreybrodt Reference Buhmann and Dreybrodt1985; Kaufmann & Dreybrodt Reference Kaufmann and Dreybrodt2007), even though we here consider a simplified dissolution kinetics that is controlled just by diffusion. Kaufmann & Dreybrodt (Reference Kaufmann and Dreybrodt2007) have in fact shown that two quasi-linear regimes (for $s<0.3$ and $s>0.3$) can be considered in the dissolution rate of calcite. These regimes have linear coefficients separated by an order of magnitude and they have been qualitatively plotted with dashed lines in figure 4(b).

4.1. Hydrodynamics

Since the hydrodynamic problem is longitudinally invariant, the continuity equation (2.2) at the linear order ($\varepsilon$) imposes $v_1'+\mathrm {i}\,kw_1=0$. Therefore, we can introduce the scalar streamfunction $\varphi$, such that

(4.2a,b)\begin{equation} v_1={-}\mathrm{i} k \varphi,\quad w_1=\varphi'. \end{equation}

After some algebraic manipulations, the Navier–Stokes equations (2.2)-(2.3) are reduced to the Orr–Sommerfeld equation for a domain longitudinally invariant (e.g. Chang & Demekhin Reference Chang and Demekhin2002; Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2011)

(4.3)\begin{equation} \varphi^{iv}-2 k^2 \varphi''+k^4 \varphi =0, \end{equation}

with boundary conditions

(4.4)\begin{gather} \varphi'=\varphi=0 \quad \text{in}\ \zeta=0, \end{gather}

(4.5)\begin{gather}\varphi ''+k^2 \varphi=\varphi ^{(3)}-3 k^2 \varphi '-\mathrm{i} k (\eta_1+h_1) (k^2 {\textit{Re}}\,{\textit{We}}+\delta)=0\quad \text{in}\ \zeta=1. \end{gather}

We recall that (4.4) are the no-slip and impermeability conditions at the rock–water interface (corresponding to (2.7a,b)), while (4.5) are the linearized dynamic conditions at the water–air interface, (2.8a,b). The advantage of the Orr–Sommerfeld approach is that the transverse hydrodynamic problem is reduced to the biharmonic equation (4.3), the solution thereof is

(4.6)\begin{equation} \varphi=\mathrm{i}\,(\delta+k^2 {\textit{Re}} {\textit{We}})(\eta_1+h_1)f(\zeta), \end{equation}

where

(4.7)\begin{equation} f(\zeta)=\frac{(k \zeta \coth (k \zeta )-1) (k \sinh{k}+\cosh{k})-k^2 \zeta \cosh{k}}{k (2 k^2+\cosh (2 k)+1)}\sinh{(k \zeta)}. \end{equation}

Combining (4.2a,b) and (4.6) readily provides the solution for $v_1$ and $w_1$ (where the latter can be shown to be much larger than the former). Because of the longitudinal invariance of the linear problem, $u_1$ is not explicitly related to $\varphi$. Thus, for its evaluation, we directly address the conservation of longitudinal momentum (2.3), which at order $\varepsilon$ reads

(4.8)\begin{equation} u_1''-k^2 u_1={-}k^2 u_0'(\eta_1+\zeta h_1)+{\textit{Re}}\,u_0'v_1 +2 u_0'' h_1, \end{equation}

with boundary conditions

(4.9a,b)\begin{equation} u_1\vert_{\zeta=0}=0,\quad u_1'\vert_{\zeta=1}=0. \end{equation}

The terms in the left side of (4.8) are the vertical and transverse diffusion of momentum, respectively. The terms in the right side arise from advection (${\textit {Re}}\,u_0'v_1$) and the change of the vertical coordinate (from $y$ to $\zeta$, see (2.6)).

By solving the (4.8) through the method of variation of parameters (e.g. Bender & Orszag Reference Bender and Orszag2013), we eventually obtain

(4.10)\begin{equation} u_1=2\textrm{sech}{k}\left[(\eta_1+h_1) \sinh (k \zeta )/k - \eta_1 \cosh (k -k \zeta )\right]+2(1-\zeta ) (\eta_1+\zeta h_1), \end{equation}

where the contribution of $v_1$ in (4.10) has been neglected as it would considerably complicate the solution without any perceptible numerical change.

4.2. Concentration

At $\mathcal {O}(\varepsilon )$, the linearized advection–diffusion equation (2.4) reads

(4.11)\begin{equation} c_1''-k^2 c_1 + \varUpsilon(\zeta,X) =u_0(\zeta)\partial_X c_1 \end{equation}

where $\varUpsilon$ is the inhomogeneous term

(4.12)\begin{equation} \varUpsilon(\zeta,X)=k^2 c_0' (\eta_1+\zeta h_1)-{\textit{Pe}}\,c_0' v_1-2 h_1 c_0'' -u_1 \partial_X c_0, \end{equation}

and the initial and boundary conditions (2.12) and (2.13a,b) are

(4.13a–c)\begin{equation} c_1\vert_{X=0} =0,\quad c_1\vert_{\zeta=0}=0,\quad c_1'\vert_{\zeta=1}=0. \end{equation}

The solution to the problem (4.11)–(4.13a–c) is (Polyanin & Nazaikinskii Reference Polyanin and Nazaikinskii2015)

(4.14)\begin{equation} c_1=\int_0^X\int_0^1 \varUpsilon(\zeta_m,X_m)\mathcal{G}(\zeta,\zeta_m,X-X_m)\,\textrm{d}\zeta_m \,\textrm{d} X, \end{equation}

where $\zeta _m$ and $X_m$ are dummy variables, and $\mathcal {G}$ is the modified Green's function

(4.15a,b)\begin{equation} \mathcal{G}(\zeta,\zeta_m,X)= \sum_{m=0}^{\infty} \frac{\mathcal{H}_m(\zeta)\mathcal{H}_m(\zeta_m)}{\|\mathcal{H}_m^2\|}\mathcal{F}_m(X),\quad \|\mathcal{H}_m^2\|=\int_0^1 u_0 \mathcal{H}_m^2 \,\textrm{d}\zeta. \end{equation}

Here, $\mathcal {H}_m(\zeta )$ and $\mathcal {F}_m(X)$ are the solutions associated with the homogeneous part of (4.11). They can be obtained, after separation of variables, following the procedure used to solve the base-state problem (3.4)–(3.5a–c) for $c_0$. The $X$-problem leads to

(4.16)\begin{equation} \mathcal{F}_m(X)=\exp{\left(-\lambda^2_mX\right)}. \end{equation}

The $\zeta$-problem defines the Sturm–Liouville problem

(4.17)\begin{gather} \mathcal{H}_m''(\zeta)+\left(\lambda^2_m u_0(\zeta)-k^2\right) \mathcal{H}_m(\zeta)=0, \end{gather}

(4.18a,b)\begin{gather}\mathcal{H}_m\vert_{\zeta=0}=0,\quad \mathcal{H}_m'\vert_{\zeta=1}=0, \end{gather}

whose solution for the eigenfunctions $\mathcal {H}_m$ (up to a constant factor) is

(4.19)\begin{equation} \mathcal{H}_m(\zeta)=\exp{\left[-\frac{\lambda_m}{2}(1-\zeta)^2\right]} {\Phi}\left(\frac{k^2+\lambda_m-\lambda_m^2}{4\lambda_m},\frac{1}{2};\lambda_m(1-\zeta)^2\right). \end{equation}

Substituting (4.19) into the first of the boundary conditions (4.18a,b) provides the transcendental equation for the eigenvalues $\lambda _m>0$

(4.20)\begin{equation} {\Phi}\left(\frac{k^2+\lambda_m-\lambda_m^2}{4\lambda_m},\frac{1}{2};\lambda_m\right)=0, \end{equation}

which highlights how the eigenvalues $\lambda _m$ are functions of the wavenumber $k$.

We recall that, to evaluate the rock dissolution, we only need the concentration gradient at the rock–water interface (see (2.14)). Using the solution (4.14) and separating the terms related to $h_1$ and $\eta _1$, the gradient at the rock–water interface reads

(4.21)\begin{equation} c_1'\vert_{\zeta=0}=\mathcal{I}_{h_1} h_1+\mathcal{I}_{\eta_1}\eta_1, \end{equation}

where the expressions for the integral terms $\mathcal {I}_{h_1}$ and $\mathcal {I}_{\eta _1}$ are reported in appendix A.

4.3. Dispersion relationship

At this point, (2.11) and (2.14) are adopted as solvability equations. These describe the temporal evolution of the rock–water and water–air interfaces and, at order $\varepsilon$, they read

(4.22)\begin{gather} \omega h_1 = \mathrm{i} k \varphi \vert_{\zeta=1} \end{gather}

(4.23)\begin{gather}\omega \eta_1 = \frac{1}{\gamma {\textit{Pe}}} c_1'\vert_{\zeta=0}-h_1V. \end{gather}

Substituting the solutions for the streamfunction (4.6) and the concentration gradient (4.21), we obtain the linear system

(4.24)\begin{equation} \omega \left(\begin{array}{c} h_1 \\ \eta_1 \end{array}\right) = \left(\begin{array}{cc} a_1 & a_2\\ a_3 & a_4 \end{array}\right)\cdot\left(\begin{array}{c} h_1 \\ \eta_1 \end{array}\right), \end{equation}

where

(4.25)\begin{gather} a_1=a_2=\frac{(k^2 {\textit{Re}}\,{\textit{We}}+\delta) (k-\cosh{k} \sinh{k})}{k (2 k^2+\cosh (2 k)+1)}, \end{gather}

(4.26a,b)\begin{gather}a_3=\frac{\mathcal{I}_{h_1}}{\gamma\,{\textit{Pe}}}-V,\quad a_4=\frac{\mathcal{I}_{\eta_1}}{\gamma\,{\textit{Pe}}}. \end{gather}

The system (4.24) admits two non-trivial solutions for the eigenvalues

(4.27)\begin{equation} \omega_{\eta,h}=\frac{(a_1+a_4)\pm\sqrt{(a_1+a_4)^2-4a_1(a_4-a_3)}}{2}, \end{equation}

which we discriminate with the subscripts $h$ and $\eta$, as it will be shown in the following that the former is related to the free surface and the latter to the rock surface.

By assuming the rigid lid approximation (RLA) on the free surface i.e. imposing $h_t=0$ in the kinematic condition (2.11), the first equation in the system (4.24) reduces to $h_1=-\eta _1$. In this case, the problem provides just one eigenvalue

(4.28)\begin{equation} \omega_{{RLA}}=V+\frac{(\mathcal{I}_{\eta_1}-\mathcal{I}_{h_1})}{\gamma{\textit{Pe}}}. \end{equation}

By substituting $h_1=-\eta _1$ into the streamfunction (4.6), it is also evident that the RLA induces zero spanwise and vertical velocities, i.e. $v_1=w_1=0$. The existence of a secondary flow in the cross-sectional plane is in fact strictly related to the free surface dynamics. Instead, the streamwise velocity perturbation (4.10) is still non-null as it is triggered by the bottom perturbation. In the next section, we will discuss on the validity and the convenience of this approximate solution in the context of karren formation.