2. Governing equations

We consider a two-dimensional viscous, incompressible thin liquid film flowing under the action of gravity on a saturated inclined anisotropic porous plate of height $d$. The porous medium is considered to be bounded on one side by a rigid wall (see figure 1). Fundamental fluid properties, i.e. density ($\rho$), surface tension ($\sigma$) and kinematic viscosity ($\nu$), are supposed to remain constant. The thickness of the entire domain, comprising the porous substrate and the liquid layer, is denoted by $H$. The origin of the vertical axis is chosen at the solid boundary, $x$ refers to the coordinate in the flow direction and $y$ to the cross-stream direction. The porous medium is uniform, with a constant porosity $\varepsilon$, but anisotropic, with different permeabilities $\kappa _{xH}$ and $\kappa _{yH}$ in the streamwise $x$ and cross-stream $y$ directions. The ‘one-domain’ approach considers a composite domain where the properties of the medium vary continuously from those of the porous medium to those of the liquid medium in a thin transitional region which replaces the liquid–porous interface. Material properties, such as porosity ($\varepsilon$) and permeability in the $y$ direction ($\kappa _y$), show significant and continuous variations in the transitional layer from constant values $\varepsilon _H$ and $\kappa _{2H}$ in the porous medium to the values $1$ and $\infty$ in the liquid layer. Apart from this zone, the porous medium is considered to be homogeneous.

Figure 1. Schematic diagram of a falling film over an anisotropic saturated porous medium.

The dimensionless governing equations for the liquid and porous medium are given by (Hirata, Goyeau & Gobin Reference Hirata, Goyeau and Gobin2009)

(2.1a)$$\begin{gather} u_x+v_y = 0, \end{gather}$$

(2.1b)$$\begin{gather}\frac{Re}{\varepsilon(y)}\left[{u_t}+\frac{1}{\varepsilon}({u u_x}+{v u_y})+{uv}\frac{\partial}{\partial y}\left(\frac{1}{\varepsilon}\right)\right] ={-}Re \,p_x +\frac{1}{\varepsilon}(u_{xx}+u_{yy}) -\frac{ u}{{\kappa}_x(y)} +\frac{Re}{Fr^2}, \end{gather}$$

(2.1c)$$\begin{gather}\frac{Re}{\varepsilon(y)}\left[{v_t}+\frac{1}{\varepsilon}({u v_x}+{v v_y})+{v^2}\frac{\partial}{\partial y}\left(\frac{1}{\varepsilon}\right)\right] ={-}Re \,p_y +\frac{1}{\varepsilon}(v_{xx}+v_{yy}) -\frac{v}{{\kappa}_y(y)} -\cot\theta\,\frac{Re}{Fr^2}. \end{gather}$$

The boundary conditions are

(2.1d)\begin{equation} \mbox{at} \quad y = 0,\quad u = 0, \quad v = 0. \end{equation}

At the free surface $y = H$ we have the conditions

(2.1e)$$\begin{gather} H_t+u H_x = v, \end{gather}$$

(2.1f)$$\begin{gather}-4u_x H_x + (u_y+v_x)(1-{H_x^2}) = 0, \end{gather}$$

(2.1g)$$\begin{gather}p-\frac{2}{Re (1+{H_x^2})} ({H_x^2} u_x - (u_y+v_x)H_x+v_y) ={-} \frac{We\,H_{xx}}{(1+{H_x^2})^{3/2}}. \end{gather}$$

The distributions of porosity, $\varepsilon (y)$, and permeabilities, $\kappa _x(y)$ and $\kappa _y(y)$, are assumed to follow hyperbolic tangent profiles, which allow one to control easily the thickness of the transition layer of typical thickness $\varDelta$:

(2.1h)$$\begin{gather} \varepsilon(y) = \frac{1+\varepsilon_H}{2}+\frac{1-\varepsilon_H}{2} \tanh\left[\frac{y-\delta}{\varDelta}\right], \end{gather}$$

(2.1i)$$\begin{gather}\frac{1}{{\kappa_x}(y)} = \frac{1}{Da}\left(\frac{1}{2}-\frac{1}{2} \tanh \left[\frac{y-\delta}{\varDelta}\right]\right), \quad \frac{1}{{\kappa_y}(y)} =\frac{\xi}{{\kappa_x}(y)}, \end{gather}$$

where $\delta$ refers to the dimensionless thickness of the porous medium. In the limit $\varepsilon \rightarrow 1$, $\kappa _x, \kappa _y \rightarrow \infty$, equations (2.1b) and (2.1c) reduce to the classical momentum balance for flow of an incompressible fluid on a rigid wall.

System (2.1) involves five dimensionless parameters. Besides the Reynolds, Froude, Weber and Darcy numbers, $Re$, $Fr$, $We$ and $Da$,

(2.2a–d)\begin{equation} Re=\frac{H_N U_N}{\nu}, \quad Fr=\frac{U_N}{\sqrt{g\sin \theta\, H_N}},\quad We=\frac{\sigma}{\rho H_N {U_N^2}},\quad Da =\frac{ {\kappa_{xH} }}{{H_N^2}}, \end{equation}

we define an anisotropy parameter $\xi$ as

(2.3)\begin{equation} \xi = {\kappa_{xH}/\kappa_{yH}}, \end{equation}

which compares the permeabilities in the bulk of the porous medium. The characteristic scales for length and velocity, $H_N$ and $U_N$, respectively, are chosen to correspond to the total thickness of the entire porous and liquid layers, and to the free-surface velocity of the uniform film which characterises the wavy motion of the free surface. The governing equations have been made dimensionless with a pressure scale equal to $\rho g H_N$, as the flow is driven by gravity. It is useful to introduce the Kapitza number $Ka = \sigma /[\rho \nu ^{4/3} (g \sin \theta )^{1/3}]= We \,Fr^{4/3} Re^{2/3}$, which compares surface tension, viscosity and gravity. This parameter depends only on the fluid properties and the inclination of the wall.

2.1. Base flow

The base flow corresponds to a film of constant thickness whose stability characteristics are investigated within the framework of linear stability analysis. Let $(U_B, P_B)$ represent the solution for the base flow, then it satisfies

(2.4a)$$\begin{gather} \frac{u_{y y}}{\varepsilon(y)}-\frac{u}{\kappa_x(y)}+\frac{Re}{Fr^2} = 0, \end{gather}$$

(2.4b)$$\begin{gather}-Re\,p_y -\frac{\cot \theta\, Re}{Fr^2} = 0, \end{gather}$$

with the boundary conditions

(2.4c)\begin{equation} u = 0 \quad \mbox{at} \ y = 0, \quad u_y = 0 \quad \mbox{at} \ y = H, \quad\mbox{and}\quad p = 0 \quad \mbox{at} \ y = H. \end{equation}

We introduce a linear differential operator $L \equiv \partial _{yy}-(\varepsilon /\kappa _x)\mathrm {Id}$ for our convenience to account for the viscous diffusion and viscous drag at the solid phase in the porous medium, where $\mathrm {Id}$ refers to the identity operator. By construction, $U_B(y=H)=1$, as the free-surface velocity of the uniform film solution is the velocity scale. Therefore (2.4a) can be written as

(2.5a–c)\begin{equation} L u ={-} \frac{Re}{Fr^2}\varepsilon, \quad u(0)= u'(H)= 0, \quad u(H)=1, \end{equation}

which is a boundary-value problem with an adjustable parameter $Re/Fr^2$ in the porous medium depending upon the properties of the medium (porosity and permeability) and the geometry of the flow such that $u(H) = 1$ is satisfied (Samanta et al. Reference Samanta, Goyeau and Ruyer-Quil2013). The anisotropy of the porous medium does not affect the base flow, which is unidirectional.

We use the AUTO07p software (Doedel et al. Reference Doedel, Champneys, Fairgrieve, Kuznetsov, Sandstede and Wang2008) to solve (2.5a–c) numerically by continuation. To start the continuation, we have considered the free-surface flow on a rigid wall ($\varepsilon = 1,\ \kappa _x^{-1} = 0$). We then gradually adjust the porosity and permeability to the desired values. AUTO07p is allowed to adapt the mesh, as it is equipped with a mesh refinement algorithm, so it can easily compute the changes of porosity and permeability at the liquid–porous interface. We have checked for the presence of a sufficient number of mesh points in the transition layer. We have also checked that the thickness of the transitional layer $\varDelta$ is sufficiently low to ensure that the solution is independent of this parameter.

Figure 2 compares the base-flow solution for different values of the Darcy number and for porous-layer thickness $\delta = 0.5$ (figure 2a) and $\delta = 0.1$ (figure 2b). For all the numerical results presented in this work, the porous region is chosen to be relatively thick, with $\delta$ set to $0.5$. As stated above, the transitional layer is thin enough with $\varDelta = 0.001$ (this value is kept constant for the remainder of this work). Even at a low value of the Darcy number ($Da=0.001$), a significant flow is observable at the top of the porous layer, even though the velocity profile in the bulk of the porous medium is flat and nearly negligible. Three regions can be considered: two momentum boundary layers (Brinkman sublayers) at the top and bottom of the porous medium, for which a significant shear is observed; and a Darcy region in between, where the velocity is nearly uniform. For higher values of the Darcy number ($Da=0.01$ and $Da=0.1$), the shear exerted by the liquid flow affects the flow in the porous medium in its entirety and the Darcy sublayer is removed. The extension of the Brinkman sublayers can be determined by balancing the viscous diffusion $\partial _{yy} u/\varepsilon _H$ and the Darcy drag $u/Da$, which gives the estimate $\delta _B = \sqrt {Da/\varepsilon _H}$. When the Brinkman sublayers do not invade the whole porous medium, i.e. $\delta _B < \delta /2$, an approximate solution to the base flow is a liquid film on a solid substrate for which the no-slip condition is displaced at $y= 1 - \delta + \delta _B$, which gives the approximate solution

(2.6)\begin{equation} {\widetilde{U}}_B = \frac{(\delta -\delta_B -y)(y -2H + \delta - \delta_B)}{(H-\delta+\delta_B)^2}, \end{equation}

and therefore

(2.7)\begin{equation} \frac{Re}{Fr^2} ={-} \frac{{\rm d}^2{\widetilde{U}}_B}{{{\rm d} y}^2} \approx \frac{2}{(H - \delta + \delta_B)^2}. \end{equation}

Figure 2. Base-flow velocity profiles $U_B(y)$ and approximated base-flow velocity profiles $\widetilde {U}(y)$ as functions of the cross-stream coordinate $y$ for different permeabilities, when $\varepsilon _H = 0.78$, $\varDelta = 0.001$, and (a) $\delta = 0.5$ and (b) $\delta = 0.1$.

The base-flow velocity, $U_B$, the solution to (2.5a–c), is compared to its approximation, ${\widetilde {U}}_B$, the solution to (2.6), in figure 2(a). An excellent agreement in the liquid region of the flow is observed. Hence, the free-surface instability of a film on a porous substrate can be reduced to the instability of a film on an effective impermeable boundary located at the bottom $y= \delta - \delta _B$ of the top Brinkman layer (Samanta et al. Reference Samanta, Goyeau and Ruyer-Quil2013). In the next section, we will consider the effect of the anisotropy of the porous substrate on the linear stability of the film.

3. Linear stability analysis

We consider a uniform layer of unit thickness $H=1$ and introduce a perturbation expansion of the base-flow solution:

(3.1a,b)\begin{gather} u(x, y, t) = U_B(y) + ub(x, y, t), \quad v(x, y, t) = vb(x, y, t), \end{gather}

(3.1c)\begin{gather} p(x, y, t) = P_B(x) + pb(x, y, t). \end{gather}

We further introduce a streamfunction for the perturbation quantities and a modal decomposition with wavenumber $k=k_r + {\rm i} k_i$ and phase speed $c=c_r + {\rm i} c_i$,

(3.2)\begin{equation} \varPsi(x, y, t) = \psi(y) \exp[{{\rm i}}k(x - ct)], \end{equation}

so that $ub = {\partial \varPsi }/{\partial y }$ and $vb = - {\partial \varPsi }/{\partial x }$ automatically satisfy the continuity equation for the perturbed flow. The Orr–Sommerfeld boundary-value problem then reads

(3.3a)\begin{align} (\mathrm{D}^2-k^2)^2 \psi &= {{\rm i}} k \,Re\left[\left(\frac{U_B}{\varepsilon}-c\right)(\mathrm{D}^2-k^2)\psi -\frac{{U_B^{\prime\prime}} \psi}{\varepsilon}+3{U_B^{\prime}} \frac{\varepsilon' \psi}{\varepsilon^2}+U_B\frac{\psi \varepsilon''}{\varepsilon^2} \right.\nonumber\\ &\quad \left. - \,3U_B\frac{\psi{\varepsilon^{\prime 2}}}{\varepsilon^3}- \left(\frac{U_B}{\varepsilon}-c\right)\frac{\varepsilon'\psi'}{\varepsilon} \right]+\frac{\varepsilon'}{\varepsilon}(\mathrm{D}^3-k^2\mathrm{D})\psi+\frac{\varepsilon \psi''}{\kappa_x(y)}-\frac{\varepsilon k^2 \xi \psi}{\kappa_x(y)}\nonumber\\ &\quad +\varepsilon\psi'\frac{\partial}{\partial y}\left(\frac{1}{\kappa_x(y)}\right)+ \frac{\varepsilon'}{\varepsilon} \psi'''+\frac{\varepsilon''}{\varepsilon}\psi'' -2\frac{{\varepsilon^{\prime 2}}}{\varepsilon^2}\psi'' -k^2\frac{\varepsilon'}{\varepsilon}\psi' , \end{align}

along with the boundary conditions

(3.3b)$$\begin{gather} {\psi\lvert_0} = 0, \quad {\psi'\lvert_0} = 0, \end{gather}$$

(3.3c)$$\begin{gather}{\psi''\lvert_1} + k^2{\psi\lvert_1} ={U_B^{\prime\prime}}\lvert_1\frac{\psi\lvert_1}{({{U_B}\lvert_1} - c)}, \end{gather}$$

(3.3d)\begin{gather} {\psi'''\lvert_1} - 3 k^2{\psi'\lvert_1} = {{\rm i}}\, Re\, k [({{U_B}\lvert_1} - c){\psi'\lvert_1} - {{U_B^{\prime}}\lvert_1} {\psi\lvert_1}] \nonumber\\ \qquad \qquad\qquad\qquad\qquad\qquad\qquad+{{\rm i}}k \left(\cot \theta \,\frac{Re}{Fr^2}+k^2 \,We\,Re\right) \frac{\psi\lvert_1}{(c-U_B\lvert_1)}, \end{gather}

where $\varepsilon (1) = 1$ and $1/\kappa (1) = 0$ have been used to simplify the continuity of the normal stress (3.3d). The symbol $\mathrm {D}$ and primes refer to derivatives with respect to the cross-stream coordinate $y$. The last four terms in (3.3a) can be referred to as the Brinkman second-order correction terms. These terms arise only in the thin transitional region and are negligibly small, so, in our calculation, we ignored these terms (Brinkman Reference Brinkman1949). The above boundary-value problem is actually an eigenvalue problem for the complex phase speed $c$. This eigenvalue problem is solved with AUTO07p for different inclination angles and different anisotropies $\xi$ of the porous plate. We only looked for the surface instability modes, starting our continuation procedure at the trivial solution for $k=0$.

We first consider the spatial stability of a film on a slightly inclined porous plane for a set of parameters corresponding to an experiment by Liu, Schneider & Gollub (Reference Liu, Schneider and Gollub1995), $\theta =4.6^\circ$, $\sigma =69\ {\rm N}\ {\rm m}^{-1}$, $\rho =1130\ {\rm kg}\ {\rm m}^{-3}$ and $\nu =5.02\times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$, which gives $Ka=769.8$. We thus consider a real angular frequency $\omega =kc$, the instability of the film being signalled by a positive spatial growth rate, i.e. $-k_i >0$. Figure 3 presents a typical marginal stability curve ($-k_i= 0$) in the Reynolds number $Re$ versus wavenumber $k_r$ plane. As expected, the threshold of the free-surface instability occurs at $k=0$ for a critical value of the Reynolds number, which depends on the inclination of the plane and the properties of the porous medium. Note that anisotropy does not affect the base flow and appears in (3.3a) only through the product $k^2\xi$. As a consequence, the threshold of the long-wave instability of the film is not affected by the anisotropy of the porous medium and remains equal to its value for the isotropic case $\xi =1$, where the instability threshold closely corresponds to the case of an effective no-slip condition achieved at the bottom of the top Brinkman momentum diffusion layer, i.e. at $y=1-\delta + \delta _B$. In this case, the threshold of instability corresponds to a liquid film of thickness $1-\delta + \delta _B$, which gives, for our set of parameters,

(3.4)\begin{equation} Re_c \approx \frac{5}{4}\,\frac{\cot \theta}{1 - \delta + \delta_B}. \end{equation}

The prediction (3.4) of the instability threshold is indicated by vertical lines on figure 3, in striking agreement with our numerical findings.

Figure 3. Marginal stability curves corresponding to the Orr–Sommerfeld eigenvalue problem for different permeabilities and different anisotropy parameters, when $\varepsilon _H = 0.78$, $\delta = 0.5$, $\varDelta = 0.001$, $\theta = 4.6^{\circ }$ and ${Ka = 769.8}$: (a) $Da = 0.001$ and (b) $Da = 0.01$.

Figure 4 presents the spatial growth rate $- k_i$ for $Re=50$ as a function of $k_r$ (recall that $\omega _i=0$) and shows the influence of the anisotropy parameter. In the isotropic case $\xi =1$, raising the mean permeability ($Da = 0.01$) enhances the range of the unstable wavenumber window and the spatial growth rate. Higher mean permeability corresponds to larger Brinkman layers and thicker effective thickness $1 -\delta + \delta _B$. As a consequence, in the vicinity of the porous–fluid interface, the liquid accelerates and the convective transport in the liquid layer is enhanced, which has a clear destabilising effect.

Figure 4. Spatial growth rate for different permeabilities and different anisotropy parameters, when $\varepsilon _H = 0.78$, $\delta = 0.5$, $\varDelta = 0.001$, $\theta = 4.6^{\circ }$, $Ka = 769.8$ and $Re = 50$: (a) $Da = 0.001$ and (b) $Da = 0.01$.

A low cross-stream permeability ($\xi \gg 1$) affects the stability of the film differently at low ($Da=0.001$) and moderate ($Da=0.01$) permeabilities, as can be observed from figures 3 and 4. The isotropic case ($\xi =1$) is compared to high cross-stream permeability ($\xi =0.001$) and low cross-stream permeability ($\xi =1000$) cases. Anisotropy affects the growth rate and phase velocity only for the larger values of the wavenumber $k_r$, again owing to the fact that $\xi$ is grouped with $k^2$ in a single term of (3.3a). This dual effect is a consequence of the organisation of the flow in the porous layer into three distinct sublayers: the top and bottom Brinkman layers, for which the shear of the base-flow velocity is significant; and the Darcy sublayer, where a plug flow is observed (cf. figure 2). Raising the anisotropy parameter $\xi$ promotes the drag of the solid matrix (proportional to $\xi /Da$), which prevents the exchange of mass between the liquid and porous regions. As a consequence, the flow in the Darcy sublayer is almost a constant plug flow, irrespective of the film-height fluctuations. Conversely, in the top Brinkman sublayer, raising $\xi$ promotes the viscous diffusion of film-height fluctuations as a consequence of the continuity of the streamwise velocity at the liquid–porous interface. The result is a reduction of viscous damping of the waves in the Darcy sublayer, whereas the opposite is true in the Brinkman top sublayer. The net balance of these two adverse effects depends on the respective extensions of the Darcy and Brinkman sublayers. At $Da=0.001$, the Brinkman sublayers are thin ($\delta _B \approx 0.04$) and the destabilising contribution of mitigating streamwise viscous diffusion in the Darcy sublayer dominates. Conversely, at $Da=0.01$, the Darcy sublayer is absent, as can be observed from figure 2 and the damping effect of an enhanced diffusion in the Brinkman top sublayer dominates.

However, a film flowing on a porous substrate with high cross-stream permeability, i.e. $\xi \ll 1$, presents stability that is equivalent to a flow on an isotropic porous medium, as the marginal stability curves, growth rates and phase speeds are identical at $\xi =1$ and $\xi =0.001$. Indeed, lowering the anisotropy parameter $\xi$ below one has no effect on the stability of the film, the viscous damping of the film-height fluctuations in the porous medium then being governed by the streamwise permeability instead of the cross-stream one.

Samanta et al. (Reference Samanta, Goyeau and Ruyer-Quil2013) stated that the stability of a film on a porous media is equivalent to the stability of a liquid film on an impermeable effective boundary located at $y=\delta - \delta _B$, where the Brinkman layer thickness is defined by $\delta _B = \sqrt {Da/\varepsilon _H}$; hereinafter we refer to this model as the ‘no-slip model’.

Figure 5 presents the marginal stability curves when recast in the plane of dimensional flow rate $q_N$ versus dimensional frequency. For a thin porous layer ($\delta = 0.1$), we recover the results obtained by Samanta et al. (Reference Samanta, Goyeau and Ruyer-Quil2013) for the isotropic case ($\xi = 1$). All the curves collapse onto a single curve irrespective of the value of the Darcy number. This curve corresponds to the stability of a liquid film on a solid plate. We conclude that, for thin porous layers, a film flowing on a porous substrate follows the same dynamics as a film on a solid impermeable substrate. However, for a thicker porous region ($\delta = 0.5$), the marginal stability curves present noticeable differences, the permeability of the porous medium significantly promoting the stability of the flow, except close to the instability threshold. Mass exchange at the liquid–porous boundary plays a role in mitigating the instability far from its threshold, as the viscous diffusion in the porous medium can then contribute to the attenuation of velocity fluctuations.

Figure 5. Marginal stability curves for different thicknesses of an anisotropic porous layer for different Darcy values and different $\xi$ values: dimensional cut-off frequency $f_c$ versus dimensional flow rate $q_N$, for $\varepsilon _H = 0.78$, $\varDelta = 0.001$, $\theta = 4.6^{\circ }$ and (a) $\delta = 0.5$, $\xi = 1$, (b) $\delta = 0.5$, $\xi = 1000$, (c) $\delta = 0.1$, $\xi = 1$ and (d) $\delta = 0.1$, ${\xi = 1000}$. ‘OS solid’ refers to an impermeable no-slip boundary.

At a large permeability ($Da=0.1$), the instability occurs at a lower flow rate than for the solid-wall situation. This slight discrepancy in the instability threshold is a consequence of the large difference observed for the velocity profiles (see figure 2). We observe again that lowering the cross-stream permeability $Da/\xi$ has a dual influence: augmenting the range of unstable frequencies for $Da$ below $0.001$, and reducing it above $0.01$. At $Da=0.1$, however, the observed reduction of unstable frequencies is very moderate. A low cross-stream permeability inhibits mass exchange between the porous and liquid layers, decouples the flows in these two layers and prevents the damping of film-height fluctuations by viscous diffusion in the Darcy sublayer. However, due to the continuity of velocity at the porous–liquid interface, the Brinkman top sublayer contributes more efficiently to the damping of the film-height fluctuations. The net balance of these two contributions depends on the relative extensions of the Darcy and Brinkman sublayers. We note that, at $Da=0.1$ and $\xi =1000$, the cross-stream permeability is not low enough to efficiently contribute to the damping of short waves, which is then controlled principally by surface tension and the viscous diffusion in the liquid layer.

4. Weighted residual modelling

In this section, a semi-analytic model is derived to capture the long-wave wavy regime of the film flow. The procedure relies on choosing a set of variables whose evolution mimics the complete dynamics of the flow. The total height ($\text {porous medium} + \text {film}$) $H = h + \delta$ and the local flow rate $q=\int _0^H u\,{{\rm d} y}$ are chosen to characterise the flow evolution. We show in §§ 4.4 and 5 that this strategy is inadequate when a low cross-stream permeability of the porous medium (high values of $\xi$) is considered. As we underline below, this deficiency results from a decoupling between the flow in the porous and liquid regions.

4.1. Three-equation model

We introduce two different flow rates to allow for the velocity distribution in the porous medium to be decoupled from that in the liquid. We define

(4.1a)\begin{equation} q_l = \int_0^H \mathrm{Id}_l u\,{{\rm d} y} \quad \hbox{and} \quad q_p = \int_0^H \mathrm{Id}_p u\,{{\rm d} y}, \end{equation}

with

(4.1b)\begin{equation} \mathrm{Id}_l = 0.5 + 0.5 \tanh[(y-\delta)/\varDelta]\quad\mbox{and}\quad \mathrm{Id}_p = 0.5 - 0.5 \tanh[(y-\delta)/\varDelta], \end{equation}

so that $\mathrm {Id}_l=1$ and $\mathrm {Id}_p=0$ in the liquid region and conversely $\mathrm {Id}_l=0$ and $\mathrm {Id}_p=1$ in the porous medium. Thus, $q_l\approx \int _\delta ^H u\,{{\rm d} y}$ is the local discharge of liquid above the porous–liquid interface $y=\delta$, and $q_p\approx \int _0^\delta u\,{{\rm d} y}$ is the rate of liquid flowing in the porous medium. Since $\mathrm {Id}_l + \mathrm {Id}_p =1$, the total flow rate $q$ is equal to the sum of $q_l$ and $q_p$. Integrating the continuity equation (2.1a) across the fluid and the porous layer, we get the exact mass balance

(4.2)\begin{equation} \partial_t H + \partial_x (q_l+ q_p) = 0. \end{equation}

Two evolution equations can be obtained for $q_l$ and $q_p$ using the weighted residual technique once a velocity profile is assumed (closure assumption). We consider slow space and time evolutions $\partial _{x,t}\sim \gamma$, where $\gamma$ is a formal small parameter that counts the order of derivation with respect to time and space of the different variables. The velocity distribution across the porous and fluid layers is assumed to remain close to the Nusselt flat-film distribution, deviations from the Nusselt profile being $O(\gamma )$ corrections induced by the free-surface deformations. We thus decompose $u$ into

(4.3)\begin{equation} u = q_l(x,t)\,f_l(y; H(x,t)) + q_p(x,t)\,f_p(y; H(x,t)) + \gamma {\tilde{u}}. \end{equation}

We use the linearity of the equations (2.5a–c) satisfied by the base flow to write the velocity profile as the superposition of two solutions to

(4.4a)\begin{equation} L\, f_i ={-} C_{ij} \mathrm{Id}_j, \end{equation}

where $L$ refers again to the operator $\partial _{yy} - (\varepsilon /\kappa _x){\mathrm {Id}}$. For convenience, we have written the above equation in a condensed format using Einstein notation, where $i,j=l$ or $p$ and repetition of indices indicates summation. The associated boundary conditions are

(4.4b)\begin{equation} f_l(0)=f_p(0)=0 \quad \hbox{and}\quad \partial_y\,f_l(H)=\partial_y\,f_p(H)= 0. \end{equation}

The constants $C_{ij}$ are adjusted to satisfy the integral constraints

(4.4c)\begin{equation} \int_0^H \mathrm{Id}_i \,f_j\,{{\rm d} y} = \delta_{ij}, \end{equation}

implied by the definition (4.1) of the variables $q_l$ and $q_p$. Here $\delta _{ij}$ refers to the Kronecker delta function. Constants $C_{ij}$ are therefore dependent on the geometry of the porous medium and therefore are functions of the film height $H$. The decomposition (4.3) is made unique with the two gauge conditions

(4.4d)\begin{equation} \int_0^H \mathrm{Id}_i {\tilde{u}} \,{{\rm d} y} = 0, \end{equation}

which guarantee that the decomposition (4.3) satisfies the definitions (4.1) of the flow rates $q_l$ and $q_p$ in the liquid and porous layers.

The momentum balances (2.1b) and (2.1c) are next simplified within the framework of the long-wave expansion to yield a Prandtl-like equation after elimination of the pressure field. From the continuity equation (2.1a) we get $v = - \int _0^y \partial _x u\,{{\rm d} y} = O(\gamma )$, so that inertial terms in the cross-stream momentum balance equation (2.1c) are $O(\gamma ^2)$ and can be omitted. Integration of (2.1c) thus gives

(4.5)\begin{align} Re \,p &= 2 \partial_y {v\vert_H} - Re \,We\, \partial_{xx} H + \int_H^y \frac{1}{\varepsilon(y)}{L} v \,{{\rm d} y} + (1-\xi) \int_H^y \frac{v}{\kappa} \,{{\rm d} y} \nonumber\\ &\quad - \frac{Re}{Fr^2} \cot\theta (y-H) + O(\gamma^2), \end{align}

where the leading-order contribution from surface tension, $-We \,Re \,\partial _{xx}H$, has been retained, although it is formally $O(\gamma ^2)$, since surface tension prevents the breaking of waves. Substituting (4.5) into (2.1b), we get the boundary-layer equation

(4.6)\begin{align} & Re \left[u_t + \frac{1}{\varepsilon}(u u_x + v u_y) + u v \frac{\partial}{\partial y} \left( \frac{1}{\varepsilon} \right) \right] \nonumber\\ &\quad = \varepsilon \left[ \frac{Re}{Fr^2} (1- \cot\theta \,H_x) + Re\, We \,H_{xxx} \right] - \varepsilon(y) \frac{\partial}{\partial x} \int_H^y \frac{1}{\varepsilon} L v\,{{\rm d} y} \nonumber\\ &\qquad - \varepsilon (1-\xi) \frac{\partial}{\partial x} \int_H^y \frac{v}{\kappa} \,{{\rm d} y} + u_{xx}+ L u + 2\varepsilon \frac{\partial}{\partial x }(u_x \vert_H), \end{align}

with the associated boundary conditions

(4.7a)$$\begin{gather} u\vert_0 = v\vert_0 = 0, \end{gather}$$

(4.7b)$$\begin{gather}\partial_y u\vert_H = 4 \partial_x u\vert_H \partial_x H - \partial_x v\vert_H, \quad v \vert_H = \partial_t H + u\vert_H \partial_x H. \end{gather}$$

Following the weighted residual method, we insert (4.3) into the truncated momentum balance (4.6) and average it with weights that satisfy

(4.8a,b)\begin{equation} L w_i ={-}\frac{Re}{Fr^2}\mathrm{Id}_i, \quad w_i(0) = \partial_y w_i(H)=0. \end{equation}

This choice of weights enables us to make use of the gauge conditions (4.4d) to get rid of the $O(\gamma )$ corrections arising from the evaluation of the viscous terms:

(4.9)\begin{equation} \int_0^H w_i L u \,{{\rm d} y} = w_i(H) \partial_y u|_H + \int_0^H L w_i u \,{{\rm d} y} = w_i(H) [ 4 u_x\vert_H H_x - v_x\vert_H ] - \frac{Re}{Fr^2} q_i. \end{equation}

We finally obtain two averaged momentum balances of the form

(4.10)\begin{align} S_{i}^{(k)} \partial_t q_{i} &={-}F_{ij}^{(k)} \frac{q_j}{H} \partial_x q_i + G_{ij}^{(k)} \frac{q_j q_i}{H^2} \partial_x H + \frac{1}{Fr^2} \left[I^{(k)} b(H) H - \frac{q_k}{H^2}\right]\nonumber\\ &\quad + \frac{1}{Re} \left[ (J_i^{(k)}+\xi {\widetilde{J}}_i^{\,(k)}) \frac{q_i}{H^2} (\partial_x H)^2 - (K_i^{(k)}+\xi {\widetilde{K}}_i^{\,(k)}) \frac{\partial_x q_i \partial_x H}{H} \right. \nonumber\\ &\quad \left.\vphantom{\frac{\partial_x q_i \partial_x H}{H} } - (L_i^{(k)}+\xi {\widetilde{L}}_i^{\,(k)}) \frac{q_i}{H} \partial_{xx} H + (M_i^{(k)}+\xi {\widetilde{M}}_i^{\,(k)}) \partial_{xx} q_i \right], \end{align}

where $b(H) = 1 - \cot \theta \,\partial _x H + We\,Fr^2 \partial _{xxx} H$ combines the gravitational acceleration and the pressure gradient that drive the flow. Equations (4.2) and (4.10) form a three-equation system of evolution equations for the three variables $H$, $q_l$ and $q_p$.

All 54 coefficients involved in (4.10) depend on the geometry and therefore on $H$. The expressions of the coefficients are given in Appendix A. We used the software AUTO07p to compute (4.4), (4.8a,b) and (A1) and tabulated the coefficients (A2) as a function of the film height $H$.

The derived system of equations is consistent up to order $\gamma$ for convective terms and up to $O(\gamma ^2)$ for diffusion terms. We underline that consistency up to $O(\gamma )$ is a requisite to ensure that the instability threshold is captured adequately. Inclusion of viscous diffusion terms enables one to recover correctly the amplitude of the capillary ripples, which precede hump waves and solitary waves and govern their dynamics.

4.2. Two-equation models

The three-equation model (4.2) and (4.10) will be contrasted with the two-equation approach adopted in Samanta et al. (Reference Samanta, Goyeau and Ruyer-Quil2013). Choosing the film elevation $H$ and the total flow rate $q=q_l+q_p$ as the set of variables, a set of evolution equations can be obtained with the above weighted residual method with the ansatz

(4.11)\begin{equation} u = q \,f(y; H) + \gamma {\tilde{u}}, \end{equation}

where $f$ corresponds to the base-flow profile $U_B$ and satisfies

(4.12a–d)\begin{equation} L \,f ={-}C, \quad f(0) =0, \quad \partial_y \,f(H)=0, \quad \int_0^H f \,{{\rm d} y} = 1. \end{equation}

The same procedure as that outlined in § 4.1 is then followed with a weight defined by

(4.13a,b)\begin{equation} L w ={-}\frac{Re}{Fr^2},\quad \partial_y w(H)=0. \end{equation}

Integrating (4.6) with the weight $w$ then leads to

(4.14)\begin{align} S(H) \partial_t q &={-}F(H)\frac{q}{H} \partial_x q +G(H) \frac{q^2}{H^2} \partial_x H + \frac{1}{Fr^2} \left[I(H) b(H) H -\frac{q}{H^2}\right]\nonumber\\ &\quad+\frac{1}{Re} \left[(J(H) + \xi \widetilde{J}(H))\frac{q}{H^2} (\partial_x H)^2 - (K(H) + \xi \widetilde{K}(H)) \frac{\partial_x q \partial_x H}{H} \right.\nonumber\\ &\quad \left.\vphantom{\frac{\partial_x q \partial_x H}{H}} -(L(H) + \xi \widetilde{L}(H))\frac{q}{H} \partial_{xx} H + (M(H) + \xi \widetilde{M}(H)) \partial_{xx} q \right]. \end{align}

The expressions of the coefficients from $S(H)$ to $\widetilde {M}(H)$ correspond to (A1) and (A2) once indices have been dropped out. For an isotropic porous medium $\xi =1$, the averaged momentum balance (4.14) is fully identical to the corresponding one obtained by Samanta et al. (Reference Samanta, Goyeau and Ruyer-Quil2013), as the differences in the expressions of the coefficients arise from a different choice of writing the ansatz (4.11) which does not affect the result.

4.3. Approximation

In order to get a more easily handled model, we propose here a drastic simplification of our modelling attempt. This simplification relies on a polynomial approximation of the velocity profiles $f_i$ and weights $w_i$ ($i$ referring to either $l$ or $p$). To proceed, we separate the flow into three regions: a Brinkman sublayer of thickness $\delta _B=\sqrt {Da/\varepsilon _H}$ at the bottom of the porous layer; a uniform-flow Darcy region in the bulk of the porous layer from $y=\delta _B$ to $y=\delta -\delta _B$; and an effective liquid layer from $y=\delta -\delta _B$ to the top of the film flow $y=H$. We thus approximate $f_i$ and $w_i$ as

(4.15)\begin{equation} \left.\begin{gathered} f_i = a_i y,\quad w_i = a_{wi} y,\quad\hbox{for}\ 0\leqslant y\leqslant \delta_B,\\ f_i = a_i \delta_B,\quad w_i = a_{wi} \delta_B,\quad\hbox{for}\ \delta_B\leqslant y\leqslant \delta - \delta_B,\\ f_i = a_i \delta_B + b_i {\widetilde{U}}_B ,\quad w_i = a_{wi} \delta_B + b_{wi} {\widetilde{U}}_B , \quad\hbox{for}\ y\geqslant \delta - \delta_B, \end{gathered}\right\} \end{equation}

where ${\widetilde {U}}_B$ is given by (2.6). This approximation implies that the operator $L$ reduces to $\partial _{yy}$ in the Brinkman sublayers $\delta - \delta _B\leqslant y \leqslant \delta$ and to $- (\varepsilon _H/Da) \mathrm {Id}$ in the Darcy layer, or bulk, of the porous layer. This approximation also assumes that the transitional layer has a negligible thickness and that the thickness of the Brinkman sublayer $\delta _B$ is small.

Approximation (4.15) yields explicit expressions for the coefficients (A2), of which we retain only the leading-order contributions or up to $O(\delta _B)$. Additionally, the flow rate in the porous layer $q_p$ can be inferred to be proportional to the Darcy number, so that $q_p = O(\delta _B^2)$ is small and the corresponding coefficients are truncated to their leading order with respect to $\delta _B$. Since the velocity distribution in the liquid layer and the top Brinkman sublayer $\delta \leqslant y \leqslant \delta -\delta _B$ is close to the velocity distribution of a liquid film on a solid interface at $y=\delta - \delta _B$, we introduce the effective film thickness $h= H - \delta + \delta _B$ and the effective thickness of the porous medium ${\tilde {\delta }} = \delta - \delta _B$. We obtain explicit expressions for the coefficients, which are listed in Appendix B.

The averaged momentum balances can be further simplified by considering ${q_p=O(Da)}$ and keeping leading-order terms for $Da \ll 1$. We thus retain terms involving the anisotropy parameter $\xi$ up to $O(\xi Da^{3/2})$:

(4.16a)\begin{align} \partial_t q_l - \frac{h}{6{\tilde{\delta}}}\partial_t q_p &=- \frac{17}{7}\frac{q_l}{h}\partial_x q_l + \frac{9}{7}\frac{q_l^2}{h^2} \partial_x h +\frac{1}{Fr^2} \left[ \frac{5h}{6} b(h) - \frac{5(1-{\tilde{\delta}})^2}{4 h^2} q_l \right] \nonumber\\ &\quad + \frac{1}{Re}\left\{ \left[4 -\frac{3}{2}\xi\left(\frac{\delta_B}{h} \right)^3\right] \frac{q_l (\partial_x h)^2}{h^2} - \left[\frac{9}{2} -\xi\left(\frac{\delta_B}{h} \right)^3\right]\frac{\partial_x q_l \partial_x h}{h} \right.\nonumber\\ &\quad - \left[ 6 -\frac{1}{2}\xi\left(\frac{\delta_B}{h} \right)^3 \right]\frac{q_l}{h}\partial_{xx} h + \left[\frac{9}{2} -\frac{1}{4}\xi\left(\frac{\delta_B}{h} \right)^3\right]\partial_{xx} q_l \nonumber\\ &\quad \left. + \, \frac{5}{12}\xi\frac{\delta_B}{h} \partial_{xx} q_p\right\}, \end{align}

and

(4.16b)\begin{align} \partial_t q_l +\partial_t q_p &={-}\frac{12}{5}\frac{q_l}{h}\partial_x q_l + \frac{6}{5}\frac{q_l^2}{h^2} \partial_x h +\frac{1}{Fr^2} \left[ (h + {\tilde{\delta}} \varepsilon_H) b(h) - \frac{(1-{\tilde{\delta}})^2\varepsilon_H}{2 Da} q_p \right] \nonumber\\ &\quad + \frac{1}{Re}\left\{ \left[ \frac{3(2 h + {\tilde{\delta}}\varepsilon_H) }{h} - 3\xi\frac{{\tilde{\delta}}^2}{h^2} \right] \frac{q_l (\partial_x h)^2}{h^2} \right.\nonumber\\ &\quad + \left[-\frac{3(2 h + {\tilde{\delta}}\varepsilon_H)}{h} + 2 \xi \frac{{\tilde{\delta}}^2}{h^2}\right]\frac{\partial_x q_l\partial_x h}{h} + \left[-\frac{3(4h+{\tilde{\delta}}\varepsilon_H)}{2h} + \xi\frac{{\tilde{\delta}}^2}{h^2} \right] \frac{q_l}{h}\partial_{xx} h\nonumber\\ &\quad \left. + \,\frac{1}{2}\left[9 + \frac{3{\tilde{\delta}} \varepsilon_H}{h} - \frac{\xi {\tilde{\delta}}^2}{h^2} \right]\partial_{xx} q_l + \frac{\xi}{3} \frac{{\tilde{\delta}}^2\varepsilon_H}{Da} \partial_{xx} q_p \right\}. \end{align}

These are completed by the mass balance rewritten in terms of $h$, $q$ and $q_p$, i.e.

(4.17)\begin{equation} \partial_t h + \partial_x (q_l + q_p) = 0. \end{equation}

The influence of anisotropy on the momentum balance in the liquid layer (4.16a) is weak, as all terms involving $\xi$ in this equation are functions of the ratio of the extensions of the Brinkman sublayer and liquid layer $\delta _B/h$. These terms originate from the solid-matrix drag occurring in the Brinkman top sublayer. More precisely, these terms arise from the integral

(4.18)\begin{equation} \int_0^H w_l \varepsilon \xi \partial _x \left(\int_H^y \frac{v}{\kappa}\,{{\rm d} y}\right){{\rm d} y} \approx \int_{\delta -\delta_B}^\delta w_l \varepsilon \xi \partial _x \left(\int_H^y \frac{v}{\kappa} \,{{\rm d} y}\right){{\rm d} y}, \end{equation}

since $1/\kappa =0$ in the liquid region and the weight function $w_l$, defined by (4.13a,b), is nearly zero in the Darcy sublayer, and set to zero in the approximation (4.15). Considering that all terms involving $\xi$ in (4.16a) are at most $O(\xi \delta _B^3)$, since $q_p=O(\delta _B^2)$, provides a criterion at which viscous diffusion in the Brinkman sublayer becomes significant, i.e. $\xi Da^{3/2} = O(1)$.

From (4.16b), we can expect an efficient diffusion of the flow rate in the porous medium in the case of a weak cross-stream dimensionless permeability $Da/\xi$. One limit is then worthy of interest, i.e. $q_p\to \mathrm {const.}$, which may be achieved when $\xi /Da\gg 1$.

Considering simultaneously $\xi Da^{3/2}\ll 1$ and $\xi /Da\gg 1$ yields

(4.19)\begin{align} \partial_t q_l &={-} \frac{17}{7}\frac{q_l}{h}\partial_x q_l + \frac{9}{7}\frac{q_l^2}{h^2} \partial_x h +\frac{1}{Fr^2} \left[ \frac{5h}{6} b(h) - \frac{5(1-{\tilde{\delta}})^2}{4 h^2} q_l \right]\nonumber\\ &\quad + \frac{1}{Re}\left\{ 4\frac{q_l (\partial_x h)^2}{h^2} - \frac{9}{2} \frac{\partial_x q_l \partial_x h}{h} -6 \frac{q_l}{h}\partial_{xx} h + \frac{9}{2} \partial_{xx} q_l \right\}, \end{align}

which corresponds to a liquid film flow on an impermeable solid wall at $y={\tilde {\delta }} = \delta - \delta _B$ with a no-slip boundary condition. In this case, the wavy motion of the liquid is unaffected by the porous medium. This is a consequence of neglecting the solid-matrix drag occurring in the Brinkman top sublayer and also neglecting the mass exchange at the liquid–porous interface.

The approximate momentum balance (4.19) can also be derived from (4.16a) for isotropic or moderately anisotropic porous media ($\xi =O(1)$) in the limit $q_p\ll 1$, which is achieved at very low permeabilities ($Da\to 0$) or for thin porous layers (small values of $\delta$).

In the remainder of this paper, we discuss and compare the solutions to the three-equation models (4.2), (4.10), and (4.17), (4.16), and to the two-equation models (4.2), (4.14), and (4.17), (4.19). For convenience, these models are referred to as three-equation full and approximate models, and two-equation full and no-slip models, respectively.

4.4. Linear stability analysis

We next turn to the stability analysis of the flat-film solution based on the three-equation full and approximate models (4.2), (4.10) and (4.17), (4.16), and the two-equation model (4.2), (4.14). We compare our results to the solutions of the Orr–Sommerfeld equations (3.3a) in an attempt to validate our modelling approach.

Figure 6 compares the marginal stability curves of the three- and two-equation models to the Orr–Sommerfeld analysis in the limit of high and low cross-stream permeability (i.e. $\xi = 1$ and $\xi =1000$, respectively). If the three-equation model satisfactorily reproduces the Orr–Sommerfeld analysis, the two-equation full model does not. This inadequacy is a consequence of the ansatz (4.11), which yields an overestimation of the contribution of the cross-stream permeability on the viscous diffusion of the waves (coefficients ${\widetilde {L}}(1)$ and ${\widetilde {M}}(1)$). We conclude that the three-equation models offer reliable alternatives to the Orr–Sommerfeld analysis and are quite accurate. In the isotropic case ($\xi =1$) and for relatively large permeability ($Da=0.01$), the three-equation approximate model is slightly less accurate than the full model, predicting reduced intervals of unstable wavenumbers than the full model and Orr–Sommerfeld analysis. However, this is expected, as the polynomial approximation leading to the approximate model assumes that the base-flow velocity profile presents a constant velocity distribution in the bulk of the porous layer (Darcy sublayer). Such a Darcy sublayer is already not observed at $Da=0.01$ (cf. figure 2).

Figure 6. Comparisons of the marginal stability curves for the two-equation full (4.2), (4.14), three-equation full (4.2), (4.10), and approximate model (4.17), (4.16), with the Orr–Sommerfeld analysis. All the curves are plotted with the parameter set $\varepsilon _H = 0.78$, $\delta = 0.5$, $\varDelta = 0.001$, $Ka = 769.8$ and $\theta = 4.6 ^{\circ }$: (a) $Da = 0.01$, $\xi = 1000$; and (b) $Da = 0.001$, $\xi = 1000$.

5. Nonlinear travelling-wave solutions

We next turn to the construction of the nonlinear TW solutions to the models. These solutions correspond to waves propagating at a constant speed and with a constant shape. The variables are then functions of the coordinate $\zeta =x - ct$ in the moving frame at speed $c$ of the wave. The system of equations then reduces to a set of ordinary differential equations, which are recast into an autonomous dynamical system. The first step is the integration of the mass balance (4.2) with respect to $\zeta$, which gives $q=q_l+q_p= cH + q_0$, where $q_0 =\int _0^H (u-c)\,{\textrm {d} y}$ is the constant flow rate in the moving frame of the wave. Substitution of $cH + q_0$ for $q$ then leads to

(5.1)\begin{equation} \frac{{\rm d}\boldsymbol{U}}{{\rm d}\zeta} = F(\boldsymbol{U}, c), \end{equation}

where $\boldsymbol {U}$ is a three-dimensional vector $(H,H',H'')^\textrm {t}$ for the two-equation models (4.2), (4.14), and (4.17), (4.19), and a four-dimensional vector $(H,H',H'',q_p)^\textrm {t}$ for the three-equation models (4.2), (4.10), and (4.17), (4.16), the primes denoting derivatives with respect to $\zeta$. The locations of the fixed-point solutions to (5.1) satisfy

(5.2)\begin{equation} H^3 \{I^{(l)}(H) + I^{(p)}(H)\} - cH -q_0 =0, \end{equation}

where $I^{(l)}$ and $I^{(p)}$ are defined in (A2m), (B1a) and (B1b). Equation (5.2) admits three solutions, of which only two may be positive. Consequently, a maximum of two fixed points is observed in the phase space for any value of the phase speed $c$.

The moving-frame flow rate $q_0$ is determined by imposing an integral constraint $\langle q \rangle =I^{(l)}(1)+I^{(p)}(1)$, which corresponds to the conservation of mass in a time-dependent simulation with open-domain boundary conditions at the outlet and a periodic forcing at the inlet (Scheid et al. Reference Scheid, Ruyer-Quil, Kalliadasis, Velarde and Zeytounian2005).

The autonomous dynamical system (5.1) constitutes an integral boundary-value problem. Branches of solutions are constructed by varying $c$, which is the principal continuation parameter of the computations. In the context of falling-film flows on an impermeable substrate, several studies (Chang, Demekhin & Kopelevitch Reference Chang, Demekhin and Kopelevitch1993; Sisoev & Shkadov Reference Sisoev and Shkadov1999; Shkadov & Sisoev Reference Shkadov and Sisoev2004) have been devoted to obtaining a detailed picture of the different TW branches of solutions to the famous Shkadov model (Shkadov Reference Shkadov1967). The first family of $\gamma _1$ waves bifurcates from the marginal stability curve through a Hopf bifurcation. These ‘slow’ waves have a speed smaller than that of spatially amplified infinitesimal waves at the same frequency. The corresponding solution branch terminates at small frequency as negative pulses with a deep trough followed by capillary waves. However, the waves experimentally observed at low frequency are ‘fast’ waves, i.e. having speed larger than that of infinitesimal waves, and have a large hump preceded by capillary ripples. This second family of waves ($\gamma _2$) bifurcates from a period-doubling bifurcation of the $\gamma _1$ branch. By considering the symmetry breaking of a wave made of $n$ identical $\gamma _1$ waves, many more branches of solutions, denoted here as $\gamma _n$, have been found. Sisoev & Shkadov (Reference Sisoev and Shkadov1999) showed that the $\gamma _n$ waves, $n \geqslant 3$, bifurcate at a low value of $Re$ in pairs of slow and fast waves that terminate at low wavenumber as solitary waves having different numbers of troughs and bumps. On increasing $Re$, Chang et al. (Reference Chang, Demekhin and Kopelevitch1993) observed a series of pinching of the $\gamma _2$ branch of solutions with the $\gamma _n$ waves.

Different branches of TWs obtained in the conditions of Liu and Gollub's experiments (Liu, Schneider & Gollub Reference Liu, Schneider and Gollub1995) are displayed in figure 7 at increasing values of the Reynolds number for a fixed inclination $\beta =4^\circ$ and Kapitza number $\varGamma =2340$. Among a wide variety of solution branches, we choose to show only those that relate to the branch of fast $\gamma _2$ waves. In particular, only the one-humped $\gamma _2$ fast waves are displayed, whereas $n$-humped fast waves can be found. However, one-humped $\gamma _2$ waves are always the fastest at a given frequency, and therefore dominate the wave dynamics.

Figure 7. Bifurcation diagrams of the TW solutions to the three-equation full model (4.2), (4.10) in the plane of velocity $c$ versus dimensional frequency $f$. The Reynolds number $R = q_N/\nu$ based on the flow rate at the inlet is: (a) $Re = 33.8\ (R = 15)$, (b) $Re = 36\ (R = 16)$, (c) $Re = 42.8\ (R = 19)$ and (d) $Re = 45.5\ (R = 20.2)$. Parameters are $Da = 0.01$, $\varepsilon _H = 0.78$, $\delta = 0.5$, $\varDelta = 0.001$, $\theta = 4^\circ$, $Ka = 2340$ and $\nu = 2.3\times 10^{-6}\ \textrm {m}^2\ \textrm {s}^{-1}$.

Close to criticality ($Re_c \approx 29.1$) and for the isotropic case ($\xi =1$), the $\gamma _2$ family emerges from a Hopf bifurcation whereas the $\gamma _1$ waves appear through a period-doubling bifurcation of the $\gamma _2$ waves. This situation is reversed at higher values of $Re$ (compare figure 7a and 7b). The $\gamma _1$ and $\gamma _2$ families form an imperfect bifurcation with a permutation of branch connections for $33.8 < Re <36$. As the Reynolds number is increased, the fast-wave $\gamma _2$ branch experiences several collisions with other branches. Figure 7(c,d) displays one such pinching, for $42.7< Re <45.5$. Each of these pinching events is reminiscent of the imperfect bifurcation affecting the $\gamma _1$ and $\gamma _2$ wave branches. They give rise to the secondary slow-wave solution branches, denoted by $\gamma _1'$ and $\gamma _1''$ branches in figure 7. The profiles of the waves corresponding to the different branches are illustrated in figure 8 for the frequency $f=1$ Hz. Figures 8(b,d) correspond to slow waves made of one or several troughs followed by capillary oscillations before returning to the level of the flat film. The shape of the high-speed end of the $\gamma _2$ branch (figure 8a) is different, with a main hump preceded by capillary ripples. The $\gamma _1$ waves correspond to a unique trough whereas $\gamma _1'$ and $\gamma _1''$ present several troughs.

Figure 8. Wave profiles of the TW solutions to the three-equation full model (4.2), (4.10) corresponding to the different wave families indicated in figure 7(d): the (a) $\gamma _2$, (b) $\gamma _1$, (c) $\gamma _1'$ and (d) $\gamma _1''$ wave families. Here $\lambda$ corresponds to the period. Parameters are $Da = 0.01$, $\varepsilon _H = 0.78$, $\delta = 0.5$, $\varDelta = 0.001$, $\theta = 4^\circ$, $Ka = 2340$, $Re = 45.5\ (R = 20.2)$, frequency $f=1$ Hz and $\nu = 2.3\times 10^{-6}\ \textrm {m}^2\ \textrm {s}^{-1}$.

Anisotropy significantly affects the sequence of bifurcations and the emergence of TWs. Notably, close to criticality, the slow $\gamma _1$ wave branch was not found. At higher values of $Re$, the sequence of pinching events giving rise to multiple trough-like $\gamma _1^{(n)}$ waves is modified. Raising $\xi$ lowers the wave speed as a result of a more efficient viscous diffusion in the porous medium. It also affects the shape of the waves, as can be observed from figure 8. Large values of $\xi$ correspond to a more efficient viscous diffusion in the Brinkman sublayer, which impacts the amplitude of short capillary waves. The effect of anisotropy on TWs strongly depends on the value of the Darcy number. At $Da=0.001$, the effect is weak, as can be observed from figure 9. This is expected, as low permeabilities imply weak flows in the porous medium.

Figure 9. The same as figure 7 except for $Da = 0.001$: (a) $Re = 42.8\ (R = 19)$ and (b) $Re = 45.5\ (R = 20.2)$.

The approximate model (4.17), (4.16) captures remarkably well the complex bifurcation diagram of TW solutions and predicts correctly the lowering of the wave speed as the anisotropy parameter is raised from $\xi =1$ to $\xi =1000$. In particular, the agreement with the full three-equation model (4.2), (4.10) is convincing even at $Da=0.01$, which is unexpected, as the derivation of the approximate model requires $Da \ll 1$ (compare figure 10 to figures 7 and 9). However, the obtained bifurcation diagrams show some small discrepancies with the result from the full three-equation model. For instance, the multiple-trough $\gamma _1''$ branch connects to a two-hump fast wave branch at $Da=0.01$ and $\xi =1$. Note that multiple-hump fast waves can be obtained with the full three-equation model (4.2), (4.10), but have not been displayed for simplification.

Figure 10. Bifurcation diagrams of the TW solutions to the three-equation approximate model (4.17), (4.16) at $Re=42.8$: (a) $Da=0.01$ and (b) $Da=0.001$. The other parameters are identical to those in figure 7.

Finally, we conclude this section by discussing $\gamma _2$ limit cycles of large extensions, which we refer to as ‘solitary waves’, as these waves tend to organise and dominate the time evolution of noise-driven falling films (Chang et al. Reference Chang, Demekhin and Kopelevitch1993). Figure 11 presents limit cycles of large extensions ($\lambda =200$) which approach homoclinicity. We compared the solutions to the three-equation full model (4.2), (4.10), and approximate model (4.17), (4.16), and to the two-equation full model (4.2), (4.14), and no-slip model (4.17), (4.19). For a small value of the Darcy number ($Da=0.001$), a remarkable agreement (not shown) is achieved between the solutions of the full and approximate three-equation models, as the polynomial approximations for the velocity profiles employed in § 4.3 hold in the limit $Da \ll 1$. As the Darcy number is raised to $Da=0.01$, some departures of the solutions from the approximate model to the full model are noticeable. However, the approximate model remains a convincing substitute for the full model.

Figure 11. Nonlinear TWs for the two-equation full (4.2), (4.14), and three-equation full (4.2), (4.10), and approximate (4.17), (4.16) models, and for the no-slip model (4.17), (4.19) (impermeable wall at $y = \delta -\delta _B$), when $\varepsilon _H = 0.78$, $\delta = 0.5$, $\varDelta = 0.01$, $Ka = 769.8$, $\theta = 4.6 ^{\circ }$ and $\lambda = 200$. Here $H$ indicates the thickness of the entire layer (porous and liquid). (a) $Da = 0.01$, $\xi = 1$; (b) $Da = 0.01$, $\xi = 1$; (c) $Da = 0.01$, $\xi = 1000$; and (d) $Da = 0.01$, $\xi = 1000$.

In the case of an isotropic porous medium ($\xi =1$), the wave profiles obtained with the two-equation full model and the no-slip model nearly collapse onto single curves. This agreement can be comprehended by computing the flow rate in the porous medium ${\tilde {q}}_p=q\int _0^H \mathrm {Id}_p \,f(y;H)\,{\textrm {d} y}$ corresponding to the single velocity profile ansatz (4.11). As can be observed from figure 11(b), ${\tilde {q}}_p$ presents moderate fluctuations so that ${\tilde {q}}_p$ can be assumed somewhat constant. Since $q_p=\mathrm {const.}$ or $q_p \ll 1$ are the conditions yielding the no-slip model, this explains the observed concordance of results from the no-slip and two-equation full model at $\xi =1$. At $\xi =1000$, solitary-wave solutions to the two-equation full model present capillary ripples whose number and amplitude are significantly lower than with the other models, to the point that, for $Da=0.01$, only one capillary ripple is observable at the front of the wave (see figure 11c). This discrepancy results from the overestimation of the contributions of the porous medium to the viscous diffusion terms in the momentum balance (4.14) implied by the one-velocity ansatz (4.11).

A comparison of the solutions of the three-equation models at $\xi =1$ and $\xi =1000$ shows that the fluctuations of the flow rate $q_p$ in the porous medium are efficiently attenuated by lowering the cross-stream permeability (thus raising $\xi$). This explains the convergence of the solutions to the three-equation and no-slip models at $\xi =1000$. Finally, let us note that the intensity of the mass exchange between the porous and liquid regions, as reflected by the fluctuations of $q_p$, is predicted to be much higher with the three-equation models than with the two-equation full model in the isotropic case ($\xi =1$). We conclude that the assumption of a complete slaving of the velocity distribution in the porous medium to that in the liquid medium, i.e. ansatz (4.11), as employed by Samanta et al. (Reference Samanta, Goyeau and Ruyer-Quil2013), is too restrictive.

6. Time-dependent simulations

In this section, we present some time-dependent numerical simulations of the three-equation full and approximate models (4.2), (4.10), and (4.17), (4.16) and the no-slip model (4.17), (4.19) with outlet open boundary conditions and forcing at the inlet (open flow). Computations have been performed using finite differences and the method of lines, which can be easily employed in the case of evolution equations (Cellier & Ruyer-Quil Reference Cellier and Ruyer-Quil2019), for which the system of partial differential equations can be written as

(6.1)\begin{equation} \boldsymbol{U}_t = F(\boldsymbol{U}, \boldsymbol{U}_x, \boldsymbol{U}_{xx}, \ldots), \end{equation}

where the indices refer to temporal or spatial derivatives; and $\boldsymbol {U}$ is the vector that contains the unknowns (in this particular case $h$, $q$ and $q_p$). The spatial derivatives of our system (6.1) are first discretised and replaced by algebraic approximations using second-order central finite differences. As a consequence, (6.1) is transformed into an initial-value problem for the vector $\widetilde {\boldsymbol {U}}$ of discrete unknowns on the chosen numerical grid:

(6.2)\begin{equation} \frac{{\rm d}\widetilde{\boldsymbol{U}}}{{\rm d}t} = {\widetilde{F}}(\widetilde{\boldsymbol{U}}). \end{equation}

An in-house solver has been developed in Julia language in order to integrate (6.2). This solver uses the libraries ForwardDiff (Revels, Lubin & Papamarkou Reference Revels, Lubin and Papamarkou2016) and DifferentialEquations (Rackauckas & Nie Reference Rackauckas and Nie2017). The ForwardDiff library implements forward-mode automatic differentiation algorithms and thus enables one to automatically and efficiently compute the Jacobian of the function ${\widetilde {F}}$. As a result, changes of the model are easily implemented in the code. DifferentialEquations offers a wide range of schemes and algorithms for integration in time. We used a Crank–Nicolson scheme with adaptive time-stepping. This scheme is recommended for highly stiff problems and limits the numerical cost.

We present some simulations of a time-periodic wavy evolution generated by a regular forcing at the inlet:

(6.3)\begin{equation} q_i = I^{(i)}(1) [1 + A \sin(\omega t)]. \end{equation}

Outlet open boundary conditions are selected in order to reproduce experimental conditions. They consist of simple no-flux conditions at the domain outlet. Owing to the intrinsic convective nature of the studied phenomena, induced oscillations only affect a few nodes of the numerical domain at the outlet. Parameters are chosen to correspond to an experiment by Liu & Gollub (Reference Liu and Gollub1994) with a dimensionless flow rate $R=q_N/\nu =19.33$, a forcing frequency $f=4.5$ Hz, angle $\theta = 6.4^\circ$, $\nu =6.27\times 10^{-6}\ \textrm {m}^2\ \textrm {s}^{-1}$ and $Ka=526$. We chose $\delta =0.5$ and a moderate Darcy number $Da=0.01$. The Reynolds number has been adjusted to $Re={3R}/[2(1 -{\tilde {\delta }})]=47.3$, which corresponds to a film on an impermeable substrate at $y = {\tilde {\delta }}$ (no-slip model). The no-slip model (4.17), (4.19) reproduces accurately the nonlinear multi-peaked wave evolution reported in the experiment (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville2002). As it develops downstream, the flow remains periodic in time at the inlet forcing frequency. The initial exponential growth gives way to a multi-peaked wave with the growth of a secondary peak. A phase locking occurs and a modulated wave-train emerges. Our simulations with the full and approximate three-equation models (4.2), (4.10), and (4.17), (4.16) in the isotropic case $\xi =1$ underline the effect of a thick porous layer on the wave evolution. In contrast to the remarkable agreement reported by Samanta et al. (Reference Samanta, Goyeau and Ruyer-Quil2013) for a thin porous medium at $\delta =0.1$, our simulations show that the flow in the porous medium significantly affects the wave evolution, mitigating the growth of a secondary peak, delaying the phase locking and the emergence of the modulated wave-train. The approximate and full models give results in close agreement.

In contrast to the no-slip model, the simulations performed with the full and approximate models at $\xi =1000$ show the emergence of waves with a lower amplitude than for the isotropic case ($\xi =1$). A quicker onset of the final modulated wave-train is observed. Spatial modulations of this wave-train are also weaker and quickly attenuated. This is clearly an effect of the viscous diffusion of the liquid momentum in the top Brinkman sublayer of the porous medium, which is promoted at low cross-stream permeabilities. In fact, we have checked that setting coefficients $\widetilde {J}^{\,l}_i$, $\widetilde {K}^{\,l}_i$, $\widetilde {L}^{\,l}_i$ and $\widetilde {M}^{l}_i$ to zero in the averaged momentum balance (4.10) for the liquid layer yields results in close agreement with the no-slip model at $\xi =1000$. Similar results are obtained for the approximate model by dropping all terms proportional to $\xi$ in (4.16a).

We next turn to a simulation of a solitary-like wave-train generated at a lower frequency $f=1.5$ Hz, for which we discuss the effect of varying the Darcy number (see figure 13). At a moderate value of the Darcy number ($Da=0.01$), raising anisotropy leads to a damping of the capillary ripples, which is consistent with the above observations for $f=4.5$ Hz. However, lowering $Da$ to $Da=0.001$, anisotropy has an opposite effect, the amplitude of capillary ripples being raised as $\xi$ is increased. These observations are consistent with our linear stability analysis (§ 3). Raising the cross-stream permeability affects the short capillary waves differently depending on the relative extension of the Darcy and Brinkman sublayer in the porous region. Results from the approximate models (not shown) are in close agreement with the full model. We note that $\xi Da^{3/2}$ equals $0.05$ for $\xi =1000$ and $Da=0.001$, and $1.5$ for $\xi =1000$ and $Da=0.01$, so that the contribution of cross-stream permeability to the viscous diffusion in the Brinkman top sublayer is weak in the former case and significant in the latter one. As expected, a close agreement is found between the solution to the no-slip model and the full three-equation model for $Da=0.001$ and $\xi =1000$. This agreement results again from the efficient damping of the mass exchange between the liquid and porous layer, so that $q_p$ is nearly constant at large values of $\xi$.

Figure 12. Simulations of a wave-train generated by a periodic inlet forcing using the no-slip model (4.17), (4.19), and the three-equation full and approximate models (4.2), (4.10), and (4.17), (4.16). Parameters correspond to an experiment conducted by Liu & Gollub (Reference Liu and Gollub1994): $Re=47.3$, $f=4.5$ Hz, $\theta = 6.4^\circ$, ${\nu =6.27}\times 10^{-6}\ \textrm {m}^2\ \textrm {s}^{-1}$ and $Ka=526$.

Figure 13. Simulations of a wave-train generated by a periodic inlet forcing using the three-equation full model (4.2), (4.10) and the no-slip model (4.17), (4.19) for different Darcy numbers. Parameters correspond to an experiment conducted by Liu & Gollub (Reference Liu and Gollub1994): $Re=47.3$, $f=1.5$ Hz, $\theta = 6.4^\circ$, $\nu =6.27\times 10^{-6}\ \textrm {m}^2\ \textrm {s}^{-1}$ and $Ka=526$. (a) $Da = 0.001$ and (b) $Da = 0.01$.

Finally, we present time-dependent simulations of the three-equation full model and the no-slip model in an extended domain. In order to simulate the natural evolution of the film as a response to ambient noise, we superimpose to a constant inlet flow rate a white-noise signal of weak amplitude, such that the inlet film flow rates $q_i(x=0,t=0)= I^{(i)}g(t)$, where $g(t)$ fluctuates in the interval $[0.9, 1.1]$. Figure 14 presents a snapshot of the film evolution at the end of the simulation. The typical coarsening dynamics of falling-film flows is clearly observable, as large-amplitude waves travel faster and capture the waves that precede them. The downstream dynamics of the flow is then organised by solitary waves in interaction (Chang & Demekhin Reference Chang and Demekhin2002). In the isotropic case ($\xi =1$), the flow rate $q_p$ in the porous medium may rise to quite large values under the large-amplitude solitary waves, corresponding to roughly three times its value for the Nusselt uniform-film solution. We thus conclude that the wavy dynamics of the film triggers an efficient exchange of mass between the porous and liquid media. Raising the cross-stream permeability ($\xi =1000$) has a strong influence on the wave dynamics. Close to the inlet, nonlinear waves emerge more quickly out of the growth of random perturbations. Solitary waves observed at the outlet of the numerical domain have a lower amplitude.

Figure 14. Simulations of the response of the film to noise at the inlet based on the three-equation full model (4.2), (4.10) and the no-slip model (4.17), (4.19). Panels ($a$) and ($b$) are snapshots of the free-surface elevation and flow rates in the liquid and porous region at the end of the simulation. Parameters are $Da = 0.01$, ${\varepsilon _H = 0.78}$, $\delta = 0.5$, $\varDelta = 0.001$, $Ka = 769.8$ and $\theta = 4.6^{\circ }$.

A mitigation of the capillary ripples preceding the solitary waves is also evident, as well as the onset of bound states made of several linked waves and multi-peaked waves. All these phenomena arise from the influence of viscous diffusion on the velocity of the waves, an effect sometimes referred to as ‘viscous dispersion’ (Pradas, Tseluiko & Kalliadasis Reference Pradas, Tseluiko and Kalliadasis2011). In contrast, the no-slip model (4.17), (4.19) presents interacting solitary waves whose shapes are similar to the isotropic case ($\xi =1$), whereas the simulated wave evolution at the inlet is closer to the anisotropic one ($\xi =1000$). This can be comprehended by considering that raising $\xi$ attenuates the fluctuations of $q_p$ and the no-slip model corresponds to $q_p=\mathrm {const.}$ Raising $\xi$ also increases the viscous dispersion of the waves prompted by the Brinkman sublayer, an effect that is not accounted for by the no-slip model.

7. Summary and conclusions

The effects of anisotropy on the stability and nonlinear evolution of a falling film on an inclined plane have been investigated in the framework of a one-domain composite description of the entire liquid and porous domains introduced by Beckermann et al. (Reference Beckermann, Viskanta and Ramadhyani1988).

Extending the weighted residual modelling approach, a three-equation semi-analytical model (4.2), (4.10) has been derived by decoupling the flows in the porous and liquid regions. Two separate unknowns have thus been introduced to characterise the flow in these two regions, i.e. the flow rates, $q_l$ and $q_p$. The coefficients in the three-equation model have been numerically computed and tabulated. In the limit $Da\ll 1$, we have shown that this model can be simplified to (4.17), (4.16) using polynomial approximations of the velocity profile.

Orr–Sommerfeld stability analysis of the base flow has revealed a non-trivial dual effect of anisotropy on the film stability depending on the permeability of the porous medium. Raising $\xi$ lowers the exchange of mass between the liquid and porous regions and therefore reduces the viscous damping of the film-height fluctuations in the Darcy sublayer at the core of the porous region. However, viscous diffusion of the film-height fluctuations is enhanced at the top Brinkman sublayer. The net balance of these two mechanisms depends on the relative size of the Brinkman and Darcy sublayers, therefore on the Darcy number.

In the isotropic case, when the thickness of the porous medium is relatively small (cf. figure 5), we observed the linear stability of the film on the porous substrate to be nearly equivalent to the stability of a film on an impermeable no-slip effective boundary located at $y = \delta - \delta _B$ (no-slip model). For a thicker porous layer, our computations show a weak amplification of the instability at the threshold at large permeabilities and an attenuation otherwise. This slight discrepancy on the instability threshold probably results from the significant differences observed for the velocity profiles in the porous region at large permeabilities (see figure 2). The attenuation of the instability further from threshold as compared to the no-slip model results from the viscous damping of the velocity fluctuations in the porous medium, which is allowed by the exchange of mass at the porous–liquid interface.

Our study of nonlinear waves has been based on the weighted residual three-equation models that we derived to account for the decoupling of the flows in the porous and liquid regions. This approach has been validated in the linear regime by comparison with the Orr–Sommerfeld analysis. The approximate model (4.17), (4.16) offers a welcome simplification of the full model (4.2), (4.10), which remains valid even at moderate values of the Darcy number. It also enables us to propose a criterion, i.e. $\xi \,Da^{3/2}=O(1)$, above which the viscous diffusion in the Brinkman top sublayer becomes significant.

We have constructed nonlinear TW solutions by a continuation method and performed time-dependent simulations with open boundary conditions at the outlet (open flow). Raising $\xi$ affects the wave shape and speed, modifying the connections of TW branches. TW solutions are close to the corresponding ones for the no-slip model (effective impermeable boundary at $y=\delta -\delta _B)$ at low permeabilities. However, for the isotropic case $\xi =1$, the intensity of the fluctuations of the flow rate $q_p$ in the porous medium, and therefore the exchange of mass at the liquid–porous interface, is significant and impacts the amplitude of the capillary ripples.

Time-dependent simulations of the film response to a periodic forcing at the inlet have revealed a sharp effect of anisotropy on the time-periodic wave-trains. Raising $\xi$ attenuates the fluctuations of $q_p$, thus limiting the mass exchange at the liquid–porous interface and mitigating the damping of film-height fluctuations in the Darcy sublayer. This mechanism amplifies the capillary ripples at low values of $Da$ (cf. figure 13). However, damping of the capillary ripples is instead observed at larger values of $Da$ when $\xi$ is raised. This is explained by the increased efficiency of the viscous diffusion of the waves, as the intensity of the flow in the Brinkman sublayer is raised along with the permeability. Raising the anisotropy parameter $\xi$ at moderate and large permeabilities also has a regularising effect on the inception of primary wave-trains generated by an inlet forcing (cf. figure 12), lowering the wave amplitude and damping the spatial modulations of the wave-trains. However, these effects on the stability and nonlinear dynamics of the liquid film are observed for fairly large values of the anisotropy parameter ($\xi =10^3$). If very anisotropic porous media can be found in Nature, for instance, for faulty rocks where $\xi$ can be as high as $10^4$ (Evans, Forster & Goddard Reference Evans, Forster and Goddard1997; Farrell, Healy & Taylor Reference Farrell, Healy and Taylor2014), it is difficult to think of a porous medium at the scale of a liquid film exhibiting such large anisotropy.

Finally, our time-dependent simulations of a noise-driven wavy film show that the fluctuations of the flow rate in the porous medium can be large due to the onset of large-amplitude solitary waves that dominate the dynamics of the film. Raising $\xi$ at $Da=0.01$ has a complex effect on the wave dynamics, modifying wave-to-wave interactions as a result of an enhanced dispersion of the waves prompted by the viscous diffusion in the Brinkman sublayer. We conclude that the influence of the flow in the porous layer on the film dynamics is much stronger than admitted in previous studies. The approximate three-equation model (4.17), (4.16) offers a reliable alternative to the more complex semi-analytical model (4.2), (4.10) and thus constitutes an interesting tool for later work.