1. Introduction
The dynamics and stability of liquid film flow have always received much attention as they directly affect the heat and mass transfer characteristics in industrial applications, especially under phase change conditions, such as water-cooled walls (Zhang et al. Reference Zhang, Wang, Jiang and Wu2024), condensers (Zhang et al. Reference Zhang, Jia, Dang and Qi2022), evaporators (Wang et al. Reference Wang, Li, Xu, Yao, Liu, Su and Wang2020), absorbers (Jenks & Narayanan Reference Jenks and Narayanan2008) and heat pipes (Zhang & Nikolayev Reference Zhang and Nikolayev2023). Consequently, throughout several decades, a large amount of theoretical and experimental work on liquid films has been conducted, which has been reviewed by Oron, Davis & Bankoff (Reference Oron, Davis and Bankoff1997), Craster & Matar (Reference Craster and Matar2009) and Blossey (Reference Blossey2012).
The driving forces of liquid film flow primarily consist of gravity and shear stress on the surface. Since the groundbreaking experiment by the Kapitza family in 1949, the gravity-driven film has intrigued many researchers due to its wide range of spatial and temporal structures (Chang Reference Chang1994; Zheng & Stone Reference Zheng and Stone2022). For small Reynolds numbers (Re = O(1)), where the flow is predominantly influenced by surface tension and the inertial force is much smaller, the Kuramoto–Sivashinsky (KS) equation can be obtained through weakly nonlinear theory (Chang, Demekhin & Kopelevich Reference Chang, Demekhin and Kopelevich1993). For moderate Reynolds numbers (
$\varpi \cdot Re = O(1)$, where 
$\varpi = d/l$, l is the wavelength of the disturbance and d is the average thickness of the thin film), where the surface tension is comparable to the inertial force, a commonly employed approach involves using the integral boundary layer (IBL) equation along with the application of weighted-residual methods (Ruyer-Quil, & Manneville Reference Ruyer-Quil and Manneville2000; Oron, Gottlieb & Novbari Reference Oron, Gottlieb and Novbari2009). These equations suggest that the falling film yields a rich spectrum of fascinating wave dynamics at different parameters, such as solitary waves, periodic waves, bifurcations and chaotic solutions (Chang & Demekhin Reference Chang and Demekhin2002). For high Reynolds numbers, 
$\varpi \cdot Re \gg 1$, the inertial force becomes predominant while the surface tension plays a secondary role, the transition to turbulence may occur and full-scale Navier–Stokes equations must be considered in numerical simulations (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2011).
The film flow driven by shear stress on the surface has also been examined extensively. Smith & Davis (Reference Smith and Davis1982, Reference Smith and Davis1983a,Reference Smith and Davisb) performed linear stability analysis on such flows and have identified two types of instabilities in thin films. The first type is the surface-wave instability, which can be induced by either wind stress or Marangoni forces. When considering two-dimensional traveling waves, a critical Reynolds number is detected in a linear velocity profile, whereas a long-wave instability arises in a parabolic velocity profile. Recently, Patne, Agnon & Oron (Reference Patne, Agnon and Oron2021) extended these findings to three-dimensional waves in a liquid film with an oblique temperature gradient. Jurman & McCready (Reference Jurman and McCready1989) examined the linear and weakly nonlinear stability of a thin, horizontal liquid film sheared by a concurrent turbulent gas flow. In another study, Frank (Reference Frank2006, Reference Frank2008) numerically investigated long nonlinear two-dimensional traveling waves on a liquid film driven by laminar flow of a gas, showing the waves with equal amplitudes but different phase speeds. The second type is the convective instability appearing in thermocapillary liquid layers, which is driven by mechanisms within the bulk of the layer. The critical mode includes the stationary rolls and travelling waves.
In many practical problems of liquid film, the driving force comprises gravity as well as surface shear stress. The instability of this flow has also received much attention. Miesen & Boersma (Reference Miesen and Boersma1995) studied the linear stability of a vertically falling film sheared by a concurrent gas flow, showing the dependence of the critical Reynolds number on the Weber number, on the curvature of the liquid velocity profile and on the properties of the gas. Miladinova et al. (Reference Miladinova, Slavtchev, Lebon and Legros2002) examined the long-wave instabilities occurring in non-uniformly heated falling films, where the flow is influenced by both gravity and thermocapillary forces. Building upon this research, Mukhopadhyay & Mukhopadhyay (Reference Mukhopadhyay and Mukhopadhyay2021) extended the investigation by incorporating the impact of odd viscosity. Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) conducted a study on the dynamics of a thin laminar liquid film flowing down an inclined plane when turbulent gas flows above it. Later, this work was generalized by Vellingiri, Tseluiko & Kalliadasis (Reference Vellingiri, Tseluiko and Kalliadasis2015) to the absolute and convective instabilities. In addition, the effects of condensation (Aktershev & Alekseenko Reference Aktershev and Alekseenko2005), evaporation (Mohamed, Dallaston & Biancofiore Reference Mohamed, Dallaston and Biancofiore2021), a vibrating inclined plane (Samanta Reference Samanta2021), surfactants (Hossain, Ghosh & Behera Reference Hossain, Ghosh and Behera2022) and thermocapillary instabilities (Choudhury & Samanta Reference Choudhury and Samanta2023) have been considered by several authors.
Unlike the three situations mentioned above, the liquid film on a horizontal plate in a gravitational field can flow under the action of surface shear forces. This scenario is common in film cooling, where a thin liquid film on a solid surface can move under friction of a forced gas and cool down the surface. Film cooling has important applications in the fields of electronic devices (Kabov, Kuznetsov & Legros Reference Kabov, Kuznetsov and Legros2004) and aero-engines (Acharya & Kanani Reference Acharya and Kanani2017). The stability of the liquid film directly relates to its cooling effect. Numerous studies have been carried out to explore the linear stability of these flows, including the two-fluid boundary layer stability (Özgen, Degrez & Sarma Reference Özgen, Degrez and Sarma1998), the convective/absolute instability in laminar and turbulent gas–liquid two-layer channel flow (Naraigh, Spelt & Shaw Reference Ó Naraigh, Spelt and Shaw2013), and the role of water layer depth on the growth rate of the Miles and rippling instabilities (Kadam, Patibandla & Roy Reference Kadam, Patibandla and Roy2023). Gravity is not the primary driving force; however, it remains a vital factor contributing to the flow instability. However, as far as we know, there has been little investigation on the nonlinear instability of this problem. Although large-amplitude disturbances and instances of film rupture in sheared liquid layers have been observed in experiments (Hirokawa, Ohta & Kabov Reference Hirokawa, Ohta and Kabov2015; Wang et al. Reference Wang, Zhang, Gong, Bai and Ma2017), the evolution of nonlinear waves in liquid films under the influences of surface tension and gravity is still unclear, and this is the central objective of the current study. Here, we use theoretical and numerical approaches to examine the nonlinear dynamics of perturbation waves in shear flows of thin films. Given that most instances of rupture and instability phenomena in thin films typically occur on longer scales (Oron et al. Reference Oron, Davis and Bankoff1997; Craster & Matar Reference Craster and Matar2009), the long-wave theory (lubrication theory) is applied in our analysis. By examining various parameters and conditions, we illustrate the diverse and complex dynamical behaviours exhibited by liquid films. Notably, our findings highlight the significant influences of gravity, in terms of both magnitude and direction, on the evolution of disturbance waves.
Due to the widespread presence of wave steepening and film rupture in shear liquid films, which exhibit significant deformation, linear and weakly nonlinear stability analyses fall short in fully understanding these phenomena. Thus, this paper investigates finite-amplitude perturbations with strong nonlinearity. The results indicate that within certain parameter ranges (wave speed, wave number, gravity, etc.), the liquid film can still sustain steady traveling wave solutions, while the relative variation of the film thickness reaches O(1). For instance, we present the relevant parameters for solitary waves in figure 10 and two wave patterns for steady-state periodic waves in figures 12 and 14. In contrast, beyond these ranges, the perturbations cannot stably exist and are likely to continuously grow, eventually leading to film rupture. Additionally, we demonstrate examples of highly nonlinear waves through our experiment in figure 2 and the numerical simulation in figure 15. Therefore, our findings contribute to the understanding of deformation and rupture in shear liquid films.
The paper is organized as follows. Section 2 presents the physical model and mathematical formulations, wherein the dimensionless governing equations and boundary conditions are derived. Section 3 is dedicated to the numerical results, wherein the solitary and periodic solutions of travelling waves are examined. The effects of gravity and surface tension are discussed. We perform an energy analysis in § 4 and conduct a numerical simulation for the gas–liquid two-phase flow in § 5. Finally, § 6 concludes the paper by summarizing the key findings.
2. Problem formulation
We consider a liquid film on an infinitely large plate placed horizontally in a gravitational field, as shown in figure 1. The gas–liquid interface experiences a constant shear stress 
$\tau $, which is caused by wind shear. In this case, the fluid within the liquid film undergoes motion due to the shear force acting on the interface, while gravity is parallel to the wall-normal direction. Here, x and z represent the streamwise and normal directions, respectively. In the following, we will examine two cases: one where the direction of gravity is opposite to the normal direction (figure 1a), and the other where they are the same (figure 1b). Because the gravity direction points to the lighter fluid, the so-called Rayleigh–Taylor instability would arise in figure 1(b).
Schematic of the liquid film flow on a horizontal plate driven by a constant shear stress at the gas–liquid interface. Two scenarios are depicted: (a) where the direction of gravity is opposite to the wall-normal direction and (b) where they are the same.
The case in figure 1(a) has been tested in the experiment of Wang et al. (Reference Wang, Zhang, Gong, Bai and Ma2017). However, the case in figure 1(b) can be realized by dispensing a certain amount of liquid onto a solid plate, then using the wind to drive the liquid film lying on the underside of the horizontal plane. We perform this experiment shown in figure 2 and a supplementary movie is available at https://doi.org/10.1017/jfm.2024.895.
Disturbance waves in a sheared liquid film lying on the underside of the horizontal plane: (a) photograph and (b) schematic of the experiment.
2.1. Governing equations
The liquid is supposed to be an incompressible Newtonian fluid with the dynamic viscosity 
$\mu $ and density 
$\rho $. The average thickness of the liquid film is d, then the scales of velocity and time can be defined as 
$\bar{U} = \tau d/\mu $ and 
${t_0} = \mu /\tau $, respectively. The surface tension coefficient is 
$\sigma $. The following dimensionless groups can be defined:
\begin{equation}Re = \frac{{\rho d\bar{U}}}{\mu } = \frac{{\rho {d^2}\tau }}{{{\mu ^2}}},\quad S = \frac{{\rho d\sigma }}{{{\mu ^2}}},\quad Ca = \frac{{Re}}{S},\quad G = \frac{{\rho g{d^2}}}{{\mu \bar{U}}} = \frac{{\rho gd}}{\tau }.\end{equation}Here, Re is the Reynolds number; S is the non-dimensional surface-tension number; Ca is the capillary number, which measures the magnitude of the surface deformation; and G is the Galileo number, which measures the gravity effect.
The dimensionless governing equations are given below, which are the continuity equation and the momentum equation, respectively,
\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0,\end{gather}
\begin{gather}Re\left( {\frac{{\partial \boldsymbol{u}}}{{\partial {\kern 1pt} t}} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}} \right) ={-} \boldsymbol{\nabla }p + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau } - G\boldsymbol{\nabla }z.\end{gather}
Here, 
$\boldsymbol{u} = (u,v,w)$, p and 
$\boldsymbol{\tau }$ stand for the velocity, pressure and stress tensor, respectively. For a Newtonian fluid,
\begin{equation}\boldsymbol{\tau } = \boldsymbol{S},\end{equation}
where 
${\textstyle{1 \over 2}}\boldsymbol{S} = {\textstyle{1 \over 2}}(\boldsymbol{\nabla }\boldsymbol{u} + {(\boldsymbol{\nabla }\boldsymbol{u})^\textrm{T}})$ is the strain-rate tensor.
To simplify the problem, we only consider two-dimensional (2-D) flows in the present work. In general, three-dimensional (3-D) long-wave dynamics is preferable to two-dimensional analysis in liquid films (Tomlin, Papageorgiou & Pavliotis Reference Tomlin, Papageorgiou and Pavliotis2017). For example, in falling films, 2-D waves with straight wave fronts are excited first. As the waves travel further downstream, the 2-D wave fronts begin to break and 3-D patterns arise. However, our experiment suggests that if the spanwise length of the liquid film is much smaller than its streamwise length, the disturbance waves can be considered approximately two-dimensional (figure 2 and supplementary movie). The fluid used in our experiment is corn oil, with its viscosity 
$\mu = 45.6 \times {10^{ - 3}}\;\textrm{kg}\;{(\textrm{m} \cdot \textrm{s})^{ - 1}}$, density 
$\rho \approx 0.92 \times {10^3}\;\textrm{kg}\;{\textrm{m}^{ - 3}}$, surface tension coefficient 
$\sigma \sim 32 \times {10^{ - 3}}\;\textrm{kg}\;{\textrm{s}^{ - 2}}$, film thickness 
$d\sim 4 \times {10^{ - 4}}\;\textrm{m}$ (measured by Keyence laser displacement sensor CL-3000), spanwise length 
${l_1}\sim 6 \times {10^{ - 3}}\;\textrm{m}$, streamwise length 
${l_2}\sim {10^{ - 1}}\;\textrm{m}$, wavelength 
$l\sim 1.4 \times {10^{ - 2}}\;\textrm{m}$ and characteristic flow velocity 
$\bar{U}\sim 4 \times {10^{ - 2}}\;\textrm{m}\;{\textrm{s}^{ - 1}}$. The dimensionless numbers are 
$Re\sim O(0.3)$, 
$G\sim O(\textrm{0}\textrm{.8})$, 
$Ca\sim O(0.\textrm{06})$, all conforming to the conditions set in this paper.
Thus, the governing equations of two-dimensional flows are
\begin{gather}\frac{{\partial u}}{{\partial x}} + \frac{{\partial w}}{{\partial z}} = 0,\end{gather}
\begin{gather}Re\left( {\frac{{\partial u}}{{\partial t}} + u\frac{{\partial u}}{{\partial x}} + w\frac{{\partial u}}{{\partial z}}} \right) ={-} \frac{{\partial p}}{{\partial x}} + \frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {z^2}}},\end{gather}
\begin{gather}Re\left( {\frac{{\partial w}}{{\partial t}} + u\frac{{\partial w}}{{\partial x}} + w\frac{{\partial w}}{{\partial z}}} \right) ={-} \frac{{\partial p}}{{\partial z}} + \frac{{{\partial ^2}w}}{{\partial {x^2}}} + \frac{{{\partial ^2}w}}{{\partial {z^2}}} - G.\end{gather}
The boundary condition at the wall 
$(z = 0)$ is
\begin{equation}\boldsymbol{u} = (u,w) = 0,\end{equation}
while the boundary conditions on the free surface 
$(z = \eta (x\textrm{,}t))$ are
\begin{gather}{\partial _t}\eta + u \cdot {\partial _x}\eta = w\textrm{,}\end{gather}
\begin{gather}- p + \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\sigma }\boldsymbol{\cdot }\boldsymbol{n} ={-} C{a^{ - 1}}({\nabla _s}\boldsymbol{\cdot }\boldsymbol{n}),\end{gather}
\begin{gather}\boldsymbol{t}\boldsymbol{\cdot }\boldsymbol{\sigma }\boldsymbol{\cdot }\boldsymbol{n} = 1,\end{gather}where
\begin{gather}\boldsymbol{\sigma } = \boldsymbol{\nabla }\boldsymbol{u} + {(\boldsymbol{\nabla }\boldsymbol{u})^\textrm{T}},\quad {\nabla _s} = (\boldsymbol{I} - \boldsymbol{nn})\boldsymbol{\cdot }\boldsymbol{\nabla },\quad {\nabla _s}\boldsymbol{\cdot }\boldsymbol{n} ={-} \frac{1}{{\sqrt {{E^3}} }}\frac{{{\partial ^2}\eta }}{{\partial {x^2}}},\end{gather}
\begin{gather}\boldsymbol{n} = \frac{1}{{\sqrt E }}\left( { - \frac{{\partial \eta }}{{\partial x}},1} \right),\quad \boldsymbol{t} = \frac{1}{{\sqrt E }}\left( {1,\frac{{\partial \eta }}{{\partial x}}} \right),\quad E = 1 + {\left( {\frac{{\partial \eta }}{{\partial x}}} \right)^2}.\end{gather}2.2. Long-wave approximation
To investigate the disturbance wave in a liquid film, we use the long-wave approximation in the analysis, assuming that the wavelength l of the disturbance is much greater than the average thickness d of the thin film, i.e. 
$\varpi = d/l \ll 1$. We choose l and d as the length scales in the x and z directions, respectively, and perform the following scale transformation on the parameters in the governing equation,
\begin{equation}x = {\varpi ^{ - 1}}X,z = Z,\quad u = U,w = \varpi W,\quad p = P,t = {\varpi ^{ - 1}}\tilde{t}.\end{equation}Thus, (2.3)–(2.6) can be written as
\begin{gather}\frac{{\partial U}}{{\partial X}} + \frac{{\partial W}}{{\partial Z}} = 0,\end{gather}
\begin{gather}\varpi Re\left( {\frac{{\partial U}}{{\partial \tilde{t}}} + U\frac{{\partial U}}{{\partial X}} + W\frac{{\partial U}}{{\partial Z}}} \right) ={-} \varpi \frac{{\partial P}}{{\partial X}} + {\varpi ^2}\frac{{{\partial ^2}U}}{{\partial {X^2}}} + \frac{{{\partial ^2}U}}{{\partial {Z^2}}},\end{gather}
\begin{gather}{\varpi ^2}Re\left( {\frac{{\partial W}}{{\partial \tilde{t}}} + U\frac{{\partial W}}{{\partial X}} + W\frac{{\partial W}}{{\partial Z}}} \right) ={-} \frac{{\partial P}}{{\partial Z}} + {\varpi ^3}\frac{{{\partial ^2}W}}{{\partial {X^2}}} + \varpi \frac{{{\partial ^2}W}}{{\partial {Z^2}}} - G.\end{gather}
The boundary condition at the wall 
$(z = 0)$ is
\begin{equation}\boldsymbol{u} = (U,W) = 0,\end{equation}
while those on the free surface 
$(z = \eta (x\textrm{,}t))$ are
\begin{gather}W = \frac{{\partial \eta }}{{\partial \tilde{t}}} + U\frac{{\partial \eta }}{{\partial X}},\end{gather}
\begin{gather}\begin{aligned}[b]& {-}\left( {1 + {\varpi^2}{{\left( {\dfrac{{\partial \eta }}{{\partial X}}} \right)}^2}} \right)P + 2{\varpi ^3}{\left( {\dfrac{{\partial \eta }}{{\partial X}}} \right)^2}\dfrac{{\partial U}}{{\partial X}} - 2\varpi \dfrac{{\partial \eta }}{{\partial X}}\left( {\dfrac{{\partial U}}{{\partial Z}} + {\varpi^2}\dfrac{{\partial W}}{{\partial X}}} \right)\\ &\quad + 2\varpi \dfrac{{\partial W}}{{\partial Z}} = C{a^{ - 1}}{\varpi ^2}\dfrac{{{\partial ^2}\eta }}{{\partial {X^2}}} \cdot {\left( {1 + {\varpi^2}{{\left( {\dfrac{{\partial \eta }}{{\partial X}}} \right)}^2}} \right)^{ - 1/2}}, \end{aligned}\end{gather}
\begin{gather}2{\varpi ^2}\frac{{\partial \eta }}{{\partial X}}\left( {\frac{{\partial W}}{{\partial Z}} - \frac{{\partial U}}{{\partial X}}} \right) + \left( {1 - {\varpi^2}{{\left( {\frac{{\partial \eta }}{{\partial X}}} \right)}^2}} \right)\left( {{\varpi^2}\frac{{\partial W}}{{\partial X}} + \frac{{\partial U}}{{\partial Z}}} \right) = \left( {1 + {\varpi^2}{{\left( {\frac{{\partial \eta }}{{\partial X}}} \right)}^2}} \right).\end{gather}2.2.1. Governing equation at small Reynolds numbers
Then, we consider the flow at small Reynolds numbers Re ≤ O(1). Meanwhile, 
$S \gg 1$. Thus, it can be seen from (2.1) that the influence of surface tension on film flow is much greater than that of inertial force. We seek the solutions of velocity and pressure as a perturbation series in powers of the small parameter 
$\varpi$,
\begin{equation}(U,W,P) = ({U_0},{W_0},{P_0}) + \varpi ({U_1},{W_1},{P_1}) + O({\varpi ^2}).\end{equation}
Substituting (2.11) into (2.8)–(2.10), the approximate equations for various orders of 
$\varpi$ can be obtained. The zero-order approximation is
\begin{equation}\frac{{\partial {U_0}}}{{\partial X}} + \frac{{\partial {W_0}}}{{\partial Z}} = 0,\quad \frac{{{\partial ^2}{U_0}}}{{\partial {Z^2}}} = 0,\quad \frac{{\partial {P_0}}}{{\partial Z}} ={-} G.\end{equation}
The boundary condition at 
$z = 0$ is
\begin{equation}\boldsymbol{u} = ({U_0},{W_0}) = 0,\end{equation}
while those at 
$z = \eta (x\textrm{,}t)$ are
\begin{equation}{P_0} ={-} {\overline {Ca} ^{ - 1}}\frac{{{\partial ^2}\eta }}{{\partial {X^2}}},\quad \frac{{\partial {U_0}}}{{\partial Z}} = 1.\end{equation}
Considering the dominant role of surface tension in the hydrodynamics of film flow, we introduce 
${\overline {Ca} ^{ - 1}} = {\varpi ^2}C{a^{ - 1}}$ in the above equation. So the solutions of (2.12)–(2.14) can be determined as follows:
\begin{equation}{U_0} = Z,\quad {W_0} = 0,\quad {P_0} = G(\eta - Z) - {\overline {Ca} ^{ - 1}}\frac{{{\partial ^2}\eta }}{{\partial {X^2}}}.\end{equation}Then, we pay attention to the first-order approximation. The governing equations are
\begin{gather}\frac{{\partial {U_1}}}{{\partial X}} + \frac{{\partial {W_1}}}{{\partial Z}} = 0,\end{gather}
\begin{gather}Re\left( {\frac{{\partial {U_0}}}{{\partial \tilde{t}}} + {U_0}\frac{{\partial {U_0}}}{{\partial X}} + {W_0}\frac{{\partial {U_0}}}{{\partial Z}}} \right) ={-} \frac{{\partial {P_0}}}{{\partial X}} + \frac{{{\partial ^2}{U_1}}}{{\partial {Z^2}}},\end{gather}
\begin{gather}\frac{{\partial {P_1}}}{{\partial Z}} = \frac{{{\partial ^2}{W_0}}}{{\partial {Z^2}}}.\end{gather}
The boundary condition at 
$z = 0$ is
\begin{equation}\boldsymbol{u} = ({U_1},{W_1}) = 0,\end{equation}
while those at 
$z = \eta (x\textrm{,}t)$ are
\begin{equation}{P_1} = 2\left( {\frac{{\partial {W_0}}}{{\partial Z}} - \frac{{\partial {U_0}}}{{\partial Z}}\frac{{\partial \eta }}{{\partial X}}} \right),\quad \frac{{\partial {U_1}}}{{\partial Z}} = 0.\end{equation}Then, the solutions of (2.16)–(2.18) can be derived as follows:
\begin{gather}{U_1} = \left( {\frac{1}{2}{Z^2} - \eta Z} \right)\left( {G\frac{{\partial \eta }}{{\partial X}} - {{\overline {Ca} }^{ - 1}}\frac{{{\partial^3}\eta }}{{\partial {X^3}}}} \right),\end{gather}
\begin{gather}{W_1} = \frac{1}{2}\frac{{\partial \eta }}{{\partial X}}{Z^2}\left( {G\frac{{\partial \eta }}{{\partial X}} - {{\overline {Ca} }^{ - 1}}\frac{{{\partial^3}\eta }}{{\partial {X^3}}}} \right) + \left( { - \frac{1}{6}{Z^3} + \frac{1}{2}\eta {Z^2}} \right)\left( {G\frac{{{\partial^2}\eta }}{{\partial {X^2}}} - {{\overline {Ca} }^{ - 1}}\frac{{{\partial^4}\eta }}{{\partial {X^4}}}} \right),\end{gather}
\begin{gather}{P_1} ={-} 2\frac{{\partial \eta }}{{\partial X}}.\end{gather}Substituting (2.15) and (2.19) into (2.10a), we can obtain
\begin{equation}\frac{{\partial \eta }}{{\partial \tilde{t}}} + \frac{{\partial \eta }}{{\partial X}}\eta - \varpi \frac{\partial }{{\partial X}}\left( {\frac{1}{3}{\eta^3}\left( {G\frac{{\partial \eta }}{{\partial X}} - {{\overline {Ca} }^{ - 1}}\frac{{{\partial^3}\eta }}{{\partial {X^3}}}} \right)} \right) = 0.\end{equation}This equation is consistent with the corresponding situation of Oron et al. (Reference Oron, Davis and Bankoff1997).
2.2.2. Linear stability analysis
We conduct the linear stability analysis for (2.20). A small normal mode disturbance is introduced to the basic state as follows:
\begin{equation}\eta = 1 + A\exp [\textrm{i}(kx - \omega t)],\end{equation}
where A is the amplitude, 
$\omega $ is the frequency and k is the wavenumber. The following dispersion relation can be derived:
\begin{equation}\omega = k - \textrm{i} \cdot {\textstyle{1 \over 3}}{k^2}(G + {k^2}{\overline {Ca} ^{ - 1}}).\end{equation}
When G > 0 (corresponding to figure 1a), 
${\mathop{\rm Im}\nolimits} (\omega ) < 0$, meaning that the flow is linearly stable. When G < 0 (corresponding to figure 1b), the flow is linearly unstable if
\begin{equation}{k^2} < k_c^2 = |G|\cdot \overline {Ca} = \rho |g|{d^2}/\sigma .\end{equation}This agrees with the result of Oron et al. (Reference Oron, Davis and Bankoff1997).
2.2.3. Weakly nonlinear stability
If we perform a weak nonlinear expansion on (2.20), the famous KS equation can be determined as follows: let 
$\eta = 1 + \varpi \eta ^{\prime}$, 
$t^{\prime} = \varpi \tilde{t}$, 
$\xi = X - \tilde{t}$, we can obtain
\begin{equation}\left. {\begin{array}{@{}c@{}} {\dfrac{{\partial \eta }}{{\partial \tilde{t}}} = \dfrac{{\partial \eta }}{{\partial \xi }} \cdot \dfrac{{\partial \xi }}{{\partial \tilde{t}}} + \dfrac{{\partial \eta }}{{\partial t^{\prime}}} \cdot \dfrac{{\partial t^{\prime}}}{{\partial \tilde{t}}} ={-} \varpi \dfrac{{\partial \eta^{\prime}}}{{\partial \xi }} + {\varpi^2}\dfrac{{\partial \eta^{\prime}}}{{\partial t^{\prime}}},}\\ {\dfrac{{\partial \eta }}{{\partial X}} = \dfrac{{\partial \eta }}{{\partial \xi }} \cdot \dfrac{{\partial \xi }}{{\partial X}} = \varpi \dfrac{{\partial \eta^{\prime}}}{{\partial \xi }},\quad \dfrac{{{\partial^2}\eta }}{{\partial {X^2}}} = \varpi \dfrac{{{\partial^2}\eta^{\prime}}}{{\partial {\xi^2}}}.} \end{array}} \right\}\end{equation}By substituting (2.24) into (2.20) and neglecting the third-order term, we can derive
\begin{equation}\frac{{\partial \eta ^{\prime}}}{{\partial t^{\prime}}} + \eta ^{\prime}\frac{{\partial \eta ^{\prime}}}{{\partial \xi }} - \frac{1}{3}\left( {G\frac{{{\partial^2}\eta^{\prime}}}{{\partial {\xi^2}}} - {{\overline {Ca} }^{ - 1}}\frac{{{\partial^4}\eta^{\prime}}}{{\partial {\xi^4}}}} \right) = 0.\end{equation}If the direction of gravity is the same as the wall-normal direction (figure 1b), then G < 0, and the following substitution is made:
\begin{equation}\xi = \sqrt { - \frac{1}{{\overline {Ca} \cdot G}}} \chi ,\quad \eta ^{\prime} = \frac{{ - 4G}}{3}\sqrt { - \overline {Ca} \cdot G} H,\quad t^{\prime} = \frac{3}{{\overline {Ca} {G^2}}}T.\end{equation}Thus, the KS equation can be obtained as follows:
\begin{equation}\frac{{\partial H}}{{\partial T}} + 4H\frac{{\partial H}}{{\partial \chi }} + \frac{{{\partial ^2}H}}{{\partial {\chi ^2}}} + \frac{{{\partial ^4}H}}{{\partial {\chi ^4}}} = 0.\end{equation}Since there has been extensive research on this equation (Chang et al. Reference Chang, Demekhin and Kopelevich1993; Chang & Demekhin Reference Chang and Demekhin2002), we will no longer discuss it separately. Instead, our primary focus will be on the solutions of (2.20) that exhibit nonlinearity of significant magnitude.
2.2.4. Traveling wave
When there is a traveling wave solution of (2.20), we can perform Galilean transformation on the equation and introduce a new coordinate system 
$x^{\prime} = X - C\tilde{t}$, 
$t^{\prime} = \tilde{t}$, so 
$\partial /\partial \tilde{t} ={-} C(\partial /\partial x^{\prime}) + \partial /\partial t^{\prime}$, 
$\partial /\partial X = \partial /\partial x^{\prime}$. Let
\begin{equation}{\textstyle{1 \over 3}}\varpi G = \hat{G},\quad {\textstyle{1 \over 3}}\varpi {\overline {Ca} ^{ - 1}} = Q,\end{equation}then the following equation can be obtained:
\begin{equation}\frac{{\partial \eta }}{{\partial t^{\prime}}} + (\eta - C)\frac{{\partial \eta }}{{\partial x^{\prime}}} - \hat{G}\frac{\partial }{{\partial x^{\prime}}}\left( {{\eta^3}\frac{{\partial \eta }}{{\partial x^{\prime}}}} \right) + Q\frac{\partial }{{\partial x^{\prime}}}\left( {{\eta^3}\frac{{{\partial^3}\eta }}{{\partial {{x^{\prime}}^3}}}} \right) = 0.\end{equation} It can be seen from (2.28b) that 
$Q \ge 0$. Thus, we consider two cases. When Q = 0, (2.29) can be simplified as follows:
\begin{equation}\frac{{\partial \eta }}{{\partial t^{\prime}}} + (\eta - C)\frac{{\partial \eta }}{{\partial x^{\prime}}} - \hat{G}\frac{\partial }{{\partial x^{\prime}}}\left( {{\eta^3}\frac{{\partial \eta }}{{\partial x^{\prime}}}} \right) = 0.\end{equation}
When Q > 0, let 
$x^{\prime} = a\varsigma $, 
$t^{\prime} = a\tilde{\tau }$, where 
$a = {Q^{1/3}}$, then (2.26) can be simplified as follows:
\begin{equation}\frac{{\partial \eta }}{{\partial \tilde{\tau }}} + (\eta - C)\frac{{\partial \eta }}{{\partial \varsigma }} - \tilde{G}\frac{\partial }{{\partial \varsigma }}\left( {{\eta^3}\frac{{\partial \eta }}{{\partial \varsigma }}} \right) + \frac{\partial }{{\partial \varsigma }}\left( {{\eta^3}\frac{{{\partial^3}\eta }}{{\partial {\varsigma^3}}}} \right) = 0,\end{equation}
where 
$\tilde{G} = \hat{G}/a$.
For the sake of writing convenience, we will unify (2.30a) and (2.30b) in the following form:
\begin{equation}\frac{{\partial \eta }}{{\partial t}} + (\eta - C)\frac{{\partial \eta }}{{\partial x}} - \tilde{G}\frac{\partial }{{\partial x}}\left( {{\eta^3}\frac{{\partial \eta }}{{\partial x}}} \right) + Q\frac{\partial }{{\partial x}}\left( {{\eta^3}\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right) = 0,\end{equation}where Q = 0,1.
We restrict our attention to two types of travelling waves: the first one is the solitary wave, which has
\begin{equation}x \to \pm \infty ,\quad \frac{{\partial \eta }}{{\partial x}} \to 0.\end{equation}The second one is the periodic wave, which has
\begin{equation}{\left. {\frac{{{\partial^n}\eta }}{{\partial {x^n}}}} \right|_{x = 0}} = {\left. {\frac{{{\partial^n}\eta }}{{\partial {x^n}}}} \right|_{x = {\mathop{T}\limits^{\frown}}}},\quad n = 0,1,2 \ldots .\end{equation}
Here, 
${\mathop{T}\limits^{\frown}}$ represents the period in the spatial direction.
The finite difference method is used in our numerical simulation. In Appendices A and B, we show the finite difference schemes for two types of waves. The validation of grid convergence is presented in Appendix C.
3. Numerical results
3.1. 
$\tilde{G} = 0$
First, we consider the case at 
$\tilde{G} = 0$. When Q = 0, (2.31) becomes the Burgers equation. Although its properties are well known, to make comparisons with other situations, we present the evolution of 
$\eta $ for the Burgers equation in figure 3. The initial conditions are 
$\eta |_{t = 0}\,{=}\,1 + 0.315\exp ( - 2.5 \cdot (x - 9)^2) \cdot \sin [2(x - 9)]$ for figure 3(a) and 
$\eta|_{t = 0}\,{=}\,1 + 0.15\sin x$ for figure 3(b). It can be seen that 
$|\partial \eta /\partial x|$ increases rapidly and tends to infinity within a finite time in some places. Then the equation will become singular.
Evolution of 
$\eta $ at Q = 0, C = 1 and 
$\tilde{G} = 0$: (a) the solitary wave; (b) the periodic wave.
For Q = 1, we can see in figure 4 that the surface tension effectively suppresses the disturbance. Here, the initial conditions of figure 4(a,b) are identical to those of figures 3(a) and 3(b), respectively. For figures 4(c) and 4(d), their initial conditions are 
$\eta {|_{t = 0}} = 1 + 0.58\exp ( - 10{(x - 9)^2}) \cdot \sin [2(x - 9)]$ and 
$\eta {|_{t = 0}} = 1 + 0.15\sin (2x)$, respectively. These results indicate that as the wavenumber increases, the disturbance decays more rapidly.
Evolution of 
$\eta $ at Q = 1, C = 1 and 
$\tilde{G} = 0$: (a,c) the solitary wave; (b,d) the periodic wave.
This can be attributed to the additional pressure exerted on liquids due to surface tension, which is described by the Young–Laplace equation: 
$\Delta p = \sigma ({k_1} + {k_2})$. Here, k 1 and k 2 are two principal curvatures of the surface. This pressure can push the liquid at the wave crest towards its equilibrium position, preventing the generation of singularities. As the wavenumber increases, the principal curvature at the wave crest also increases, leading to a faster decay of perturbations. The same trend holds true for disturbances with larger amplitudes 
$\textrm{(max|}\eta - 1\textrm{|} = 0.3)$.
3.2. $\tilde{G} > 0,\;Q = 0$
Then, we examine the effect of gravity at Q = 0 and 
$\tilde{G} > 0$, meaning that the direction of gravity is opposite to the wall-normal direction. Figure 5 shows the evolution of 
$\eta $ at Q = 0, C = 1 and 
$\tilde{G} = 1$. The initial conditions of figures 5(a) and 5(b) are identical to those of figures 3(a) and 3(b), respectively. It can be found that the perturbation is stabilized by the gravity at 
$\tilde{G} > 0$. This result agrees with the experiment reported by Wang et al. (Reference Wang, Zhang, Gong, Bai and Ma2017), where a thin liquid film is sheared by the co-flowing air from above and heated from below in a horizontal aluminium channel.
Evolution of 
$\eta $ at 
$Q = 0$, C = 1 and 
$\tilde{G} = 1$: (a) the solitary wave; (b) the periodic wave.
Further research indicates that the greater the gravity, the stronger the suppressive effect on disturbances. In addition, an increase in wavenumber results in a faster decay of disturbances, a trend similar to that observed in figure 4.
To a certain extent, the problem we examined is similar to the cases of Kelvin–Helmholtz instability and Rayleigh–Taylor instability. In these scenarios, two uniform fluids are separated by a horizontal boundary at 
$\eta = 1$. The lower fluid has a density of 
${\rho _1}$, while the upper fluid has a density of 
${\rho _2}$. Gravity acts vertically downwards in this setup. The two fluids are streaming with the constant velocities 
${U_1}$ and 
${U_2}$.
In the linear stability analysis (2.21), the result of interfacial instability for two inviscid fluids is (Chandrasekhar Reference Chandrasekhar2013)
\begin{align}\omega = \frac{k}{{{\rho _1} + {\rho _2}}}\left[ {({\rho_1}{U_1} + {\rho_2}{U_2}) \pm \sqrt {\frac{{{\rho_1} + {\rho_2}}}{k}[\sigma {k^2} + ({\rho_1} - {\rho_2})g] - {\rho_1}{\rho_2}{{({U_1} - {U_2})}^2}} } \right].\end{align}Thus, instabilities appear when
\begin{equation}\frac{{{\rho _1} + {\rho _2}}}{k}[\sigma {k^2} + ({\rho _1} - {\rho _2})g] - {\rho _1}{\rho _2}{({U_1} - {U_2})^2} < 0.\end{equation}
The instability arising from differences in velocity 
$({U_1} - {U_2})$ is known as the Kelvin–Helmholtz instability, while that resulting from differences in density is called the Rayleigh–Taylor instability. It is evident that surface tension always serves to stabilize, whereas gravity acts as a destabilizing force only when the density of the upper fluid exceeds the density of the lower fluid: 
$({\rho _1} - {\rho _2}) < 0$.
Since these conclusions are derived from linear stability analysis for inviscid fluids, they may not be directly applicable to nonlinear waves in a sheared liquid film for viscous fluids. Nonetheless, the aforementioned outcomes suggest that the influence of 
$Q,\tilde{G}$ on the instability bears similarities.
3.3. $\tilde{G} < 0$
The results above indicate that disturbances are suppressed by the surface tension (Q = 1), whereas the gravity at 
$\tilde{G} < 0$ is destabilizing. Therefore, we can speculate that if a balance is reached between these two factors, the disturbance waves may persist. Further analysis confirms this possibility.
Suppose (2.31) has a steady-state solution, meaning 
$\partial \eta /\partial t = 0$, then the equation becomes
\begin{equation}(\eta - C)\frac{{\partial \eta }}{{\partial x}} - \tilde{G}\frac{\partial }{{\partial x}}\left( {{\eta^3}\frac{{\partial \eta }}{{\partial x}}} \right) + \frac{\partial }{{\partial x}}\left( {{\eta^3}\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right) = 0.\end{equation}Integrating (3.3) yields
\begin{equation}\frac{1}{2}{(\eta - C)^2} - \tilde{G}{\eta ^3}\frac{{\partial \eta }}{{\partial x}} + {\eta ^3}\frac{{{\partial ^3}\eta }}{{\partial {x^3}}} = {C_1},\end{equation}
where 
${C_1}$ is a constant.
3.3.1. Solitary waves
First, we consider solitary waves. Let
\begin{equation}\eta ^{\prime} = \frac{{\partial \eta }}{{\partial x}},\quad \eta ^{\prime\prime} = \frac{{{\partial ^\textrm{2}}\eta }}{{\partial {x^\textrm{2}}}},\quad \eta ^{\prime\prime\prime} = \frac{{{\partial ^3}\eta }}{{\partial {x^3}}}.\end{equation}
When 
$x \to - \infty $, 
$\eta \to {\eta _\infty }$, 
$\eta ^{\prime},\eta ^{\prime\prime\prime} \to 0$, then
\begin{equation}{\eta ^3}( - \tilde{G}\eta ^{\prime} + \eta ^{\prime\prime\prime}) = {\textstyle{1 \over 2}}[{({\eta _\infty } - C)^2} - {(\eta - C)^2}].\end{equation}There are two stationary points in (3.6), which are
\begin{equation}{\kappa _1} = {\eta _\infty },\quad {\kappa _2} = 2C - {\eta _\infty }.\end{equation}
Now we discuss the stability of these stationary points. By introducing the phase space coordinates
$({\eta _1},{\eta _2},{\eta _3})$, where 
${\eta _1} = \eta ,{\eta _2} = \eta ^{\prime},{\eta _3} = \eta ^{\prime\prime}$, we can derive
\begin{equation}{\eta ^{\prime}_1} = {\eta _2},\quad {\eta ^{\prime}_2} = {\eta _3},\quad {\eta ^{\prime}_3} = \frac{{[{{({\eta _\infty } - C)}^2} - {{({\eta _1} - C)}^2}]}}{{2\eta _1^3}} + \tilde{G}{\eta _2}.\end{equation}
We regard 
${\eta ^{\prime}_i}(i = 1,2,3)$ as a velocity field in the phase space, then the incompressibility condition is satisfied,
\begin{equation}
\boldsymbol{\nabla }\boldsymbol{\cdot }\eta ^{\prime} = \frac{{\partial{\eta^{\prime}_1}}}{{\partial {\eta _1}}} + \frac{{\partial
{\eta ^{\prime}_2}}}{{\partial {\eta _2}}} + \frac{{\partial {\eta ^{\prime}_3}}}{{\partial{\eta_3}}} = 0.\end{equation}
The two stationary points in the phase space are 
${O_1}({\kappa _1},0,0)$ and 
${O_2}({\kappa _2},0,0)$. Then, we consider an asymptotic solution near a stationary point as follows:
\begin{equation}\eta = \bar{\eta } + \tilde{\eta },\end{equation}
where 
$\bar{\eta } = {\kappa _1},{\kappa _2}$, 
$\tilde{\eta } = A\exp (\lambda x)$, A is a small amplitude and 
$\lambda $ is an eigenvalue. Substituting (3.10) into (3.8) and retaining the first-order small quantity, the eigenvalue equation can be derived.
For 
${O_1}({\kappa _1},0,0)$, its eigenvalue equation is
\begin{equation}{\lambda ^3} - \tilde{G}\lambda + \frac{{{\eta _\infty } - C}}{{\eta _\infty ^3}} = 0.\end{equation}The roots of (3.11) are
\begin{align}{\lambda _1} = m + n,\quad {\lambda _2} ={-} \frac{1}{2}(m + n) + \frac{{\sqrt 3 }}{2}i(m - n),\quad {\lambda _3} ={-} \frac{1}{2}(m + n) + \frac{{\sqrt 3 }}{2}i(n - m),\end{align}where
\begin{gather}m = \sqrt[3]{{ - \frac{{{\eta _\infty } - C}}{{2\eta _\infty ^3}} + \sqrt {{{\left( {\frac{{{\eta_\infty } - C}}{{2\eta_\infty^3}}} \right)}^2} - {{\left( {\frac{{\tilde{G}}}{3}} \right)}^3}} }},\end{gather}
\begin{gather}n = \sqrt[3]{{ - \frac{{{\eta _\infty } - C}}{{2\eta _\infty ^3}} - \sqrt {{{\left( {\frac{{{\eta_\infty } - C}}{{2\eta_\infty^3}}} \right)}^2} - {{\left( {\frac{{\tilde{G}}}{3}} \right)}^3}} }}.\end{gather}
As 
$\tilde{G} < 0$, we can derive that 
$m > 0,n < 0$, 
${\mathop{\rm Im}\nolimits} ({\lambda _2}) > 0$ and 
${\mathop{\rm Im}\nolimits} ({\lambda _3}) < 0$. When 
${\eta _\infty } - C > 0$, we have 
${\lambda _1} < 0$ and 
$Re ({\lambda _2}) = Re ({\lambda _3}) > 0$; while 
${\eta _\infty } - C < 0$, we have 
${\lambda _1} > 0$ and 
$Re ({\lambda _2}) = Re ({\lambda _3}) < 0$.
For 
${O_2}({\kappa _2},0,0)$, its eigenvalue equation is
\begin{equation}{\lambda ^3} - \tilde{G}\lambda + \frac{{{\eta _\infty } - C}}{{{{({\eta _\infty } - 2C)}^3}}} = 0.\end{equation}The roots of (3.14) are
\begin{align}{\mu _1} = a + b,\quad {\mu _2} ={-} \frac{1}{2}(a + b) + \frac{{\sqrt 3 }}{2}i(a - b),\quad {\mu _3} ={-} \frac{1}{2}(a + b) + \frac{{\sqrt 3 }}{2}i(b - a),\end{align}where
\begin{gather}a = \sqrt[3]{{\frac{{C - {\eta _\infty }}}{{2{{({\eta _\infty } - 2C)}^3}}} + \sqrt {{{\left( {\frac{{{\eta_\infty } - C}}{{2{{({\eta_\infty } - 2C)}^3}}}} \right)}^2} - {{\left( {\frac{{\tilde{G}}}{3}} \right)}^3}} }},\end{gather}
\begin{gather}b = \sqrt[3]{{\frac{{C - {\eta _\infty }}}{{2{{({\eta _\infty } - 2C)}^3}}} - \sqrt {{{\left( {\frac{{{\eta_\infty } - C}}{{2{{({\eta_\infty } - 2C)}^3}}}} \right)}^2} - {{\left( {\frac{{\tilde{G}}}{3}} \right)}^3}} }}.\end{gather}
The existence of O 2 requires that 
$2C - {\eta _\infty } > 0$. Using the condition 
$\tilde{G} < 0$, we can derive that 
$a > 0,b < 0$, 
${\mathop{\rm Im}\nolimits} ({\mu _2}) > 0$ and 
${\mathop{\rm Im}\nolimits} ({\mu _3}) < 0$. When 
${\eta _\infty } - C > 0$, we have 
${\mu _1} > 0$ and 
$Re ({\mu _2}) = Re ({\mu _3}) < 0$; while when 
${\eta _\infty } - C < 0$, we have 
${\mu _1} < 0$ and 
$Re ({\mu _2}) = Re ({\mu _3}) > 0$.
When 
${\eta _\infty } = C$, 
${\lambda _\textrm{1}} = Re ({\lambda _2}) = Re ({\lambda _3}) = 0$, and two stationary points overlap. The trajectory near the stationary point will not be attracted or repelled. Therefore, we will not consider this case. In addition, when 
$2C - {\eta _\infty } \le 0$, there is only one stationary point. However, solitary waves are not found in this condition. Therefore, we restrict our attention to the region 
$C > {\textstyle{1 \over 2}}{\eta _\infty },\;C \ne {\eta _\infty }$ in the following. In this case, both stationary points have eigenvalues with positive real parts, indicating their absolute instability.
We provide an example at 
$C = {\textstyle{5 \over 6}}{\eta _\infty },\;{\eta _\infty } = 1$ to illustrate the behaviour of trajectories near stable points, while the case at 
$C > {\eta _\infty }$ is similar. In the vicinity of the stationary point O 1, the disturbance exhibits exponential convergence,
\begin{equation}\eta = \bar{\eta } + {A_1}\exp ({\lambda _1}x),\end{equation}
where 
${\lambda _1} < 0$. However, the manner in which the disturbance diverges exhibits oscillatory behaviour,
\begin{equation}\eta = \bar{\eta } + {A_2}\exp (Re ({\lambda _2})x)\cos ({\mathop{\rm Im}\nolimits} ({\lambda _2})x + \varphi ),\end{equation}
where 
$Re {\lambda _2} > 0$ and 
$\varphi $ is the phase angle. In contrast, in the vicinity of the stationary point O 2, the disturbance exhibits oscillatory convergence,
\begin{equation}\eta = \bar{\eta } + {A_3}\exp(Re ({\mu _\textrm{2}})x\textrm{)}\cos ({\mathop{\rm Im}\nolimits} ({\mu _\textrm{2}})x + \varphi ),\end{equation}
where 
$Re ({\mu _2}) < 0$, while it diverges exponentially,
\begin{equation}\eta = \bar{\eta } + {A_4}\exp ({\mu _1}x),\end{equation}
where 
${\mu _1} > 0$. It should be noted that 
${A_i}(i = 1\sim 4)$ are real numbers.
We select points near the stationary point O 1 as initial points and employ numerical integration methods (refer to Appendix D) to investigate the phase trajectory of (3.9). The results indicate that the evolution of 
$\eta $ shows a quasi-periodic oscillation pattern between O 1 and O 2. The specific process is displayed in figure 6.
Evolution of 
$\eta $, suggesting heteroclinic trajectories between two stationary points O 1 and O 2: (a) the sketch; (b) the numerical result at 
$\tilde{G} ={-} 3$, Q = 1, 
${\eta _\infty } = 1$ and 
$C = {\textstyle{5 \over 6}}{\eta _\infty }$.
Initially, 
$\eta $ oscillates near O 1 with an increasing amplitude, demonstrating the presence of the unstable two-dimensional manifold 
$W_1^u$ of O 1. Subsequently, it oscillates towards O 2 with a decreasing amplitude, representing the stable two-dimensional manifold 
$W_\textrm{2}^s$ of O 2. After reaching a certain distance, 
$\eta $ exponentially deviates upwards from O 2, corresponding to the unstable one-dimensional manifold 
$W_\textrm{2}^u$ of O 2. Later, it exponentially converges towards O 1 from below, signifying the stable one-dimensional manifold 
$W_\textrm{1}^s$ of O 1, and after a certain distance, this pattern repeats itself as the system oscillates away from O 1 again. Consequently, the aforementioned situation suggests the existence of heteroclinic trajectories between the two stationary points O 1 and O 2, as illustrated in figure 7.
Heteroclinic contour of the stationary points O 1 and O 2.
In the following, we discuss the impact of 
$\tilde{G}$ on the evolution of 
$\eta $. The computation reveals that there are quasi-periodic solutions (see figure 8) when 
$\tilde{G}$ falls within a suitable region. However, once 
$\tilde{G}$ surpasses the allowable range, the liquid film will break 
$(\eta \to 0)$, as shown in figure 9. The approximate range of 
$\tilde{G}$ at 
$C = {\textstyle{5 \over 6}},{\eta _\infty } = 1$ is 
$[ - 10, - 3]$, while the specific region is contingent upon the initial disturbance.
Quasi-periodic solution of 
$\eta $ at 
$\tilde{G} ={-} 5$, Q = 1, 
${\eta _\infty } = 1$ and 
$C = {\textstyle{5 \over 6}}$.
Evolution of 
$\eta $ at 
$\tilde{G} ={-} 2$, Q = 1, 
${\eta _\infty } = 1$ and 
$C = {\textstyle{5 \over 6}}$.
The allowable range of 
$\tilde{G}$ at different C is displayed in figure 10 where 
${\eta _\infty } = 1$. As 
$2C - {\eta _\infty } > 0$, we have 
$C > 0.5$. For fixed C, stationary solitary waves could only appear when 
$( - \tilde{G}) \in [Mi,Mx]$. Once 
$\tilde{G} < 0$, we can find an allowable range of wave speed C.
Allowable range of 
$\tilde{G}$ for stationary solitary waves at different C. Here, 
${\eta _\infty } = 1$.
To examine the properties of quasi-periodic solutions, we compute 3000 cycles at 
$\tilde{G} ={-} \textrm{5}$ for different 
${\eta _\infty }$. The results in figure 11 indicate that the cycle length (Cl) steadily decreases and converges to a constant value, while the maximum value of 
$\eta $ within one cycle 
$({\eta _{max}})$ remains nearly constant.
Variations of (a) cycle length and (b) 
${\eta _{max}}$ in 3000 cycles at 
$\tilde{G} ={-} 5$, Q = 1, 
${\eta _\infty } = 1$ and 
$C = {\textstyle{5 \over 6}}{\eta _\infty }$.
However, we do not find any other structures of solitary waves that appeared in the KS equation, such as the homoclinic trajectory and the heteroclinic trajectory between a stationary point and a limit cycle (Chang & Demekhin Reference Chang and Demekhin2002). The reason may be related to the nonlinearity of the disturbance. The KS equation describes the weak nonlinear instability near the equilibrium position 
$\eta = 1$ with a disturbance amplitude of 
$|\varpi \eta ^{\prime}|\ll 1$. However, the amplitude of the disturbance reaches O(1) in this paper. The strong nonlinearity may lead to rapid and sustained growth of disturbances. In our computation, disturbances that continuously grow near an equilibrium point, if not promptly attracted by another equilibrium point, are prone to persistent divergence, ultimately leading to film rupture. Therefore, steady-state solitary waves are most likely to be found in the heteroclinic orbits between two equilibrium points.
3.3.2. Periodic waves
For periodic waves, we begin by conducting a linear stability analysis as (3.1). Substituting (3.1) into (2.31) and considering 
$|A|\ll 1$, we can determine that a steady-state solution exists at
\begin{equation}C = 1,\quad {k^2} ={-} \tilde{G}.\end{equation}
Despite the linear nature of the aforementioned analysis, numerical simulations demonstrate that these results remain qualitatively valid even at 
$|A|= \textrm{0}\textrm{.5}$.
We examine the evolution of a periodic wave with the initial condition 
$\eta {|_{t = 0}} = \textrm{1} + A\sin x$. The results show that although nonlinearity is significant when A is large, the distributions of 
$\eta $ still resemble the sine function (see figure 12).
Distributions of 
$\eta $ for periodic waves at t = 30: (I) 
$\tilde{G} ={-} 1$, C = 1.009; (II) 
$\tilde{G} ={-} 1.001$, C = 1.033; (III) 
$\tilde{G} ={-} 1$, C = 1.085. At initial time, 
$\eta {|_{t = 0}} = 1 + A\sin x$, while A = 0.15, 0.3 and 0.5 for the three respective cases.
However, the amplitude A varies with time t. Figure 13 illustrates that the amplitude increases at 
$\tilde{G} ={-} \textrm{1}\textrm{.1}$, but decreases at 
$\tilde{G} ={-} \textrm{1}\textrm{.0005}$. Therefore, a steady-state solution may exist when 
$\tilde{G}$ falls between these two values.
Variation of amplitude A with time t. The initial condition is 
$\eta {|_{t = 0}} = 1 + 0.3\sin x$.
Additionally, we find other wave patterns for steady-state periodic waves. Figure 14 illustrates one such pattern, where the curves are observed in a stationary coordinate system. The wave speed of these periodic waves is C = 0.76 and the period is T = 8.56.
Another pattern of steady-state periodic wave at 
$\tilde{G} ={-} 2$.
Further calculations indicate that the evolutions of periodic waves are sensitive to the initial conditions (including wave number, amplitude, superposition of multiple waves, etc.) as well as the parameter 
$\tilde{G}$. In some cases, the amplitude of the disturbance may continue to increase, thus the film thickness becomes very thin in some places, tending towards film rupture (for example, when 
$\tilde{G} ={-} 10$ and 
$\eta {|_{t = 0}} = \textrm{1} + \textrm{0}\textrm{.3}\sin x$). We will conduct a detailed study on this in our subsequent work to determine the parameter ranges for the existence of various patterns of steady-state periodic waves, and further assess the parameter conditions for film rupture.
4. Energy analysis
To analyse the evolution of disturbances from an energy perspective, we use 
$h = \eta - 1$ to represent the deviation of the liquid surface from its equilibrium position, and perform a Fourier expansion on it,
\begin{equation}h = \sum\limits_{n = 1}^\infty {[{A_n}(t)\cos (nkx) + {B_n}(t)\sin (nkx)]} .\end{equation}We set C = 1 in (2.31), meaning that we observe in a reference frame moving at a speed of 1, then (2.31) can be rewritten as
\begin{equation}\frac{{\partial h}}{{\partial t}} + h\frac{{\partial h}}{{\partial x}} - \tilde{G}\frac{\partial }{{\partial x}}\left( {{{(1 + h)}^3}\frac{{\partial h}}{{\partial x}}} \right) + Q\frac{\partial }{{\partial x}}\left( {{{(1 + h)}^3}\frac{{{\partial^3}h}}{{\partial {x^3}}}} \right) = 0.\end{equation}Multiplying (4.2) by h and integrating over a wavelength range yields
\begin{align}\frac{\partial }{{\partial t}}\int_{ - {\rm \pi}/k}^{{\rm \pi}/k} {\left( {\frac{1}{2}{h^2}} \right) \cdot \textrm{d}\kern0.7pt x} ={-} \int_{ - {\rm \pi}/k}^{{\rm \pi} /k} {\left( {{h^2}\frac{{\partial h}}{{\partial x}} + h\frac{\partial }{{\partial x}}\left( {{{(1 + h)}^3}\left( {Q\frac{{{\partial^3}h}}{{\partial {x^3}}} - \tilde{G}\frac{{\partial h}}{{\partial x}}} \right)} \right)} \right) \cdot \textrm{d}\kern0.7pt x} .\end{align}By performing integration by parts and using the periodic conditions, we can obtain
\begin{equation}\frac{\partial }{{\partial t}}\int_{ - {\rm \pi}/k}^{{\rm \pi}/k} {\left( {\frac{1}{2}{h^2}} \right) \cdot \textrm{d}\kern0.7pt x} \textrm{ } = \int_{ - {\rm \pi}/k}^{{\rm \pi}/k} {{{(1 + h)}^3}\left( {Q\frac{{\partial h}}{{\partial x}}\frac{{{\partial^3}h}}{{\partial {x^3}}} - \tilde{G}{{\left( {\frac{{\partial h}}{{\partial x}}} \right)}^2}} \right) \cdot \textrm{d}\kern0.7pt x} = {P_1} + {P_2},\end{equation}where
\begin{gather}{P_1} = Q\int_{ - {\rm \pi}/k}^{{\rm \pi}/k} {{{(1 + h)}^3}\left( {\frac{{\partial h}}{{\partial x}}\frac{{{\partial^3}h}}{{\partial {x^3}}}} \right) \cdot \textrm{d}\kern0.7pt{x_2}} ,\end{gather}
\begin{gather}{P_2} ={-} \tilde{G}\int_{ - {\rm \pi}/k}^{{\rm \pi}/k} {{{(1 + h)}^3}{{\left( {\frac{{\partial h}}{{\partial x}}} \right)}^2} \cdot \textrm{d}\kern0.7pt x} ,\end{gather}
representing the contributions of surface tension and gravity to the energy of the disturbance, respectively. We only consider the case 
$1 + h = \eta > 0$, where the liquid film does not rupture.
For 
${P_2}$, as 
${(\partial h/\partial x)^2} \ge 0$, the sign of 
${P_2}$ is entirely dependent on 
$\tilde{G}$. Thus, when 
$\tilde{G} > 0(\tilde{G} < 0)$, gravity causes the disturbance to decay (grow).
For 
${P_1}$, if 
$|h|\ll 1$, meaning in the case of small disturbances, then
\begin{equation}{P_1} \approx Q\int_{ - {\rm \pi}/k}^{{\rm \pi}/k} {\left( {\frac{{\partial h}}{{\partial x}}\frac{{{\partial^3}h}}{{\partial {x^3}}}} \right) \cdot \textrm{d}\kern0.7pt {x_2}} ={-} Q\int_{ - {\rm \pi}/k}^{{\rm \pi}/k} {{{\left( {\frac{{{\partial^2}h}}{{\partial {x^2}}}} \right)}^2} \cdot \textrm{d}\kern0.7pt {x_2}} < 0.\end{equation}
If the disturbance is a monochromatic wave 
$h = {A_1}(t)\cos (kx)$, then
\begin{equation}{P_1} ={-} QA_1^2{k^4}\int_{ - {\rm \pi}/k}^{{\rm \pi}/k} {{{(1 + h)}^3}{{\sin }^2}kx \cdot \textrm{d}\kern0.7pt {x_2}} < 0.\end{equation}
Therefore, in both of these scenarios, surface tension must cause the disturbance to decay. For more general cases, it is not possible to directly determine the sign of 
${P_1}$. Our calculations do not find situations where the surface tension causes the disturbance to grow.
In summary, in the small-Reynolds-number regime studied in this paper, the only source of energy that can contribute to the growth of disturbances is gravity, and in this case, the liquid film must be positioned below a horizontal plate.
5. Numerical simulation
To compare with our findings, we perform a numerical simulation for the gas–liquid two-phase flow using the software Fluent, as presented in figure 15. The liquid (corn oil) film is sheared by an incompressible gas (air) flow. The system is assumed to be two-dimensional. All parameters are listed in table 1. We employ the PISO (pressure-implicit with splitting of operators) scheme to solve the flow equations, and use the VoF (volume of fluid) method for tracking moving interfaces. In figure 15(a), the top and bottom boundaries are solid walls, while the left and right sides employ periodic boundary conditions. In the horizontal direction, 3000 grids are used, while in the vertical direction, we set 38 grids for the gas section and 20 grids for the liquid section.
Numerical simulation of the gas–liquid two-phase flow: (a) the schematic and the evolution of the liquid film at (b) t = 1.13 s, (c) t = 2.79 s and (d) t = 3.79 s.
Parameters in the numerical simulation for the gas–liquid two-phase flow.
The simulation results show that the evolution process of the liquid film is sensitive to initial disturbances and the flow is unsteady. In figure 15(b–d), we display the evolution within a certain region at three different times. We can see a small disturbance at t = 1.13 s. The amplitude increases at t = 2.79 s and the long-wave approximation is satisfied. Additionally, the relative variation of the film thickness exceeds 30 % when t = 3.79 s. These results are qualitatively consistent with the theoretical analysis in figures 12 and 14, as well as the experimental finding presented in figure 2.
6. Conclusion
The current study focuses on investigating the nonlinear dynamics of a sheared liquid film on a horizontal plate using theoretical and numerical methods. Due to the low Reynolds number, the flow is primarily governed by surface tension while inertial forces playing a much smaller role. By applying the long-wave theory, we derive the governing equation for film thickness, and we employ finite difference schemes in our numerical simulations. The effects of gravity and surface tension on disturbances are examined, measured by two dimensionless parameters 
$\tilde{G}$ and Q, respectively. The evolutions of solitary and periodic waves are presented. While the shear stress on the surface drives the film flow, the magnitude and direction of gravity are still crucial for the instability.
When effects of gravity and surface tension are neglected 
$(\tilde{G} = Q = 0)$, the governing equation simplifies to the Burgers equation. In this case, singularities may occur at certain points within a finite time. However, by increasing the surface tension, disturbances can be suppressed, preventing the occurrence of singularities. Additionally, disturbances decay more rapidly at higher wavenumbers. When the direction of gravity is opposite to the wall-normal direction 
$(\tilde{G} > 0)$, the perturbation is stabilized by gravity.
In contrast, when the direction of gravity aligns with the wall-normal direction 
$(\tilde{G} < 0)$, gravity becomes destabilizing. Therefore, stationary travelling waves can exist when the effects of gravity and surface tension reach equilibrium. In the case of weakly nonlinearity, the Kuramoto–Sivashinsky equation is derived. For the steady-state solution of solitary waves, we find two stationary points, both of which are unstable. When 
$\tilde{G}$ falls within an appropriate range, quasi-periodic oscillations occur between these two points, suggesting the presence of heteroclinic trajectories. However, if 
$\tilde{G}$ exceeds the permissible range, the liquid film will break in some places. The precise boundaries of 
$\tilde{G}$ depend on the initial disturbance. The evolutions of periodic waves are sensitive to 
$\tilde{G}$ and initial disturbances, including wave number, amplitude, superposition of multiple waves, etc. Two patterns of the steady-state periodic waves are presented, one of which resembles a sine function even when nonlinearity is significant.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2024.895.
Funding
This work has been supported by the National Natural Science Foundation of China (No. 12372247) and Ningbo Municipality Key Research and Development Program (No. 2022Z213).
Declaration of interests
The authors report no conflict of interest.
Author contributions
K.-X.H. made substantial contributions to the conception of the work, wrote the paper for important intellectual content and approved the final version to be published. He is accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. K.D. made contributions to the numerical computation and validation. Q.-S.C. provided editing and writing assistance.
Appendix A. Finite difference schemes for solitary waves
In the finite difference calculation of (2.31), the temporal direction only has a first-order derivative, so we adopt a simple first-order accuracy scheme. In the calculation, reducing the time step can improve the time accuracy without affecting the computational stability. However, in the spatial direction, there are derivatives from first to fourth orders. Meanwhile, reducing the space step may lead to divergence in the calculation, thus a higher-order accuracy scheme is used. We adopt the following central difference schemes that are fourth-order accurate in the spatial direction:
\begin{gather}\left( {\frac{{\partial \eta }}{{\partial x}}} \right)_j^n = \frac{{ - \eta _{j + 2}^n + 8\eta _{j + 1}^n - 8\eta _{j - 1}^n + \eta _{j - 2}^n}}{{12\Delta x}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^2}\eta }}{{\partial {x^2}}}} \right)_j^n = \frac{{ - \eta _{j + 2}^n + 16\eta _{j + 1}^n - 30\eta _j^n + 16\eta _{j - 1}^n - \eta _{j - 2}^n}}{{12{{(\Delta x)}^2}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_j^n = \frac{{ - \eta _{j + 3}^n + 8\eta _{j + 2}^n - 13\eta _{j + 1}^n + 13\eta _{j - 1}^n - 8\eta _{j - 2}^n + \eta _{j - 3}^n}}{{8{{(\Delta x)}^3}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_j^n = \frac{{ - \eta _{j + 3}^n + 12\eta _{j + 2}^n - 39\eta _{j + 1}^n + 56\eta _j^n - 39\eta _{j - 1}^n + 12\eta _{j - 2}^n - \eta _{j - 3}^n}}{{6{{(\Delta x)}^4}}} + O(\Delta {x^4}).\end{gather}The superscript n = 1–M and subscript j = 1–N represent the grid numbers in the temporal and spatial directions, respectively. Meanwhile, we employ the following first-order accurate explicit difference scheme in the temporal direction:
\begin{gather}\eta _j^{n + 1} = \eta _j^n - \Delta t\left[ \begin{array}{@{}l@{}} (\eta_j^n - C)\left( {\dfrac{{\partial \eta }}{{\partial x}}} \right)_j^n - \tilde{G}\left( {{{(\eta_j^n)}^3}\left( {\dfrac{{{\partial^2}\eta }}{{\partial {x^2}}}} \right)_j^n + 3{{(\eta_j^n)}^2}{{\left( {\left( {\dfrac{{\partial \eta }}{{\partial x}}} \right)_j^n} \right)}^2}} \right)\\ \quad +\,Q\left( {{{(\eta_j^n)}^3}\left( {\dfrac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_j^n + 3{{(\eta_j^n)}^2}\left( {\dfrac{{\partial \eta }}{{\partial x}}} \right)_j^n\left( {\dfrac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_j^n} \right) \end{array} \right].\end{gather}The difference schemes for boundary conditions are given as follows.
At the left boundary, it can be seen from (A1)–(A4) that we need to provide the schemes at j = 1,2 for 
$(\partial \eta /\partial x)_j^n$ and 
$({\partial ^2}\eta /\partial {x^2})_j^n$, and the schemes at j = 1–3 for 
$({\partial ^3}\eta /\partial {x^3})_j^n$ and 
$({\partial ^4}\eta /\partial {x^4})_j^n$. Therefore, we adopt the following second-order accurate difference schemes for 
$(\partial \eta /\partial x)_j^n$ and 
$({\partial ^2}\eta /\partial {x^2})_j^n$:
\begin{gather}\left( {\frac{{\partial \eta }}{{\partial x}}} \right)_1^n = \frac{{ - 3\eta _1^n + 4\eta _2^n - \eta _3^n}}{{2\Delta x}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{\partial \eta }}{{\partial x}}} \right)_\textrm{2}^n = \frac{{\eta _\textrm{3}^n - \eta _1^n}}{{2\Delta x}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^2}\eta }}{{\partial {x^2}}}} \right)_1^n = \frac{{2\eta _1^n - 5\eta _2^n + 4\eta _3^n - \eta _4^n}}{{{{(\Delta x)}^2}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^2}\eta }}{{\partial {x^2}}}} \right)_\textrm{2}^n = \frac{{\eta _\textrm{3}^n - 2\eta _\textrm{2}^n + \eta _1^n}}{{{{(\Delta x)}^2}}} + O(\Delta {x^2}).\end{gather}
For 
$({\partial ^3}\eta /\partial {x^3})_j^n$ and 
$({\partial ^4}\eta /\partial {x^4})_j^n$, the second-order accurate difference schemes are given as
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_1^n = \frac{{ - 2.5\eta _1^n + 9\eta _2^n - 12\eta _3^n + 7\eta _4^n - 1.5\eta _5^n}}{{{{(\Delta x)}^3}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_2^n = \frac{{ - 1.5\eta _1^n + 5\eta _2^n - 6\eta _3^n + 3\eta _4^n - 0.5\eta _5^n}}{{{{(\Delta x)}^3}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_3^n = \frac{{ - 0.5\eta _1^n + \eta _2^n - \eta _4^n + 0.5\eta _5^n}}{{{{(\Delta x)}^3}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_1^n = \frac{{3\eta _1^n - 14\eta _2^n + 26\eta _3^n - 24\eta _4^n + 11\eta _5^n - 2\eta _6^n}}{{{{(\Delta x)}^4}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_2^n = \frac{{2\eta _1^n - 9\eta _2^n + 16\eta _3^n - 14\eta _4^n + 6\eta _5^n - \eta _6^n}}{{{{(\Delta x)}^4}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_3^n = \frac{{\eta _1^n - 4\eta _2^n + 6\eta _3^n - 4\eta _4^n + \eta _5^n}}{{{{(\Delta x)}^4}}} + O(\Delta {x^2}).\end{gather}Similarly, at the right boundary, we adopt the following schemes:
\begin{gather}\left( {\frac{{\partial \eta }}{{\partial x}}} \right)_N^n = \frac{{3\eta _N^n - 4\eta _{N - 1}^n + \eta _{N - 2}^n}}{{2\Delta x}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{\partial \eta }}{{\partial x}}} \right)_{N - 1}^n = \frac{{\eta _N^n - \eta _{N - 2}^n}}{{2\Delta x}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^2}\eta }}{{\partial {x^2}}}} \right)_N^n = \frac{{2\eta _N^n - 5\eta _{N - 1}^n + 4\eta _{N - 2}^n - \eta _{N - 3}^n}}{{{{(\Delta x)}^2}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^2}\eta }}{{\partial {x^2}}}} \right)_{N - 1}^n = \frac{{\eta _N^n - 2\eta _{N - 1}^n + \eta _{N - 2}^n}}{{{{(\Delta x)}^2}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_N^n ={-} \frac{{ - 2.5\eta _N^n + 9\eta _{N - 1}^n - 12\eta _{N - 2}^n + 7\eta _{N - 3}^n - 1.5\eta _{N - 4}^n}}{{{{(\Delta x)}^3}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_{N - 1}^n ={-} \frac{{ - 1.5\eta _N^n + 5\eta _{N - 1}^n - 6\eta _{N - 2}^n + 3\eta _{N - 3}^n - 0.5\eta _{N - 4}^n}}{{{{(\Delta x)}^3}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_{N - 2}^n ={-} \frac{{ - 0.5\eta _N^n + \eta _{N - 1}^n - \eta _{N - 3}^n + 0.5\eta _{N - 4}^n}}{{{{(\Delta x)}^3}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_N^n = \frac{{3\eta _N^n - 14\eta _{N - 1}^n + 26\eta _{N - 2}^n - 24\eta _{N - 3}^n + 11\eta _{N - 4}^n - 2\eta _{N - 5}^n}}{{{{(\Delta x)}^4}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_{N - 1}^n = \frac{{2\eta _N^n - 9\eta _{N - 1}^n + 16\eta _{N - 2}^n - 14\eta _{N - 3}^n + 6\eta _{N - 4}^n - \eta _{N - 5}^n}}{{{{(\Delta x)}^4}}} + O(\Delta {x^2}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_{N - 2}^n = \frac{{\eta _N^n - 4\eta _{N - 1}^n + 6\eta _{N - 2}^n - 4\eta _{N - 3}^n + \eta _{N - 4}^n}}{{{{(\Delta x)}^4}}} + O(\Delta {x^2}).\end{gather}Appendix B. Finite difference schemes for periodic solutions
The schemes of boundary condition for periodic solutions can be derived using a fourth-order accurate central difference scheme consistent with Appendix A. This approach takes advantage of periodicity 
${\eta _i} = {\eta _{N + i}}$. The results are presented as follows:
\begin{gather}\left( {\frac{{\partial \eta }}{{\partial x}}} \right)_1^n = \frac{{ - \eta _3^n + 8\eta _2^n - 8\eta _N^n + \eta _{N - 1}^n}}{{12\Delta x}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{\partial \eta }}{{\partial x}}} \right)_2^n = \frac{{ - \eta _4^n + 8\eta _3^n - 8\eta _1^n + \eta _N^n}}{{12\Delta x}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^2}\eta }}{{\partial {x^2}}}} \right)_1^n = \frac{{ - \eta _3^n + 16\eta _2^n - 30\eta _1^n + 16\eta _N^n - \eta _{N - 1}^n}}{{12{{(\Delta x)}^2}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^2}\eta }}{{\partial {x^2}}}} \right)_2^n = \frac{{ - \eta _4^n + 16\eta _3^n - 30\eta _2^n + 16\eta _1^n - \eta _N^n}}{{12{{(\Delta x)}^2}}} + O(\Delta {x^4}), \end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_1^n = \frac{{ - \eta _4^n + 8\eta _3^n - 13\eta _2^n + 13\eta _N^n - 8\eta _{N - 1}^n + \eta _{N - 2}^n}}{{8{{(\Delta x)}^3}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_2^n = \frac{{ - \eta _5^n + 8\eta _4^n - 13\eta _3^n + 13\eta _1^n - 8\eta _N^n + \eta _{N - 1}^n}}{{8{{(\Delta x)}^3}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_3^n = \frac{{ - \eta _6^n + 8\eta _5^n - 13\eta _4^n + 13\eta _2^n - 8\eta _1^n + \eta _N^n}}{{8{{(\Delta x)}^3}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_1^n = \frac{{ - \eta _4^n + 12\eta _3^n - 39\eta _2^n + 56\eta _1^n - 39\eta _N^n + 12\eta _{N - 1}^n - \eta _{N - 2}^n}}{{6{{(\Delta x)}^4}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_2^n = \frac{{ - \eta _5^n + 12\eta _4^n - 39\eta _3^n + 56\eta _2^n - 39\eta _1^n + 12\eta _N^n - \eta _{N - 1}^n}}{{6{{(\Delta x)}^4}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_3^n = \frac{{ - \eta _6^n + 12\eta _5^n - 39\eta _4^n + 56\eta _3^n - 39\eta _2^n + 12\eta _1^n - \eta _N^n}}{{6{{(\Delta x)}^4}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{\partial \eta }}{{\partial x}}} \right)_N^n = \frac{{ - \eta _2^n + 8\eta _1^n - 8\eta _{N - 1}^n + \eta _{N - 2}^n}}{{12\Delta x}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{\partial \eta }}{{\partial x}}} \right)_{N - 1}^n = \frac{{ - \eta _1^n + 8\eta _N^n - 8\eta _{N - 2}^n + \eta _{N - 3}^n}}{{12\Delta x}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^2}\eta }}{{\partial {x^2}}}} \right)_N^n = \frac{{ - \eta _2^n + 16\eta _1^n - 30\eta _N^n + 16\eta _{N - 1}^n - \eta _{N - 2}^n}}{{12{{(\Delta x)}^2}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^2}\eta }}{{\partial {x^2}}}} \right)_{N - 1}^n = \frac{{ - \eta _1^n + 16\eta _N^n - 30\eta _{N - 1}^n + 16\eta _{N - 2}^n - \eta _{N - 3}^n}}{{12{{(\Delta x)}^2}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_N^n = \frac{{ - \eta _3^n + 8\eta _2^n - 13\eta _1^n + 13\eta _{N - 1}^n - 8\eta _{N - 2}^n + \eta _{N - 3}^n}}{{8{{(\Delta x)}^3}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_{N - 1}^n = \frac{{ - \eta _2^n + 8\eta _1^n - 13\eta _N^n + 13\eta _{N - 2}^n - 8\eta _{N - 3}^n + \eta _{N - 4}^n}}{{8{{(\Delta x)}^3}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^3}\eta }}{{\partial {x^3}}}} \right)_{N - 2}^n = \frac{{ - \eta _1^n + 8\eta _N^n - 13\eta _{N - 1}^n + 13\eta _{N - 3}^n - 8\eta _{N - 4}^n + \eta _{N - 5}^n}}{{8{{(\Delta x)}^3}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_N^n = \frac{{ - \eta _3^n + 12\eta _2^n - 39\eta _1^n + 56\eta _N^n - 39\eta _{N - 1}^n + 12\eta _{N - 2}^n - \eta _{N - 3}^n}}{{6{{(\Delta x)}^4}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_{N - 1}^n = \frac{{ - \eta _2^n + 12\eta _1^n - 39\eta _N^n + 56\eta _{N - 1}^n - 39\eta _{N - 2}^n + 12\eta _{N - 3}^n - \eta _{N - 4}^n}}{{6{{(\Delta x)}^4}}} + O(\Delta {x^4}),\end{gather}
\begin{gather}\left( {\frac{{{\partial^4}\eta }}{{\partial {x^4}}}} \right)_{N - 2}^n = \frac{{ - \eta _1^n + 12\eta _N^n - 39\eta _{N - 1}^n + 56\eta _{N - 2}^n - 39\eta _{N - 3}^n + 12\eta _{N - 4}^n - \eta _{N - 5}^n}}{{6{{(\Delta x)}^4}}} + O(\Delta {x^4}).\end{gather}Appendix C. The validation of grid convergence
We examine the evolution of 
$\eta $ at Q = 1, C = 1 and 
$\tilde{G} ={-} 0.5$ for a solitary wave. The initial conditions are the same as those in figure 3(a). The distribution of 
$\eta $ at t = 0.03 is presented in figure 16. The computation results are the same for different grids.
Distribution of 
$\eta $ at Q = 1, C = 1, 
$\tilde{G} ={-} 0.5$ and t = 0.003 for a solitary wave: (I) the space step 
$\Delta x = 0.06$ and the time step 
$\Delta t = 3 \times {10^{ - 7}}$; (II) 
$\Delta x = 0.04$ and 
$\Delta t = {10^{ - 7}}$; (III) 
$\Delta x = 0.03$ and 
$\Delta t = {10^{ - 8}}$.
Appendix D. Numerical integration method
The first-order accurate difference scheme is adopted for the initial position,
\begin{gather}{\eta ^{\prime\prime}_2} = {\eta ^{\prime\prime}_1} + {\eta ^{\prime\prime\prime}_1} \cdot \Delta x.\end{gather}For other positions, we use the following second-order accurate difference scheme:
\begin{gather}{\eta ^{\prime\prime}_i} = {\eta ^{\prime\prime}_{i - 1}} + \frac{{(3{\eta^{\prime\prime\prime}_{i - 1}} - {\eta ^{\prime\prime\prime}_{i - 2}})}}{2}\Delta x + O(\Delta {x^2}),\end{gather}
\begin{gather}{\eta ^{\prime}_i} = {\eta ^{\prime}_{i - 1}} + \frac{{({\eta ^{\prime\prime}_{i - 1}} + {\eta ^{\prime\prime}_i})}}{2}\Delta x + O(\Delta {x^2}),\end{gather}
\begin{gather}{\eta _i} = {\eta _{i - 1}} + \frac{{({\eta ^{\prime}_{i - 1}} + {\eta ^{\prime}_i})}}{2}\Delta x + O(\Delta {x^2}),\end{gather}
\begin{gather}{\eta ^{\prime\prime\prime}_i} = \frac{{{{({\eta _\infty } - C)}^2} - {{({\eta _i} - C)}^2}}}{{2\eta _i^3}} + \tilde{G}{\eta ^{\prime}_i} + O(\Delta {x^2}).\end{gather}
