1. Introduction
The occurrence of falling films in a wide range of industrial and daily-life applications has driven significant interest in the literature and led to comprehensive reviews (see e.g. Chang Reference Chang1994; Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Craster & Matar Reference Craster and Matar2009; Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012). The complex topological features on the surface of such films have fascinated the scientific community since the ground-breaking experiments by Kapitza (Reference Kapitza1948). The desire to isolate the fundamental mechanisms underlying the genesis and development of two-dimensional and three-dimensional waves has led to numerous experimental studies – see, for instance, Kapitza (Reference Kapitza1948), Tailby & Portalski (Reference Tailby and Portalski1962), Liu, Paul & Gollub (Reference Liu, Paul and Gollub1993), Alekseenko et al. (Reference Alekseenko, Nokariakov, Pokusaev and Fukano1994), Liu, Schneider & Gollub (Reference Liu, Schneider and Gollub1995), Oron et al. (Reference Oron, Davis and Bankoff1997), Park & Nosoko (Reference Park and Nosoko2003), Craster & Matar (Reference Craster and Matar2009) and Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012) and references therein. These works have uncovered the generation of ‘families’ of waves and the transition from two- to three-dimensional waves.
Large two-dimensional wave humps dominate the early stages of film development, after which two kinds of secondary transitions help create the spatio-temporal chaos associated with solitary wave structures in falling films (Liu et al. Reference Liu, Paul and Gollub1993; Cheng & Chang Reference Cheng and Chang1995). Chang et al. (Reference Chang, Demekhin, Kalaidin and Ye1996) discussed the presence of streamwise two-dimensional secondary instability leading to the coalescence and coarsening of the initially saturated two-dimensional waves. Additionally, a secondary three-dimensional instability initiates the spanwise transformation of the two-dimensional waves (Joo & Davis Reference Joo and Davis1992). Two avenues for the transition from two- to three-dimensional waves exist depending on the ratio of the spanwise to streamwise noise level at the inlet: an out-of-phase three-dimensional chequerboard evolution of the two-dimensional wave front is observed at sufficiently large cross-stream noise level (Chang et al. Reference Chang, Cheng, Demekhin and Kopelevich1994); and a synchronous horseshoe-shaped modulation of the wave front is seen at weak spanwise noise levels (Liu et al. Reference Liu, Schneider and Gollub1995).
Drawing inspiration from the work of Liu et al. (Reference Liu, Schneider and Gollub1995), Scheid, Ruyer-Quil & Manneville (Reference Scheid, Ruyer-Quil and Manneville2006) developed a low-dimensional weighted residual integral boundary layer model to study the two- to three-dimensional transition and found that the herringbone pattern is largely dependent on the initial conditions. Knowledge of the effect of the synchronous spanwise instability has led to the design of the experiments conducted by Park & Nosoko (Reference Park and Nosoko2003), who were able to isolate the horseshoe-shaped solitary waves of prescribed spanwise and streamwise wavelength, while bypassing the secondary two-dimensional wave dynamics. Using a similar modulation approach, Dietze et al. (Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014) performed numerical simulations of three-dimensional waves, which they then used to provide a comprehensive study of flow structures present within the inertia-dominated large wave hump region as well as within the visco-capillary region.
Surfactants are surface-active species that act to decrease surface tension, additionally introducing variations of this quantity that give rise to Marangoni stresses, which drive fluid away from regions of high surfactant concentration (low surface tension). In the context of falling film flows, surfactants have a stabilising effect, a concept perceived in the early studies on this topic (Tailby & Portalski Reference Tailby and Portalski1962; Benjamin Reference Benjamin1964; Whitaker Reference Whitaker1964). In recent years, linear stability studies have been a primary tool to study the stabilising effect of surfactants on the falling film wave dynamics (Blyth & Pozrikidis Reference Blyth and Pozrikidis2004; Shkadov, Velarde & Shkadova Reference Shkadov, Velarde and Shkadova2004; Wei Reference Wei2005; Pereira & Kalliadasis Reference Pereira and Kalliadasis2008; Karapetsas & Bontozoglou Reference Karapetsas and Bontozoglou2013, Reference Karapetsas and Bontozoglou2014; Bhat & Samanta Reference Bhat and Samanta2018; Hu, Fu & Yang Reference Hu, Fu and Yang2020). More specifically, the role of insoluble surfactants was discussed in greater detail by Blyth & Pozrikidis (Reference Blyth and Pozrikidis2004), who via numerical solutions of the Orr–Sommerfeld equation have shown that surfactants act to stabilise the interface. The authors also comment on the competition of the ever-present interfacial Yih mode and surfactant-concentration-dependent Marangoni mode: at zero Reynolds numbers, the latter is more dominant, but, as the Reynolds number departs from zero, the growth rate of the Yih mode increases, eventually overtaking the Marangoni mode. Additionally, the critical Reynolds number at which the Yih mode produces an instability was shown to increase in the presence of surfactants.
Using a very similar formulation, Pereira & Kalliadasis (Reference Pereira and Kalliadasis2008) found that Marangoni effects dampen the interfacial Kapitza mode and reduce the critical Reynolds number for onset of instability. The authors also acknowledged that surfactants act to reduce the amplitude of free-surface solitary waves and commented on the fact that the emergent dynamics are such that the surfactant accumulates below the primary hump of free-surface solitary waves, especially at low Marangoni numbers. More recently, the effect of insoluble surfactants on the long-wave instability of falling films was investigated by Hu et al. (Reference Hu, Fu and Yang2020), who, using a weighted residual model, found that surfactant surface elasticity decreases the temporal growth rate and increases the critical Reynolds number. The effect of surfactant solubility on the long-wave instability was investigated in detail by Karapetsas & Bontozoglou (Reference Karapetsas and Bontozoglou2014), who found the interfacial concentration gradient to be the strongest for insoluble surfactants and to weaken with increasing surfactant solubility.
Experimentally, the damping of wave activity has been shown by the works of Georgantaki, Vlachogiannis & Bontozoglou (Reference Georgantaki, Vlachogiannis and Bontozoglou2012) and Georgantaki, Vlachogiannis & Bontozoglou (Reference Georgantaki, Vlachogiannis and Bontozoglou2016). More recently, Bobylev et al. (Reference Bobylev, Guzanov, Kvon and Kharlamov2019) investigated the effect of varying concentration of surfactants and observed that at high concentrations the damping effect is reversed, with wave structures beginning to grow again.
In this paper, we study for the first time the effect of insoluble surfactants on vertical falling films in a three-dimensional, nonlinear framework. We implement the same initialisation approach as Dietze et al. (Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014) and demonstrate the emergence of oscillatory behaviour in the developing wave fronts for intermediate values of a parameter that characterises the relative significance of the Marangoni stresses. We use our detailed numerical results to elucidate the mechanism underlying this phenomenon. In § 2, we provide details of the problem formulation along with the numerical technique used to carry out the computation. We discuss our numerical results in § 3 and include those from a validation study wherein our predictions are benchmarked against the experimental observations of Park & Nosoko (Reference Park and Nosoko2003). Finally, concluding remarks are provided in § 4.
2. Problem formulation and numerical method
With the purpose of studying the three-dimensional wave development of vertically falling films in the presence of insoluble surfactant species, we utilise a front-tracking/level-set method – also known as the level contour reconstruction method as in Shin & Juric (Reference Shin and Juric2002, Reference Shin and Juric2009) and Shin, Chergui & Juric (Reference Shin, Chergui and Juric2017). The continuity and momentum equations are solved in a three-dimensional Cartesian domain, 
${\boldsymbol {x}}=(x,y,z)$, as shown in figure 1, with a coupling for the transport of interfacial surfactant species. Bulk surfactant transport is not considered in this work. This is a common assumption in the context of several previous falling liquid film studies (Blyth & Pozrikidis Reference Blyth and Pozrikidis2004; Wei Reference Wei2005; Pereira & Kalliadasis Reference Pereira and Kalliadasis2008; Bhat & Samanta Reference Bhat and Samanta2018; Hu et al. Reference Hu, Fu and Yang2020). Insoluble surfactants are also naturally occurring and we believe that the results from this paper could be of great interest to any experimental set-up that incorporates the use of insoluble surfactants such as NBD-PC (1-palmitoyl-2-
$ \{ $12-[(7-nitro-2-1,3-benzoxadiazol-4-yl)amino]dodecanoyl
$ \} $-sn-glycero-3 -phosphocholine) (Strickland, Shearer & Daniels Reference Strickland, Shearer and Daniels2015; Fallest et al. Reference Fallest, Lichtenberger, Fox and Daniels2010). The gas and liquid are assumed to be immiscible, incompressible Newtonian fluids. The full set of dimensional equations for this method can be found in Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018).
(a) Initial (
$\tilde t=0$) three-dimensional wave profile; (b) schematic representation of the problem in the 
$x$–
$z$ (
$y=0$) plane showing the initial film thickness distribution, 
$\tilde \delta$, and initial streamwise velocity profile, 
$\tilde {u}_x$.
To render the equations dimensionless, we use the following scalings:
\begin{equation} \tilde{\boldsymbol{x}}=\frac{\boldsymbol{x}}{\delta_0}, \quad \tilde{\boldsymbol{u}}=\frac{\boldsymbol{u}} {U_0}, \quad \tilde{t}=\frac{t}{\delta_0/U_0}, \quad \tilde{p}=\frac{p}{\rho_l U_0^2}, \quad \tilde{\sigma}=\frac{\sigma}{\sigma_s}, \quad \tilde{\varGamma}=\frac{\varGamma}{\varGamma_\infty}, \end{equation}
where the tildes designate dimensionless quantities. Here, 
$t$, 
$\boldsymbol {u}$ and 
$p$ denote time, velocity, and pressure, respectively, and the density of the liquid is 
$\rho _l$. The mean velocity and thickness of a waveless falling film, as presented theoretically by Nusselt (Reference Nusselt1923), are designated by 
$U_0$ and 
$\delta _0$, respectively. The surfactant-free surface tension is 
$\sigma _s$, while the surface tension coefficient varying with the interfacial surfactant concentration 
$\varGamma$ is given by an equation of state 
$\sigma =\sigma (\varGamma )$ given below; the concentration at saturation is given by 
$\varGamma _\infty$.
Using the relations in (2.1a–f), the dimensionless forms of the continuity and momentum equations are respectively expressed as
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}=0, \end{equation}and
\begin{align} \tilde{\rho}\left(\frac{\partial \tilde{\boldsymbol{u}}}{\partial \tilde{t}}+\tilde{\boldsymbol{u}} \boldsymbol{\cdot}\boldsymbol{\nabla} \tilde{\boldsymbol{u}}\right) & = - \boldsymbol{\nabla} \tilde{p} + \frac{1}{Re}\boldsymbol{\nabla}\boldsymbol{\cdot} [ \tilde{\mu} (\boldsymbol{\nabla} \tilde{\boldsymbol{u}} +\boldsymbol{\nabla} \tilde{\boldsymbol{u}}^\textrm{T})] +\frac{\boldsymbol{e}_x}{Fr^2} \nonumber\\ &\quad +\frac{1}{We} \int_{\tilde{A}(\tilde{t})} (\tilde{\sigma} \tilde{\kappa} \boldsymbol{n} + \boldsymbol{\nabla}_{\!s\,} \tilde{\sigma} ) \delta(\tilde{\boldsymbol{x}} - \tilde{\boldsymbol{x}}_f)\,\textrm{d}\tilde{A} , \end{align}
where 
$\tilde \kappa$ denotes the interface curvature, 
$\boldsymbol {\nabla }_{\!s}$ the surface gradient operator, 
$\boldsymbol {n}$ the outward-pointing unit normal to the interface, and 
$\boldsymbol {e}_x$ a unit vector in the direction of gravity (i.e. 
$x$-direction). Here 
$\tilde {\boldsymbol {x}}_f$ is the parametrisation of the time-dependent interface area 
$\tilde {A}(\tilde {t})$, where 
$\delta (\tilde {\boldsymbol {x}}-\tilde {\boldsymbol {x}}_f)$ is the three-dimensional Dirac delta function, which vanishes everywhere except at the interface localised at 
$\tilde {\boldsymbol {x}}=\tilde {\boldsymbol {x}}_f$. The density 
$\tilde {\rho }$ and viscosity 
$\tilde {\mu }$ are given by
\begin{equation} \tilde{\rho} (\tilde{\boldsymbol{x}},\tilde{t} ) =\frac{\rho_g}{\rho_l} + \left(1-\frac{\rho_g}{\rho_l}\right) \textrm{H}(\tilde{\boldsymbol{x}},\tilde{t} ), \quad \tilde{\mu} (\tilde{\boldsymbol{x}},\tilde{t} ) =\frac{\mu_g}{\mu_l} + \left(1-\frac{\mu_g}{\mu_l}\right) \textrm{H}(\tilde{\boldsymbol{x}},\tilde{t} ). \end{equation}
Here, 
$\textrm {H}(\tilde {\boldsymbol {x}},\tilde {t})$ represents a smoothed Heaviside function, which is zero in the gas phase and unity in the liquid phase, while the subscripts 
$l$ and 
$g$ designate the individual liquid and gas phases, respectively.
Surfactant transport in the context of insoluble surfactants is governed by
\begin{equation} \frac{\partial \tilde{\varGamma}}{\partial \tilde{t}} + \boldsymbol{\nabla}_{\!s} \boldsymbol{\cdot} (\tilde{\varGamma}\tilde{\boldsymbol{u}}_{\text{t}}) =\frac{1}{Pe_s} {\nabla}^2_{\!s} \tilde{\varGamma}, \end{equation}
where 
$\tilde {\boldsymbol {u}}_{{t}} = (\tilde {\boldsymbol {u}}_{{s}}\boldsymbol {\cdot }\boldsymbol {t})\boldsymbol {t}$ is the tangential velocity vector, in which 
$\tilde {\boldsymbol {u}}_{{s}}$ is the surface velocity and 
${\boldsymbol {t}}$ is the unit tangent to the interface. This type of interfacial transport equation (2.5) follows a similar approach as in Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018). The dimensionless parameters appearing in these equations are given by
\begin{equation} Re=\frac{\rho_lU_0\delta_0}{\mu_l},\quad We= \frac{\rho_l U^2_0 \delta_0}{\sigma_s},\quad Fr=\frac{U_0}{g^{1/2} \delta^{1/2}}, \quad Pe_s=\frac{U_0 \delta_0}{D_s}, \end{equation}
where 
$Re$, 
$We$ and 
$Fr$ are the Reynolds, Weber and Froude numbers, respectively. The gravitational acceleration, assumed to be acting only in 
$x$-direction, is given by 
$g$. The surface Péclet number, 
$Pe_s$, describes the relative importance of surface diffusion to convection, where 
$D_s$ is the diffusion coefficient.
The decrease of 
$\sigma$ in relation to 
$\varGamma$ is modelled using a Langmuir relation (Shin et al. Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018):
\begin{equation} \tilde{\sigma}=1 + \beta_s \ln(1 -\tilde{\varGamma}). \end{equation}
The surfactant elasticity parameter is defined as 
$\beta _s = RT\varGamma _\infty /\sigma _s$, where 
$R$ is the ideal gas constant and 
$T$ is the absolute temperature. Marangoni stresses, which arise from gradients of surface tension, can be expressed in terms of 
$\tilde \varGamma$ by
\begin{equation} \frac{1}{We} \boldsymbol{\nabla}_{\!s\,}\tilde{\sigma} \boldsymbol{\cdot}{\boldsymbol{t}} \equiv \frac{\tilde{\tau}}{We} =-Ma\frac{1}{(1-\tilde{\varGamma})}\boldsymbol{\nabla}_{\!s} \tilde{\varGamma} \boldsymbol{\cdot} {\boldsymbol{t}}, \end{equation}
where 
$Ma\equiv \beta _s/We = RT\varGamma _\infty /\rho _l U_0^2 \delta _0$ is a Marangoni parameter.
The procedure for solving the above equations follows the same approach as in Shin et al. (Reference Shin, Chergui and Juric2017, Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018). In summary, the Navier–Stokes equations (2.3) are solved by a projection method (Chorin Reference Chorin1968), where the temporal discretisation follows a second-order Gear scheme (Wang & Wen Reference Wang and Wen2006) with implicit time integration for the viscous terms. As for the spatial discretisation, a staggered-mesh marker-and-cell (MAC) method (Harlow & Welch Reference Harlow and Welch1965) has been used on a uniform finite-difference grid with a second-order essentially non-oscillatory (ENO) advection scheme (Shu & Osher Reference Shu and Osher1988). The pressure 
$\tilde {p}$ is located at the centres of the cells while the three components of the velocity 
$\tilde {\boldsymbol {u}}$ are located at the faces. All spatial derivatives are approximated using second-order centred differences. The interfacial surfactant transport, (2.5), is treated in a manner consistent with the Lagrangian front-tracking approach and is discretised using triangular interface elements. The detailed numerical solution scheme is presented in Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018).
The numerical set-up of the problem closely follows the previous construct of Dietze et al. (Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014). We impose a periodic boundary condition in both streamwise and spanwise directions of the domain shown in figure 1(a). The bottom wall of the domain is treated as a no-slip boundary, whereas a no-penetration boundary condition is prescribed for the top domain boundary. Following the same formulations for the initial film thickness and streamwise velocity profile as Dietze et al. (Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014), we set
\begin{equation} \tilde \delta|_{t=0}=\bar{\tilde \delta}\left[1+0.2\cos\left(\frac{2{\rm \pi} \tilde x}{\tilde{\lambda}_x}\right)+0.05\cos\left(\frac{2{\rm \pi} \tilde y}{\tilde{\lambda}_y}\right)\right],\quad \tilde u_l=\frac{Re}{Fr^2} (\tilde\delta\tilde{z}-0.5\tilde{z}^2 ), \end{equation}
where 
$\tilde {\lambda }_x$ and 
$\tilde {\lambda }_y$ are the dimensionless domain length and width, respectively, and are motivated by the experimental set-up of Park & Nosoko (Reference Park and Nosoko2003). In their work, Dietze et al. (Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014) acknowledge the importance of 
$\bar {\tilde \delta }$ as a control parameter for liquid volume (since 
$\tilde {V}_l=\bar {\tilde \delta } \tilde {\lambda }_x\tilde {\lambda }_y$), Reynolds number and streamwise wave frequency in numerical cases where periodic boundary conditions are employed. In this work, we have set up the base surfactant-free validation case using the height, 
$\tilde {H}$, and the mean film thickness, 
$\tilde {\bar \delta }$, estimated in the work of Dietze et al. (Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014). As a result, our domain size has the following dimensions: 
$0.025\ \textrm {m} \times 0.02\ \textrm {m} \times 0.0012\ \textrm {{m}}$.
The targeted flow conditions for the surfactant-free base case are from Dietze et al. (Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014): 
$Re=59.3$, 
$We=0.159$, 
$Fr=4.45$ and 
$f=17\ \textrm {Hz}$. The selected fluid properties are representative of an air–water system at 
$25\,^{\circ }\textrm {C}$: 
$\rho _l/\rho _g=841.41$, 
$\mu _l/\mu _g=48.22$ and 
$\sigma _s=0.072\ \textrm {N}\,\textrm {m}^{-1}$. For surfactant-laden cases, the initial surfactant distribution is uniform across the interface (e.g. 
$\tilde {\varGamma }=0.1$) as shown in figure 1(b). The surface Péclet number also remains unchanged from 
$Pe_s=10$ for all surfactant simulations. The Marangoni parameter is the key focus of this work and is varied in the range 0.63–1.88. For the chosen initial surfactant concentration and Marangoni parameters, the maximum overall surface tension reduction is 
$3.2\,\%$, which is sufficiently low to allow us to analyse the effect of Marangoni stresses without the need for running additional simulations where 
$\tilde \tau =0$. To assess the grid dependence of the results, we have used two uniform Cartesian grids for the surfactant-free base case: 
$M_1=768\times 576\times 64$ and 
$M_2=1536\times 1152\times 128$. All surfactant simulations were run using mesh 
$M_1$, unless stated otherwise in the text.
3. Results
We start the discussion of the results by comparing our numerical predictions for the surfactant-free case against the experimental results of Park & Nosoko (Reference Park and Nosoko2003) (see figure 2a,b). The numerical results are given as snapshots of the wave fronts at times corresponding to subsequent periods, whereas the experimental shadowgraph shows the evolution of the waves in the streamwise direction. During the initial stages, the dynamics are dominated by small-amplitude sinusoidal undulations that develop into well-defined ‘horseshoe’ shapes separated by flat regions. We also capture the development of the capillary waves, which precede the horseshoes, and their interactions. It is evident from figure 2(b) that our numerical technique permits the simulation of the wave evolution with good accuracy. The disparity with the experimental observations of Park & Nosoko (Reference Park and Nosoko2003) in terms of overestimation of the horseshoe spacing is similar to that in the work of Dietze et al. (Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014); these authors argued that this is an artefact of the imposed constant spanwise wavelength, 
$\tilde {\lambda }_y$, whereas in the experiments this wavelength varies as the wave fronts travel downstream.
(a) Experimental snapshots from Park & Nosoko (Reference Park and Nosoko2003) and simulation results for (b) surfactant-free and (c–e) surfactant-laden falling films. The times associated with panels (b–e) are: (b) 
$\tilde {t}=(60,116,170,223,277,330)$; (c) 
$Ma=0.63$ at 
$\tilde {t}=(60,116,170,223,625,1111)$; (d) 
$Ma=1.25$ at 
$\tilde {t}=(60,116,170,223,636,794)$; and (e) 
$Ma=1.88$ at 
$\tilde {t}=(60,116,170,223,277,330)$. The surfactant-free and surfactant-laden Reynolds, Weber and Froude numbers are 
$59.3$, 
$0.159$ and 
$4.45$, respectively. The surface Péclet number is kept constant for all surfactant cases at 
$Pe_s=10$.
In addition to the validation against experimental data, we have also carried out mesh sensitivity studies. The two key monitored quantities, the kinetic energy 
$E_k=\int _{V}({\boldsymbol {u}}^2/2)\,\textrm {d}V$ (see figure 7a) and interfacial surface area 
$\tilde {A}$ (see figure 7b), respectively, associated with meshes 
$M_1$ and 
$M_2$ are essentially indistinguishable, highlighting the mesh independence of our results.
Following the validation step, we then proceed to add surfactant species to the flow. It should be noted that all surfactant simulations were run until no further topological changes in the interface shape were detected. Although sinusoidal wave segments dominate the first stages of wave development for all studied cases, significant differences in the individual wave evolution stages can be observed in figure 2(c–e) as we increase the Marangoni parameter. For 
$Ma=0.63$ (see figure 2d), a horseshoe-shaped wave develops similarly to the surfactant-free case; however, its curvature is smaller and the arc connecting its legs, which bulges upwards for the surfactant-free case, is almost completely flattened. A decrease in the number and increase in wavelength of the capillary wave structures preceding the horseshoe-shaped wave can be seen at 
$\tilde {t}=170$ for the surfactant-free (see figure 2b) and surfactant-laden (see figure 2c) cases, respectively. The last two panels for 
$Ma=0.63$ show a divergence of wave development in comparison to what is observed experimentally for the surfactant-free case. At 
$\tilde {t}=625$, we observe that the locations of the horseshoe-shaped and horizontal wave segments are flipped around before the wave stabilises into a quasi-two-dimensional shape at 
$\tilde {t}=1111$. An explanation of this wave development, and the underlying role of the surfactants, is provided below.
Attention is now turned towards the 
$Ma=1.25$ case presented in figure 2(d), where we see that the rise in 
$Ma$ preserves the initial sinusoidal wave pattern at the early stages of wave development (i.e. up to 
$\tilde {t}=223$) until surfactants act to change the mode of the trailing wave segment from three-dimensional to quasi-two-dimensional, bypassing the intermediary step of spanwise oscillation observed for 
$Ma=0.63$. For 
$Ma=1.25$, we also observe further suppression of capillary wave structures. Finally, we can see the dominant effect of the highest Marangoni parameter on the flow development in figure 2(e), where the initially developed sinusoidal shape of the wave is preserved for the entire duration of the simulation, and the capillary wave development is nearly suppressed.
In figure 3, we examine in greater detail the effect of increasing 
$Ma$ on the leading capillary structures and trailing wave fronts. The two-dimensional cut performed in the 
$y=0$ plane reveals that, for all cases investigated, the presence of surfactants suppresses the peak of the trailing wave humps (see figure 3a). This behaviour can be explained further via inspection of the interfacial surfactant concentration in figure 3(b) where we see that the peak of 
$\tilde \varGamma$ observed ahead of the wave crest gives rise to a bi-directional Marangoni stress, which drives flow away from the region of high 
$\tilde \varGamma$ and acts to curb the amplitude of the waves. In figure 3(c), the ‘rigidification’ effect induced by the Marangoni stresses is evident by the reduction of 
$\tilde {u}_t$ of the travelling wave front. We also observe that the suppression of the wave peaks becomes more effective with increasing 
$Ma$.
Effect of varying the Marangoni parameter on two-dimensional projection in the 
$x$–
$z$ (
$\tilde y / \tilde W =0$) plane of film thickness (a), interfacial surfactant concentration (b) and the streamwise component of the interfacial velocity (c), all parameters remaining unchanged from figure 2.
Next our attention is turned towards the effect of surfactants on the capillary wave structures, where magnification of the region in terms of 
$\tilde \delta$ and 
$\tilde \varGamma$ is given in panels 3(a) and 3(b), respectively. For 
$Ma=0.63$ and 
$Ma=1.25$, the capillary wave structures are still present, however, their amplitude is suppressed significantly. This Marangoni-driven damping is caused by a local 
$\tilde \varGamma$ maximum at the peak of each oscillation, which drives the fluid away from the crest. Further evidence in support of the rigidification effect is seen in the decrease of the peak amplitude 
$\tilde {u}_t$ of each capillary structure (see figure 3c). For 
$Ma=1.88$, we observe significant thickening of the capillary region and near-complete elimination of the oscillatory structures (viz. the enlarged insets in figure 3). Further examination of the 
$\tilde \varGamma$ field in figure 3(b) reveals the presence of a concentration gradient that gives rise to a Marangoni stress that drives fluid towards the trough of the trailing hump structure, resulting in the near-complete flattening of the overall wave topology.
The presence and understanding of vortical structures in the trailing wave hump is a critical element when designing efficient heat and mass transfer processes. First we consider the vortical structure in the 
$x$–
$z$ plane shown in the reference frame of the wave crest (see figure 4). The presence of a recirculation zone for the surfactant-free case is associated with increased heat and mass transfer capability (see figure 4a). On the other hand, surfactant addition actively suppresses the recirculation zone of the trailing hump, where similar effects were seen for all 
$Ma$ values examined. The rigidification of the interface brought on by the presence of Marangoni stresses is also the main reason for the uniformity of the vorticity field.
Vortical structure evolution in the reference frame of the wave crest viewed in the 
$x$–
$z$ (
$\tilde y/\tilde W=0$) plane for surfactant-free (a) and surfactant-laden (b) films, respectively, with the darkness of the shading indicating the magnitude of the dimensionless vorticity, 
$|\tilde \varOmega |$. In (b), 
$Ma=1.88$ and all other parameters remain unchanged from figure 2.
In figure 5 the main vortical structures in the 
$x$–
$y$ slice of the large wave humps are displayed in the wave reference frame for two wall distances for the surfactant-free case (see figure 5a,b). Several vortical structures persist over the height of the film. Firstly, two large co-rotating vortices are observed in the horseshoe leg region of the flow. In figure 5(b), the centres of these structures move towards the horizontal section of the large wave. Meanwhile, at a location closer to the wall, the horizontal wave portion is represented by two counter-rotating vortices, which are subsequently disrupted by the movement of the horseshoe vortices. Similar observations concerning the trailing hump vortical structures were made by Dietze et al. (Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014), except for the fourth counter-rotating vortex located at the horizontal wave segment. The authors only observe this vortex to appear as a consequence of the disruptive movement of the large horseshoe co-rotating vortices. Upon addition of surfactants, the three-dimensional flow structures are entirely suppressed. This observation persists across all studied 
$Ma$ values and all wall distances. This result is detrimental for heat and mass transfer operations, as they have been shown to receive a wave-induced boost in three-dimensional films (Alekseenko et al. Reference Alekseenko, Nokariakov, Pokusaev and Fukano1994). The apparent drawback of surfactant-laden flows is the elimination of spanwise flow structures and subsequent suppression of the convective transport in this direction.
Vortical structure evolution in the wave crest reference frame viewed in the 
$x$–
$y$ plane at 
$\tilde z/\tilde H=0.2$ and 
$\tilde z/\tilde H=0.28$ for surfactant-free case (a,b) and surfactant-laden case at 
$\tilde z/\tilde H=0.2$ (c). Vortical structure shown in wave crest reference frame. In (c), 
$Ma=1.88$ and all other parameters remain unchanged from figure 2.
We now examine the spanwise oscillatory motion of the wave structures observed for 
$Ma=0.63$. In figure 6, we present snapshots of the three-dimensional wave shape at 
$\tilde {t}=223$, 625 and 1111 complemented by spanwise two-dimensional representations of the interface, streamwise velocity, 
$\tilde {u}_x$ (where the velocity is given in the reference frame of 
$\tilde {u}_x$ at 
$\tilde {y}=0$), interfacial concentration, 
$\tilde \varGamma$, and arising spanwise Marangoni stresses, 
$\tilde {\tau }$. We observe that the non-uniform distribution of 
$\tilde \varGamma$ at 
$\tilde t=223$ (see figure 6a–c) gives rise to spanwise Marangoni stresses, which drive fluid flow from the horizontal part of the wave hump towards the legs of the horseshoe-shaped wave and also from the tip of the horseshoe towards its legs. The combined effect of the Marangoni stress is to bridge the gap between the tip of the horseshoe and the horizontal portion of the wave. This effect, however, is sufficiently strong so as to promote the development of the middle portion of the wave segment, causing it to accelerate in relation to the adjoining regions, giving rise to a spanwise bulge (see figure 6d–f). A new local peak in 
$\tilde \varGamma$, which coincides with the spanwise peak of the bulge, leads to a 
$\tilde \tau$ structure that induces the final stabilisation of the wave topology (see figure 6g–i). Here, the nearly uniform distribution of 
$\tilde \varGamma$ results in the elimination of all Marangoni stresses.
Three-dimensional wave dynamics, with the darkness of the shading indicating the magnitude of surfactant interfacial concentration, 
$\tilde {\varGamma }$, shown in (a), (d) and (g); two-dimensional projection of interface and streamwise velocity component in the 
$x$–
$y$ plane (
$\tilde z / \tilde H =0.2$) shown in (b), (e) and (h); and 
$\tilde {\varGamma }$ and Marangoni stress in the 
$x$–
$y$ plane (
$\tilde z / \tilde H =0.2$) shown in (c), (f) and (i). Panels (a–c), (d–f) and (g–i) correspond to 
$\tilde {t}=223$, 
$625$ and 
$1111$, respectively, with 
$Ma=0.63$ and the rest of the parameters remain unchanged from figure 2.
Finally, in figure 7, we show the influence of surfactants on the kinetic energy, defined as 
$E_k=\int _V (\rho \boldsymbol {u}^2/2) \,\textrm {d}V$, and the surface area, normalised by their initial values, for the same parameters as in figure 2. Inspection of the kinetic energy plot in figure 7(a) reveals that increasing 
$Ma$ acts to decrease the overall value of 
$\tilde {E}_k$. The amplitude of the oscillations in 
$\tilde {E}_k$ observed at early times for 
$Ma=0.63$ is also all but suppressed with increasing 
$Ma$. A further increase in 
$Ma$ to 
$Ma=1.88$ rigidifies the flow and eliminates completely any oscillation in 
$\tilde {E}_k$. In figure 7(b), we see that the presence of surfactant reduces the initial, linear growth rates in interfacial area for all cases, with this effect becoming particularly pronounced at high 
$Ma$, in line with the recent observations of Hu et al. (Reference Hu, Fu and Yang2020).
Temporal evolution of (a) kinetic energy, 
$\tilde E_k$, and (b) surface area, 
$\tilde A$, scaled on the initial kinetic energy, 
$E_{k0}$, and interface area, 
$A_0$, respectively. The rest of the parameters remain unchanged from figure 2.
4. Concluding remarks
Three-dimensional numerical simulations of vertically falling liquid films in the presence of insoluble surfactants were carried out for the first time. The numerical predictions for the surfactant-free case were benchmarked against the experimental observations of Park & Nosoko (Reference Park and Nosoko2003). For the surfactant-laden case, emphasis was placed on isolating the effect of the Marangoni stresses on the dynamics. The results demonstrate the emergence of oscillations at the wave fronts at low values of the Marangoni parameter, 
$Ma$, mediated by the Marangoni stresses, brought about by spanwise surfactant concentration gradients; the wave fronts eventually evolve into quasi-two-dimensional structures. With increasing 
$Ma$, the Marangoni stresses led to the progressive elimination of the capillary wave structures where near-complete rigidification, and flattening of the liquid film, were observed for sufficiently large 
$Ma$. An increase in 
$Ma$ also resulted in the elimination of vortical structures within the wave crests, and significant reduction in interfacial area and system kinetic energy.
Acknowledgements
This work is supported by the Engineering and Physical Sciences Research Council, UK, through the MEMPHIS (EP/K003976/1) and PREMIERE (EP/T000414/1) Programme Grants, and by computing time at HPC facilities provided by the Research Computing Service of Imperial College London (ICL). D.J. and J.C. acknowledge support through computing time at the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Equipement National de Calcul Intensif) grant no. 2020A0082B06721. The numerical simulations were performed with code BLUE (Shin et al. Reference Shin, Chergui and Juric2017) and the visualisations have been generated using ParaView. The authors also thank U. Farooq (ICL) for fruitful discussions.
Declaration of interests
The authors report no conflict of interest.
