2.1. Surface tension dominated regime

Inspection of the normal stress balance condition (2.16) shows that a small capillary number $D=\delta ^3\overline {D}$ (here and in what follows bars denote order one quantities) induces a pressure in the film of order $1/\delta$ . To retain the external flow effects and recalling the expansion (2.21), we conclude that $C=\delta ^{-2}\overline {C}$ . The large capillary pressure drives an order one flow in the film and the appropriate expansions are $u=u_0+O(\delta )$ , $v=v_0+O(\delta )$ , $p={\delta ^{-1}}p_0+O(1)$ , $h=h_0+O(\delta )$ . To retain gravitational effects we need small inclination angles $\alpha =O(\delta )$ and write $\cot \alpha =G/\delta$ where $G$ is a constant. Substitution into (2.10)–(2.12) gives

(2.26) \begin{align} -p_{0x}+\frac {1}{R}u_{0yy}-\frac {2}{R} &=0, \end{align}

(2.27) \begin{align} p_{0y}+\frac {2}{R}G &=0, \end{align}

(2.28) \begin{align} u_{0x}+v_{0y} &=0. \\[9pt] \nonumber \end{align}

The no-slip conditions (2.13), (2.14) are unchanged besides dependent variables gaining a subscript $0$ . To retain shear stresses in (2.15) we need $C C_D=O(1)$ since the left-hand side is of $O(1)$ . Hence, we write $CC_D=\overline {\tau }$ with $\overline {\tau }$ an order one constant. The kinematic condition (2.14) and the stress balances (2.15)–(2.16) give to leading order

(2.29) \begin{align} v_0=h_{0t}+u_0 h_{0x},\quad \text{at}\quad y=h_0, \end{align}

(2.30) \begin{align} u_{0y}=\overline {\tau },\quad \text{at}\quad y=h_0, \end{align}

(2.31) \begin{align} -\overline {C}{\mathcal{H}}[h_{0x}]-\frac {1}{2}R p_0=\frac {1}{2\overline {D}}h_{0xx},\quad \text{at}\quad y=h_0, \end{align}

where the solution (2.25) has been used. A single equation for $h_0(x,t)$ is obtained as follows. Integration of (2.27) gives $p_0=\overline {p}_0(x,t)-2Gy/R$ , while integration of (2.26) with respect to $y$ , use of no slip $u_0(x,0,t)=0$ and the tangential stress balance (2.30) gives $u_0=y(y-2h_0) [({1}/{2})R\overline {p}_{0x}+1 ]+\overline {\tau } h_0$ . Mass conservation (2.28) immediately yields $v_{0}$ after the condition $v_0(x,0,t)=0$ is used. The normal stress balance (2.31) gives $R\overline {p}_0=-2\overline {C}{\mathcal{H}}[h_{0x}]+2Gh_0- ({1}/{\overline {D}})h_{0xx}$ . Using this in the solutions for $u_0$ and $v_0$ and substituting the velocities into (2.29) provides the desired evolution equation

(2.32) \begin{equation} h_{0t}+\left [-\dfrac {2h_0^3}{3}+\dfrac {\overline {\tau }}{2}h_0^2-\dfrac {2G}{3}h_0^3h_{0x}+\dfrac {\overline {2C}}{3}{\mathcal{H}}[{h_{0xx}}]h_0^3+\dfrac {1}{{3\overline {D}}}h_{0xxx}h_0^3\right ]_{x} = 0. \end{equation}

Equation (2.32) is solved on periodic domains of dimensionless size $L$ , i.e. $h(x,t)=h(x+L,t)$ . The size $L$ is a crucial dimensionless parameter with the number of linearly unstable modes increasing with $L$ . In order to compare with the results of Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024) we rescale to spatial domains of unit length. This is achieved by writing $x=L\bar {x}$ , $h=(\overline {\tau }/2)\bar {h}$ , and $t=(12 L/\overline {\tau }^2)\bar {t}$ , to cast (2.32) into (after dropping bars)

(2.33) \begin{equation} h_{t}+2\left [\left (\dfrac {\nu ^2\pi }{d_1}{\mathcal{H}}[h_{xx}]+\dfrac {\nu ^3}{d_2}h_{xxx}-\nu {h_{x}}-1\right )h^3+\dfrac {3}{2}h^2\right ]_{x}=0, \end{equation}

where

(2.34) \begin{equation} d_1=\dfrac {\pi G^2\overline {\tau }}{\overline {C}},\qquad d_2=\dfrac {\overline {D}\overline {\tau }^2 G^3}{2},\qquad \nu =\dfrac {G\overline {\tau }}{2L}, \end{equation}

where $d_{1}$ measures the induced drag, $d_2$ is a scaled Bond number (small $d_2$ implies large surface tension) and $\nu$ is inversely proportional to the size of the system (small $\nu$ corresponds to large dimensionless wavelengths). Equation (2.33) is identical to the problem studied in Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024) (a factor of $2$ is removable by the additional rescaling $\overline {\tau }\to 2\overline {\tau }$ which is due to a factor of 2 difference in the non-dimensionalisation of the shear stress in the present study), and is to be solved on spatially periodic domains of unit length. Version (2.33) is very useful because it enables us to directly study the stability and dynamics of the nonlinear bifurcation branches computed in Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024). This is undertaken later.

2.2. The inertio-capillary regime

This regime arises when the capillary number is not as small as in § 2.1 so that pressure perturbations due to surface tension and the air stream are ${\mathcal{O}}(1)$ ; this happens when $D=\delta ^2\overline {D}$ in (2.16), the angle $\cot \alpha ={\mathcal{O}}(1)$ , $C=\delta ^{-1}\overline {C}$ and $R$ , $\overline {\tau }$ being order one constants. The leading-order pressure $P_{a0}$ in the air region is given by (2.25). Expanding variables in the liquid film, $u=u_0+\delta u_1+\cdots$ , $v=v_0+\delta v_1+\cdots$ , $p=p_0+\delta p_1+\cdots$ , $h=h_0+\delta h_1+\cdots$ and substituting into (2.10)–(2.12) gives to leading order

(2.35) \begin{align} u_{0yy}=2,\qquad p_{0y}=-\frac {2}{R}\cot \alpha ,\qquad u_{0x}+v_{0y}=0, \end{align}

which gives, using the no-slip and tangential shear-stress conditions

(2.36) \begin{align} u_0=y(y-2h_0)+\overline {\tau }y,\qquad v_0=h_{0x}y^2,\qquad p_0=-\frac {2\cot \alpha }{R}y+\overline {p}_0(x,t). \end{align}

The pressure $\overline {p}_0$ follows from the normal stress balance condition (2.16) at leading order

(2.37) \begin{align} \frac {1}{2}R\overline {p}_0(x,t)=\cot \alpha \,h_0-\overline {C}{\mathcal{H}}[h_{0x}]-\frac {1}{2\overline {D}}h_{0xx}. \end{align}

The kinematic condition (2.14) gives, to leading order

(2.38) \begin{align} h_{0t}+\left (\frac {\overline {\tau }}{2}h_0^2-\frac {2}{3}h_0^3\right )_x=0. \end{align}

Similar leading-order behaviour was found in Tseluiko & Papageorgiou (Reference Tseluiko and Papageorgiou2006). Solutions of (2.38) for fairly general initial conditions (e.g. smooth spatially periodic initial data) encounter infinite slopes in finite time heralding the breakdown of the long-wave asymptotic expansion. The singularity can be regularised by incorporating high-order terms including surface tension, as done by Tseluiko & Papageorgiou (Reference Tseluiko and Papageorgiou2006) in a related problem. The analysis leading to the desired evolution equation (2.39) below, is provided in Appendix A

(2.39) \begin{align} h_{t} &+ \left [-\frac {2}{3}h^3+\frac {\overline {\tau }}{2}h^2+\delta \left (\frac {8R}{15}h^6-\frac {4R\overline {\tau }}{15}h^5-\frac {2\cot \alpha }{3}{h^3}\right )h_{x}+\frac {2\delta \overline {C}}{3}h^3{\mathcal{H}}[h_{xx}] \right. \nonumber\\&\left. +\, \frac {\delta }{3\overline {D}}h^3h_{xxx}\right ]_{x}=0. \end{align}

Equation (2.39) agrees with that in Tseluiko & Papageorgiou (Reference Tseluiko and Papageorgiou2006) in the absence of the air stream here and no electric field in that study. Equation (2.39) is solved on $L\hbox{-}$ periodic domains, i.e. $h(x,t)=h(x+L,t)$ . In the computations presented in § 5 we take $L=1$ and also increase it to $L=2\pi$ in order to evaluate the effect of domain size on the solution.

3.1. Numerical methods and validation

To fix matters, we will describe our numerical methods for the inertio-capillary (2.39). Totally analogous schemes are implemented to solve (2.33). Either equation is written as

(3.1) \begin{equation} h_{t}+{\mathcal{A}}(h) = {\mathcal{B}}(h), \end{equation}

where the nonlinear operator ${\mathcal{A}}(h)$ contains terms involving $h^3 h_{xxxx}, h^3 {\mathcal{H}}[h_{xxx}]$ and $h^3 h_{xx}$ , and all remaining terms are grouped into ${\mathcal{B}}(h)$ . For (2.39) these are

(3.2) \begin{align} {\mathcal{A}}(h) =& \left (A_{1}{\mathcal{H}}[{h_{xxx}}]+A_{2}h_{xxxx}-A_{3}h_{xx}\right )h^3, \end{align}

(3.3) \begin{align} {\mathcal{B}}(h) =& \left [B_{1}h^3+B_{2}h^2+\delta {B_{3}}\left (\dfrac {8R}{15}h^6-\dfrac {4R}{15}CC_{D}h^5\right )h_{x}\right ]_{x}\nonumber \\ &+3h^2h_{x}\left (A_{1}{\mathcal{H}}[{h_{xx}}]+A_{2}h_{xxx}-A_{3}h_{x}\right )\!, \\[9pt] \nonumber \end{align}

and analogous expressions with different coefficients $A_{i}$ and $B_{i}$ arise for (2.33). In our numerical schemes the function ${\mathcal{A}}(h)$ is treated implicitly whilst the part ${\mathcal{B}}(h)$ is evaluated explicitly. Spatially periodic domains $[0, L]$ are considered (all computations in the present work use $L=1$ unless stated otherwise). On implementing the schemes, the spatial domain $[0,L]$ is split into $N$ (odd) equidistant points and the time step is defined by ${\textrm{d}}{t}$ ; typically we take ${\textrm{d}}{t}\approx 10^{-4}$ or smaller to guarantee the conservation of the mean film thickness $\int _0^L h(x,t){\textrm d}x$ as dictated by the equation. Second-order backward differentiation formulas (BDF) are used for the time discretisation (see Akrivis et al. (Reference Akrivis, Papageorgiou and Smyrlis2011) for all BDF formulas up to sixth-order accuracy)

(3.4) \begin{equation} \dfrac {3}{2}h^{n+2}+{\textrm{d}}{t}{\mathcal{A}}(h^{n+2})=2h^{n+1}+2{\textrm{d}}t{\mathcal{B}}(h^{n+1})-\dfrac {1}{2}h^{n}-{\textrm{d}}t{\mathcal{B}}(h^{n}). \end{equation}

To start the time integration (3.4), we need $h^{1}$ starting from the initial condition $h^{0}$ . This is achieved by one step of the implicit Euler method. After computation of $h^{n}$ and $h^{n+1}$ , the implicit time stepping requires solution of a nonlinear system to obtain $h^{n+2}$ . There are $N$ unknowns, $h^{n+2}_{j}$ , at the collocation points $j = 1,2,3 \cdots N$ . Newton’s method is applied to solve the resultant nonlinear algebraic systems at the collocation points. The Jacobian matrix is computed analytically prior to each iteration by taking advantage of spectral matrices. Iterations are terminated when the error is less than $10^{-8}$ , thus guaranteeing accurate results. The square matrix of the discrete Hilbert transform, denoted by $\mathbb{H}$ , is given by

(3.5) \begin{equation} \mathbb{H}(i,j) = \begin{cases} \frac {1-(-1)^{i-j}}{2N}\cot \left (\frac {\pi ({i-j})}{2N}\right )-\frac {1+(-1)^{i-j}}{2N}\tan \left (\frac {\pi ({i-j})}{2N}\right )\!, \quad &\mathrm{mod}(i-j,N)\not \equiv {0}\\ 0, \quad & \mathrm{mod}(i-j,N) \equiv 0,\\ \end{cases} \end{equation}

and the first derivative matrix is given by

(3.6) \begin{equation} \mathbb{D}(i,j) = \begin{cases} 0, \quad & \mathrm{mod}(i-j,2) \equiv 0,\\ \frac {(-1)^{(i-j)}\pi }{L}\csc \left (\frac {\pi (i-j)}{N}\right )\!, \quad &\mathrm{mod}(i-j,2) \not \equiv 0, \end{cases} \end{equation}

where $L$ is the spatial period mentioned above (see Trefethen (Reference Trefethen2000) for detailed derivations). For the inertia-less model (2.33) we take $L=1$ and vary $\nu$ to study waves of different unscaled wavelengths, while for the inertio-capillary regime, the effect of computation domain sizes is directly examined by changing $L$ .

To verify the numerical scheme, we perform the time integration for models with and without inertia, and numerical parameters $N = 501$ and ${\textrm{d}}t = 10^{-4}$ . In the first numerical experiment, we select a small amplitude initial condition $h = 1+0.02\cos (2\pi {x})$ , so that the growth rate at an early stage can be predicted by the linear theory in Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024). Figure 2(a) shows the evolution of the wave amplitude (solid blue curve) superimposed with the prediction of linear theory; agreement is seen to be excellent prior to nonlinear saturation effects that enter around $t\approx 0.8$ . For the second computation, we compute a steady solution of (2.39) and use it as the initial condition in our code. Figure 2(b) shows the time evolution and as expected the solution remains unchanged (to within truncation errors) for a relatively large computed time $t=40$ , verifying once again the accuracy of our algorithms. To classify different types of unsteady solutions as parameters vary we follow the energy $E=\int _{0}^{1} h^2(x,t){\textrm{d}}x$ of the solutions The maxima and minima of $E(t)$ are tracked to construct numerically Poincaré maps in order to probe the complexity of the dynamics – see for example Smyrlis & Papageorgiou (Reference Smyrlis and Papageorgiou1991).

Figure 2. (a) Comparison of the growth rate for $d_1 = 0.8$ and $d_2 = 0.35$ and $\nu = 0.1$ . Here, $A$ is the wave amplitude. The numerical results is based on (2.33). (b) Comparison of transient solutions of (2.39) at $t = 0$ and $t = 40$ for $C = 23$ , $D = 0.01$ , $C_{D} = 0.1$ , $R = 2.5$ , $\alpha = \pi /8$ and $\delta = 0.3365$ .

4.1. Stability of non-uniform steady states

In Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024) we computed steady periodic wave solutions to (2.33) (equivalently (2.32)) and analysed them to develop insight into the rich bifurcation diagrams at play. Unless otherwise stated, all computations presented in this section pertain to (2.33) in order to facilitate direct comparisons with Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024). The solutions that branch off from a flat interface as infinitesimal waves (linear waves) form the primary branch, whereas solutions that bifurcate from other finite-amplitude branches are denoted as secondary branches. For primary branches, based on the values of $\nu$ obtained from weakly nonlinear analysis at the bifurcating point, we classify these further into upper and lower branches. A natural question is the stability of the steady solutions mapped out in Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024), and in cases where they are unstable the investigation of the new dynamics they get attracted to. To examine this, we choose the unimodal primary branches (upper branch in blue, lower branch in orange) computed for $d_1=0.8$ and $d_2 = 0.31$ in Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024), and reproduced here in figure 3(a). Several solutions on these two primary branches are picked and perturbed according to $h(x; t=0) = h_{0}(x)+0.1\cos {2\pi {x}}$ , where $h_{0}(x)$ is the steady non-uniform periodic solution. Figure 3(a) also shows the location of the chosen solutions (stars for the upper branch and a diamond for the solution selected on the lower branch).

Figure 3. (a) The primary upper and lower branches of steady solutions of (2.33) for $d_1=0.8$ and $d_2 = 0.31$ . The stars represent chosen solutions on upper branches, the diamond represents the chosen solution on the lower branch. (b–d) The energy and wave profiles of the perturbed solutions and the undisturbed solutions. Panel (b) corresponds to point 1, panel (c) to point 2 and (d) to point 3. For the energy curves, the blue solid lines represent the evolution of energy of perturbed solutions while the red dashed lines is the energy of undisturbed steady solutions.

We begin with the stability of the solutions on the upper branch. The parameter $q$ in what follows is the constant flux introduced by Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024) and arises by looking for steady solutions of (2.33) and integrating once in $x$ . For the state corresponding to $q=0.43$ , $\nu =0.1336$ , figure 3(b) shows the evolution of the energy $E(t)$ of the perturbed solution indicating that it decays to a constant value corresponding to the energy of the steady solution – a comparison of $h_0(x)$ and $h(x,t=19)$ as an inset indicates this, confirming that the solution is stable. In fact increasing the perturbation amplitudes of the initial condition to $0.2$ , $0.4$ and $0.6$ also produces the steady state we started with. The stability of another solution on the upper branch denoted as point 3 in figure 3(a) and corresponding to $q = 0.48$ , $\nu = 0.058$ , is shown in figure 3(d). It can be seen that the energy of the waves oscillates initially but then increase to another larger constant value (the energy of the initial steady state is depicted by the dashed line in figure 3(d)). This illustrates that the steady solution is unstable and the instabilities evolve into a bimodal travelling wave (see the wave profiles in the inset to figure 3(d) with the dashed curve denoting the initial steady state). Next, we consider the stability of the lower branch solution and in particular point 2 in figure 3(a) at $q=0.4775$ , $\nu =0.068$ . The results are given in figure 3(c) and it is seen that the energy undergoes some large amplitude oscillations before arriving to a larger constant state after a longer time. The evolution of the perturbed wave leads to the emergence of another unimodal travelling wave shown by the solid curve in the inset to figure 3(c); the initial steady state is depicted with a dashed curve. These representative results show that the stability of steady states is complicated. It is not easy to examine the nonlinear stability of every solution branch in the bifurcation diagrams and we do not pursue this further.

Figure 4. Schematics of the various attractors for $d_1=0.8$ . UT: unimodal travelling waves, TP: time-periodic travelling waves, BT: bimodal travelling waves, QP: quasi-periodic attractors, NBT: near-bimodal travelling waves, NMT: near-multimodal travelling waves (including near bimodal and near trimodal).

4.2. Initial value problems, large-time dynamics and dominant attractors

A more direct way of probing the dynamics is to solve the initial value problem subject to sinusoidal initial conditions of the form

(4.1) \begin{align} h(x; t=0) = 1+0.1[\cos (2\pi {x}/L)+\sin (2\pi {x}/L)], \end{align}

and analyse the data from a large number of numerical experiments. Different initial conditions were tested for all results presented by increasing the amplitude to values of 0.2 and adding phase shifts to the sinusoids; all reported results are attractors in the sense that they emerge for such ranges of initial conditions. In all computations that follow we choose the air-speed parameter $d_1=0.8$ and vary the surface tension parameter $d_2$ . As $d_2$ increases surface tension becomes smaller and air-stream effects become dominant. For each fixed $d_2$ we carry out extensive computations to determine the dynamics and attractors as $\nu$ varies noting that as $\nu$ decreases the unscaled wave period becomes bigger. The results of these computations are presented schematically in the solution phase space in figure 4. The diagram depicts the most attracting state that emerges from the initial value problem. These have been categorised into (i) unimodal travelling waves (UT), (ii) time-periodic travelling waves (TP), (iii) bimodal travelling waves (BT), (iv) waves with quasiperiodic oscillations in time (QP), (v) near-bimodal travelling waves (NBT) and (vi) near-multimodal travelling waves (NMT). Window boundaries, where there is a transition from one state to another, have been computed within the accuracy of the reported decimals. In addition we have used the convention that the right boundary of each window represents the largest value of $\nu$ considered belonging to the window on the left. In all cases reported trivial solutions develop if $\nu$ is larger than a critical value $\nu _0(d_1, d_2)$ - these are the values of rightmost windows in figure 4. In fact we can predict this by linearising (2.33) about $h=1$ and looking for solutions of the form $\exp (2i\pi k x+st)$ ; setting $k=1$ readily produces the value of $\nu _0$ above which only trivial solutions exist

(4.2) \begin{align} \nu _0(d_1, d_2)=\frac {d_2}{4\pi }\left [\frac {\pi }{d_1}+\left (\frac {\pi ^2}{d_1^2}-\frac {4}{d_2}\right )^{1/2}\right ]\!. \end{align}

The values given by (4.2) are in complete agreement with our simulations. In addition, all solutions bifurcating from trivial states are unimodal travelling waves, as indicated in the figure. As $\nu$ decreases further the attractors depend on the surface tension parameter $d_2$ . In what follows we describe some of these solutions to provide a quantitative picture of the dynamics.

The general trends of the dynamics are depicted in the schematic of figure 4. Recall that a decrease in $\nu$ corresponds to an increase of the spatial domain size and hence an increase in the number of unstable modes entering into the nonlinear dynamics. Starting with the first window where non-trivial states are found, we see that as $\nu$ decreases, the first window of unimodal travelling waves (UT) typically gives rise to time-periodic travelling waves (TP) for the smaller values of $d_2$ (i.e. larger surface tension forces), or bimodal and near bimodal travelling waves for larger $d_2$ . Further decrease in $\nu$ mostly produces bimodal travelling waves, and for the smaller values of $d_2=0.31, 0.33$ we observe complex dynamics that appear to be quasi-periodic in time, or time-periodic travelling waves for larger $d_2$ .

Figure 5. Unimodal travelling waves for $d_1 =0.8$ , $d_2 = 0.33$ and $\nu = 0.1$ . (a) The evolution of the wave profile; (b) the evolution of the energy.

We provide a more detailed view of the flows for the case $d_2=0.33$ ( $d_1=0.8$ has been fixed also). In figures 5–8(a) we show representative solutions from each of the attracting windows labelled UT, TP, BT and QT in figure 4. Figure 5 has $\nu =0.1$ and shows the evolution of the solution as it strongly gets attracted to a unimodal travelling wave – figure 5(a) shows the profiles as time increases and figure 5(b) shows the evolution of the energy, clearly indicating that all transients have disappeared for $t$ as small as 2, approximately. It should be noted that for values of $\nu$ just larger than $\nu =0.072$ that denotes the window boundary of the TP attractor there are very long-lived transients with oscillations that eventually decay at large times to produce a unimodal travelling wave. Figure 6 presents the solution at $\nu =0.071$ from the time-periodic travelling wave attractor. Panel 6(a) shows the evolution for $t\gt 8.25$ to ensure that transients have disappeared as confirmed by the energy plot in panel 6(b) that also clearly shows the time-periodic oscillations. Decreasing $\nu$ further to $0.062$ we enter the bimodal travelling wave window as seen by the results of figure 7; once again there is strong attraction to these states as seen from the evolution of the energy in panel 7(b) with transients gone by $t\approx 2$ . A further reduction to $\nu =0.047$ puts us in a new solution window. Here, the solution undergoes spatio-temporal oscillations that we have analysed to determine that the time oscillations are quasi-periodic. The results are provided in figure 8 which shows the spatio-temporal evolution of the interface in panel (a) and the corresponding energy signal that strongly indicates a complex dynamics. Further data analysis of the latter suggests that the flow is quasi-periodic (for example, we did not observe any fractal folding behaviour in Poincaré maps). It is noted that when $\nu \leqslant {0.041}$ another bifurcation takes place to produce, initially at least, trimodal travelling waves. The computation can be continued by further decreasing $\nu$ and tracking subsequent attractors but for brevity we restrict our results by obtaining the value of $\nu$ for the onset of the trimodal travelling waves. Therefore, the lowest values of $\nu$ reported in figure 5 for different $d_2$ , correspond to the value where solutions first develop into trimodal travelling waves. It is worth noting that we can obtain qualitatively similar schematics of the various attractors when $d_2$ is fixed and $d_1$ decreases, but for brevity the results are not included here.

Figure 6. Time-periodic attractors for $d_1 =0.8$ , $d_2 = 0.33$ and $\nu = 0.071$ . (a) The evolution of the wave profile; (b) the evolution of the energy.

Figure 7. Bimodal travelling waves for $d_1 =0.8$ , $d_2 = 0.33$ and $\nu = 0.062$ . (a) The evolution of the wave profile; (b) the evolution of the energy.

Figure 8. The quasi-periodic solution for $d_1 =0.8$ , $d_2 = 0.33$ and $\nu = 0.047$ . (a) The evolution of the wave profile; (b) The evolution of energy signals.

From figure 4, we observe that an increase of the values of $d_2$ leads to a qualitatively different dynamics. For example, increasing $d_2$ to $0.36$ , appears to remove the quasi-periodic solutions window found for $d_2=0.31$ and $d_2=0.33$ . It is also noted that a type of near-bimodal waves first appear for $d_{2} = 0.38$ . When $d_2$ further increases to $d_2 = 0.4$ , the window of time-periodic states between the unimodal and bimodal travelling waves disappears; however, there still exist time-periodic solutions between the windows of bimodal and trimodal or near-trimodal travelling waves. When $d_2 = 0.42$ and as $\mu$ is decreased, only travelling waves can be observed preceding the onset of trimodal travelling waves. Larger values of $d_2$ were also explored in order to establish dominant trends. We find that, as $d_2$ increases, travelling waves attain larger amplitudes and have a tendency to touch down onto the substrate. A typical result is given in figure 9(a) for $d_1 = 0.8$ , $d_2 =0.7$ and $\nu =0.23$ . This behaviour is also observed for the upper primary branch in Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024) when the limiting configuration of waves with fixed $d_1$ is considered. We note that there is a fundamental connection between the steady-state solutions of Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024) and the travelling wave solutions computed here. In King et al. (Reference King, Tuck and Vanden-Broeck1993) and Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024) the film thickness is picked so as to give zero phase speeds for linear waves, and these are then followed into the nonlinear regime with the mean film thickness being part of the solution. In the present initial value problems it is more convenient to fix the mean film thickness – linear waves in general have non-zero phase speeds as expected. We expect, therefore, that the structure of the branches is similar and the numerical solutions given in figure 9(a) attest to this. The travelling wave solutions emerging in long-time computations are located on the upper branches if we computed the bifurcation diagrams for the travelling waves.

Figure 9. (a) The wave profile for $d_1 = 0.8$ , $d_2 = 0.7$ and $\nu = 0.23$ . (b) The wave profile for $d_1 = 0.8$ , $d_2 = 0.4$ and $\nu = 0.057$ . The wave speed is approximately $-0.6$ . (c) The bifurcation diagrams of travelling wave solutions with fixed mean thickness of 1 for $d_1 = 0.8$ , $d_2 = 0.4$ . The wave speeds are $-0.6$ and $0$ . (d) Wave profiles for $\nu = 0.057$ on the bifurcation curves in figure 9(c).

4.3. Near-multimodal travelling waves

Next we turn to the appearance of the near-bimodal travelling waves denoted by NBT in the schematic of figure 4. Figure 9(b) shows a typical wave profile of this type for $d_1 = 0.8$ , $d_2 = 0.4$ and $\nu = 0.057$ (note that this results from the large-time evolution of the waves starting with initial condition (4.1)). The wave speed of this solution is approximately $-0.6$ . Similar wave solutions were also found on vertical liquid films in Salamon, Armstrong & Brown (Reference Salamon, Armstrong and Brown1994), where a period-doubling symmetry was used to find the secondary solutions. In order to trace the appearance of these waves, we consider the bifurcation of travelling waves having the same wave speeds in the $\nu -q$ plane, where $q$ is the flux. More specifically, looking for travelling waves of (2.33) of the form $h(x-ct)$ and integrating once, we find

(4.3) \begin{align} -\frac {c}{2h^2}+\frac {3}{2h}+\frac {\nu ^2\pi }{d_1}{\mathcal{H}}[h_{xx}]+\frac {\nu ^3}{d_2}h_{xxx}-\nu h_x-1=\frac {q}{h^3}. \end{align}

Equation (4.3) is the same as that studied by Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024) for steady waves if we set $c=0$ and allow the mean film thickness to emerge as part of the solution. For a more direct comparison with our time-dependent solutions we fix the mean thickness to be $1$ and compute the solutions with wave speed $c =0$ and $c=-0.6$ , by constructing $\nu -q$ bifurcation curves, as shown in figure 9(c) for the case $d_1 = 0.8$ and $d_2 = 0.4$ . The bifurcation curves are similar to the steady primary upper branches in our previous study Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024), while all solutions here have mean thickness $1$ as mentioned above. Since the bifurcation curves in figure 9(c) branch off from infinitesimal linear waves, they are denoted as primary upper branch solutions in this section. It should be noted that the bifurcation computed here are similar to the primary branch in Salamon et al. (Reference Salamon, Armstrong and Brown1994) since the solutions in both studies start from a flat free surface. The difference is that the flux of the primary solution branch in Salamon et al. (Reference Salamon, Armstrong and Brown1994) is fixed. We see that branches with different wave speeds have similar shapes, and the solutions corresponding to $\nu = 0.057$ on these branches are shown collectively in figure 9(d)). By comparison of the wave profiles we preclude the possibility that the near-bimodal travelling waves are located on these primary branches. However, these type of near-bimodal solutions were found in Meng et al. (Reference Meng, Papageorgiou and Vanden-Broeck2024) when they calculated the secondary bifurcation branches. It is reasonable to conclude, therefore, that the near-bimodal travelling waves are on the secondary bifurcation branch. This in turn implies that the unimodal, near-bimodal and bimodal travelling waves emerging in the time-dependent computations presented in this study belong to different bifurcation diagrams. Additional evidence for this follows by examination of the energy of these travelling waves as $\nu$ varies, and in figure 10 we present the cases $d_1=0.8$ , $d_2 = 0.42$ and $d_{2} = 0.38$ . In figure 10(a), the jumps in the energy show that there exist three distinct solution branches in the range $0.08\leqslant \mu \leqslant 0.12$ , corresponding to the bimodal, near-bimodal and unimodal travelling waves, respectively. We also examine the energy variation for $d_{2} = 0.38$ . It can be seen from figure 10(b), that the jump also exists for this case. For $d_{2} = 0.38$ , the window of the near-bimodal waves is very narrow. The energies of the solutions in this window are very close, therefore, in the plot we only include the energy of one solution (others are graphically indistinguishable). The fact that near-multimodal travelling waves appear when the value of $d_2$ is above a threshold reveals the change of stability of solutions on the secondary bifurcation branches as $d_2$ increases. When $d_2$ is relatively small (smaller than $0.38$ , roughly), the solutions on the secondary bifurcation branches are all unstable and hence travelling wave solutions of this type do not appear as large-time solutions of our time-dependent computations. The stable travelling wave solutions appear to be exclusively located on the primary bifurcation branches. As the value of $d_2$ increases above $0.38$ , some solutions on secondary bifurcation branches become stable and the near-bimodal waves (analogously for the near-trimodal waves – confirmed numerically but results not included here) appear as global attractors in the time-dependent computations.

Figure 10. The variation of the energy of travelling waves for $d_1 = 0.8$ and different $d_{2}$ . Panel (a) shows $d_{2} = 0.42$ . From left to right, the three branches correspond to waves of the bimodal, near-bimodal and unimodal type. Panel (b) shows $d_{2} = 0.38$ . From left to right, the three branches correspond to waves of the near-bimodal, bimodal and unimodal type.

Figure 11. The boundary for the onset of bimodal (the dashed red line and dotted blue line), trimodal (the dashed-dotted line) and tetramodal(the solid line) travelling wave solutions in the $d_2-\mu $ plane for $d_1 = 0.8$ . The boundary for the onset of tetramodal travelling waves is predicted by the boundary of bimodal travelling waves by the scaling $\nu \rightarrow {\nu /2}$ . The data marked by crosses are obtained numerically.

4.4. Bimodal, trimodal and tetramodal travelling waves

In figure 11, we depict the boundaries for the onset of both bimodal and trimodal waves in $(d_2,\nu )$ space, for fixed $d_1=0.8$ ; as surmised from the results of the schematic (figure 4), bimodal and trimodal travelling waves appear as $\nu$ is decreased for a fixed value of $d_2$ . The dashed (red) line in figure 11 shows the boundary of $\nu$ below which bimodal travelling waves emerge (the open circles are results from the computations and the straight line is the result of a least squares fit of the open circle data). With a decrease in $\nu$ trimodal travelling waves appear and the boundary is now depicted by the dash-dotted line (purple), with diamonds representing computations and the straight line is a least squares fit of the diamond data. Further decrease in $\nu$ produces tetramodal travelling waves as shown – the crosses come from nonlinear computations and the solid (magenta) line is derived from the bimodal curve using the symmetry scaling $\nu \rightarrow \nu /2$ , see below. Interestingly, similar phenomena were also found in the weakly nonlinear dynamics of perfectly conducting electrified falling liquid films by Tseluiko & Papageorgiou (Reference Tseluiko and Papageorgiou2006), where constant slopes are found for a specific boundary between attractors despite the change of the computation domain parameter (analogous to our $\nu )$ . For the fully nonlinear model studied here, the solution boundaries are more complicated. For example, the boundary of the bimodal travelling waves for $d_2 \lt 0.41$ has larger slope than that for $d_2 \geqslant {0.41}$ , which is shown in figure 11 with a dotted green line and square symbols coming from nonlinear computations. This break in the curve indicates that, for sufficiently large $d_2$ , the onset of bimodal waves happens at a smaller $\nu$ than what would be predicted by the extrapolated linear behaviour of the range of smaller $d_2$ values (e.g. $d_2\leqslant 0.41$ in this example). As seen in figure 11, there is a gap between these two boundaries in the interval $0.4\lt d_2\lt 0.41$ . This fact is consistent with what we described above in § 4.3: the near-bimodal travelling waves on the secondary bifurcation branches change their stability and become stable when $d_2$ is large enough. For the values of $\nu$ in this gap, the solutions develop into the near-bimodal waves instead of the bimodal waves – this has been confirmed numerically. Comparing the bimodal and trimodal boundaries in the range $d_2\leqslant 0.41$ , we find that the values of $\nu$ on the former boundary has values of $\nu$ that are approximately $1.5$ times of those on the trimodal wave boundary. As mentioned, the boundaries (the lines) for the bimodal and trimodal travelling waves are obtained by the least square fits of the numerical data (the circles, squares and diamonds). Using the scaling laws $\nu \to \nu /2$ of the boundaries, we can predict the boundary for the onset of tetramodal travelling waves from the boundary of bimodal waves by the scaling. To verify this predicted boundary, we carried out numerical computations to find the onset values of $\nu$ for two values of $d_2 = 0.31$ and $d_2 = 0.33$ (marked by crosses). Good agreement is found with the predicted boundary.

5.1. Low Reynolds numbers: time-periodic attractors and nonlinear travelling waves

We begin with the case $\delta =0.03$ and carry out numerical computations on periodic domains of unit length ( $\beta =2\pi$ in (5.1)), with $\bar {C}=0.69$ , $\bar {D}=13.33$ , $\bar {\tau }=1.84$ and $\alpha =\pi /8$ . The range of Reynolds numbers in this numerical experiment is $0 \lt R \leqslant 0.7$ . Note that even for the largest Reynolds number $R=0.7$ , we have ${8R}/{15}-{2}\cot {\alpha }/3\approx -1.24$ placing us firmly in the subcritical stable regime in the absence of an outer air flow. Using the chosen parameters in (5.1) we find that, for all $0\lt R\leqslant 0.7$ , there are two linearly unstable modes, $k=1$ and $k=2$ . The evolution becomes highly nonlinear even on such relatively short domains due to the energy supplied by the air stream to rapidly increase the wave amplitude via an interaction with the inertial terms – we provide quantitative details later. Our computations show that in the range $0.5\leqslant R \leqslant 0.7$ the solutions develop into nonlinear travelling waves, while at smaller $R$ in the range $0 \lt R \lt 0.5$ , no travelling waves appear and the solutions saturate to time-periodic attractors with different temporal oscillation periods. Figure 12(a) and figure 12(b) show the variation of wave speed and amplitudes of the travelling waves in the range $0.5 \leqslant R \leqslant 0.7$ . It can be seen that with the increase of the Reynolds number, the wave speeds and amplitudes increase, a feature that was also found in falling liquid film flows without the air stream, see Pumir et al. (Reference Pumir, Manneville and Pomeau1983). In fact, when we further increase $R$ slightly to $0.73$ , the solutions become unbounded (discussed in detail later). Details of the time-periodic solutions found for $0\lt R\lt 0.5$ are given in figure 12(c) which plots the time period $T$ of the oscillations against $R$ . It is seen that as $R$ increases the time period decreases monotonically at an apparently linear rate until the attractor looses stability to travelling waves at $R\geqslant 0.5$ .

Figure 12. Effect of the Reynolds number on (a) travelling wave speeds, (b) travelling wave amplitudes and (c) periods of oscillations of time-periodic solutions, for $\bar {C} = 0.69$ , $\bar {\tau } = 1.84$ , $\bar {D} = 13.33$ , $\alpha = \pi /8$ and $\delta = 0.03$ .

5.2. Reynolds number above critical: blow up in finite time

Now we turn to the phenomenon of unbounded solutions of the Benney equation documented in the absence of an air stream when the Reynolds number exceeds a critical value. Such singular behaviour of the Benney equation was first reported numerically by Pumir et al. (Reference Pumir, Manneville and Pomeau1983). The stability structure, including both bounded and unbounded solutions, was later investigated by Rosenau et al. (Reference Rosenau, Oron and Hyman1992) and Oron & Gottlieb (Reference Oron and Gottlieb2002) described the phenomenon quantitatively by numerically calculating the boundary between the regions of bounded and unbounded solutions in $(\epsilon , R)$ space, where $\epsilon$ is equivalent to our slenderness parameter $\delta$ . In later work Scheid et al. (Reference Scheid, Ruyer-Quil, Thiele, Kabov, Legros and Colinet2005) connected dynamical systems theory to the numerical results obtained by Oron & Gottlieb (Reference Oron and Gottlieb2002). They concluded that the absence of one-hump travelling wave solutions is related to the finite-time blow up of the Benney equations. Note that the results mentioned above are based on the classical Benney equation with the instability arising from the $h_{xx}$ term when $R\gt 5\cot \alpha /4$ . For the present counter-current air/liquid problem, the presence of the air stream enhances the instabilities and leads to rich phenomena when we consider the boundary between the bounded and unbounded solution regions in parameter space.

We have carried out extensive numerical computations to categorise the dynamics and different attractors in the $R-\delta$ parameter plane, in order to evaluate the novel effects of the air stream on subcritical unforced Benney solutions. All computations having $\bar {C}=\bar {\tau }=0$ would be stable. The way we construct this solution space is to fix the unscaled air-stream parameter $C$ and surface tension parameter $D$ and study the effect of varying $\delta$ . To fix matters we considered $C=23$ , $C_D=0.08$ , $D=0.012$ and $\alpha =\pi /8$ . In terms of these unscaled parameters (recall from § 2.2 that $\overline {C}=\delta C$ , $\overline {D}=D/\delta ^2$ and $\overline {\tau }=C C_D={\mathcal{O}}(1)$ ) (2.39) reads

(5.2) \begin{align} h_{t} &+ \left[-\frac {2}{3}h^3+\frac {{C C_D}}{2}h^2+\delta \left (\frac {8R}{15}h^6-\frac {4R{C C_D}}{15}h^5-\frac {2\cot \alpha }{3}{h^3}\right )h_{x} \right. \nonumber\\ & \left. +\, \frac {2\delta ^2{C}}{3}h^3{\mathcal{H}}[h_{xx}]+\frac {\delta ^3}{3{D}}h^3h_{xxx} \right ]_{x}=0. \end{align}

The spatial period is set to $L=1$ and the initial value problem of (5.2) is computed as $R$ and $\delta$ are varied. This is more convenient than working with the scaled (2.39) which requires $\overline {C}$ and $\overline {D}$ to change according to the scalings given above every time $\delta$ is changed. The results are presented collectively in figure 13. For every value of the slenderness parameter $\delta$ , there exists a critical value of the Reynolds number below which the solutions are bounded at large times and above it a finite-time blow up is found. Below this critical value the dynamics of the bounded solutions is rich and depends on the value of $\delta$ . Figure 13 identifies three regions, delineated by the horizontal dashed lines. For $\delta$ between $0.02$ and $0.023$ we find bounded solutions to the left of the dashed-green curve that are time-periodic travelling waves. When $\delta$ is approximately between $0.024$ and $0.028$ the bounded solutions to the left of the purple-dashed curve now evolve to give bimodal travelling waves, and for $\delta$ larger than approximately $0.029$ and $0.043$ we obtain bounded unimodal travelling waves to the left of the red-dashed curve. It is interesting to note that, in Oron & Gottlieb (Reference Oron and Gottlieb2002), the boundaries in the $R-\delta$ plane are given for what they term $n = 1$ and $n=2$ solutions that represent the one-hump and two-hump travelling waves. As outlined above the situation is different here and bounded attractors below critical values of $R$ are more complicated and analogous to those for the inertialess solutions presented in figure 4 in § 4.2.

Figure 13. The boundary between the domains of bounded and unbounded solutions for the Benney equation computed numerically for $C = 23$ , $C_{D} = 0.08$ , $D = 0.012$ , $\alpha = \pi /8$ .

Figure 14. (a) The shape of the surface at three different times with the parameters $C = 22$ , $D = 0.011$ , $\delta = 0.0232$ , $C_{D} = 0.08$ , $\alpha = \pi /8$ and $R = 2$ just preceding the blow-up. The computation domain $L = 1$ . (b) The size of different terms at $t = 9.7867$ . (c) The shape of the surface at three different times with the parameters $C = 22$ , $D = 0.011$ , $\delta = 0.0232$ , $C_{D} = 0.08$ , $\alpha = \pi /8$ and $R = 3$ just preceding the blow-up. The computation domain $L = 2\pi$ . (d) Growth rate for $L = 1$ (above) and $L = 2\pi$ (below) for the two cases discussed here. The wavenumbers with positive growth rate are marked by circles.

Some typical results showing blow up are depicted in figures 14(a) and 14(b). The parameters for this run are $\delta =0.0232$ and $R=2$ . We are slightly above the critical Reynolds boundary and two unstable modes are present according to the dispersion relation (5.1) with $L=1$ the growth rate curve is shown in the upper part of figure 14(d) with the two unstable modes $k=1$ and $k=2$ indicated by a circle. Figure 14(a) superposes the interfacial profiles at three different times just before blow up – the times are given in the inset to figure 14(a). A narrowing of the wave profiles can be observed coupled with a rapid growth in the wave amplitude. In the adjacent panel 14(b) we show the spatial behaviour of the three competing nonlinear terms in (2.39) at the final computed time $t=9.7867$ just prior to the blow up. to better understand the singularity. The spatial dependence of the inertial term, surface tension term and air-induced pressure term at $t= 9.7867$ are plotted in figure 14(b). These terms are the surface tension term shown in blue, the inertial term shown in red and the air-induced pressure term shown in orange. The data indicate that there is a dominant balance between surface tension and inertia, with the air-induced pressure being asymptotically subdominant. The role of the air stream is to induce instability and drive the system into a highly nonlinear regime that produces a balance between surface tension and inertia. It is expected, therefore, that the singular structure will be the same as for the Benney equation as described by Pumir et al. (Reference Pumir, Manneville and Pomeau1983). If the singularity takes place at $t=t_s$ and at position $x=x_s$ , then the self-similar solution near the singularity is

(5.3) \begin{equation} h \sim (t_{s}-t)^{-1/9}G(\xi ),\qquad \xi =\frac {x_{s}-x}{(t_{s}-t)^{1/6}}, \end{equation}

where the shape function $G$ satisfies a nonlinear ordinary differential equation in $-\infty \lt \xi \lt \infty$ . We do not pursue this analysis further here since it has been addressed for the Benney problem, however, we have confirmed that scalings in (5.3) are consistent with our numerical data. We have also carried out computations on larger domains that support more linearly unstable modes. A blow up singularity is encountered once again, but with additional nonlinear wavy structures in the final profile. The case $L=2\pi$ is included in figure 14(c) that shows the solution at three consecutive times close to the singularity. The wave structures arise due to the additional unstable modes (there are now 10 unstable modes) as shown in the lower panel of figure 14(d). Note that both growth rate curves in figure 14(d) are of Turing type and support short-wave unstable modes due to the presence of the air stream. These new modes are exclusively due to the air stream since in its absence the flow would be stable. In addition, instabilities without the air stream are only possible above a critical Reynolds number and are of the long-wave type.

In the last part of this section, we explore the effect of the air stream $\overline {C}$ on the critical value of $R$ that heralds a transition from bounded to unbounded solutions. We do this for a fixed value $\delta =0.04$ , a fixed value of $D=0.012$ (equivalently $\overline {D}=7.5$ ), $\alpha = \pi /8$ , and also fix $C_D=0.08$ , i.e. $\overline {\tau }=CC_D=2\overline {C}$ , given the scaling $C=\overline {C}/\delta$ . Hence, we can study the blow up behaviour in terms of the two parameters $\overline {C}$ and $R$ . Solution of the initial value problem allowed us to numerically estimate the boundary between bounded and unbounded solutions in the $(\overline {C},R_c)$ space, where $R_c$ is our estimate for the critical Reynolds number. The results are given in figure 15. It can be seen that the presence of the air stream enhances singularity formation of the air-stream Benney equation in the sense that blow up occurs at smaller values of the Reynolds number as the air speed increases. The decrease in $R_c$ is monotonic and quite strong – an increase of $\overline {C}$ from approximately $0.9$ to $1.1$ , causes a tenfold decrease in $R_c$ from $2.2$ to $0.2$ . These results suggest that the air stream can be used systematically to trigger large amplitude waves in the flow in the system.

Figure 15. The effect of the scaled air-stream parameter $\overline {C}$ on the critical value $R_c$ of the Reynolds number where the transition from bounded to unbounded solutions occurs. Here, $\delta = 0.04$ , $C_{D} = 0.08$ ( $\overline {\tau }=CC_D=\overline {C}C_D/\delta =2\overline {C}$ ), $\bar {D} = 7.5$ ( $D=0.012$ ), $\alpha = \pi /8$ . Diamonds represent computational data.

5.3. Effects of non-uniform shear stress

In this section the effects of non-uniform tangential stress that can be introduced by the air stream are investigated for both the surface tension dominated and the inertio-capillary regimes. Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) treated the interface as a rigid wall, to leading order, and modelled the tangential and normal stresses induced by a counter-flowing turbulent gas. The analytical results were also compared with the experiments of Thorsness, Morrisroe & Hanratty (Reference Thorsness, Morrisroe and Hanratty1978) that measured the spatial variation of shear stress by the turbulent flow along a solid wall with small sinusoidal indentations. The theory shows good agreement with the experiments. In light of their results, we now discuss the effect of the variation of the tangential stress on the models developed here under the assumption of a constant shear stress. Our main aim is to evaluate the effect of non-uniform shear derived from theoretical and experimental studies, on the nonlinear terminal dynamics of the inertio-capillary regime.

Following Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011), the tangential shear can be represented by the order-one constant term ( $\tau _{w0}$ ) and higher-order variations ( $\tau _{w1}\cdots$ ). The surface tension dominated model (2.32) arises at leading order and only order-one terms are retained, i.e. $\tau _{w0}$ . The non-uniform shear stress has no effect on the surface tension dominated regime. For the inertio-capillary regime, the higher-order variation $\tau _{w1}$ will enter into the next-order equation in the formulation. It is straightforward to modify our derivation to include this term. The modified equation reads

(5.4) \begin{eqnarray} h_{t}+\left [-\frac {2}{3}h^3+\frac {\overline {\tau }}{2}h^2+\delta \left (\frac {8R}{15}h^6-\frac {4R\overline {\tau }}{15}h^5-\frac {2\cot \alpha }{3}{h^3}\right )h_{x}+\frac {2\delta \overline {C}}{3}h^3{\mathcal{H}}[h_{xx}]\right .\nonumber \\ \left .+\frac {\delta }{3\overline {D}}h^3h_{xxx}+\delta \dfrac {\bar {\tau }}{2}\tau _{w1}[h]h^2\right ]_{x}=0, \end{eqnarray}

where $\tau _{w1}[h]$ is a linear non-local operator involving $h$ and has been calculated by Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) (its symbol $\widehat {\tau _{w1}}(k)$ say, in Fourier space is known). From the dispersion relation (figure 12 in Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011)), it can be seen that the inclusion of this term introduces instabilities and alters the linear wave speed. It should be noted that, in their paper, the positive $x$ direction is opposite to that in the present paper, implying that for our coordinate system the real part of $\widehat {\tau _{w1}}(k)$ is positive. We expect therefore, that for sufficiently low Reynolds numbers (i.e. prior to any blow-up events), it is more likely that the nonlinear waves travel upwards. We do not pursue travelling waves further and concentrate more on blow-up events.

We now focus on the blow-up behaviour of solutions to (5.4) incorporating spatial variations depending on the interfacial position in the tangential shear stress. We examine the cases with parameters $C = 23$ , $D = 0.012$ , $\alpha = \pi /8$ and $C_{d} = 0.08$ (the same parameters as in figure 13). Two values of $\delta$ are chosen, $\delta = 0.038$ and $\delta = 0.042$ . The initial conditions are the same as those used in figure 13. For the constant tangential stress, the critical Reynolds numbers above which blow-up occurs for these two values of $\delta$ are approximately $1.1$ and $1.6$ , respectively. Solution of (5.4) incorporating shear-stress variations are given in figure 16. The figure shows the energy evolution of solutions corresponding to $\delta =0.038$ (a) and $\delta =0.042$ (b), for different Reynolds numbers. It is strongly suggested that the critical Reynolds numbers where the solutions go unbounded do not change noticeably, indicating that inclusion of the next-order variation of shear stress has little effect on the blow-up boundary. Note that, if we do not assume that the air speed is large enough to induce the modulated pressure term, the equation obtained reduces to the model in Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011). In such a regime, the unstable modes introduced by the air into the system arise from the variation of shear stress alone (the mean part of the shear stress contributes an advective term linearly). It is expected that the critical Reynolds numbers where solutions go unbounded will change, but the self-similar solutions near the singular time remain same and produce a balance between surface tension and inertia. For the solutions with low Reynolds numbers, the inclusion of the non-uniform shear stress potentially alters the speed of nonlinear travelling waves.

Figure 16. The energy evolution of transient solutions of different values of $\delta$ : (a) $\delta = 0.038$ ; (b) $\delta = 0.042$ .