Roll waves on laminar sheet flow of Newtonian fluid with negligible surface tension

Abstract

We conduct direct numerical simulations (DNS) to study the temporal and spatial developments of the roll waves on a laminar sheet flow of Newtonian fluid. The DNS unveil the physics of the wavefront and show the limitation of the widely used shallow-layer approximations. The most prominent wave, the front runner, is determined by the DNS for the first time in studying the spatial development of the laminar sheet flow with negligible surface tension. Depending on the Froude and Reynolds numbers, the front runner can be a multi-peaked undular bore or a single-peaked non-breaking or breaking wave. The simulation has uncovered an extended region behind the wavefront, where the bed-friction stress is much higher than the corresponding friction in the undisturbed uniform flow. It also produces an uplift velocity needed in the description of wave breaking. For comparison, we also examine the nonlinear development of the instability using two-equation and four-equation shallow-layer models. The two-equation shallow-layer model has produced the bulk of the wave profile but is deficient because it fails to predict the uplift velocity and the substantial increase in bed friction in the frontal region. The four-equation shallow-layer model correctly predicts the bed friction but cannot produce the breaking wave. The simulations also determine the celerity and amplitude of the front runner to follow a linear relationship, qualitatively similar to the roll waves in a turbulent flow.

1. Introduction

Roll waves exist in thin films of millimetre thickness and at depths of many metres in landslide mudflows and volcanic lava flows. This paper focuses on the large-scale roll wave instability affected mainly by inertia and gravity. We conduct direct numerical simulations (DNS) in this paper, solving the full Navier–Stokes equations for the nonlinear instability of the sheet flow of Newtonian fluid where the surface tension is negligible.

The aim of the present investigation differs from previous studies conducted for roll waves on film with a small thickness in the millimetre range affected mainly by surface tension. A teardrop smooth hump with short-wavelength capillary bow waves in front of the hump are characteristics of roll waves dominated by surface tension and viscosity (Kapitza 1948; Julien & Hartley 1986; Liu, Paul & Gollub 1993; Alekseenko, Nakoryakov & Pokusaev 1994; Liu & Gollub 1994; Nosoko et al. 1996; Ramaswamy, Chippada & Joo 1996; Malamataris, Vlachogiannis & Bontozoglou 2002; Gao, Morley & Dhir 2003; Park & Nosoko 2003).

On the other hand, gravity and inertia are the dominant influences in mudflow and debris-flow avalanches. Hungr, Morgan & Kellerhals (1984) and Hungr (1988) found debris flow at 1–5 metre depths in western Canada and Japan to be laminar, with apparent viscosity approximately $3 \times 10^6$ times that of water. Hunt (1994) analysed field and laboratory data, and concluded that the debris flows are Newtonian fluid flow in the laminar regime.

Giordano, Russell & Dingwell (2008), Hui & Zhang (2007), Giordano & Dingwell (2003) and Giordano et al. (2006) modelled volcanic magma and lava also as viscous Newtonian fluids. Giordano et al. (2008) stated that the viscosity of natural silicate magmas generally spans from $10^{-1}$ to $10^{14}$ Pa s, depending on the melt temperature and composition. The wave phenomena, such as regular hydraulic jump, undulating hydraulic jump, oblique hydraulic jump and standing waves, were observed in channelized lava flows during natural volcanic events (Finch & Macdonald 1955; Takahashi & Griggs 1987; Wolfe et al. 1988; Heslop et al. 1989). Appendix A shows how and why debris torrents and volcanic magma can be modelled as laminar sheet flows with negligible surface tension.

Shkadov (1967) examined the laminar flow instability using a depth-averaged model including the surface tension. Ng & Mei (1994) formulated the theory of mudflow instability by extending Shkadov's model to power-law fluids, and found a ‘permanent roll waves’ solution between discontinuities using their shallow-layer formulation.

Other studies of roll waves were devoted to turbulent flow with quadratic friction (Jeffreys 1925; Dressler 1949; Brock 1967; Zanuttigh & Lamberti 2002; Balmforth & Mandre 2004; Que & Xu 2006; Richard & Gavrilyuk 2012; Ivanova et al. 2017; Yu & Chu 2022, 2023). Jeffreys (1925) analysed the linear stability of turbulent sheet flow using a shallow-layer model, and found the critical Froude number ${\textit {Fr}}=2$ – a critical Froude number very different from the value ${\textit {Fr}} = 0.527$ obtained by Yih (1955, 1963) and Benjamin (1957) for laminar sheet flow. The hyperbolic equations of the shallow-layer model for the turbulent flow also produce solutions with mathematical discontinuities. Dressler (1949) found the ‘permanent roll waves’ by matching the smooth profiles between discontinuities. Richard & Gavrilyuk (2012) and Richard (2024) proposed supplemental energy and enstrophy equations in modelling the front of roll waves in turbulent flow, and found improved agreement with laboratory observation of the roll waves by Brock (1967).

One focus of the present study of the roll waves in the laminar sheet flow is the physics in a region immediately behind the wavefront. The other focus is examining the relationship between temporal and spatial instability. The temporal instability depends on a fixed wavelength specification. In spatial instability development, the prominent wave emerges as a front runner, which advances with increasing wavelength, wave amplitude and celerity. Yu & Chu (2022, 2023) studied these problems in turbulent flow. This paper is about laminar flow of Newtonian fluid with negligible surface tension. In the DNS, we solve the Navier–Stokes equations using the volume of fluid (VOF) method. For comparison, we also conduct the simulation using the two-equation shallow-layer simulations (SLS) model of Ng & Mei (1994) and the four-equation shallow-layer weighted residual integral boundary layer (WRIBL) model of Ruyer-Quil & Manneville (2000).

This paper has seven sections and four appendices, including this Introduction. Section 2 describes the formulation for DNS, the two-equation SLS model and the four-equation shallow-layer WRIBL model. Section 3 reviews the existing linear instability theories. Section 4 presents the simulations of the roll waves developed from temporal instability of a specified wavelength using a periodic boundary condition. The difference between the DNS and SLS is presented in this section using the analytical solution of Ng & Mei (1994) as a reference. Section 5 presents the spatial development of the roll waves. In this section, we determine the dependence of the celerity on the roll wave's amplitude. Section 6 examines a more accurate but still inadequate shallow-layer WRIBL model. Section 7 further explains the roll wave's spatial development versus its temporal development, and discusses the limitations of the two-equation and four-equation shallow-layer models. Appendix A explains why surface tension is negligible in geophysical debris and lava flows. Appendix B describes the adaptive quadtree mesh used in the DNS, and its refinement for numerical accuracy. Appendix C provides the computation details of the four-equation shallow-layer model. Appendix D evaluates the computation resources requirement for DNS compared with two-equation and four-equation shallow-layer approximations.

2. Formulation and numerical methods

We conduct DNS of the roll wave instability by solving the full Navier–Stokes equations using a two-phase VOF solver. The simulation starts from the initial instability of the sheet flow to the full nonlinear development of the roll waves. For comparison, we also conduct simulation using a two-equation SLS model of Ng & Mei (1994) and a four-equation WRIBL model by Ruyer-Quil & Manneville (2000). Figure 1 shows two roll wave instability problems examined in this paper. Figure 1(a) shows the temporal development of the roll waves produced by a periodic disturbance. Figure 1(b) shows the spatial development of the roll waves produced by a localized disturbance.

Figure 1. (a) Periodic wave trains fully developed nonlinearly from temporal instability on an incline of angle $\theta$ relative to the horizontal. (b) Wave packet characterized by the front runner (FR) developed from spatial instability. The dashed lines show the perturbation to the initial steady and uniform flow, which has depth $H$ and average depth velocity $U$. The four identical waves of wavelength $\lambda$ in (a) are initiated by periodic disturbance simulated using the periodic boundary conditions. On the other hand, the ‘wave packet’ and the front runner in (b) are developed from a localized disturbance near the inlet. The colour map shows the magnitude of the velocity $\sqrt {u^2+w^2}/U$, where $(u, w)$ are velocity components in the $x$- and $z$-directions, respectively.

2.1. Two-dimensional DNS

The governing equations for the DNS are

(2.1)$$\begin{gather} \frac{\partial u_i }{\partial x_i} =0, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \rho u_i}{\partial t}+ \frac{ \partial \rho u_i u_j}{\partial x_j} =- \frac{\partial p}{\partial x_i} + \frac{\partial \tau_{ij}}{\partial x_j} + \rho g_i, \end{gather}$$

(2.3)$$\begin{gather}\frac{\partial \alpha}{\partial t}+ u_i\,\frac{\partial \alpha}{\partial x_i}=0, \end{gather}$$

where the $i$ and $j$ subscripts follow Einstein's summation convention, $u_i$ is the velocity vector in the $x_i$-direction, and $g_i$ is the gravity component in the $x_i$-direction. The deviatoric stress tensor is related to the strain-rate tensor as follows:

(2.4)\begin{equation} \tau_{ij} = \mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right), \end{equation}

where $\mu$ is the dynamic viscosity. For a two-dimensional simulation of the flow down an incline shown in figure 1, the components of the velocity vector are $u_i = (u,w)$, and the components of gravity are $g_i= g (\sin \theta, \cos \theta )$ in the $x$- and $z$-directions, respectively. The volume fraction of the liquid $\alpha$, varying from 0 to 1, determines the volumetric average of density and viscosity as follows:

(2.5)$$\begin{gather} \rho = \alpha \rho_f +(1-\alpha){\rho_a}, \end{gather}$$

(2.6)$$\begin{gather}\mu =\left[ \frac{\alpha}{\mu_f} + \frac{1-\alpha}{\mu_a} \right]^{-1}, \end{gather}$$

where the subscripts $f$ and $a$ refer to the liquid and air phases, respectively. Tryggvason, Scardovelli & Zaleski (2011) recommended the harmonic average as specified by (2.6) for large viscosity ratios. The piecewise linear interface calculation geometric VOF solver in the Basilisk platform (Popinet 2009; Popinet & Collaborators 2013–2024) determines the location of the air–liquid interface. The VOF scheme is momentum-conserving, enabling the momentum flux to be advected consistently with the mass flux (Fuster & Popinet 2018; Zhang, Popinet & Ling 2020), thus improving the computational robustness for cases with large density and viscosity ratios. The DNS are conducted using an adaptive quadtree mesh. The refined mesh has maintained within a less than 1 % computational error before wave breaking, according to the mesh refinement study in Appendix B. The accuracy of the simulation is not determined when the wave breaks randomly, entrapping air below the waves.

2.2. Two-equation SLS model

The two equations for the SLS are

(2.7)$$\begin{gather} \frac{\partial h}{\partial t} + \frac{\partial \bar{u} h}{\partial x} = 0, \end{gather}$$

(2.8)$$\begin{gather}\frac{\partial \bar{u}h}{\partial t} + \frac{\partial \left(\tfrac{6}{5} {\bar{u}^2}{h}+\tfrac{1}{2}g'h^2\right)}{\partial x} = g'hS_o - 3\,\frac{\mu}{\rho_f}\, \frac{\bar{u}}{h}, \end{gather}$$

where $g^\prime =g\cos \theta$ is the reduced gravity on the slope, and $S_o= \tan \theta$ is the channel slope. The formulation for the two-equation SLS model assumes hydrostatic pressure distribution and self-similar parabolic velocity profiles:

(2.9)\begin{equation} {u} (z)= \frac{3}{2} \left[1- \left(1-\frac{z}{h}\right)^{2}\right] \bar{u} , \end{equation}

where $\bar {u}=({1}/{h})\int _0^h {u}(z) \,{\rm d}z$ is the depth-averaged velocity. Ng & Mei (1994) developed the two-equation SLS model for power-law fluid to study landslide mudflow. The power-law index $n=1$ reduces the model to Newtonian laminar flow. Equations (2.7) and (2.8) form a hyperbolic equation set. They lead to mathematical discontinuities. The ‘permanent roll wave’ solution of Ng & Mei (1994) comprises the smooth-wave profiles connecting the discontinuities.

We conduct the SLS solving (2.7) and (2.8) on the Basilisk platform. The central-upwind scheme by Kurganov & Levy (2002) is used for spatial discretization. The method by Kurganov & Lin (2007) is adopted to reduce the numerical dissipation. The generalized MINMOD TVD limiter is employed for proper shock capturing, avoiding spurious numerical oscillation. The time integration scheme is a second-order predictor–corrector scheme. The details of the solver are elaborated in Popinet (2011).

The two-equation SLS model is one-dimensional, involving only the $\bar {u}$ velocity component. But the sheet flow can have $u$ and $w$ velocity components. Once disturbed, the pressure variation over the depth of the flow will no longer be hydrostatic, and velocity variation will not be parabolic in the region behind the wavefront.

2.3. Four-equation WRIBL model

More accurate one-dimensional shallow-layer models exist to account for the deviation of the velocity profile from the parabolic profile. One such model is the WRIBL model by Ruyer-Quil & Manneville (2000). The model involves four equations, solving for depth $h$, flow rate $q=\bar {u} h$, and two variables $s_1$ and $s_2$ introduced to account for deviation from the parabolic velocity profile. We conduct the WRIBL model simulation to compare with our DNS of roll waves using the open-source code WaveMaker by Rohlfs, Rietz & Scheid (2018). The governing equations for the full second-order WRIBL model are complicated, involving lengthy equations with terms of up to third-order derivatives, and are not reproduced here. Interested readers could refer to (11), (38), (39) and (40) in Ruyer-Quil & Manneville (2000).

Ruyer-Quil & Manneville (2002) performed linear stability analysis and numerical simulation of thin film dominated by surface tension, and found good agreement of the WRIBL model with experiment, analytical solution and DNS up to $Re=200$. However, to our knowledge, the WRIBL model has not been employed to study the sheet flow without surface tension.

2.4. The compatibility condition for steady and uniform undisturbed flow

The numerical simulation for the instability starts with an undisturbed sheet flow, which must be steady and uniform with a parabolic velocity profile given by (2.9). The depth and depth-averaged velocity of the undisturbed flow are $H$ and $U$, respectively. The balance of the gravity driving force and the bed friction in the steady and uniform flow gives the relation

(2.10)\begin{equation} U= \frac{1}{3}\left[\frac{\rho g \sin{\theta}}{\mu}\,H^{2}\right], \end{equation}

which in dimensionless form is

(2.11)\begin{equation} {{\textit{Fr}}^2}=\tfrac{1}{3}\,{Re}\,S_o, \end{equation}

where ${\textit {Fr}} = U/\sqrt {g'H}$ is the Froude number, and ${Re}=\rho UH/\mu$ is the Reynolds number. Equation (2.11) is a condition to ensure a steady and uniform undisturbed flow. Meaningful stability analysis must satisfy this compatibility condition. The result of DNS depends on any two of the three dimensionless parameters ${\textit {Fr}}$, $Re$ and $S_o$. However, the result of SLS depends only on ${\textit {Fr}}$ (Ng & Mei 1994).

3. Linear instability theories

The instability of laminar sheet flow has been the subject of many studies. We summarize the neutral stability boundaries obtained from various approximations for their dependence on the Froude and Reynolds numbers, ignoring the surface tension effect.

3.1. The first-order long-wavelength approximation

Yih (1955, 1963) and Benjamin (1957) analysed the linear instability of the sheet flow using a long-wavelength approximation of the Orr–Sommerfeld equations, and found that the Froude number needs to be smaller than the critical value ${\textit {Fr}} = \sqrt {5/18} \approx 0.527$ as a requirement for stability of the laminar flow. Ng & Mei (1994) found the critical Froude number to be ${\textit {Fr}}= \sqrt {1/3} \approx 0.577$ using the two-equation shallow-layer model. The thick blue dashed line and blue dotted lines in figure 2 delineate the neutral stability boundaries of ${\textit {Fr}} = 0.527$ and 0.577, respectively.

Figure 2. The dependence of the neutral stability boundaries on Froude number and Reynolds number expressed in brown, green and purple colour according to (a) Chang & Demekhin (2002), (b) Lee & Mei (1996) and (c) Ruyer-Quil & Manneville (2002), respectively. The thick curve is for the perturbation wavelength $S_o\lambda /H=4.72$, while the thin curve is for a much shorter wavelength, $S_o\lambda /H=$ 0.472. The thick blue dashed line is the neutral stability boundary obtained by Yih (1955, 1963), Benjamin (1957) and Ruyer-Quil & Manneville (2002). The blue dotted line is the neutral stability boundary according to the two-equation shallow-layer model of Ng & Mei (1994). The thin grey dotted lines are the relation between ${\textit {Fr}}$ and $Re$ for a given $S_o$, according to the compatibility condition (2.11). The black open squares, red open circles and red solid circles denote the stable flows, non-breaking roll waves and breaking roll waves obtained from DNS for wavelength $S_o\lambda /H=4.72$.

The stability analysis of Yih (1955, 1963) and Benjamin (1957), and the two-equation model of Ng & Mei (1994), were both first-order approximations in long-wavelength perturbation expansions. The neutral stability boundaries obtained from these first-order approximations depend only on the Froude number. However, the critical Froude number 0.577 obtained by the two-equation model is slightly greater than the correct value 0.527 because of the two-equation model assumption of hydrostatic pressure and unchanged parabolic velocity distribution over the depth of the laminar flow.

3.2. Higher-order long-wavelength approximation

Subsequent stability analyses by Chang & Demekhin (2002), Lee & Mei (1996) and Ruyer-Quil & Manneville (2002) to higher orders of approximation have shown that stability also depends on the Reynolds number and the wavelength of the perturbation.

The neutral stability condition derived from the long-wavelength expansion of the Orr–Sommerfeld equation by Chang & Demekhin (2002) is

(3.1)\begin{equation} 4\,{Re}\,{ \sqrt{9009(-5+18\,{{\textit{Fr}}}^2) \left(\frac{S_o \lambda}{H}\right)^2+5{\rm \pi}^2} } = { 9 {\rm \pi}\sqrt{1\,107\,249} \, {{\textit{Fr}}}^2 } . \end{equation}

Lee & Mei (1996) captured the non-hydrostatic effect in their shallow-layer formulation, and obtained the following neutral stability condition:

(3.2)\begin{equation} {Re}\,{\sqrt{(-1+3\,{{\textit{Fr}}}^2) \left(\frac{S_o \lambda}{H}\right)^2}} = {6 {\rm \pi}\, \sqrt{\frac{14}{5}} \, {{\textit{Fr}}}^2 }. \end{equation}

The neutral stability condition by the simplified second-order WRIBL model of Ruyer-Quil & Manneville (2002) is

(3.3)\begin{equation} {Re}\left(\frac{S_o\lambda}{H}\right) = 6{\rm \pi}\,\sqrt{\frac{3}{35}}\,{\textit{Fr}}^2\, \sqrt{\frac{-3\,{\textit{Fr}}^2+5 \sqrt{333\,{\textit{Fr}}^2+1470}\,{\textit{Fr}}+105}{18\,{\textit{Fr}}^2-5}}. \end{equation}

Equations (3.1), (3.2) and (3.3) are different in form from the original neutral stability conditions given in Chang & Demekhin (2002), Lee & Mei (1996) and Ruyer-Quil & Manneville (2002). We eliminated the dependence of the original stability conditions on $S_o$ using compatibility condition (2.11).

Figure 2 shows the neutral stability boundaries on the ${\textit {Fr}}$–${Re}$ plane. The grey dotted lines show the channel slope $S_o$ that satisfies the compatibility condition. The brown, green and purple colour lines in figures 2(a–c) represent (3.1), (3.2) and (3.3), respectively. The neutral stability condition depends on the wavelength of the perturbation. The thick and thin lines in the figure are neutral stability boundaries for the wavelength parameter $S_o \lambda /H= 4.72$ and 0.472, respectively. The selected wavelength $S_o \lambda /H= 4.72$ represents those roll waves observed in the laboratory by Julien & Hartley (1986). On the other hand, the wavelength $S_o \lambda /H= 0.472$ represents the much shorter wavelength cases, with greater dependency on the wavelength and Reynolds number effects.

3.3. Stability dependence on wavelength and Reynolds number

The higher-order approximations have led to the roll wave instability dependence on the wavelength parameter $S_o \lambda /H$ and the Reynolds number $Re$. As shown by the thick and thin curves in figure 2 of the neutral stability boundaries for the perturbation wavelengths $S_o\lambda /H=4.72$ and 0.472, the critical Reynolds numbers are approximately $Re_c \approx 4.50$ and 45.0, respectively. Therefore, the laminar sheet flow with a small Reynolds number $Re < 4.50$ would be stable to perturbation wavelengths shorter than $S_o\lambda /H=4.72$; even the Froude number is greater than the critical value 0.527 determined by the first-order approximation.

The laminar sheet flow is more unstable when perturbed by the longer wavelength of $S_o\lambda /H=4.72$. The stability of the sheet flow depends only on the Froude number in the limiting case when the perturbation wavelength is infinitely long.

4. Temporal development of periodic roll waves

The hydrodynamic stability of sheet flow has been examined using linear and nonlinear theories. Assuming periodic disturbance, the eigenvalues of the linearized system of equations determine the initial exponential growth rate. The growth rate diminishes with the increase in the waves’ amplitude due to the nonlinear effect. The waves become fully developed, in equilibrium with the nonlinear effect, as the growth diminishes to zero.

4.1. The DNS for roll waves with $S_o \lambda /H=4.72$

We conduct the DNS from the initial instability to the fully developed nonlinear roll waves, then compare the results with the permanent roll waves of Ng & Mei (1994). The simulation starts with a wavelength $S_o \lambda /H=4.72$ perturbation to find the fully developed nonlinear waves when the growth rate diminishes to zero. Figure 3 shows the DNS compared with the SLS. The blue dashed lines in the figure outline the wave profiles obtained by SLS. The horizontal thin lines are the crest and trough of the permanent roll wave solution by Ng & Mei (1994).

Figure 3. The periodic roll waves in a quasi-steady state developed from temporal instability. The colour maps show the velocity contours of the roll waves obtained from DNS. The blue dashed lines outline the wave profiles obtained by SLS. The horizontal thin lines are the crest and trough of the permanent roll wave solution by Ng & Mei (1994). The iso-lines for the $w$ velocity component, denoted by the green curves on the right-hand side of the figure, show the significant uplifting velocity at the wavefront. Settings are (a i,ii) ${\textit {Fr}}=0.840$, $Re=5.88$, $S_o=0.360$; (b i,ii) ${\textit {Fr}}=1.10$, $Re=7.26$, $S_o=0.500$; (c i,ii) ${\textit {Fr}}=0.840$, $Re=70.6$, $S_o=0.0300$; (d i,ii) ${\textit {Fr}}=0.600$, $Re=72.0$, $S_o=0.0150$; (e i,ii) ${\textit {Fr}}=0.840$, $Re=141$, $S_o=0.0150$; (f i,ii) ${\textit {Fr}}=1.14$, $Re=260$, $S_o=0.0150$; (g i,ii) ${\textit {Fr}}=1.40$, $Re=392$, $S_o=0.0150$.

4.1.1. Stable cases

The temporal development of the waves depends on the wavelength parameter and the Reynolds number $Re = \rho UH/\mu$ and the Froude number ${\textit {Fr}}=U/\sqrt {g^\prime H}$ of the base flow. The DNS amplitude diminishes to zero for the two cases with $({\textit {Fr}}, Re) = (0.840, 4.23)$ and $(1.10, 4.85)$. These two cases, denoted by the black open square symbols in the stability diagram figure 2, are on the stable side of the neutral stability boundary with $S_o \lambda /H=4.72$. Therefore, the diminishing DNS amplitude for these two cases agrees with the stability theories by Chang & Demekhin (2002), Lee & Mei (1996) and Ruyer-Quil & Manneville (2002).

4.1.2. Unstable cases a, b, c, d, e, f and g

We also conduct DNS for cases on the unstable side of the neutral stability boundaries with $S_o \lambda /H=4.72$. The Reynolds number and Froude number of these unstable cases labelled as a, b, c, d, e, f and g in figure 2 are $({\textit {Fr}}, Re) = (0.840, 5.88)$, $(1.10, 7.26)$, $(0.840, 70.6)$, $(0.600, 72.0)$, $(0.840, 141)$, $(1.14, 260)$ and $(1.40, 392)$, respectively. The perturbation amplitude of these unstable periodic disturbances increases with time and eventually approaches a fully nonlinear quasi-steady state when growth diminishes to zero. Figure 3 shows these fully nonlinear roll wave profiles in their quasi-steady state. The plots on the left-hand side of the figure show the roll waves over the entire wavelength from $S_ox/H=0$ to 4.72. The right-hand plots show the close-up views of the wavefronts in undistorted scales.

The wavefront in the left-hand plots of figure 3, expressed in terms of the compressed coordinate $S_o x/H$, appears like a discontinuity, the amplitude of which closely agrees with the permanent roll wave solution of Ng & Mei (1994). However, in the close-up undistorted scales in the right-hand plots, the wavefront is either a multi-peaked undular bore or a single-peaked non-breaking or breaking wave.

4.2. The SLS and permanent roll waves

We also conducted the SLS to compare them with the DNS. The SLS profiles obtained from the hyperbolic model equations (2.7) and (2.8) shown as dashed lines in figure 3 are infinitely steep at the wavefront. The amplitudes of the wavefronts obtained by the SLS are in perfect agreement with the ‘permanent roll wave’ solution of Ng & Mei (1994) delimited between the blue horizontal thin lines in the figure.

On the other hand, the SLS for this case a, however, has produced a wave profile delineated by the blue dashed line, with depth $\hat {h}_{SLS}=1.29$ at the wave crest, and depth $\check {h}_{SLS}=0.824$ at the wave trough, which are in perfect agreement with the permanent roll wave solution of Ng & Mei (1994) delimited between the blue horizontal thin lines, because the Froude number ${\textit {Fr}} = 0.840$ in this case is greater than the critical value 0.577 according to the one-dimensional two-equation model of Ng & Mei (1994).

4.3. Smooth wavefronts and undular bore for marginally unstable cases a, b and d

The difference between the DNS and SLS profiles occurs mainly at the wavefront, where the long wavelength approximation breaks down. The wavefront is smooth and gradually varying for the marginally unstable cases a and b, where the Reynolds number is slightly greater than the critical value, which is slightly above the neutral stability boundaries as shown in figure 2. Figures 3(a) and 3(b) show the smooth wavefronts with the roll wavefront amplitude slightly smaller than in the SLS model.

The other significant difference between the DNS and SLS is the profile for the marginally unstable case d. The Froude number of this case is ${\textit {Fr}} = 0.600$, which is higher but very close to the critical value ${\textit {Fr}}_c=$ 0.527. However, the Reynolds number for this case, $Re = 72$, is quite large. The profile behind the wavefront in this marginally unstable case, with a large Reynolds number, is the in the multi-peaked undular bore profile as shown in figure 2(d).

4.4. Breaking wavefront for unstable cases c, e, f, and g

The sheet flows in cases d, e, f and g are on the same slope $S_o=0.0150$. The frontal region of case d is a multi-peaked undular bore. As the Froude number increases from ${\textit {Fr}} = 0.600$ to 0.840, 1.14 and 1.40, away from the critical value 0.527, the wavefronts in cases c, e, f and g become progressively more chaotic, as shown in figures 2(c) and 2(e–g). We identify the smooth waves using open circle symbols, and breaking waves using solid circle symbols, in figure 2.

The breaking wavefront shows the entrainment of air bubbles. The process of air entrapment is chaotic. The present level of mesh refinement predicts the initiation of breaking and the bulk of the roll wave but may not correctly simulate the detailed chaotic process of air entrapment.

4.5. The perturbation wavelength associated with the temporal analysis

The simulation cases from a to g presented in figure 3 are for a specified perturbation wavelength parameter $S_o \lambda /H = 4.72$. The results would be different, and the neutral stability boundary would be the thin curves in figure 2, if we had selected a much shorter wavelength $S_o \lambda /H =0.472$. Cases a and b would become stable, and cases c and d become marginally unstable because cases a and b would be on the stable side of the neutrally stable boundaries for $S_o \lambda /H = 0.472$, which are the thin curves in figure 2.

On the other hand, if the perturbation wavelength is infinitely long, then the stability of laminar sheet flow would depend only on the critical Froude number. The flow will become unstable due to the disturbance of infinitely long wavelength as long as the Froude number exceeds the critical value ${\textit {Fr}}_c=0.527$ – unaffected by the Reynolds number.

5. Spatial development of the front runner produced by localized disturbances

The classical hydrodynamic stability theory derived from the normal-mode approach assumes periodic disturbance, which leads to periodic roll waves, whose amplitudes are not determined until the wavelength is specified. The permanent roll waves by Ng & Mei (1994) also depend on the wavelength.

To overcome this difficulty of indeterminacy, Yu & Chu (2022, 2023) studied the spatial development of roll waves produced by a localized disturbance, focusing on the dominant leading waves – the front runner. The works of Yu and Chu were for the turbulent flow of clear waters.

We conduct DNS for the spatial development of the roll waves on the laminar sheet flow, initiated by localized disturbance at the inlet, then follow the growth of the disturbance for a series of unstable sheet flows of undisturbed Froude number ${\textit {Fr}} = 0.70$, 0.84 and 1.00. The corresponding undisturbed flow Reynolds numbers of these unstable sheet flows are $Re = 98.0$, 35.3 and 50.0, and the channel slopes are $S_o=0.015$, 0.06, 0.06, respectively. The Froude number ${\textit {Fr}}$, the Reynolds number $Re$, and the channel slope $S_o$ are selected for the simulation to satisfy the compatibility equation (2.11) so that gravity is in balance with the bed friction stress in the undisturbed base flow, which is steady and uniform flow.

The most prominent roll wave in spatial development is the front runner, which leads the waves in the packet developed from instability. Yu & Chu (2022) studied the front runner in turbulent sheet flow. The front runner in a wave packet was also observed in a mudflow simulation by Liu & Mei (1994) and by Meza & Balakotaiah (2008), who studied the nonlinear waves on vertically falling films dominated by viscosity and surface tension. Meza & Balakotaiah (2008) used the word ‘tsunami’ in exaggeration to describe the front runner in their study of the falling film.

Figures 4, 5, 6 and 7 show the development of the front runner. The depth $h_{FR(DNS)}$ and the depth-averaged velocity $\bar {u}_{FR(DNS)}$ of the front runner obtained from DNS are compared with the corresponding results $h_{FR(SLS)}$ and $\bar {u}_{FR(SLS)}$ obtained from SLS in a table at the bottom of each figure. On the left-hand sides of the figures, the depth and velocity profiles are on the compressed coordinate $S_ox/H$. The right-hand sides show the close-up profiles of the front runner on the undistorted scales.

Figure 4. The development of the front runner produced by a localized disturbance at the inlet on an unstable sheet flow with ${\textit {Fr}}=0.560$, ${Re}=15.7$ and $S_o=0.060$. The colour maps show the DNS results, while the blue dashed lines denote the free-surface profiles obtained by the SLS. The green curves are iso-lines for the $w$ velocity component. Here: (a i,ii) $S_ot_{DNS}U/H=4.25$, $S_ot_{SLS}U/H=3.94$; (b i,ii) $S_ot_{DNS}U/H=8.50$, $S_ot_{SLS}U/H=8.02$; (c i,ii) $S_ot_{DNS}U/H=12.8$, $S_ot_{SLS}U/H=12.2$; (d i,ii) $S_ot_{DNS}U/H=17.0$, $S_ot_{SLS}U/H=16.4$; (e i,ii) $S_ot_{DNS}U/H=21.3$, $S_ot_{SLS}U/H=20.6$.

Figure 5. The development of the front runner produced by a localized disturbance at the inlet on an unstable sheet flow with ${\textit {Fr}}=0.700$, ${Re}=98.0$ and $S_o=0.0150$. The colour maps show the DNS results, while the blue dashed lines denote the free-surface profiles obtained by the SLS. The green curves are iso-lines for the $w$ velocity component. Here: (a i,ii) $S_ot_{DNS}U/H=2.00$, $S_ot_{SLS}U/H=1.87$; (b i,ii) $S_ot_{DNS}U/H=4.00$, $S_ot_{SLS}U/H=3.77$; (c i,ii) $S_ot_{DNS}U/H=6.00$, $S_ot_{SLS}U/H=5.68$; (d i,ii) $S_ot_{DNS}U/H=8.00$, $S_ot_{SLS}U/H=7.61$; (e i,ii) $S_ot_{DNS}U/H=10.0$, $S_ot_{SLS}U/H=9.54$.

Figure 6. The development of the front runner produced by a localized disturbance at the inlet on an unstable sheet flow with ${\textit {Fr}}=0.840$, ${Re}=35.3$ and $S_o=0.0600$. The colour maps show the DNS results, while the blue dashed lines denote the free-surface profiles obtained by the SLS. The green curves are iso-lines for the $w$ velocity component. Here: (a i,ii) $S_ot_{DNS}U/H=4.60$, $S_ot_{SLS}U/H=4.32$; (b i,ii) $S_ot_{DNS}U/H=9.20$, $S_ot_{SLS}U/H=8.74$; (c i,ii) $S_ot_{DNS}U/H=13.8$, $S_ot_{SLS}U/H=13.2$; (d i,ii) $S_ot_{DNS}U/H=18.4$, $S_ot_{SLS}U/H=17.8$; (e i,ii) $S_ot_{DNS}U/H=23.0$, $S_ot_{SLS}U/H=22.3$.

Figure 7. The development of the front runner produced by a localized disturbance at the inlet on an unstable sheet flow with ${\textit {Fr}}=1.00$, ${Re}=50.0$ and $S_o=0.0600$. The colour maps show the DNS results, while the blue dashed lines denote the free-surface profiles obtained by the SLS. The green curves are iso-lines for the $w$ velocity component. Here: (a i,ii) $S_ot_{DNS}U/H=4.80$, $S_ot_{SLS}U/H=4.53$; (b i,ii) $S_ot_{DNS}U/H=9.60$, $S_ot_{SLS}U/H=9.13$; (c i,ii) $S_ot_{DNS}U/H=14.4$, $S_ot_{SLS}U/H=13.8$; (d i,ii) $S_ot_{DNS}U/H=19.2$, $S_ot_{SLS}U/H=18.4$; (e i,ii) $S_ot_{DNS}U/H=24.0$, $S_ot_{SLS}U/H=23.2$.

The front runner development is not dependent on the form of the inlet disturbance. Yu & Chu (2022) examine the instability due to six types of localized disturbances. For the present simulation, the disturbance is ‘type c of constant ${\textit {Fr}}$’, which, according to Yu & Chu (2022), is a small rise in depth in the form of a half-sinusoid, and at the same time a slight increase in velocity to keep the disturbance Froude number ${\textit {Fr}}$ unchanged.

To reproduce the wave packet from the instability, we introduced a localized perturbation to the uniform flow. The depth-averaged initial perturbation at $t=0$ is

(5.1a)$$\begin{gather} h(x, t=0)= \begin{cases} H\left[1+\epsilon\sin\left(\dfrac{2{\rm \pi} x}{\ell}\right)\right] & {\rm for}\ 0< x\leq\dfrac{\ell}{2},\\ H & {\rm for} \ x > \dfrac{\ell}{2}, \end{cases} \end{gather}$$

(5.1b)$$\begin{gather}\bar{u} (x, t=0 ) = {\textit{Fr}}\,\sqrt{g' [h(x,t=0) ]^3}. \end{gather}$$

The depth of the initial disturbance $h(x, t=0)$, is the upper half of a sinusoid of length $\ell$ raised to an amplitude $\epsilon$ over the extent $0< x\leq {\ell }/{2}$, while $\bar {u}$ is selected to keep the Froude number of the disturbance equal to the Froude number of the undisturbed flow. This initial perturbation defined by (5.1) is similar to the ‘type c disturbance of constant ${\textit {Fr}}$’, which is one of the six inlet disturbances examined by Yu & Chu (2022) and Yu & Chu (2023), who found that all six produced the same front runner.

For DNS, we assume parabolic variation of velocity over the depth. The initial profiles for $\alpha$, $u$ and $w$ are

(5.2a)$$\begin{gather} \alpha(x, z, t=0)= \begin{cases} 1 & {\rm for} \ 0< z\leq h(x, t=0),\\ 0 & {\rm for} \ z>h(x, t=0), \end{cases} \end{gather}$$

(5.2b)$$\begin{gather}u (x, z, t=0 ) = \begin{cases} \dfrac{3}{2} \, \bar{u} (x, t=0 ) \left[1-\left(1-\dfrac{z}{h(x, t=0)}\right)^2 \right] & {\rm for} \ 0< z\leq h(x,t=0),\\[10pt] \dfrac{3}{2} \, \bar{u} (x, t=0 ) & {\rm for} \ z>h(x, t=0), \end{cases} \end{gather}$$

(5.2c)$$\begin{gather}w (x, z, t=0 ) = 0. \end{gather}$$

The front runner is not sensitive to the amplitude and length of the initial inlet disturbance, as demonstrated in the previous papers by Yu & Chu (2022, 2023). For the simulations shown in figures 4, 5, 6 and 7, $\epsilon =0.2$ and $S_o\ell /H =2$.

5.1. Profile of the front runner

Depending on the Froude and Reynolds numbers, the front runner can be a multi-peaked undular bore or a single-peaked non-breaking or breaking wave. However, according to the SLS model, the wavefront is nothing but a discontinuity.

The front runner for the case $({\textit {Fr}}, Re) = (0.56, 15.7)$ shown in figure 4 has a smooth wavefront. The case $({\textit {Fr}}, Re) = (0.70, 98.0)$ shown in figure 5 has a smooth frontal region in the form of an undular bore.

For the higher Froude number cases with $({\textit {Fr}}, Re) = (0.84, 35.3)$ and $(1.00, 50.0)$ shown in figures 6 and 7, the steep wavefront eventually breaks, entrapping air in the form of a plunging wave breaker. The process of breaking is somewhat chaotic. The Reynolds number of the front runner is as large as $Re_{FR(DNS)}={\rho h_{FR}\bar {u}_{FR}}/{\mu }=218$ nearly in the transition from laminar to turbulent flow in these high Froude number cases. The laboratory observation of laminar falling films by Adomeit & Renz (2000) has found isolated turbulent spots interacting with surface waves at Reynolds numbers varying from 75 to 200. The wave breaking and entrapment of air in the front runner are chaotic, as shown in figures 6(e) and 7(e). Further mesh refinement would be required to resolve the turbulent spots (if any) in our DNS.

The table at the bottom of each figure provides the critical data obtained from the simulations. The front runner velocity obtained by DNS, ${u}_{FR(DNS)}$, is approximately 2 % larger than the velocity ${u}_{FR(SLS)}$ obtained by SLS.

For the case shown in figure 6, the front runner Froude numbers are ${\textit {Fr}}_{FR(DNS)} =\bar {u}_{FR}/\sqrt {g' h_{FR}}=0.909$, 1.04, 1.15, 1.32 and 1.54, and the front runner Reynolds numbers are $Re_{FR(DNS)} ={\rho h_{FR}\bar {u}_{FR}}/{\mu }=64.0$, 79.9, 95.5, 121 and 164, at front runner positions $S_ox_{FR}/H=13.8$, 26.8, 40.6, 55.3 and 71.7, respectively. A breaking wave is observed to occur as the front runner advances to a position $S_ox_{FR}/H=40$, beyond which the motion progressively becomes chaotic.

For the case shown in figure 7, the front runner Froude numbers are ${\textit {Fr}}_{FR(DNS)} =1.20$, 1.26, 1.32, 1.42 and 1.54, and the front runner Reynolds numbers are $Re_{FR(DNS)} = 107$, 122, 143, 173 and 218, at front runner positions $S_ox_{FR}/H=13.5$, 26.2, 39.8, 54.4 and 70.3, respectively. At position $S_ox_{FR}/H=70.3$, the front runner Froude number in this case has value ${\textit {Fr}}_{FR(DNS)} = 1.54$, which is getting close to the critical Froude number ${\textit {Fr}} = 2.0$ for turbulent roll waves according to the stability analysis of Jeffreys (1925).

5.2. The front runner's depth and velocity

The velocity of the front runner $\bar {u}_{FR}/U$ increases with its depth $h_{FR}/U$ following closely a linear relationship as shown in figure 8. The strong linear dependence of velocity on depth is observed in debris flow, and it has been used by Hungr et al. (1984) and Wu et al. (1990) to show debris flow to be in a laminar regime.

Figure 8. The front runner's depth $h_{FR}/H$ and its relation to the velocity $\bar {u}_{FR}/U$ for (a) ${\textit {Fr}} = 0.56$ ($Re = 15.7$, $S_o=0.060$), (b) ${\textit {Fr}} = 0.70$ ($Re = 98.0$, $S_o=0.015$), (c) ${\textit {Fr}} = 0.84$ ($Re = 35.3$, $S_o=0.060$) and (d) ${\textit {Fr}} = 1.00$ ($Re = 50.0$, $S_o= 0.060$).

5.3. Front runner's celerity–amplitude relationship

The celerity of the front runner depends on wave amplitude. Figure 9 shows the celerity–amplitude relationship. We use the front runner's depth-averaged speed $\bar {u}_{FR}/U$ to represent its amplitude, and the celerity $c_{FR}/U$ is related to the front runner's amplitude $\bar {u}_{FR}/U$. The DNS results all follow approximately a linear relationship denoted by the dashed lines in the figure. The SLS produces slightly different trends. The front runner is the most prominent wave because it advances faster than all trailing waves. It is destructive because the wave advances at a celerity speed greater than the base flow.

Figure 9. The front runner's celerity $c_{FR}/U$ and its dependence on the speed of the front runner $\bar {u}_{FR}/U$ for (a) ${\textit {Fr}} = 0.56$ ($Re = 15.7$, $S_o=0.060$), (b) ${\textit {Fr}} = 0.70$ ($Re = 98.0$, $S_o= 0.015$), (c) ${\textit {Fr}} = 0.84$ ($Re = 35.3$, $S_o=0.060$) and (d) ${\textit {Fr}} = 1.00$ ($Re = 50.0$, $S_o=0.060$).

The square symbol in the figure denotes the limiting condition when $c_{FR}=U+\sqrt {g'H}$ and $\bar {u}_{FR}/U = 1$. The dashed line traces back to the square symbol. For ${\textit {Fr}} = 0.56$, 0.70, 0.84 and 1.00, the slopes of the dashed lines are ${\rm d} c_{FR}/{\rm d} \bar {u}_{FR}=1.10$, 1.21, 1.18 and 1.15, which are nearly constant. The corresponding limiting conditions are ($c_{FR}/U, \bar {u}_{FR}/U) = (1+{{\textit {Fr}}}^{-1}, 1)=(2.79, 1)$, $(2.43, 1)$, $(2.19, 1)$ and $(2.00, 1)$.

The dependence of celerity on amplitude explains the existence of the front runner. The dominant wave with the most significant amplitude advances with the maximum celerity among the waves, and emerges from the pack of waves as the front runner.

Yu & Chu (2022) found that the turbulent roll wave has similar celerity–amplitude relations. Meza & Balakotaiah (2008) studied the dynamics of nonlinear large-amplitude waves on vertically falling films for viscosity and surface tension effects. The celerity–amplitude relationships of the vertical falling films also follow approximately linear relationships, independent of the forcing frequency and amplitude of the disturbance introduced at the inlet. The similarity in the production of waves due to nonlinear instability between our laminar sheet flow, the turbulent sheet flow studied by Yu & Chu (2022) and the vertical falling film problem by Meza & Balakotaiah (2008) is remarkable – despite the difference in inertia, gravity, viscosity and surface tension effects between these problems.

5.4. Bed friction and uplifting velocity at the wavefronts

The difference in results between DNS and SLS occurs mainly in the frontal region of the roll wave where, strictly speaking, the long-wavelength approximation fails. One parameter of significance is the bed friction $\tau _b=[\mu (\partial u/\partial z)]|_{z=0}$. The other is the $w$ velocity component, which also is significant in the frontal region. The $w$ velocity occurs mainly near the free surface. Searching for the maximum $w$ velocity over the depth gives the ‘uplifting velocity’ $\hat {w}$. The significant increases in bed friction and uplifting velocity are also the characteristics of breaking shallow-water waves in a tank (Perlin, He & Bernal 1996; Mostert & Deike 2020; Liu et al. 2023).

Figure 10 shows the variation of the bed-friction shear stress $\tau _{b}$ relative to the undisturbed $\overline {\tau _b}$, $\tau _{b}/\overline {\tau _b}$, in the frontal region over a distance from $(x-x_{FR})/H=-30$ to 10. The blue solid line in the figure shows the bed-friction shear stress obtained by DNS. The red dashed line is the bed friction significantly underestimated by SLS. The SLS model may produce a profile for the bulk body of the roll wave, but is inadequate in describing the roll wave's frontal region. The inability to predict bed-friction shear stress correctly is due to the assumed self-similar parabolic velocity profile in formulating the SLS model.

Figure 10. The bed-shear stress $\tau _b/\overline {\tau _b}$ profiles near the wavefronts for the simulation cases (a) ${\textit {Fr}}=0.700$, ${Re}=98.0$ ($S_o=0.0150$), (b) ${\textit {Fr}}=0.840$, ${Re}=35.3$ ($S_o=0.0600$), (c) ${\textit {Fr}}=1.00$, ${Re}=50.0$ ($S_o=0.0600$) shown in figures 5, 6 and 7, respectively. The blue solid lines and red dashed lines denote the DNS and SLS results, respectively.

The tables at the bottom of figure 10 provide the peak values of the bed friction obtained from DNS and SLS. The DNS give the relative bed friction varying from $\hat {\tau }_{b(DNS)}/\overline {\tau _b} = 1.37$ to 1.95. The substantial increase of bed friction in the frontal region suggests the need for an accurate model to study soil erosion in debris flow and avalanche.

Figure 11 shows the uplifting velocity $\hat {w}/U$ profile in the frontal region obtained from DNS. The uplifting velocity rises sharply before the breaking of the waves. It then becomes slightly chaotic with the air entrapment by the plunging breaker. The uplifting velocity can be upwards (positive) or downwards (negative) after the breaking. The SLS model could not predict the breaking of the wave without the $w$ velocity component in the model.

Figure 11. The peak uplifting velocity $\hat {w}/U$ profile near the wavefronts for the simulation cases (a) ${\textit {Fr}}=0.700$, ${Re}=98.0$ ($S_o=0.0150$), (b) ${\textit {Fr}}=0.840$, ${Re}=35.3$ ($S_o=0.0600$), (c) ${\textit {Fr}}=1.00$, ${Re}=50.0$ ($S_o=0.0600$) shown in figures 5, 6 and 7, respectively.

5.5. Length of the frontal region

The wavefront is a region of excessive bed-friction stress and uplifting velocity. The profiles of ${\tau }_b^{(DNS)}/\overline {\tau _b}$ and $\hat {w}/U$ denoted by the blue lines in figures 10 and 11 show that the excessive bed friction and uplifting velocity occur over the length varying from $(x-x_{FR})/H \approx -10$ to 0. Therefore, the length of the frontal region is approximately $10 H$, ten times the depth of the undisturbed flow.

5.6. Front runner velocity and pressure variation over the depth

The velocity variation over the depth in the frontal region is not precisely parabolic, and the pressure is not hydrostatic. Figures 12(a)–12(c) show the velocity profiles, while figures 12(d–f) show the corresponding DNS pressure profiles at the front runner obtained by the DNS. The dashed lines are the profiles of the undisturbed flow.

Figure 12. The front runner velocity $u_{FR}/\bar {u}_{FR}$ as a function of $z/h_{FR}$, and front runner pressure $p_{FR}/(\rho g^\prime H)$ as a function of $z/h_{FR}$ obtained by DNS. The dashed line delineates the reference profile of the undisturbed flow. (a–c) The velocity profiles; (d–f) the pressure profiles, for (a,d) ${\textit {Fr}}=0.70$, ${Re}=98$, $S_o=0.015$, (b,e) ${\textit {Fr}}=0.84$, ${Re}=35$, $S_o=0.060$, and (c,f) ${\textit {Fr}}=1.0$, ${Re}=50$, $S_o=0.060$.

The front runner region's DNS velocity deviates slightly from the parabolic profile given by (2.9). The velocity gradient near the channel bed exceeds the undisturbed flow, leading to higher bed-friction stress in the region near the wavefront.

The vertical pressure distribution obtained by DNS is approximately linear. However, the value $p_{FR}/(r g^\prime h_{FR})$ at the bottom is slightly smaller than unity, which reflects the effect of the vertical uplift velocity and air entrapment at the wavefront of the front runner.

6. Four-equation WRIBL model

Ruyer-Quil & Manneville (2000) introduced an accurate four-equation WRIBL model. We perform the WRIBL model simulation using the open-source code WaveMaker by Rohlfs et al. (2018). The WRIBL model equations have many high-order derivative terms. Therefore, the spectral method is used for spatial discretization to guarantee simulation convergence, while a fifth-order Runge–Kutta adaptive scheme is used for time integration.

Figures 13 and 14 show the close-up profiles near the wavefronts for the velocity and bed-shear stress in the region near the wavefront, obtained by DNS, two-equation SLS and four-equation WRIBL model simulation for the cases ${\textit {Fr}}=0.600$, ${Re}=72.0$ ($S_o=0.0150$) and ${\textit {Fr}}=0.840$, ${Re}=141$ ($S_o=0.0150$), respectively. The simulation was conducted for the temporal development with a periodic boundary condition over one wavelength $S_o \lambda /H=4.72$. Appendix C provides the profiles over the entire wavelength obtained from the simulations.

Figure 13. The close-up profiles in the region near the wavefront obtained from (a,d) DNS, (b,e) full second-order SLS, and (c,f) WRIBL simulation, for ${\textit {Fr}}=0.600$, ${Re}=72.0$ ($S_o=0.0150$). (a–c) Velocity contours; (d–f) bed-friction stress profiles.

Figure 14. The close-up profiles in the region near the wavefront obtained from (a,d) DNS, (b,e) two-equation SLS, and (c,f) the four-equation WRIBL model simulation, for ${\textit {Fr}}=0.840$, ${Re}=141$ ($S_o=0.0150$). (a–c) Velocity contours; (d–f) bed-friction stress profiles.

Despite the improvement in predicting the bed-friction stress at the wavefront as shown in figures 13(f) and 14(f), the four-equation WRIBL model still does not produce the undular bores and wave breaking predicted by the DNS shown in figures 13(a) and 14(a), respectively. The WRIBL model agrees better with DNS than with SLS. However, the WRIBL simulation wavefront is too steep, unlike any profiles predicted by the DNS.

7. Summary and conclusion

We conduct numerical simulations (DNS) to study roll waves’ temporal and spatial development on a laminar sheet flow with negligible surface tension. The DNS results are compared with the two-equation shallow-layer simulation (SLS) and the four-equation WRIBL model to determine the limitations of these approximation models.

The DNS for temporal development are to show consistency with existing theories of roll waves with different orders of approximations. The existing linear stability theories developed using the normal mode approach are for disturbance of a specified wavelength. The ‘permanent roll wave’ solution obtained by Ng & Mei (1994) also depends on the wavelength.

Roll waves’ spatial development is the primary focus of our DNS study. The spatial development depends not on wavelength specification but on Froude number and distance away from the localized disturbance. Our DNS form the first-ever study of the spatial development of the roll waves on laminar sheet flow with negligible surface tension. The front runner in the spatial development can be a multi-peaked undular bore or a single-peaked non-breaking or breaking wave. The bed friction is significant in the frontal region of the front runner. The frontal length with significant bed friction and uplifting velocity is approximately ten times the undisturbed flow depth.

The two-equation SLS model has produced the bulk profile of roll waves with the least computational cost among the three models. The model fails to predict the significant increase in bed friction in the frontal region of the wave. It slightly underestimates the roll wave amplitude and velocity compared with DNS. Much effort has been made to improve the two-equation shallow-layer model. Ruyer-Quil & Manneville (2000) have proposed a hierarchy of shallow-layer models with different degrees of simplification. Their four-equation shallow-layer model can better predict the frontal region's bed friction. The two-equation SLS model by Ng & Mei (1994), and the four-equation WRIBL model, are both one-dimensional models, solving the $\bar {u}$ component of velocity. These one-dimensional models, however accurate they may be, cannot produce the undular bore or wave breaking phenomena correctly.

The present simulation of the laminar sheet flow ignores surface tension. The results are therefore not comparable to small-scale laboratory experiments by Julien & Hartley (1986) and Liu & Gollub (1994). However, surface tension should be negligible in large-scale geophysical flow where gravity and inertia are the dominant influences. Hunt (1994), for example, has suggested modelling debris flow using a Newtonian model of laminar sheet flow. The modelling of mudflow with negligible surface tension, such as in the works of Ng & Mei (1994), Liu & Mei (1994) and Balmforth & Liu (2004), are closely related to the laminar sheet examined in this paper.

We found that the celerity of the roll waves in laminar sheet flow depends on the waves’ amplitude following a linear relationship. Yu & Chu (2022) found that the turbulent roll wave has similar celerity–amplitude relations. Meza & Balakotaiah (2008) studied the dynamics of nonlinear large-amplitude waves on vertically falling films for viscosity and surface tension effects. The celerity–amplitude relationships of the vertical falling films also follow linear relationships independent of the forcing frequency and amplitude of the disturbance introduced at the inlet. The similarity in the production of waves due to nonlinear instability between our laminar sheet flow, the turbulent sheet flow studied by Yu & Chu (2022) and the vertical falling film problem by Meza & Balakotaiah (2008) is remarkable – despite the difference in inertia, gravity, viscosity and surface tension effects between these problems.

The DNS require significant computation resources, as discussed in Appendix D. On the other hand, the SLS are practical for some applications. The extension of the present simulation study for Newtonian fluid flow to non-Newtonian mudflow is straightforward, and it is our plan for a future investigation.