1. Introduction
Surface wave breaking plays a significant role in the exchange of mass, momentum and energy across the air–sea boundary (Melville Reference Mostert and Deike1996; Iafrati Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009; Garbe et al. Reference Garbe2014; Deike Reference Deike2022). When waves break, they generate turbulence, entrain air and enhance mixing, processes that significantly affect oceanic and atmospheric dynamics (Melville Reference Mostert and Deike1996; Deike Reference Deike2022). The underlying physics of wave breaking, particularly spilling and plunging breakers, has been extensively studied using both experiments and numerical simulations (Rapp & Melville Reference Richard, Duran and Fabrèges1990; Chen et al. Reference Chen, Kharif, Zaleski and Li1999; Duncan et al. Reference Duncan1999; Duncan Reference Duncan2001; Qiao & Duncan Reference Ramberg and Griffin2001; Liu & Duncan Reference Liu and Duncan2003, Reference Liu and Duncan2006; Iafrati Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009; Deike et al. Reference Deike, Popinet and Melville2015, Reference Deike, Melville and Popinet2016).
Early studies by Chen et al. (Reference Chen, Kharif, Zaleski and Li1999) showed that vorticity generation and energy dissipation are strongly localised near the wave crest, corroborating experimental observations of wave-induced turbulence. These findings, obtained through two-dimensional simulations, captured plunging and spilling breaker dynamics. Building on this, Iafrati (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009) provided a comprehensive numerical framework for analysing the mechanics of both plunging and spilling breakers. That work classified different wave regimes as a function of the initial wave steepness
$\epsilon$
for fixed Reynolds and Weber numbers (
$\textit{Re}=10^4$
,
$\textit{We}=100$
). Waves remain regular for
$\epsilon \lt 0.33$
, transition to spilling breaker when
$0.33\leqslant \epsilon \lt 0.37$
, and become plunging breakers for
$\epsilon \geqslant 0.37$
. These results agree with the theoretical analysis by Grue & Fructus (Reference Grue and Fructus2010). Building on earlier studies, Deike et al. (Reference Deike, Popinet and Melville2015, Reference Deike, Melville and Popinet2016) used direct numerical simulations to investigate the mechanisms of energy dissipation in breaking waves, highlighting differences between plunging and spilling regimes. Their work emphasised the transition to turbulence, vorticity generation and air entrainment as central features of breaking dynamics. It is important to note that the relationship between initial wave steepness and breaking type is not universal. Experimental and numerical studies show that the limiting steepness for breaking spans a wide range, depending on the generation mechanism, spectral bandwidth, and conditions (Duncan Reference Duncan1981, Reference Duncan, Qiao, Philomin and Wenz1983; Ramberg & Griffin Reference Rapp and Melville1987; Rapp & Melville Reference Richard, Duran and Fabrèges1990; Wu & Yao Reference Wu and Yao2004; Babanin et al. Reference Babanin, Chalikov, Young and Savelyev2010). As reviewed comprehensively by Perlin et al. (Reference Peskin2013), a gravity-capillary wave has been observed over
$\epsilon \approx 0.15$
–
$0.48$
, and in some configurations, spilling breakers may exhibit equal or even greater steepness compared to plunging breakers.
Breaking waves undergo a cascade of processes, from jet impact and vortex formation, to bubble entrainment and droplet ejection, that enhance mixing and momentum transfer across the interface (Liu & Duncan Reference Liu and Duncan2006; Deike Reference Deike2022). High-fidelity simulations by Wang et al. (Reference Wee, Wagoner, Kamat and Basaran2016) captured the formation of bubbles, droplets and spray, offering insight into the multiphase dynamics following wave impact. To model energy dissipation, Mostert & Deike (Reference Mostert, Popinet and Deike2020) developed a theoretical framework relating dissipation rates to wave height, water depth and beach slope, which was validated numerically and extended empirically to offshore conditions. For deep-water plunging breakers, however, inertial dissipation is more appropriately described by the crest-length scaling of Drazen et al. (Reference Drazen, Melville and Lenain2008), which links dissipation to wave kinematics and breaker geometry. In a follow-up work, Mostert et al. (Reference Muradoglu and Tryggvason2022) examined the transition from laminar to turbulent flow, showing that the energy dissipation rates and the size distributions of bubbles and droplets become independent of the Reynolds number once it exceeds a sufficiently high threshold. Recent studies have further examined the dynamics of plunging breakers, such as the work by Di Giorgio et al. (Reference Di Giorgio, Pirozzoli and Iafrati2022), who performed both two- and three-dimensional simulations across different Reynolds numbers, resolving the evolution of the air–water interface and associated bubble structures. Similarly, Hu et al. (Reference Iafrati2024) investigated the influence of wave slope and Reynolds number on spray formation and air entrainment during breaking.
Most of the recent studies on breaking waves have focused on clean, uncontaminated interfaces. However, natural surface-active agents (surfactants), which are ubiquitous in marine environments due to biological activity and pollution, are often present in the ocean surface (Lapham Reference Laxague, Zappa, Soumya and Wurl1998; Laxague et al. Reference Liao, Franses and Basaran2024). Surfactants may be either insoluble (residing primarily at the air–water interface) or soluble (dissolved in the bulk fluid and able to adsorb to and desorb from the interface). In both cases, spatial variations in concentration in the interface result in surface tension gradients, leading to Marangoni-stress-driven flow from regions of low to high surface tension (Manikantan & Squires Reference Melville2020). Surfactant-driven Marangoni stresses can modify interfacial phenomena. For example, in capillary singularities, surfactant-induced flows can retard thread thinning and promote the formation of microthreads (Ambravaneswaran & Basaran Reference Ambravaneswaran and Basaran1999; Craster et al. Reference Craster, Matar and Papageorgiou2002; Liao et al. Reference Liu, Erinin, Liu and Duncan2006; Kamat et al. Reference Kwak and Lee2018; Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020; Wee et al. Reference Whitham2020). In small-scale systems related to wave breaking, such as bubbles and droplets, surfactants have been found to resist jet break-up, slow down the droplet ejection after the bursting of a bubble, and alter aerosol production (Takagi & Matsumoto Reference Teigen, Song, Lowengrub and Voigt2011; Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021c ; Ma et al. Reference Manikantan and Squires2023).
Surfactants alter every stage of wave breaking, such as capillary wave formation, and the structure of the turbulent flow beneath the crest (Ceniceros Reference Ceniceros2003; Liu & Duncan Reference Liu and Duncan2006, Reference Longuet-Higgins2007). Two-dimensional simulations of capillary waves with insoluble surfactant by Ceniceros (Reference Ceniceros2003) showed that Marangoni stresses reduce both the bulge size and the wave amplitude while sharpening the toe curvature. Liu & Duncan (Reference Liu and Duncan2003, Reference Liu and Duncan2006) revealed similar trends with a suppression of capillary ripples, smoother crest evolution, and a smaller, blunter bulge even at modest surfactant concentrations in wave-tank experiments on weak spilling breakers with the soluble surfactant (sodium dodecyl sulphate). By tracking the time-varying bulge geometry, they linked these macroscopic changes directly to surface-rheological properties, concluding that surfactant-induced surface tension gradients delay toe shedding and reduce energy dissipation. Plunging breakers respond to the presence of surfactants even more dramatically. Erinin et al. (Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023) did experiments with Triton X-100 and surfactant-depleted water, and showed that surfactants deform the jet tip, fragment the entrained air cavity, and weaken the downward plunge. They also presented numerical simulations that modelled surfactants as tracer particles along the interface. While these confirmed the qualitative trends, they could not capture the fully coupled surfactant-induced Marangoni stresses. Extending this work, Liu et al. (Reference Liu and Duncan2025) demonstrated that soluble surfactants reduce both droplet count and ejection velocity while shifting the size distribution towards larger diameters. Together, these studies establish that even trace amounts of surfactants suppress small-scale capillary activity, delay the onset of turbulence beneath the crest, and reshape the air entrainment and spray signature of breaking waves.
As revealed by the foregoing review, considerable progress has been made in understanding the role of surfactants in breaking wave dynamics. However, to the best of our knowledge, the specific influence of surfactant-induced Marangoni stresses remains poorly understood. In this study, we investigate the influence of insoluble surfactants on breaking wave dynamics, with a particular focus on regular and spilling regimes, using direct numerical simulations. The rest of this paper is structured as follows. In § 2, the numerical method, governing dimensionless parameters, problem configuration, and validation, is introduced. Section 3 presents the results, and concluding remarks are given in § 4.
2. Problem formulation and numerical method
High-resolution simulations were performed by solving the two-phase incompressible Navier–Stokes equations with surface tension in a three-dimensional Cartesian domain
$\boldsymbol{x} = (x, y, z )$
(see figure 1). The interface between the gas and the liquid is described by a hybrid front-tracking/level-set method, where (insoluble) surfactant transport is explicitly resolved at the interface (Shin et al. Reference Shin and Juric2018). Here, and in what follows, all variables are made dimensionless (represented by tildes) using the characteristic scales
\begin{align} \quad \tilde {\boldsymbol{x}}=\frac {\boldsymbol{x}}{\lambda }, \quad \tilde {t}=\frac {t}{t_{r}}, \quad \tilde {\boldsymbol{u}}=\frac {\boldsymbol{u}} {c}, \quad \tilde {p}=\frac {p}{\rho _{l}c^2} , \quad \tilde {\sigma }=\frac {\sigma }{\sigma _s}, \quad \tilde {\varGamma }=\frac {\varGamma }{\varGamma _\infty }, \end{align}
Schematic representation of the flow configuration. Initial configuration showing the computational domain in a three-dimensional Cartesian coordinate system,
$\boldsymbol{x} = (x, y, z)$
.
where,
$t$
,
$\boldsymbol {u}$
and
$p$
stand for time, velocity and pressure, respectively. The physical parameters correspond to the liquid density
$\rho _l$
, liquid viscosity
$\mu _l$
, and surfactant-free surface tension
$\sigma _s$
. The characteristic velocity is
$c=\sqrt {g\lambda }$
, with
$\lambda$
the wavelength of the wave; hence the characteristic time scale is
$t_{r}= \lambda /c$
. The interfacial surfactant concentration
$\varGamma$
is scaled by the saturation interfacial concentration
$\varGamma _{\infty }$
. As a result of this scaling, the dimensionless equations are
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\tilde {\boldsymbol{u}}=0, \end{align}
\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 {\boldsymbol{i}_z }{\textit{Fr}^2} + \frac {1}{\textit{Re}}\, \boldsymbol{\nabla }\boldsymbol{\cdot }\left [ \tilde {\mu } (\boldsymbol{\nabla }\tilde {\boldsymbol{u}} +\boldsymbol{\nabla }\tilde {\boldsymbol{u}}^{\rm T}) \right ] \notag \end{align*}
\begin{align} {}+\frac {1}{\textit{We}} \int _{\tilde {A}\tilde {(t)}} (\tilde {\sigma } \tilde {\kappa } \boldsymbol {n} + \boldsymbol{\nabla}_{\!s} \tilde {\sigma })\, \delta \left (\tilde {\boldsymbol{x}} - \tilde {\boldsymbol {x}}_{_f} \right )\mbox{d}\tilde {A}, \end{align}
\begin{align} \frac {\partial \tilde {\varGamma }}{\partial \tilde {t}}+\boldsymbol{\nabla}_{\!s} \boldsymbol{\cdot }(\tilde {\varGamma }\tilde {\boldsymbol {u}}_{\text{t}})=\frac {1}{\textit{Pe}_s} {\nabla} ^2_s \tilde {\varGamma }, \end{align}
where the density and viscosity are given by
$\tilde {\rho }=\rho _g/\rho _{l} + (1 -\rho _g/\rho _{l} )\, \mathcal{H} (\tilde {\boldsymbol {x}},\tilde {t} )$
and
$\tilde {\mu }=\mu _g/\mu _{l}+ (1 -\mu _g/\mu _{l} )\, \mathcal{H} ( \tilde {\boldsymbol {x}},\tilde {t} )$
, wherein
$\mathcal{H} ( \tilde {\boldsymbol {x}},\tilde {t} )$
represents a Heaviside function, which is zero in the gas phase and unity in the liquid phase, while the subscripts
$g$
and
$l$
designate the gas and liquid phases, respectively, and
$\tilde {\boldsymbol {u}}_{{t}}= ( \tilde {\boldsymbol {u}}_{{s}} \boldsymbol{\cdot }\boldsymbol {t} ) \boldsymbol {t}$
is the tangential velocity at the interface in which
$\tilde {\boldsymbol {u}}_{{s}}$
represents the interfacial velocity,
$\kappa$
is the curvature, and
$\boldsymbol{t}$
is the unit tangent vector to the interface. The interfacial gradient is given by
$\boldsymbol{\nabla}_{\!s}= ({\boldsymbol{I}}-\boldsymbol{n}\boldsymbol{n} )\boldsymbol{\cdot }\boldsymbol{\nabla}$
, wherein
$\boldsymbol{I}$
is the identity tensor, and
$\boldsymbol{n}$
is the outward-pointing unit normal. In addition,
$\delta$
is a Dirac delta function, equal to unity at the interface and zero otherwise, and
$\tilde {A} (\tilde {t})$
is the time-dependent interface area. As justified by the final paragraph of Stone (Reference Sussman, Smereka and Osher1990), in our frame of reference,
$\boldsymbol {u} \boldsymbol{\cdot }\boldsymbol {n}=0$
, which gives rise to (2.4).
The dimensionless groups that appear in the governing equations are defined as
\begin{align} \textit{Re} =\frac {\rho _l c\lambda} {{\mu _l}}, \quad \textit{We} =g^{1/2}\lambda \sqrt {\frac {\rho _l}{\sigma _s}}, \quad \textit{Fr} =\frac { c }{\sqrt {g \lambda }}, \quad \textit{Pe}_s=\frac { c \lambda }{\mathcal{D}_s},\quad \beta _s= \frac {\mathcal{R} \mathcal{T} \varGamma _\infty }{\sigma _s}, \end{align}
where
$\textit{Re}$
,
$\textit{We}$
,
$\textit{Fr}$
and
$\textit{Pe}_s$
stand for the Reynolds, Weber, Froude and (interfacial) Péclet numbers. As a consequence of our non-dimensionalisation, the phase speed is
$c=\sqrt {g\lambda }$
, so the Froude number becomes
$\textit{Fr}=c/\sqrt {g\lambda }=1$
, consistent with the formulation of Iafrati (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009). We note that a depth-based definition would yield
$\textit{Fr}=\sqrt {2}$
for the present configuration, as described by Treske (Reference Wang, Yang and Stern1994), Chanson (Reference Chanson2009), Gavrilyuk et al. (Reference Gavrilyuk, Liapidevskii and Chesnokov2016) and Richard et al. (Reference Shin, Chergui and Juric2019). The parameter
$\beta _s$
is the surfactant elasticity number that is a measure of the sensitivity of the surface tension
$\sigma$
to the interface surfactant concentration
$\varGamma$
. Here,
$\mathcal{R}$
is the ideal gas constant value (
$\mathcal{R} = 8.314$
J K
$^{-1}$
mol
$^{-1}$
),
$\mathcal{T}$
denotes temperature, and
$\mathcal{D}_s$
stands for the interfacial diffusion coefficient.
The nonlinear Langmuir equation is used to describe
$\sigma$
in terms of
$\varGamma$
:
\begin{align} \tilde {\sigma }=\text{max} [0.05, 1 + \beta _s \ln { (1 -\tilde {\varGamma })]}, \end{align}
where a lower limit of
$\sigma$
has been set to
$0.05$
, below which the Langmuir equation of state may diverge. This formulation is common in previous studies (Muradoglu & Tryggvason Reference Peck and Sigurdson2008; Teigen et al. Reference Temam2011; Atasi et al. Reference Atasi, Haut, Pedrono, Scheid and Legendre2018; Shin et al. Reference Shin and Juric2018; Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021c
). Finally, the Marangoni stress
$\tilde {\tau }$
is expressed as a function of
$\tilde {\varGamma }$
as
\begin{align} \frac {1}{\textit{We}} \boldsymbol{\nabla}_{\!s} \tilde {\sigma }\boldsymbol{\cdot }{\boldsymbol{t}} \equiv \frac {\tilde {\tau }}{\textit{We}} = - \frac {\textit{Ma} }{1-\tilde {\varGamma }}\, \boldsymbol{\nabla}_{\!s}\tilde {\varGamma }\boldsymbol{\cdot }{\boldsymbol{t}}, \end{align}
where
$\textit{Ma}=\beta _s/\textit{We}=\mathcal{R}\, \mathcal{T} \varGamma _\infty /\rho _l c^2\lambda$
is the Marangoni parameter. Tildes are dropped henceforth.
The incompressible Navier–Stokes equation is integrated numerically using a second-order finite difference formulation on a staggered grid (Harlow & Welch Reference Hu, Liu, Yang and Hu1965; Temam Reference Treske1968). The flow field is represented on a fixed, uniform Eulerian mesh. All spatial derivatives are discretised with second-order centred differences, with the exception of the nonlinear advective terms for which a second-order essentially non-oscillatory scheme is employed to ensure numerical stability and limit spurious oscillations (Sussman et al. Reference Takagi and Matsumoto1994). Viscous stresses are likewise evaluated using a second-order centred scheme. The incompressibility constraint is enforced through a fractional-step projection method, in which the pressure is obtained by solving a Poisson equation using a multigrid iterative solver (Chorin Reference Chorin1968; Kwak & Lee Reference Landau and Lifshitz2004). To accurately resolve interfacial dynamics, the numerical framework incorporates an adaptive Lagrangian mesh within a hybrid front-tracking/level-set methodology (Shin & Juric Reference Stone2007; Shin et al. Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2017). Coupling between the Lagrangian interface and the Eulerian flow solver is achieved via an immersed boundary formulation (Peskin Reference Pico, Kahouadji, Shin, Chergui, Juric and Matar1977).
The Lagrangian interface is advanced in time according to
${\rm d}x/{\rm d}t=V$
, where the interface velocity
$V$
is obtained by second-order bilinear interpolation of the Eulerian velocity field. Time integration of the interface position is performed using a second-order Runge–Kutta scheme. Parallel computation is implemented through a domain-decomposition approach, with data communication across subdomains handled using the message passing interface (Shin et al. Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2017). The surfactant conservation equation on the evolving interface is discretised within the hybrid Lagrangian–Eulerian level contour reconstruction method framework, following Muradoglu & Tryggvason (Reference Peck and Sigurdson2008) for the transient and convective terms, which are integrated conservatively over each triangular interface element using the Leibniz transport theorem. The surface tension coefficient
$\sigma$
is evaluated locally on the Lagrangian interface from the surfactant concentration. We refer the reader to Shin et al. (Reference Shin and Juric2018) for more details of the surfactant implementation, including validation against standard benchmark test cases.
2.1. Numerical set-up and parameters
In this work, we study the evolution of a single breaking wave with wavelength
$\lambda$
. The computational domain is defined as
$(\lambda \times \lambda /4 \times \lambda )$
, as shown in figure 1. While previous studies have commonly employed cubic domains of size
$\lambda ^3$
, the reduced spanwise extent adopted here is intentional, as we are not focusing on the strongly three-dimensional dynamics associated with the plunging regime. By reducing the spanwise dimension, we limit the computational cost while retaining three-dimensional dynamics. This approach has also been adopted by Di Giorgio et al. (Reference Di Giorgio, Pirozzoli and Iafrati2025). No variations in surfactant concentration are observed along the spanwise direction, indicating that the dynamics is effectively two-dimensional. The wave propagates in the
$x$
-direction, with periodic boundary conditions in both the streamwise (
$x$
) and spanwise (
$y$
) directions. No-slip and no-penetration conditions are applied at the bottom boundary. We emphasise that the chosen spanwise extent is sufficient to capture the essential dynamics of the breaking wave while maintaining computational efficiency.
The initial free-surface profile
$\eta (x,y)$
follows previous work by Chen et al. (Reference Chen, Kharif, Zaleski and Li1999) and Iafrati (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009), which corresponds to a third-order Stokes wave solution that neglects depth effects and third-order corrections to the fundamental frequency component, computed up to third order in the wave amplitude
$a$
. The surface elevation is prescribed as
\begin{align} \eta (x,y) = \frac {a}{\lambda } \left ( \cos (kx) + \frac {1}{2} \varepsilon \cos (2kx) + \frac {3}{8} \varepsilon ^2 \cos (3kx) \right ), \end{align}
where
$k = 2\pi /\lambda$
is the wavenumber, and
$\varepsilon = ak$
is the initial wave steepness. As noted by Iafrati (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009), the omission of the third-order correction term has a negligible effect on the resulting dynamics; we therefore adopt this initialisation, and focus on surfactant-induced effects. The initial velocity field is obtained under the assumption of a free surface (e.g. negligible air pressure). The velocity profile in the water is derived from the velocity potential associated with the free-surface profile, and is expressed as
\begin{align} u = \varOmega a{\rm e}^{kz} \cos (kx), \quad v = \varOmega a{\rm e}^{kz} \sin (kx), \end{align}
where
$\varOmega = \sqrt {gk(1 + \varepsilon ^2)}$
accounts for the nonlinear correction by removing higher-order terms (Iafrati Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009; Whitham Reference Wu and Yao2011). We note that other studies have considered more physically detailed conditions, in which a velocity profile in the air is also defined (Deike et al. Reference Deike, Popinet and Melville2015). Finally, we note that surface tension forces are not included in the initialisations, but their influence will become relevant once the wave propagates in the streamwise direction.
The dimensionless parameters used in this study are consistent with experimentally realisable systems. The density and viscosity ratios are set to
$\rho _g / \rho _l = 1/850$
and
$\mu _g / \mu _l = (17.4 \times 10^{-6}) / (8.9 \times 10^{-4}) = 1.96 \times 10^{-2}$
, respectively, based on standard properties of air and water. The water depth is set to
$h = \lambda /2$
, consistent with the configuration used by Deike et al. (Reference Deike, Popinet and Melville2015), who showed that this ratio has minimal influence on wave evolution in the regimes of interest. Following Iafrati (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009), we fix the Weber and Reynolds numbers at
$\textit{We} = 100$
and
$\textit{Re} = 10^4$
. The chosen Weber number corresponds to water waves with wavelength approximately
$27$
cm, corresponding to Reynolds number
$\textit{Re} \approx 4.4 \times 10^5$
. Resolving flows at such high Reynolds would require prohibitively fine resolution, which is not feasible within the constraints of our fixed-grid approach. While adaptive mesh refinement has enabled higher-
$\textit{Re}$
simulations in other studies (Deike et al. Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018; Mostert et al. Reference Muradoglu and Tryggvason2022), our focus is limited to regular (non-breaking) and weakly spilling breakers. The chosen parameters therefore strike a balance between computational feasibility and physical realism, while capturing the dominant surfactant-laden dynamics. In this work, we have considered insoluble surfactants whose critical micelle concentration (CMC), i.e.
$\varGamma _\infty \sim \mathcal{O}(10^{-6})$
mol m
$^{-2}$
for NBD-PC (1-palmitoyl-2–12-[(7-nitro-2–1,3-benzoxadiazol-4-yl)amino]dodecanoyl-sn-glycero-3-phosphocholine); thus we have explored the range
$0.1 \lt \beta _s\lt 0.5$
, which corresponds to
$2.9 \times 10^{-6} \lt \mathrm{CMC} \lt 1.4 \times 10^{-5}\,\mathrm{mol\,m^{-2}}$
. The values
$\beta _s = 0.3$
and
$0.5$
were selected to ensure that the interface remains strictly below the critical value of the Langmuir equation of state, so that it never enters the saturation (cut-off) regime. This guarantees that the surface tension varies monotonically with surfactant concentration, and that Marangoni stresses remain active throughout the simulation. To further avoid operating near the saturation limit of the Langmuir equation of state, where Marangoni stresses are suppressed, the initial surfactant concentration is set to
$\varGamma _0 = \varGamma _\infty /4$
in all simulations, ensuring that both surface tension reduction and Marangoni stresses remain active. Finally, the initial concentration of the surfactant-laden cases is initialised uniformly along the interface as in previous work by Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020, Reference Constante-Amores, Kahouadji, Shin, Chergui, Juric, Castrejón-Pita, Matar and Castrejón-Pita2023b
) and Kamat et al. (Reference Kamat, Wagoner, Thete and Basaran2020). We note that a real Stokes wave at the onset of breaking may possess pre-existing surfactant gradients; our choice of a uniform initial concentration is an idealisation. We have set
$\textit{Pe}_s=10^3$
following Batchvarov et al. (Reference Batchvarov, Kahouadji, Magnini, Constante-Amores, Craster, Shin, Chergui, Juric and Matar2020) and Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020), who showed that the interfacial dynamics is weakly dependent on
$\textit{Pe}_s$
beyond this value.
The simulations presented here are mesh-independent. A resolution
$\lambda /\Delta \boldsymbol{x}=512$
was found sufficient to ensure that further mesh refinement does not alter the numerical results. This resolution is consistent with that used in previous studies at the same dimensionless parameters (Iafrati Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009). We also note that the simulations in the plunging regime are terminated upon the onset of plunging. Although the numerical framework is capable of resolving post-breaking behaviour, this is not the focus of the current study. Further investigation would require refining the simulations, which would significantly increase computational cost, and consequently we do not explore this avenue here.
The numerical validation of the surfactant equations has been previously presented in Shin et al. (Reference Shin and Juric2018). Extensive mesh studies for surface-tension-driven phenomena using the same computational method can be found elsewhere (e.g. Batchvarov et al. Reference Batchvarov, Kahouadji, Magnini, Constante-Amores, Craster, Shin, Chergui, Juric and Matar2020, Reference Batchvarov, Kahouadji, Constante-Amores, Norões Gonçalves, Shin, Chergui, Juric, Craster and Matar2021; Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020, Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021b
, Reference Constante-Amores, Chergui, Shin, Juric, Castrejón-Pita and Castrejón-Pita2022, Reference Constante-Amores, Abadie, Kahouadji, Shin, Chergui, Juric, Castrejón-Pita and Matar2023a
). Additionally, liquid volume and surfactant mass conservation are met under errors of less than
$10^{-1}\,\%$
and
$10^{-2}\,\%$
, respectively.
3. Results
Table 1 summarises the observed wave behaviour as a function of the steepness parameter
$\epsilon$
with
$\textit{Re}=10^4$
and
$\textit{We}=100$
. We observe that for surfactant-free cases, the wave remains regular for
$\epsilon = 0.3$
, transitions to a spilling-type break-up for
$\epsilon = 0.324$
, and finally exhibits plunging-type break-up for
$\epsilon = 0.37$
. These findings agree with previous work by Iafrati (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009) and Grue & Fructus (Reference Harlow and Welch2010). We have included simulations with
$\epsilon =0.324$
, which aligns with the onset from regular to spilling-type breaking predicted when accounting for the secular growth term
$\epsilon ^3/8$
by Grue & Fructus (Reference Harlow and Welch2010). It should be noted that for
$\epsilon = 0.37$
, the simulation was terminated just prior to the break-up. This choice was made to focus the analysis on the wave behaviour leading up to the event. A detailed study of the post-break-up dynamics in the plunging regime is beyond the scope of this study.
Wave behaviour as a function of wave steepness
$\epsilon$
and surfactant parameter
$\beta _s$
. The steepness
$\epsilon$
denotes the initial Stokes-wave steepness prescribed in (2.8). The surfactant parameter
$\beta _s$
denotes surfactant elasticity.
Table 1. Long description
The table has four columns and twelve rows. The columns are labeled with the steepness parameter, the surfactant parameter, and three types of wave behavior: Regular, Spilling, and Plunging. The rows list different values for the steepness parameter and surfactant parameter, indicating the observed wave behavior with an X mark. Row 1: Epsilon 0.3, Beta S 0.3, Surfactant-free, Regular X. Row 2: Epsilon 0.3, Beta S 0.3, Spilling −. Row 3: Epsilon 0.3, Beta S 0.5, Regular X. Row 4: Epsilon 0.3, Beta S 0.5, Spilling −. Row 5: Epsilon 0.324, Beta S Surfactant-free, Regular −. Row 6: Epsilon 0.324, Beta S 0.3, Spilling X. Row 7: Epsilon 0.324, Beta S 0.5, Spilling −. Row 8: Epsilon 0.33, Beta S Surfactant-free, Spilling X. Row 9: Epsilon 0.33, Beta S 0.3, Spilling −. Row 10: Epsilon 0.33, Beta S 0.5, Spilling −. Row 11: Epsilon 0.35, Beta S Surfactant-free, Spilling X. Row 12: Epsilon 0.35, Beta S 0.3, Spilling −. Row 13: Epsilon 0.35, Beta S 0.5, Spilling −. Row 14: Epsilon 0.37, Beta S Surfactant-free, Plunging X.
To investigate the role of wave breaking in energy dissipation, we track the temporal evolution of the kinetic energy. Figure 2 presents the kinetic energy
$E_k=\rho\int_V(\boldsymbol{u}^2/2)\,\mathrm{d}V$
, normalised by its initial value
$E_{k0}$
, as a function of the initial wave steepness
$\epsilon$
. The figure reveals distinct trends corresponding to the different breaking regimes, highlighting the increasing intensity of energy decay with more intense wave breaking. In the regular wave regime, a nearly constant decay rate is observed. The theoretical decay of mechanical energy due to viscous damping is given by
$E_{mech}(t)=E_{mech}(0)\,{\rm e}^{-2\gamma t}$
, where
$E_{mech}$
is the mechanical energy (e.g. equal to twice the kinetic energy for small-amplitude waves), and
$\gamma =2k^2\mu _l/\rho _l$
is the damping coefficient (Landau & Lifshitz Reference Lapham1987); see Appendix D for more details. For the wave with steepness
$\epsilon =0.3$
, the numerical simulation results in
$\gamma =0.0091$
, which is in good agreement with the theoretical estimate
$\gamma \approx 0.0079$
. At higher steepness
$\epsilon =0.37$
, a monotonic decrease in kinetic energy is observed, preceded by a phase of increasing amplitude and enhanced dissipation. Although the simulation is terminated just before the onset of plunging-type break-up, the flow already exhibits intensified velocities and steep velocity gradients near the interface, resulting in large kinetic energy and enhanced energy dissipation. Nonlinear energy exchange in the finite-amplitude wave produces oscillations in
$E_k/E_{k0}$
, so the linear dissipation theory applies only to the average decay rate rather than the instantaneous dynamics.
Time evolution of kinetic energy as a function of the wave steepness
$\epsilon$
for surfactant-free cases.
The changes in wave morphology that accompany these energy dissipation trends are illustrated in figures 3(a–c), which show the spatio-temporal evolution of the interface in the
$(x{-}z)$
plane at
$y = \lambda /8$
for initial steepness values
$\epsilon = 0.3, 0.33, 0.35$
. At
$\epsilon = 0.3$
, the wave evolution remains regular over time. As the steepness increases to
$\epsilon = 0.33$
, weak signs of spilling emerge in the surfactant-free case, which become more pronounced at
$\epsilon = 0.35$
, where the formation of a forward-facing bulge near the wave crest marks the onset of spilling. The leading edge of this bulge, known as the toe, curves inwards due to surface tension effects, consistent with the mechanism described by Duncan (Reference Duncan2001). Capillary waves are also formed upstream of the toe, as reported by Duncan et al. (Reference Duncan1999). For gravity-capillary waves, Deike et al. (Reference Deike, Popinet and Melville2015) reported that the breaking threshold depends primarily on the Bond (or Weber) number and the steepness. In the large Weber number regime relevant to our simulations (
$\textit{We} =100$
,
$Bo \sim 253$
), their scaling predicts a critical steepness in the range
$\epsilon _{\textit{crit}} \approx 0.32$
. The transition observed here in the clean-wave case at
$\epsilon _{\textit{crit}} \approx 0.324$
therefore lies within the gravity-capillary breaking regime reported by Deike et al. (Reference Deike, Popinet and Melville2015). We note, however, that their simulations were performed at a significantly higher Reynolds number (
$\textit{Re} \approx 4\times 10^{4}$
) compared with the present
$\textit{Re} = 10^{4}$
, so that the agreement should not be interpreted as a strict quantitative match, but rather as a consistency of physical regimes.
Time–space plots of the interface in the
$(x{-}z)$
plane (
$y=\lambda /8$
) at equal time intervals for varying wave steepnesses
$\epsilon$
, and surfactant elasticity
$\beta _s$
. The vertical coordinate
$z$
are stretched by a factor of 4 and shifted upwards by the corresponding time to better illustrate the interface evolution. Coloured lines mark selected times analysed in detail in the text.
Figures 3(d–i) show the effect of insoluble surfactants on the wave dynamics. At
$\epsilon = 0.3$
(see figures 3
d,g), the addition of surfactant produces minimal changes compared to the surfactant-free case. This is expected, as the steepness remains below the critical value for Stokes waves (
$\epsilon _{crit} \approx 0.32$
); although surfactants reduce the surface tension, it is insufficient to generate a noticeable wave deformation. When the wave enters the spilling regime, there are no visually obvious differences between surfactant cases and the surfactant-free case since the spilling signatures remain weak at this stage. The enhancement of spilling signs can be observed for higher steepness (
$\epsilon = 0.35$
), with the effect being more pronounced at higher values of
$\beta _s$
(see figures 3
f,i).
To elucidate the coupling between the interface dynamics and surfactant concentration, figure 4 presents the spatio-temporal evolution of
$\varGamma$
at the same times as in figure 3. A clear correlation emerges between the interfacial shape and the distribution of
$\varGamma$
. As the wave evolves, surfactant initially accumulates at troughs due to surface convergence, then shifts towards the crests (where the interface exhibits divergent flow). Once the crest is formed, the local flow divergence redistributes
$\varGamma$
, causing it to diffuse away gradually. Comparing cases with
$\beta _s = 0.3$
and
$\beta _s = 0.5$
, the crest in the higher
$\beta _s$
leans forwards more prominently, indicative of surfactant-induced modifications to the interfacial dynamics and wave deformation. The largest accumulation of surfactant concentration occurs in the
$\epsilon = 0.35$
case, where enhanced crest steepening and incipient spilling intensify interfacial compression, leading to sharper concentration gradients and increased accumulation of
$\varGamma$
. Noticeably, a double-peaked distribution of
$\varGamma$
emerges (see figures 4
c, f), originating from interfacial overturning (see figures 3
f,i).
Time–space plots of normalised surfactant concentration
$\varGamma$
at time intervals as a function of the steepness
$\epsilon$
, and
$\beta _s$
. The data are plotted against arc length
$s \in [0,1]$
, with vertical coordinates representing the normalised surfactant concentration
$\varGamma$
shifted by the corresponding time to clearly show the evolution of surfactant distribution. For the surfactant-laden cases,
$\textit{Pe}_s=10^3$
and
$\varGamma_0=\varGamma _\infty /4$
. Coloured lines mark selected times analysed in detail in the text.
We now examine the coupling between the interface location, surfactant concentration
$\varGamma$
, Marangoni stress, and tangential velocity in the streamwise direction
$u_{t,x}$
. To elucidate the underlying mechanisms, we select representative time instances at which wave features are most pronounced. These times, corresponding to different steepness regimes, are indicated in different colours in figures 3 and 4. We begin with the regular wave regime. For
$\epsilon = 0.3$
, figure 5(a) shows that the presence of surfactant has a minimal effect on the interface location at late times (
$t = 5.6$
). However, figure 5(b) reveals that the surfactant concentration peaks at the wave crest, where surface tension is locally minimised. The resulting Marangoni stresses, driven by surface tension gradients, generate outward tangential flow away from the crest, i.e.
$\tau \gt 0$
on the forward part of the crest, and
$\tau \lt 0$
on the rear part of it. These Marangoni stresses act to suppress further steepening of the wave by opposing interfacial deformation, thereby reducing the wave amplitude at late times (not shown). The influence of Marangoni stresses on the surface flow is evident in figure 5(d), which shows that
$u_{t,x}$
is significantly enhanced near the crest in the presence of surfactants, compared to the clean case, particularly for the high elasticity case (
$\beta _s = 0.5$
). Despite the enhanced tangential velocity and the non-uniform surfactant distribution, the interface displacement remains small at this particular time. At
$\epsilon = 0.3$
, the inertial forces associated with wave propagation are insufficient to drive significant interface deformation or wave overturn, and the flow remains in a weakly nonlinear regime, where the balance between gravity and surface tension dominates. Moreover, surface tension remains relatively strong and effectively resists interface steepening. As a result, the Marangoni-induced tangential flow primarily redistributes momentum along the interface without producing noticeable changes in the wave shape.
Effect of the elasticity parameter
$\beta _s$
for the regular wave regime (
$\epsilon =0.3$
) at
$t=5.6$
. Two-dimensional projections of the interface,
$\varGamma$
,
$\tau$
and
$u_{t,x}$
in the
$(x{-}z)$
plane (
$y=\lambda /8$
) are shown in (a–d), respectively. Note that the abscissa corresponds in (a) to the
$x$
coordinate, and in (b–d) to the arc length
$s$
.
Figure 5. Long description
The image contains four line graphs labeled Panel A, Panel B, Panel C, and Panel D. Panel A shows a line graph with the horizontal axis labeled x and the vertical axis labeled z. The graph displays a sinusoidal wave pattern with a zoomed-in inset highlighting a specific range of x values. Panel B presents a line graph with the horizontal axis labeled s and the vertical axis labeled Γ. It shows two curves, one solid blue and one dashed red, both exhibiting a parabolic trend. Panel C features a line graph with the horizontal axis labeled s and the vertical axis labeled τ. It includes two curves, one solid blue and one dashed red, showing oscillatory behavior. Panel D illustrates a line graph with the horizontal axis labeled s and the vertical axis labeled u_t_x. It includes three curves: a solid black line representing surfactant-free conditions, a dashed red line for β = 0.3, and a solid blue line for β = 0.5. An inset zooms in on a specific range of s values, highlighting differences among the curves.
We next consider the spilling regime; the interfacial location,
$\varGamma$
,
$\tau$
and
$u_{t,x}$
are shown in figure 6. For
$\epsilon = 0.324$
, the addition of surfactant significantly modifies the shape and position of the wave crest (see figure 6
a). Surfactants also accumulate near the wave crest, leading to a slightly sharper, and more forward-leaning profile (figure 6
b). This effect becomes more pronounced with increasing surfactant elasticity, which enhances interfacial gradients and promotes stronger Marangoni stresses. The resulting non-uniform surfactant distribution leads to localised peaks in Marangoni stresses near the crest (figure 6
c), and larger tangential velocities
$u_{t,x}$
(figure 6
d). At slightly higher steepness,
$\epsilon = 0.33$
, the addition of surfactant induces significant changes in the interfacial dynamics at
$t=5.2$
. As shown in figure 6(e), a small forward jet forms near the crest. Liu & Duncan (Reference Liu and Duncan2006) observed experimentally that in surfactant-laden spilling breakers, the bulge shape flattens and becomes elongated along the slope of the front face of the wave. We observe similar qualitative behaviour in our simulations, and acknowledge that quantitative agreement is not expected due to differences in Reynolds number and surfactant solubility; the geometric trends are consistent with the experimental observations (more details are given in Appendix E). We also note that Erinin et al. (Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023) observed experimentally that the addition of surfactants leads to a jet that curls inwards for plunging breakers. The formation of this forward jet significantly reshapes both the surfactant concentration and tangential velocity fields. For
$\beta _s = 0.5$
, the surfactant distribution profile (figure 6
f) peaks at the wave crest, with larger Marangoni stresses on the front side of the crest. The resulting Marangoni stress (figure 6
g) increases with
$\beta _s$
and promotes jet elongation. A negative Marangoni-stress peak is observed, which further promotes this elongation. Figure 6(h) shows the tangential velocity field
$u_{t,x}$
along the interface. In the surfactant-free case, the tangential surface velocity first accelerates along the front face of the wave, but undergoes a localised deceleration near the jet base, where momentum is redirected into the forming jet. Surfactants redistribute surface momentum by promoting upstream acceleration along the wave front while inducing a sharply localised deceleration near the jet base through Marangoni stresses acting over a short arc length. By partially rigidifying the interface near the jet base, surfactants promote upstream surface stretching, which favours elongation of the bulge along the wave slope, and suppresses its vertical growth, consistent with experimental observations.
Effect of the elasticity parameter for the spilling regime for
$\epsilon =0.324$
at
$t=4.8$
, and
$\epsilon =0.33$
at
$t=5.2$
. Two-dimensional projections of the interface,
$\varGamma$
,
$\tau$
and
$u_{t,x}$
in the
$(x{-}z)$
plane (
$y=\lambda /8$
) are shown, for (a–d)
$\epsilon = 0.324$
and (e–h)
$\epsilon = 0.33$
.
Finally, in figure 7, we examine the effect of surfactant in greater detail by analysing the tangential velocity in the reference frame of the wave crest for both surfactant-free and surfactant-laden cases with
$\epsilon =0.33$
and
$\beta _s=0.5$
at
$t = 5.2, 5.6, 6.0$
. During the onset of spilling, a local minimum develops in the tangential velocity profile near the base of the crest – a feature that is enhanced in the presence of surfactants for spilling breakers. In the surfactant-free case, this secondary peak decays rapidly, and has disappeared by
$t=6.0$
(see figure 7
c). In contrast, for the surfactant-laden case, the secondary peak remains clearly visible at the same instant (see figure 7
f), indicating a prolonged effect of surfactant-driven Marangoni stresses near the crest. Although the crest speed is constant, the persistence and increased magnitude in the surfactant case reflect enhanced tangential shear along the interface. The extended lifetime and strengthened secondary peak therefore provide direct evidence that surfactants promote and intensify spilling by sustaining near-crest shear through Marangoni stresses.
Two-dimensional projections of the interface location and tangential velocity in the reference frame of the wave crest for (a–c) surfactant-free with
$\epsilon =0.33$
, (d–f) surfactant-laden with
$\beta _s=0.5$
, and (g–i) surfactant-laden without Marangoni stresses (i.e.
$\boldsymbol{\nabla}_{\!s} \tilde {\sigma }=0$
in (2.3)) with
$\beta _s=0.5$
, at times (a,d,g)
$t=5.2$
, (b,e,h)
$t=5.6$
, (c, f,i)
$t=6.0$
. Here,
$c_t$
denotes the tangential component of the average interfacial velocity.
Figure 7. Long description
The image contains nine line graphs arranged in a 3x3 grid, each depicting the interface location and tangential velocity in the reference frame of the wave crest for different conditions and times. Panel A, D, and G show data for surfactant-free conditions at time t = 0. Panel B, E, and H show data at time t = 0.5. Panel C, F, and I show data at time t = 1. Each panel has two y-axes: the left y-axis represents the interface location (z) ranging from −0.05 to 0.10, and the right y-axis represents the tangential velocity (u_t,x) and the difference between the tangential velocity and a constant (u_t,x − c_l) ranging from −0.5 to 0.5. The x-axis represents the horizontal position (x, s) ranging from −0.5 to 0.5 in Panels A, B, and C, and from 0 to 1 in Panels D, E, F, G, H, and I. The graphs show the evolution of the interface location and tangential velocity over time for surfactant-free conditions, surfactant-laden conditions, and surfactant-laden conditions without Marangoni stresses.
To isolate the role of Marangoni stresses in the interfacial dynamics, we perform additional simulations in which Marangoni stresses are suppressed by imposing
$\boldsymbol{\nabla}_{\!s} \sigma = 0$
in (2.3). This constraint allows surfactant to remain non-uniformly distributed along the interface while preventing the generation of Marangoni stresses; a similar approach has been used by Kamat et al. (Reference Kamat, Wagoner, Thete and Basaran2020) and Constante-Amores et al. (Reference Constante-Amores, Chergui, Shin, Juric, Castrejón-Pita and Castrejón-Pita2022, Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021c
). Figures 7(g–i) show the cases for
$\boldsymbol{\nabla}_{\!s} \sigma = 0$
with
$\epsilon =0.33$
and
$\beta _s=0.5$
at times
$t=5.2, 5.6, 6.0$
, allowing the distinct roles of Marangoni stresses and surface tension variations to be inferred in the dynamics. Comparing figures 7(d–f) (Marangoni stresses allowed) with figures 7(g–i) (no Marangoni stresses) shows that the presence of Marangoni stresses enhances the second peak in the tangential velocity distribution, which is associated with the onset of spilling. In contrast, comparison with the surfactant-free cases (figures 7
a–c) shows that a local reduction in surface tension, in the absence of Marangoni stresses, suppresses the secondary tangential velocity peak relative to the surfactant-free case. This demonstrates that a decrease in the local surface tension alone does not promote spilling. Instead, the enhancement of spilling observed in surfactant-laden simulations arises from Marangoni stresses, which generate tangential stress gradients
$\boldsymbol{\nabla}_{\!s} \sigma (\varGamma )$
that locally accelerate the interface and amplify the secondary velocity peak.
Figure 8 shows the spanwise vorticity
$\omega _y = (\boldsymbol{\nabla }\times \boldsymbol{u})_y$
in an
$(x{-}z)$
plane overlaid with the interface shape for steepness
$\epsilon = 0.324$
(figures 8
a–c) and
$\epsilon =0.33$
(figures 8
d–h). For each case, results are shown for the surfactant-free and the surfactant-laden cases with
$\beta_s=0.3$
and
$\beta _s=0.5$
(figures 8(a,d), 8(b,e,g) and 8(c, f,h), respectively). For
$\epsilon = 0.324$
, the surfactant-free case exhibits a smooth flow separation near the wave crest, forming a broad region of negative
$\omega _y$
loosely coupled to the interface (see figure 8
a). At moderate elasticity number
$\beta _s=0.3$
, a distinct interfacial shear layer forms beneath the interface due to Marangoni-driven flow, characterised by adjacent bands of positive and negative
$\omega _y$
straddling the interface, a signature of Marangoni-driven vorticity generation (see figure 8
b). At higher elasticity
$\beta _s=0.5$
, this vorticity layer becomes more localised, and a secondary vorticity layer of opposite sign develops beneath the interface, a feature absent in the surfactant-free case. For the slightly larger steepness (
$\epsilon = 0.33$
), separation is more abrupt, and early vortex roll-up occurs near the crest (see figure 8
d). The presence of surfactant significantly amplifies these effects. The interfacial shear layer strengthens, and the vorticity becomes more intense and confined in regions of high
$\varGamma$
(see figures 8
e, f). At
$\beta _s=0.5$
, we observe a larger positive vorticity near the base of the wave, coinciding with the location of strongest surfactant accumulation. This correlation is consistent with enhanced interfacial shear and sharper velocity gradients induced by Marangoni stresses acting along the interface. In addition, figures 8(g,h) show the corresponding flow fields when Marangoni stresses are suppressed,
$\boldsymbol{\nabla}_{\!s} \sigma = 0$
. In this case, the positive vorticity generated near the crest region is significantly reduced compared to the surfactant-free case. This reduction becomes more pronounced as the elasticity parameter increases from
$\beta _s = 0.3$
to
$\beta _s = 0.5$
, reflecting the stronger decrease in local surface tension at the crest for a given surfactant concentration. These results indicate that a reduction in surface tension alone does not enhance the shear layer at the crest; instead, it weakens vorticity generation. The amplification of vorticity observed in the full surfactant cases therefore requires the presence of Marangoni stresses, which introduce tangential stress gradients along the interface.
Two-dimensional representation of the interfacial location (coloured by
$\varGamma$
) overlaid with the spanwise vorticity component
$\omega _y$
in the
$(x{-}z)$
plane, with (a–c)
$\epsilon = 0.324$
at
$t = 4.8$
, (d–f)
$\epsilon = 0.33$
at
$t = 5.2$
, and (g,h)
$\epsilon =0.33$
at
$t = 5.2$
, and
$\boldsymbol{\nabla}_{\!s} \tilde {\sigma }=0$
in (2.3), for (a,d) the surfactant-free case, (b,e,g)
$\beta _s = 0.3$
and (c, f,h)
$\beta _s = 0.5$
.
Finally, we examine the role of Marangoni stresses in the generation of vorticity at the gas–liquid interface. Vorticity production at free surfaces has been the subject of several studies (Longuet-Higgins Reference Lundgren and Koumoutsakos1992; Cresswell & Morton Reference Cresswell and Morton1995; Peck & Sigurdson Reference Perlin, Choi and Tian1998; Lundgren & Koumoutsakos Reference Ma, Hessenkemper, Lucas and Bragg1999; Brøns et al. Reference Brøns, Thompson, Leweke and Hourigan2014; Constante-Amores et al. Reference Constante-Amores, Batchvarov, Kahouadji, Shin, Chergui, Juric and Matar2021a ), and is particularly relevant in air–water systems where the free-surface approximation is appropriate. In this work, we adopt a two-dimensional version of the generalised formulation derived by Constante-Amores et al. (Reference Constante-Amores, Abadie, Kahouadji, Shin, Chergui, Juric, Castrejón-Pita and Matar2023a ), who described the evolution of circulation at a deformable interface within a three-dimensional control volume in the presence of surfactants. For brevity, we present only the final form of the governing equation used in the analysis in the main text; a full derivation of the two-dimensional form is provided in Appendix A. The circulation is given by
\begin{align} \frac {{\rm D}}{{\rm D}t}\left (\int _A \omega \,{\rm d}A + \int _a^b [[\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t}]]\,{\rm d}s\right ) &= \oint _C \nu \,\boldsymbol{\nabla }\omega \boldsymbol{\cdot }\boldsymbol{n}_c\,{\rm d}s - \int _a^b \kappa \, [[(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t})]]\,{\rm d}s \nonumber \\ &\quad {}+ \frac {1}{2}\int _a^b \frac {\partial }{\partial s} [[(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})^{2}]]\,{\rm d}s - \int _a^b \frac {\partial }{\partial s} \left [ \left [\frac {p}{\rho }\right ] \right ] {\rm d}s. \end{align}
Equation (3.1) describes the rate of change of circulation within a two-dimensional control volume intersecting an arbitrary fluid interface of area
$A$
separating two fluids with different material properties. The left-hand side consists of two terms: the integral of vorticity over
$A$
, and the circulation jump along the curve from
$a$
to
$b$
, represented by the jump in tangential velocity
$[[\boldsymbol{u} \boldsymbol{\cdot }\hat {\boldsymbol{t}}]]$
. However, under the assumption of velocity continuity across the interface in both the normal and tangential directions, and given that the unit normal
$\boldsymbol{n}$
and tangential
$\boldsymbol{t}$
vectors are identical on either side of the interface, the jump in tangential velocity
$[[\boldsymbol{u} \boldsymbol{\cdot }{\boldsymbol{t}}]]$
vanishes. Similarly, the second term on the right-hand side involving the jump of the product
$(\boldsymbol{u} \boldsymbol{\cdot }{\boldsymbol{n}})(\boldsymbol{u} \boldsymbol{\cdot }{\boldsymbol{t}})$
, and the third term involving the jump of
$(\boldsymbol{u} \boldsymbol{\cdot }{\boldsymbol{n}})^2$
, also vanish. Therefore, the second term on the left-hand side and the second and third terms on the right-hand side of (3.1) are identically zero. The last right-hand-side term of (3.1) involves a jump in the pressure across the interface. Taking into account the jump in the normal stress condition at the interface, we express
$\partial [[p/\rho ]]/\partial s$
as
\begin{align} \frac {\partial }{\partial s}\left[ \left[\frac {p}{\rho }\right] \right]=-\frac {1}{\rho _1}\frac {\partial }{\partial s}(\kappa \sigma )-\frac {2}{\rho _1}\left[ \left[\mu \left (\frac {\partial ^2}{\partial s^2}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t})+\frac {\partial }{\partial s}(\kappa (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}))\right )\right] \right]+\left[ \left[\frac {1}{\rho }\right] \right]\frac {\partial p_2}{\partial s}. \end{align}
Therefore, the final expression for the circulation is
\begin{align} \frac {{\rm D}}{{\rm D}t}\left(\int _A \omega \, {\rm d}A \right) &= \oint _C \nu\, \boldsymbol{\nabla }\omega \boldsymbol{\cdot }\boldsymbol{n}\, {\rm d}s + \int _a^b \frac {\sigma }{\rho _2}\frac {\partial \kappa }{\partial s} {\rm d}s + \int _a^b \frac {\kappa }{\rho _2}\frac {\partial \sigma }{\partial s} {\rm d}s \nonumber \\ &\quad {}+\int _a^b \left ( \frac {2}{\rho _1}[[\mu ]]\left (\frac {\partial ^2}{\partial s^2}(\boldsymbol{u}\boldsymbol{\cdot } \hat {\boldsymbol{t}})+\frac {\partial }{\partial s}( \kappa (\boldsymbol{u}\boldsymbol{\cdot } \hat {\boldsymbol{n}}))\right )-\left[ \left[\frac {1}{\rho }\right] \right]\frac {\partial p_2}{\partial s}\right ) {\rm d}s , \end{align}
where the terms of the right-hand side represent mechanisms that drive the net evolution of circulation. These include the viscous diffusion, curvature gradients, surface tension gradients along the interface (including Marangoni effects), and jumps in viscosity and pressure across the interface. We note that Marangoni stresses and variations in the surface tension are distinct in the previous equation; Marangoni stresses arise from the surface tension gradient term
$\partial _s \sigma$
, which generates tangential vorticity, whereas the magnitude of surface tension enters through the curvature term
$\sigma\, \partial _s \kappa$
, and modifies the flow via capillary pressure. The results presented below indicate that the observed enhancement is associated with the former mechanism.
Figure 9 shows the right-hand-side terms of (3.3), evaluated along the arc length
$s$
in an
$(x{-}z)$
plane located at the centreline of the domain, for the two spilling-like regimes considered in this work,
$\epsilon = 0.324$
(figures 9
a–c) and
$\epsilon =0.33$
(figures 9
d–f). For each case, results are shown for the surfactant-free case and surfactant-laden cases with
$\beta_s=0.3$
and
$\beta _s=0.5$
(figures 9(a,d), 9(b,e) and 9(c,f), respectively). This allows identification of the dominant physical mechanisms contributing to circulation generation. The viscous diffusion term is negligible in all conditions, confirming that interfacial mechanisms dominate circulation production. Across all cases, the primary contributor is the term associated with jumps in viscosity and density, which is consistent with the large density ratio in the system. Notably, the magnitude of this term increases in the presence of surfactants, indicating more intense interfacial activity. Contributions from surface tension gradients (Marangoni stresses) become significant only at the largest elasticity number for
$\epsilon =0.33$
, and they are spatially localised near regions of strong gradients in
$\varGamma$
, particularly along the wave crest.
Two-dimensional representation of circulation generation mechanisms along the interface in the
$(x{-}z)$
plane, with (a–c)
$\epsilon = 0.324$
at
$t = 4.8$
, and (d–f)
$\epsilon = 0.33$
at
$t = 5.2$
, for (a,d) the surfactant-free case, (b,e)
$\beta _s = 0.3$
, and (c,f)
$\beta _s = 0.5$
.
Figure 9. Long description
Panel A: A line graph depicts the surfactant-free case. The x-axis is labeled s ranging from 0 to 1. The y-axis on the left is labeled with an integral symbol and ranges from −0.4 to 0.2. The y-axis on the right is labeled I sub sigma and ranges from 0 to 15 times 10 to the power of −4. The blue line shows fluctuations around zero, while the red line shows a peak around 0.6. Panel B: A line graph depicts the case at beta sub s equals 0.3. The x-axis is labeled s ranging from 0 to 1. The y-axis on the left is labeled with an integral symbol and ranges from -0.2 to 0.2. The y-axis on the right is labeled I sub sigma and ranges from 0 to 15 times 10 to the power of −4. The blue line shows fluctuations around zero, while the red line shows a peak around 0.6. Panel C: A line graph depicts the case at beta sub s equals 0.5. The x-axis is labeled s ranging from 0 to 1. The y-axis on the left is labeled with an integral symbol and ranges from −0.2 to 0.2. The y-axis on the right is labeled I sub sigma and ranges from 0 to 15 times 10 to the power of −4. The blue line shows fluctuations around zero, while the red line shows a peak around 0.6. Panel D: A line graph depicts the surfactant-free case. The x-axis is labeled s ranging from 0 to 1. The y-axis is labeled with an integral symbol and ranges from −2 to 1. The blue line shows fluctuations around zero. Panel E: A line graph depicts the case at beta sub s equals 0.3. The x-axis is labeled s ranging from 0 to 1. The y-axis is labeled with an integral symbol and ranges from −2 to 1. The y-axis on the right is labeled I sub sigma and ranges from 0 to 10 times 10 to the power of −3. The blue line shows fluctuations around zero, while the red line shows a peak around 0.6. Panel F: A line graph depicts the case at beta sub s equals 0.5. The x-axis is labeled s ranging from 0 to 1. The y-axis is labeled with an integral symbol and ranges from −2 to 1. The y-axis on the right is labeled I sub sigma and ranges from 0 to 10 times 10 to the power of −3. The blue line shows fluctuations around zero, while the red line shows a peak around 0.6.
4. Conclusions
We present a detailed analysis of the effect of insoluble surfactants on wave breaking dynamics in both regular and spilling regimes, using high-resolution three-dimensional simulations. At
$\textit{We}=100$
and
$\textit{Re}=10^4$
, we explore how wave steepness modulates the flow in the presence of surfactants. The hybrid front-tracking/level-set method employed enables accurate coupling of inertia, capillarity, interfacial diffusion and Marangoni stresses arising from surfactant-induced surface tension gradients. We find that surfactants have minimal impact on regular breakers, which retain similar shapes to the surfactant-free case, though a gradual interfacial rigidification is observed at longer times. In contrast, for spilling breakers, increasing elasticity induces a marked shift in dynamics, which agrees with previously reported experimental work in the same spilling regime. Surfactants enhance spilling breakers through Marangoni stresses, while surface tension reduction alone does not reproduce this effect. This is demonstrated by simulations in which Marangoni stresses are suppressed while allowing local variations of surface tension along the interface. We also show that surfactants generate interfacial vorticity through Marangoni stresses driven by surfactant concentration gradients, as quantified by the circulation along the interface.
To our knowledge, this work provides the first direct characterisation of the spatio-temporal evolution of surfactant concentration profiles during the spilling transition. Our results demonstrate that surfactant redistribution and the resulting Marangoni stresses play a key role in shaping the interfacial dynamics in the parameter regime considered here. How these mechanisms extend to higher Reynolds and Weber numbers, where inertia becomes increasingly important, remains an open question and warrants further investigation.
Future studies should extend the present analysis beyond the regular and spilling regimes into the fully plunging regime, where bubble entrainment and droplet ejection result from the dynamics. Capturing these processes requires that all relevant length and time scales be resolved, making adaptive mesh refinement indispensable. Higher Reynolds numbers and greater wave steepness must also be explored to reproduce realistic ocean conditions. In addition, incorporating soluble surfactants, with full adsorption/desorption kinetics and bulk transport, would introduce chemical time scales that could either dampen or amplify wave breaking events.
Acknowledgements
The authors acknowledge use of the Delta advanced computing and data resource. C.R.C.-A. gratefully acknowledges O.K. Matar for insightful discussions on vorticity production and L. Kahouadji for initial discussions. The numerical simulations were performed using the code BLUE (Shin et al. Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2017).
Funding
This work was supported by the National Science Foundation (award OAC 2005572) and the State of Illinois. D.J. and J.C. acknowledge HPC/AI 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 Équipement National de Calcul Intensif) under grant 2026 A0202B06721. This work used NCSA Delta CPU resources at the National Center for Supercomputing Applications through ACCESS allocations PHY250289 and MCH250041 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Vorticity generation in two-dimensional flows
A.1. Circulation
Before we examine the sources of vorticity flux near an interface, we must first consider the case of a single fluid and the vorticity equation that is given by
\begin{align} \frac {\partial \boldsymbol{\omega }}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\omega }=\boldsymbol{\omega }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}+\nu\, {\nabla} ^2\boldsymbol{\omega }. \end{align}
Now, for a two-dimensional flow,
$\boldsymbol{\omega }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}=0$
since
$\boldsymbol{\nabla }\boldsymbol{u}$
exists only in the
$(x{-}y)$
plane, whereas
$\boldsymbol{\omega }=(0,0,\omega )$
. Thus (A1) can be expressed as a two-dimensional equation for the scalar quantity
$\omega$
:
\begin{align} \frac {{\rm D}\omega }{{\rm D}t}=\nu\, {\nabla} ^2\omega . \end{align}
From Stokes’ theorem, we can write
$\oint _C \boldsymbol{v}\boldsymbol{\cdot }{\rm d}\boldsymbol{s}=\int _A (\boldsymbol{\nabla }\times \boldsymbol{v})\boldsymbol{\cdot }{\rm d}\boldsymbol{a}$
for any vector
$\boldsymbol{v}$
, where
${\rm d}\boldsymbol{s}$
is an infinitesimal displacement along the arc length on a contour
$C$
that encloses a simply connected (note that this becomes important below when we consider the case of two fluids separated by an interface) area
$A$
in the fluid, and
${\rm d}\boldsymbol{a}=\boldsymbol{e}_k\, {\rm d}A$
represents an infinitesimal region in this area;
$\boldsymbol{e}_{k}$
points into the page perpendicularly to the
$(x{-}y)$
plane. Thus if we set
$\boldsymbol{v}=\boldsymbol{u}$
, then
\begin{align} \oint _C\boldsymbol{u}\boldsymbol{\cdot }{\rm d}\boldsymbol{s}=\int _A (\boldsymbol{\nabla }\times \boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{e}_{k}\,{\rm d}A=\int _A \boldsymbol{\omega }\boldsymbol{\cdot }\boldsymbol{e}_{k}\,{\rm d}A=\int _A\omega\, {\rm d}A\equiv \varGamma _A, \end{align}
which relates the total circulation
$\varGamma _A$
to
$\omega$
. By taking a total derivative of
$\varGamma _A$
and making use of (A2), we can write
\begin{align} \frac {{\rm D}\varGamma _A}{{\rm D}t}=\frac {{\rm D}}{{\rm D}t}\int _A\omega\, {\rm d}A=\int _A\frac {{\rm D}\omega }{{\rm D}t}\, {\rm d}A=\int _A \nu\, {\nabla} ^2\omega\, {\rm d}A. \end{align}
Using the divergence theorem, the final term on the left-hand side of this equation can be re-expressed as
\begin{align} \int _A\nu\, {\nabla} ^2\omega\, {\rm d}A=\int _A \boldsymbol{\nabla }\boldsymbol{\cdot }(\nu\, \boldsymbol{\nabla }\omega )\,{\rm d}A=\oint _C\boldsymbol{n}_c\boldsymbol{\cdot }(\nu\, \boldsymbol{\nabla }\omega )\,{\rm d}s, \end{align}
where
$\boldsymbol{n}_c$
is the outward-pointing unit normal to the contour
$C$
. As a result, the total rate of change of circulation in area
$A$
,
$\varGamma _A$
, becomes
\begin{align} \frac {{\rm D}\varGamma _A}{{\rm D}t}=\frac {{\rm D}}{{\rm D}t}\int _A\omega\, {\rm d}A=\oint _C \nu\, \boldsymbol{\nabla }\omega \boldsymbol{\cdot }\boldsymbol{n}_c\, {\rm d}s. \end{align}
This is an important equation that shows that the change in circulation results directly from the flux of vorticity from the boundary (characterised by the contour
$C$
).
A.2. Near-interface vorticity jump
We consider the case of two fluids separated by an interface. In order to evaluate the left-hand side of (A6), we follow an approach that relies on a construction similar to that represented by figure 2 in Brøns et al. (Reference Brøns, Thompson, Leweke and Hourigan2014). This is necessary because the area
$A$
is no longer simply connected but is instead divided because of the presence of the interface. With reference to figure 2 in Brøns et al. (Reference Brøns, Thompson, Leweke and Hourigan2014), we write
\begin{align} \frac {{\rm D}}{{\rm D}t}\int _{A_1 \cup A_2}\omega\, {\rm d}A&=\int _{C_1}\nu _1\boldsymbol{\nabla }\omega _1\boldsymbol{\cdot }\boldsymbol{n}_c\, {\rm d}s + \int _{C_2}\nu _2\boldsymbol{\nabla }\omega _2\boldsymbol{\cdot }\boldsymbol{n}_c\,{\rm d}s +\int _{C_1'}\nu _1\boldsymbol{\nabla }\omega _1\boldsymbol{\cdot }\boldsymbol{n}_c\, {\rm d}s\nonumber\\[0pt] &\quad + \int _{C_2'}\nu _2\boldsymbol{\nabla }\omega _2\boldsymbol{\cdot }\boldsymbol{n}_c\,{\rm d}s. \end{align}
Here, we distinguish between the normal to the contour
$C$
,
$\boldsymbol{n}_c$
, and that to the interface
$I$
,
$\boldsymbol{n}$
.
We now let
$\boldsymbol{n}_c \rightarrow -\boldsymbol{n}$
from the fluid 2 side, and
$\boldsymbol{n}_c \rightarrow \boldsymbol{n}$
from the fluid 1 side, which results in
\begin{align} \frac {{\rm D}}{{\rm D}t}\int _{A_1 \cup A_2}\omega\, {\rm d}A&=\int _{C_1}\nu _1\boldsymbol{\nabla }\omega _1\boldsymbol{\cdot }\boldsymbol{n}_c\, {\rm d}s + \int _{C_2}\nu _2\boldsymbol{\nabla }\omega _2\boldsymbol{\cdot }\boldsymbol{n}_c\,{\rm d}s +\int _{C_1'}\nu _1\boldsymbol{\nabla }\omega _1\boldsymbol{\cdot }\boldsymbol{n}\, {\rm d}s \nonumber\\[0pt] &\quad - \int _{C_2'}\nu _2\boldsymbol{\nabla }\omega _2\boldsymbol{\cdot }\boldsymbol{n}\,{\rm d}s. \end{align}
We note that in region
$A_1$
, we march in an anticlockwise direction along contour
$C_1'$
where
$\boldsymbol{n}_c$
is aligned with
$\boldsymbol{n}$
and then along contour
$C_1$
. In region
$A_2$
, we march along contour
$C_2$
in an anticlockwise direction and then along contour
$C_2'$
, where
$\boldsymbol{n}_c$
points out of fluid 2 into fluid 1 in a direction opposite to that of
$\boldsymbol{n}$
; the latter explains the minus sign in the fourth term of (A8).
We now let the contours
$C_1'$
and
$C_2'$
approach the interface
$I$
so that the union of
$A_1$
and
$A_2$
approaches
$A$
, i.e.
$A_1 \cup A_2 \rightarrow A$
, which leads to
\begin{align} \frac {{\rm D}}{{\rm D}t}\int _A\omega\, {\rm d}A=\oint _C\nu\, \boldsymbol{\nabla }\omega \boldsymbol{\cdot }\boldsymbol{n}_c\, {\rm d}s -\left (\int _I\nu _2\boldsymbol{\nabla }\omega _2\boldsymbol{\cdot }\boldsymbol{n}\,{\rm d}s-\int _I\nu _1\boldsymbol{\nabla }\omega _1\boldsymbol{\cdot }\boldsymbol{n}\,{\rm d}s\right ). \end{align}
By introducing the notation
$[[q]]=q_2-q_1$
for the jump across the interface of a quantity
$q$
, this equation can be re-expressed as
\begin{align} \frac {{\rm D}}{{\rm D}t}\int _A\omega\, {\rm d}A=\oint _C\nu\, \boldsymbol{\nabla }\omega \boldsymbol{\cdot }\boldsymbol{n}_c\, {\rm d}s-\int ^b_a\left [\left [\nu\, \boldsymbol{\nabla }\omega \boldsymbol{\cdot }\boldsymbol{n}\right ]\right ]{\rm d}s. \end{align}
Equation (A10) illustrates the contribution of the jump in the vorticity flux across the interface to the rate of change of the circulation. We can calculate this jump by using
\begin{align} \frac {{\rm D}}{{\rm D}t}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t})+\kappa (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t})-\frac {1}{2}\frac {\partial }{\partial s}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})^2+\frac {\boldsymbol{\nabla }p}{\rho }\boldsymbol{\cdot }\boldsymbol{t}=\nu\, \boldsymbol{\nabla }\omega \boldsymbol{\cdot }\boldsymbol{n}, \end{align}
which comes from (23) in Lundgren & Koumoutsakos (Reference Ma, Hessenkemper, Lucas and Bragg1999) and (2.9) in Brøns et al. (Reference Brøns, Thompson, Leweke and Hourigan2014). Then we have
\begin{align} \frac {{\rm D}}{{\rm D}t}[[\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t}]]+\kappa [[(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t})]] +\frac {\partial }{\partial s}\left[ \left[\frac {p}{\rho }\right] \right]-\frac {1}{2}\frac {\partial }{\partial s}[[(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})^2]]=[[\nu\, \boldsymbol{\nabla }\omega \boldsymbol{\cdot }\boldsymbol{n}]], \end{align}
where we have used the fact that
$\boldsymbol{t}\boldsymbol{\cdot }\boldsymbol{\nabla }=\partial /\partial s$
. Substitution of (A12) into (A10) leads to
\begin{align} \frac {{\rm D}}{{\rm D}t}\left (\int _A \omega \,{\rm d}A + \int _a^b [[\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t}]]\,{\rm d}s\right ) &= \oint _C \nu \,\boldsymbol{\nabla }\omega \boldsymbol{\cdot }\boldsymbol{n}_c\,{\rm d}s - \int _a^b \kappa \, [[(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t})]]\,{\rm d}s \nonumber \\ &\quad {}+ \frac {1}{2}\int _a^b \frac {\partial }{\partial s} [[(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})^{2}]]\,{\rm d}s - \int _a^b \frac {\partial }{\partial s} \left [ \left [\frac {p}{\rho }\right ] \right ] {\rm d}s. \end{align}
Equation (A13) involves a jump in the pressure across the interface. To make progress, we must first consider the jump in the normal stress condition at the interface expressed by
\begin{align} -p_1+\mu _1 \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{D}_1\boldsymbol{\cdot }\boldsymbol{n}+\kappa \sigma =-p_2+\mu _2 \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{D}_2\boldsymbol{\cdot }\boldsymbol{n}, \end{align}
in which
$\sigma$
represents the interfacial tension, and
$\mu _1$
and
$\mu _2$
the dynamic viscosities of fluids 1 and 2, respectively. Rearranging (A14) yields
\begin{align} [[p]]=-\kappa \sigma -2\left[ \left[\mu \left (\frac {\partial }{\partial s}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t})+\kappa (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})\right )\right] \right], \end{align}
where we have made use of
\begin{align} \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{D}\boldsymbol{\cdot }\boldsymbol{n}= -2\frac {\partial }{\partial s}(\boldsymbol{u\boldsymbol{\cdot }\boldsymbol{t}})-2(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})\kappa . \end{align}
Now we re-express
$[[p/\rho ]]$
as
\begin{eqnarray} \left[ \left[\frac {p}{\rho }\right] \right]=\frac {p_2}{\rho _2}-\frac {p_1}{\rho _1}&=&\frac {p_2}{\rho _2}-\frac {p_2}{\rho _1}+\frac {p_2}{\rho _1}-\frac {p_1}{\rho _1}\nonumber \\ &=&p_2\left (\frac {1}{\rho _2}-\frac {1}{\rho _1}\right )\nonumber +\frac {1}{\rho _1}(p_2-p_1)\\ &=&p_2\left[ \left[\frac {1}{\rho }\right] \right]+\frac {[[p]]}{\rho _1}, \end{eqnarray}
from which we get that
\begin{align} \frac {\partial }{\partial s}\left[ \left[\frac {p}{\rho }\right] \right]=-\frac {1}{\rho _1}\frac {\partial }{\partial s}(\kappa \sigma )-\frac {2}{\rho _1}\left[ \left[\mu \left (\frac {\partial ^2}{\partial s^2}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{t})+\frac {\partial }{\partial s}(\kappa [\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}])\right )\right] \right]+\left[ \left[\frac {1}{\rho }\right] \right]\frac {\partial p_2}{\partial s}. \end{align}
Appendix B. Equation of state in insoluble surfactants
In this appendix, we assess the interfacial dynamics at elevated surfactant concentrations, which can lead to artificial numerical results. We emphasise the role of the equation of state relating surface tension
$\sigma$
and surfactant concentration
$\varGamma$
, given by
${\sigma }=\text{max} [0.05, 1 + \beta _s \ln { (1 -{\varGamma })]}$
. While the logarithmic equation of state
$ {\sigma }=1+\beta _s\ln (1-{\varGamma })$
is widely used (Kamat et al. Reference Kamat, Wagoner, Thete and Basaran2020; Manikantan & Squires Reference Melville2020; Constante-Amores et al. Reference Constante-Amores, Kahouadji, Shin, Chergui, Juric, Castrejón-Pita, Matar and Castrejón-Pita2023b
), and is appropriate over a broad range of concentrations, the imposed lower bound
${\sigma }=0.05$
constitutes an abrupt and non-physical modification of the system, a regularisation also adopted in studies such as Muradoglu & Tryggvason (Reference Peck and Sigurdson2008), Asadi et al. (Reference Asadi, Asgharzadeh and Borazjani2018) and Pico et al. (Reference Qiao and Duncan2024). This cut-off becomes active when the concentration exceeds the critical value
${\varGamma }_c=1-\exp (-0.95/\beta _s)$
, beyond which the surface tension is artificially pinned. Although such a regularisation would be inconsequential if
${\varGamma }\lt {\varGamma }_c$
everywhere or only in dynamically passive regions, the simulations discussed in this appendix exhibit local concentrations
${\varGamma }\gt {\varGamma }_c$
, particularly near wave crests. We reiterate that the cut-off is active in these appendix simulations, but is not active in the simulations discussed in the main text. These results are presented to demonstrate that this numerical formulation can yield qualitatively different and non-physical dynamics, underscoring the need for caution in its use and interpretation.
Figures 10 and 11 show the interfacial dynamics in the regular and spilling-like regimes for two values of the surfactant elasticity,
$\beta _s=0.5$
and
$0.9$
, corresponding to critical concentrations
$\varGamma _c \approx 0.85$
and
$\varGamma _c \approx 0.65$
, respectively, for an initial surfactant concentration
$\varGamma _0=0.5$
. In the spilling-like regime, the dynamics exhibits an abrupt qualitative change, with a transition from plunging to spilling-type breaking as
$\beta _s$
is increased. In key snapshots such as figure 11(f), the surfactant concentration in the vicinity of the wave crest clearly exceeds
$\varGamma _c$
for both values of
$\beta _s$
. Since this crest region governs the dominant interfacial dynamics, the imposed surface tension cut-off becomes active precisely where gradients in surface tension, and subsequently Marangoni stresses, are expected to play a central role. As a result, the reported relationship between surface tension reduction, Marangoni forcing, and the transition between plunging and spilling regimes may be influenced by the numerical regularisation inherent in the constitutive model. We emphasise that these limitations are not specific to the present configuration but are inherent to inertia-dominated flows with insoluble surfactants, where strong interfacial compression can drive the concentration into regimes in which commonly used logarithmic equations of state require numerical regularisation. A physically consistent route to modelling such regimes is to extract the equation of state directly from experiments, and to fit a smooth constitutive relation that remains well behaved at high concentrations. This approach has recently been adopted by Eshima et al. (Reference Eshima, Aurégan, Farsoiya, Popinet, Stone and Deike2025), leading to improved physical fidelity in simulations. In the absence of such an experimental equation of state, artificial cut-offs may lead to non-physical effects when high surfactant concentrations are reached. An alternative modelling approach is to allow surfactant solubility, enabling adsorption–desorption dynamics that naturally limits surface accumulation and alleviates the need for ad hoc regularisation.
Effect of the elasticity parameter
$\beta _s$
for the regular wave regime (
$\epsilon =0.3$
) at
$t=5.6$
. Two-dimensional projections of the interface,
$\varGamma$
,
$\tau$
and
$u_{t,x}$
in the
$(x{-}z)$
plane (
$y=\lambda /8$
) are shown in (a–d), respectively. Note that the abscissa corresponds in (a) to the
$x$
coordinate, and in (b–d) to the arc length
$s$
.
Effect of the elasticity parameter for the spilling regime for (a–d)
$\epsilon =0.324$
at
$t=4.8$
, and (e–h)
$\epsilon =0.33$
at
$t=4.4$
. Two-dimensional projections of the interface,
$\varGamma$
,
$\tau$
and
$u_{t,x}$
in the
$(x{-}z)$
plane (
$y=\lambda /8$
) are shown.
Figure 11. Long description
Panel A: A line graph shows the variable z plotted against x. The x-axis ranges from −0.5 to 0.5, and the z-axis ranges from −0.1 to 0.1. The graph includes an inset zooming in on the range from 0 to 0.1 on the x-axis and 0.06 to 0.08 on the z-axis. Panel B: A line graph shows the variable Γ plotted against s. The s-axis ranges from 0 to 1, and the Γ-axis ranges from 0.4 to 1.0. Panel C: A line graph shows the variable τ plotted against s. The s-axis ranges from 0 to 1, and the τ-axis ranges from −0.5 to 1.0. Panel D: A line graph shows the variable ux plotted against s. The s-axis ranges from 0 to 1, and the ux-axis ranges from 0 to 0.4. The graph includes a legend indicating different conditions: Surfactant-free, β = 0.5, and β = 0.9. An inset zooms in on the range from 0.6 to 0.7 on the s-axis and 0.2 to 0.4 on the ux-axis. Panel E: A line graph shows the variable z plotted against x. The x-axis ranges from −0.5 to 0.5, and the z-axis ranges from −0.1 to 0.1. The graph includes an inset zooming in on the range from −0.2 to 0 on the x-axis and 0.04 to 0.08 on the z-axis. Panel F: A line graph shows the variable Γ plotted against s. The s-axis ranges from 0 to 1, and the Γ-axis ranges from 0.2 to 0.8. Panel G: A line graph shows the variable τ plotted against s. The s-axis ranges from 0 to 1, and the τ-axis ranges from −0.5 to 0.5. The graph includes an inset zooming in on the range from 0.4 to 0.5 on the s-axis and −0.5 to 0.5 on the τ-axis. Panel H: A line graph shows the variable ux plotted against s. The s-axis ranges from 0 to 1, and the ux-axis ranges from −0.5 to 0.5. The graph includes a legend indicating different conditions: Surfactant-free, β = 0.5, and β = 0.9. An inset zooms in on the range from 0.45 to 0.50 on the s-axis and −0.5 to 0.5 on the ux-axis.
Appendix C. Mesh resolution and spanwise domain length
To demonstrate mesh-independent results, we have monitored the temporal variation of the kinetic energy of the system for resolutions
$\lambda /\Delta \boldsymbol{x}= 256,512,1024$
for waves with
$\epsilon = 0.30$
(regular) and
$\epsilon = 0.33$
(spilling), displayed in figures 12(a) and 12(b), respectively. For both cases, the coarsest grid, with
$\lambda /\Delta \boldsymbol{x} = 256$
, exhibits noticeable discrepancies, particularly in the valley regions of the kinetic energy evolution. In contrast, the cases with
$\lambda /\Delta \boldsymbol{x} = 512$
and
$1024$
collapse almost exactly onto one another, with discrepancies remaining below
$1\,\%$
over the entire simulation window. The
$\lambda /\Delta \boldsymbol{x} =1024$
resolution was therefore run only for a short time owing to the prohibitive computational cost associated with the much smaller time step. The results for
$\lambda /\Delta \boldsymbol{x} = 512$
and
$1024$
are essentially indistinguishable, demonstrating the mesh independence of our results for regular and spilling breaking cases shown in this work. All results presented in this paper correspond to simulations performed at
$\lambda /\Delta \boldsymbol{x} = 512$
. Finally, figure 13 presents a comparison of key interfacial quantities at
$t=3.0$
, including the interface profile,
$\varGamma$
,
$\tau$
and
$u_{t,x}$
for different grid resolutions with
$\epsilon = 0.33$
and
$\beta _s = 0.5$
. Excellent agreement is observed between the
$\lambda /\Delta x = 512$
and
$1024$
cases across all quantities considered, with relative differences of
$\mathcal{O}(10^{-3})$
. In addition to the instantaneous fields at
$t=3.0$
, we also compare the temporal evolution of the interface shape and interfacial concentration at several times (see figure 14). The interface profiles obtained with
$\lambda /\Delta x=512$
and
$1024$
remain nearly indistinguishable throughout the time interval considered. These results indicate that the interfacial dynamics are well resolved at
$\lambda /\Delta x = 512$
, and the results presented in the paper are mesh-independent.
Grid and spanwise domain-size independence of the kinetic energy evolution. The temporal evolution of the normalised kinetic energy
$E_k/E_{k0}$
is shown for (a)
$\epsilon =0.30$
(regular wave), and (b)
$\epsilon = 0.33$
(spilling breaker). In (a), spanwise domain widths
$\lambda$
,
$\lambda /2$
,
$\lambda /4$
and
$\lambda /8$
are also added at fixed spatial resolution
$\lambda /\Delta \boldsymbol{x}=512$
.
Grid convergence for a representative surfactant-laden spilling breaker (
$\epsilon = 0.33$
,
$\beta _s = 0.5$
) at time
$t=3.0$
. Comparisons of
$\lambda /\Delta x=256, 512, 1024$
for (a) interface location, (b)
$\varGamma$
, (c)
$\tau$
and (d)
$u_{t,x}$
.
Temporal grid convergence for the spilling breaker case
$\epsilon =0.33$
,
$\beta _s=0.5$
up to
$t=9.0$
: (a)
$z(x,t)$
and (b)
$\varGamma (s,t)$
for
$\lambda /\Delta x=512$
and
$1024$
.
Figure 14. Long description
Panel A: A line graph depicts the relationship between x and 4z plus t. The horizontal axis is labeled x with a range from −0.5 to 0.5. The vertical axis is labeled 4z plus t with a range from 0 to 8. Two data series are shown: one in black and one in red, representing different grid sizes of 512 and 1024 respectively. Panel B: Another line graph shows the relationship between s and 2z plus t. The horizontal axis is labeled s with a range from 0 to 1. The vertical axis is labeled 2z plus t with a range from 0 to 8. Similar to Panel A, two data series are displayed in black and red, corresponding to grid sizes of 512 and 1024. A legend in Panel B indicates the color coding for the grid sizes.
To verify that the chosen computational width (
$\lambda /4$
) does not influence the dynamics, we performed a width-independence study for the clean-wave case with
$\epsilon = 0.3$
, measuring the kinetic energy evolution while repeating the simulation using domain widths
$\lambda$
,
$\lambda /2$
,
$\lambda /4$
and
$\lambda /8$
, as shown in figure 12(a). The resulting kinetic energy plots collapse onto one another, demonstrating that spanwise confinement does not affect the wave dynamics in the regime studied.
Appendix D. Mechanical energy evolution
In this work, the instantaneous kinetic and gravitational potential energies are defined as
\begin{align} E_K(t) = \frac {1}{2}\int _{d\geqslant 0} \rho \left (u^2+v^2\right )\mathrm{d}x\,\mathrm{d}y, \quad E_P(t) = \int _{d\geqslant 0} \rho y \,\mathrm{d}x\,\mathrm{d}y, \end{align}
and the mechanical energy is given by
\begin{align} E_{\textit{mech}}(t)=E_K(t)+E_P(t). \end{align}
Following the approach of Iafrati (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009), we note that for weakly damped regular gravity waves,
$E_K$
and
$E_P$
exchange periodically over a wave period, while their sum
$E_{\textit{mech}}$
decays monotonically due to viscous dissipation. In this regime, the decay of the mechanical energy is given by
\begin{align} E_{\textit{mech}}(t)\sim E_{\textit{mech}}(0)\, \mathrm{e}^{-2\gamma t}, \end{align}
and the exponential envelope of
$E_K(t)$
provides an accurate estimate of the decay rate.
Figure 15 shows the computed potential energy
$E_P(t)$
evaluated at discrete times, reflecting the fact that snapshots were saved only at selected instants. For comparison, we also include the predicted mechanical energy
$E_{\textit{mech}}(t)$
. We also find that the decay of
$E_{\textit{mech}}$
is in excellent agreement with the exponential decay of the kinetic energy. Furthermore, the contribution of surface tension energy is negligible in this regime, as the Weber number is large; the same assumption was adopted by Iafrati (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2009).
Temporal evolution of the mechanical, kinetic and gravitational potential energies, normalised by their initial values for
$\epsilon =0.3$
.
Appendix E. Comparison of bulge shapes with experimental results
Figure 16 presents a qualitative comparison between the present simulations and the experimental profiles reported by Liu & Duncan (Reference Liu and Duncan2006). We emphasise that this comparison is not intended as a quantitative validation. The experiments involve soluble surfactants with adsorption–desorption kinetics, whereas the present model assumes an insoluble surfactant with no exchange with the bulk. In addition, the characteristic Reynolds and Bond numbers in the experiments are substantially larger than those considered here. Using the crest height
$H_c$
as the length scale, Liu & Duncan (Reference Liu and Duncan2006) reported
$\textit{Re} \approx 9.3\times 10^4$
,
$Bo \approx 9.2\times 10^2$
, and wave steepness
$\epsilon \approx 0.44$
. When expressed in terms of the wavelength
$\lambda$
, these correspond to
$\textit{Re} \approx 1.3\times 10^6$
and
$Bo \approx 1.9\times 10^5$
. By contrast, the present simulations are performed at
$\textit{Re} = 10^4$
and
$Bo = 253.3$
(
$\textit{We} = 100$
), indicating differences of 1–2 orders of magnitude, depending on the choice of length scale. Despite these differences, the overall bulge morphology, including crest elevation, forward leaning of the crest, and the curvature transition near the toe, shows similar qualitative trends across the datasets. This suggests that these large-scale interfacial features are reproduced consistently, even though the surfactant model and operating conditions differ substantially. The comparison should therefore be interpreted as evidence of qualitative consistency in interface shape, rather than as a strict quantitative agreement or a validation of the detailed dynamics.
Comparison of bulge shapes for clean and surfactant-laden conditions between (a) the experiments of Liu & Duncan (Reference Liu and Duncan2006) and (b) the present simulations at
$\epsilon = 0.35$
.
x=(x,y,z)
ϵ
βs
ϵ
βs
ϵ
(x−z)
y=λ/8
ϵ
βs
z
Γ
ϵ
βs
s∈[0,1]
Γ
Pes=103
Γ0=Γ∞/4
βs
ϵ=0.3
t=5.6
Γ
τ
ut,x
(x−z)
y=λ/8
x
s
ϵ=0.324
t=4.8
ϵ=0.33
t=5.2
Γ
τ
ut,x
(x−z)
y=λ/8
ϵ=0.324
ϵ=0.33
ϵ=0.33
βs=0.5
∇sσ~=0
βs=0.5
t=5.2
t=5.6
t=6.0
ct
Γ
ωy
(x−z)
ϵ=0.324
t=4.8
ϵ=0.33
t=5.2
ϵ=0.33
t=5.2
∇sσ~=0
βs=0.3
βs=0.5
(x−z)
ϵ=0.324
t=4.8
ϵ=0.33
t=5.2
βs=0.3
βs=0.5
βs
ϵ=0.3
t=5.6
Γ
τ
ut,x
(x−z)
y=λ/8
x
s
ϵ=0.324
t=4.8
ϵ=0.33
t=4.4
Γ
τ
ut,x
(x−z)
y=λ/8
Ek/Ek0
ϵ=0.30
ϵ=0.33
λ
λ/2
λ/4
λ/8
λ/Δx=512
ϵ=0.33
βs=0.5
t=3.0
λ/Δx=256,512,1024
Γ
τ
ut,x
ϵ=0.33
βs=0.5
t=9.0
z(x,t)
Γ(s,t)
λ/Δx=512
1024
ϵ=0.3
ϵ=0.35
Open access
