2.1. Governing equations

We solve the two-dimensional, incompressible, two-phase Navier–Stokes equations with interfacial effects, including Marangoni stresses from surface tension gradients. The interface between water and air is reconstructed by a coupled level-set/VOF method. The equations can be written as

(2.1) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}&=0, \end{align}

(2.2) \begin{align} \rho\! \left (\frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\right )&=-\boldsymbol{\nabla }\! p+\boldsymbol{\nabla }\boldsymbol{\cdot }(\mu (\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{\nabla }\boldsymbol{u}^{\rm T}))+\delta _s\sigma \kappa \boldsymbol{n}+\delta _s\,\boldsymbol{\nabla} _{\!s}\sigma , \end{align}

where $\boldsymbol{u}$ is the fluid velocity, $\rho$ is the fluid density, $p$ is the pressure, $\mu$ is the dynamic viscosity, $\sigma$ is the surface tension, and $\kappa$ is the interface curvature. The unit normal vector to the interface is denoted by $\boldsymbol{n}$ . The function $\delta _s$ is a Dirac delta function that localises the surface tension term at the interface, and $\boldsymbol{\nabla} _{\!s}=(\boldsymbol {I}-\boldsymbol {nn})\boldsymbol{\cdot }\boldsymbol{\nabla}$ is the surface gradient operator (Stone Reference Stone1990; Farsoiya et al. Reference Farsoiya, Popinet, Stone and Deike2024). The last term represents the tangential Marangoni stress due to gradients in $\sigma$ along the interface.

The Marangoni force is computed using an integral formulation for surface tension (Abu-Al-Saud, Popinet & Tchelepi Reference Abu-Al-Saud, Popinet and Tchelepi2018), and surfactant transport is modelled using a coupled phase-field/VOF method developed by Jain (Reference Jain2024) and implemented in Farsoiya et al. (Reference Farsoiya, Popinet, Stone and Deike2024). The evolution of the insoluble surfactant concentration is given by

(2.3) \begin{equation} \frac {\partial c}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u}c) =\boldsymbol{\nabla }\boldsymbol{\cdot }\left [D_s\left \{\boldsymbol{\nabla }c-\frac {2(0.5-\phi )\boldsymbol{n}c}{\varepsilon }\right \}\right ]\!, \end{equation}

where $c=\varGamma \phi (1-\phi )/\varepsilon$ is the volumetric surfactant concentration, and $\varGamma$ is the local surface surfactant concentration used to calculate $\sigma$ through the EOS. The second term on the right-hand side is an anti-diffusion term that keeps the scalar $c$ concentrated at the interface, which ensures that surfactant remains restricted to the interface, thereby numerically enforcing the insolubility condition. Here, $\varepsilon$ is the phase field interface thickness, which is set to $0.75\,\Delta x_{\textit{min}}$ , where $\Delta x_{\textit{min}}$ is the minimum grid spacing, following Farsoiya et al. (Reference Farsoiya, Popinet, Stone and Deike2024). Sensitivity tests on the values of $\varepsilon$ have been performed, and the chosen value $\varepsilon =0.75\,\Delta x_{min}$ provides a good balance between accuracy and numerical stability. More details are given in Appendix A. The scalar field $\phi$ defines the interface location and is governed by a conservative phase field equation

(2.4) \begin{equation} \frac {\partial \phi }{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u}\phi )=\boldsymbol{\nabla }\boldsymbol{\cdot }\left \{\zeta \left \{\varepsilon\, \boldsymbol{\nabla }\phi -\frac {1}{4}{\left [1-\tanh ^{2}\left (\frac {\psi }{2\varepsilon }\right )\right ]}\frac {\boldsymbol{\nabla }\psi }{\left |\boldsymbol{\nabla }\psi \right |}\right \}\right \}\!, \end{equation}

where $\zeta$ is a velocity scale parameter set to $\zeta = 1.1\, |\boldsymbol{u}|_{\textit{max}}$ , and $\psi = \varepsilon \log [ (\phi + \delta )/(1 - \phi + \delta ) ]$ is an auxiliary signed-distance-like function with a small regularisation parameter $\delta = 10^{-6}$ , following Farsoiya et al. (Reference Farsoiya, Popinet, Stone and Deike2024). Sensitivity tests to the values of $\zeta$ are performed in Appendix A, which show negligible effects on the results. Therefore, the surfactant is introduced by coupling (2.3) and (2.4) to solve for $\varGamma$ and then contribute to $\sigma$ by the EOS, which specific formulation is given in the next subsection. By introducing a phase field model, we avoid the explicit tracking of surfactant transport on a sharp interface with large deformations, which improves the solver’s numerical stability. The trade-off comes at the cost of an increased computational complexity and the introduction of additional numerical parameters for the phase field equation, which we test in Appendix A. Note that as verified in Farsoiya et al. (Reference Farsoiya, Popinet, Stone and Deike2024), the present numerical scheme has good surfactant mass conservation properties. This was also verified in the present configuration, with errors below 4 %.

This coupled framework is implemented in the open-source solver Basilisk (https://basilisk.fr) (Popinet Reference Popinet2009, Reference Popinet2018), which supports adaptive mesh refinement and parallel scalability, and has been extensively applied to multiphase flow problems. The spatial discretisation uses a quadtree (in two dimensions), where each computational cell can be subdivided into four child cells. The base of the tree, known as the root cell, is assigned level zero, and each refinement step increases the level by one. The maximum refinement level, denoted by $L$ , determines the finest resolution available in the simulation. Accurate simulation of surface waves requires high spatial resolution near the interface and in the boundary layers, where vorticity develops and energy dissipation mainly occurs. The wavelet adaptive refinement criteria are based on the volume fraction field, velocity field, phase field and concentration field, with corresponding thresholds $\epsilon _{\!f}=10^{-3}$ , $\epsilon _u=5\times 10^{-3}$ , $\epsilon _\phi =10^{-3}$ , $\epsilon _c=10^{-3}$ , ensuring that the key physical processes are well resolved.

2.2. The EOS

The EOS is typically nonlinear with surfactant concentration, and different types of surfactants can lead to different EOS (Liao et al. Reference Liao, Franses and Basaran2006; Muradoglu & Tryggvason Reference Muradoglu and Tryggvason2008). For surface waves, with the assumption of insolubility, a relationship between surface tension and surfactant concentration was experimentally measured in the case of insoluble surfactant Triton X using a Langmuir trough (Erinin et al. Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023). The measured relationship presents some key features: it shows that surface tension decreases nonlinearly with increasing surfactant concentration, and approaches a saturation limit corresponding to the maximum surface packing density. Therefore, we consider a smooth EOS with a two-parameter tanh model (Liao et al. Reference Liao, Franses and Basaran2006; Farsoiya et al. Reference Farsoiya, Popinet, Stone and Deike2024; Eshima et al. Reference Eshima, Aurégan, Farsoiya, Popinet, Stone and Deike2025):

(2.5) \begin{equation} \frac {\sigma }{\sigma _c}=1-\Delta \sigma _\infty \tanh\! \left (\frac {\beta }{\Delta \sigma _\infty }\frac {\varGamma }{\varGamma _0}\right )\!, \end{equation}

where $\sigma _c$ is the surface tension of a clean (surfactant-free) interface, $\varGamma$ is the local surfactant concentration, and $\varGamma _0$ is the initial (uniform) surfactant concentration. This EOS originates from Henry’s isotherm, which describes a linear dependence of surface tension on surfactant concentration (Manikantan & Squires Reference Manikantan and Squires2020). The use of a hyperbolic tangent function provides a smooth transition of surface tension at different surfactant concentrations. This form recovers the linear relation at low concentrations while smoothly approaching a finite minimum value at high concentrations, so as to avoid artificial cut-offs or discontinuities that could cause numerical instability. This EOS also captures the basic nonlinear trend observed in experiments (Erinin et al. Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023) while remaining mathematically simple and computationally efficient. Although it is not a one-to-one comparison, and the specific choice of parameters ( $\beta , \varDelta$ ) modifies the magnitude of the Marangoni number, the underlying physics governed by surface tension gradients should remain consistent.

We would also like to point out that there is no obvious good choice in terms of EOS for insoluble surfactant in terms of experimental measurements in the context of surfactant effects on surface waves. Most experimental studies do not report the EOS, or do not measure it. Part of the puzzle is to find out how much the specific surfactant, or mixture of surfactant, will affect the physics of surface waves. Detailed comparison between experiments and simulations would require knowledge of the measured EOS. These questions have motivated our choice of a simple EOS with a small number of parameters. We also used this EOS in a separate paper on jet drops formed by bubble bursting, demonstrating close agreement to experiments where we had measured the Langmuir isotherm as a way to access the EOS (Eshima et al. Reference Eshima, Aurégan, Farsoiya, Popinet, Stone and Deike2025).

Two additional dimensionless parameters are therefore introduced as we consider insoluble surfactant, $\beta$ and $\Delta \sigma _\infty$ , which is physically related to the surfactant concentration and strength, and controls the slope of $\sigma (\beta )$ and the maximum reduction in surface tension, respectively. This formulation ensures that $\sigma \to \sigma _c(1 - \Delta \sigma _\infty )$ as $\varGamma \gg \varGamma _0$ , and recovers the linear behaviour for $\varGamma \ll \varGamma _0$ , i.e. $\sigma \approx \sigma _c(1 - \beta \varGamma / \varGamma _0)$ (Farsoiya et al. Reference Farsoiya, Popinet, Stone and Deike2024). The parameter $\beta$ is tightly linked to the Marangoni number. During the surfactant redistribution process, the total interfacial surfactant mass is conserved. Adsorption and desorption processes are neglected in the present model.

Figure 1 presents $\sigma /\sigma _c$ as a function of $\varGamma /\varGamma _0$ from (2.5) for different $\beta$ with fixed $\Delta \sigma _\infty = 0.5$ . This shows how the sensitivity of surface tension to surface compression evolves with increasing $\beta$ . For small values of $\beta$ , as $\varGamma /\varGamma _0$ increases from $1$ , the normalised surface tension $\sigma /\sigma _c$ initially decreases slowly, followed by an approximately linear decrease with $\varGamma /\varGamma _0$ . For larger $\beta$ , $\sigma /\sigma _c$ decreases continuously as $\varGamma /\varGamma _0$ increases, with the slope of $\sigma (\varGamma )$ also decreasing. Parameter $\beta$ influences not only the slope but also the initial surface tension at uniform concentration ( $\varGamma /\varGamma _0 = 1$ ), which decreases with increasing $\beta$ . We note again that the fact that $\beta$ controls ${\rm d}\sigma /{\rm d}l$ and the initial surface tension is a feature common to any nonlinear EOS.

Figure 1. Illustration of the nonlinear EOS used in this study: Normalised surface tension $\sigma /\sigma _c$ as a function of $\varGamma /\varGamma _0$ for increasing values of the surfactant parameter $\beta$ (representing surfactant strength and concentration) with $\Delta \sigma _\infty =0.5$ .

In the following, the value of the surfactant parameter $\Delta \sigma _\infty$ is fixed at $0.5$ , which gives a $\sigma (\varGamma )$ EOS similar to those measured experimentally by Erinin et al. (Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023) for Triton-X using a Langmuir trough (and assuming insoluble surfactant so that the area of compression can be related to the surfactant concentration; see also Eshima et al. Reference Eshima, Aurégan, Farsoiya, Popinet, Stone and Deike2025).

2.3. Control parameters and initial conditions

The numerical configuration used here follows the one in Deike et al. (Reference Deike, Popinet and Melville2015) used to study capillary effects on breaking waves. As discussed in the Introduction, the problem can be defined by a set of dimensionless parameters: the initial wave slope $ak$ (going from non-breaking to breaking waves), the wave Reynolds and Bond numbers, and the air–water density and viscosity ratios ( $\rho _a / \rho _w = 1/850$ and $\mu _a / \mu _w = 1.96 \times 10^{-2}$ ). We will investigate the influence of $\textit{Re}$ , $\textit{Bo}$ and $ak$ on the wave patterns and dynamics for various surfactant contamination measured by the parameter $\beta$ .

Following Deike et al. (Reference Deike, Popinet and Melville2015), the wave is initialised using a nonlinear potential flow solution $\varPhi _0$ with corresponding free-surface elevation $\eta _0$ , based on a third-order Stokes wave approximation (Lamb Reference Lamb1924). This set-up provides both the initial wave profile and the associated water velocity field, so that all cases are initialised consistently, ensuring that the comparative trends across different $\textit{Bo}$ and $\beta$ remain robust. We note the work from Shelton & Trinh (Reference Shelton and Trinh2022) that included the effect of surface tension to compute Stokes waves with surface tension at high steepness. The wave propagates in a periodic domain with free-slip boundary conditions at the top and bottom walls. The wavelength is equal to the domain length, and the water depth is set to $h/\lambda = 1/2$ , ensuring that depth effects on wave evolution are negligible, consistent with Deike et al. (Reference Deike, Popinet and Melville2015).

In the presence of surfactants, two main non-dimensional parameters $\beta$ and $\Delta \sigma _\infty$ are introduced from the EOS (2.5). We explore a broad parameter space with $0.005 \leqslant \beta \leqslant 1$ for $10^1 \leqslant \textit{Bo} \leqslant 10^3$ , $ak = 0.2, 0.25, 0.3, 0.35$ , and $\textit{Re} = 10^4, 4 \times 10^4, 10^5$ , with $\Delta \sigma _\infty$ fixed as $0.5$ . Since the surfactant transport in our study is dominated by advection, $\textit{Pe}$ should be sufficiently high so that further increases have a negligible influence on the results; $\textit{Pe}$ is fixed at $40$ based on the convergence test for increasing $\textit{Pe}$ shown in Appendix A.

All simulations run up to $t = 2T_0$ , where $T_0 = 2\pi / \sqrt {gk}$ is the linear wave period derived from the gravity wave linear dispersion relation. The wave period based on the full gravity-capillary wave dispersion relation varies slightly with $\textit{Bo}$ , but the variation is minimal for the range of $\textit{Bo}$ considered here. Variations due to nonlinear effects are also small (see Deike et al. Reference Deike, Popinet and Melville2015).

Over this time, the formation of parasitic capillaries or spilling breaking waves can be observed and characterised, together with the associated energy dissipation.

To ensure numerical accuracy, grid convergence tests are performed for various combinations of $\textit{Bo}$ and $\beta$ , as shown in Appendix A. We examine both the instantaneous wave shape (after one wave period) and the time evolution of the total energy budget at grid resolutions $L = 8$ , $9$ and $10$ , where $L$ is the maximum level of adaptive refinement. The results show convergence in both interface dynamics and energy evolution. Based on this analysis, we adopt $L = 9$ (corresponding to $2^9 = 512$ grid cells per wavelength) as the standard resolution for all cases, except for $\textit{Re} = 10^5$ , where $L = 10$ is used to resolve the thinner boundary layers associated with higher Reynolds numbers.

We focus in this paper on two-dimensional simulations of surfactant effects on surface waves in order to span a wide parameter range, allowing us to rigorously investigate the effect of the various non-dimensional parameters (wave parameters $\textit{Bo}$ , $\textit{Re}$ and $ak$ , as well as surfactant parameter $\beta$ ), together with grid resolution verification. These simulations will permit us to discuss the transition in dynamics from gravity waves to parasitic capillary waves, as well as the transition from breaking to non-breaking. Energy dissipation in these various regimes is then quantified. We performed in total 400 simulations for a total cost of 40 000 CPU hours, summarised in table 1. Three-dimensional simulations would be a natural next step, allowing us to explore transverse Marangoni effects, bubble entrainment and break-up, and drop formation, to name a few, in targeted regimes of parameters, informed by the present study. The code for this work is available online (https://basilisk.fr/sandbox/ryang/surfactant/).

Table 1. Ranges of control parameters used in the simulations.

3.1. Non-monotonic suppression of parasitic capillary waves

We first examine the evolution of the wave interface at low Bond number $\textit{Bo} = 10$ and slope below the breaking threshold $ak = 0.3$ , with different values of the surfactant parameter $\beta$ . Figure 2 shows the time evolution of the interface, the surfactant concentration at the interface, and the local Marangoni number (defined in (1.2)).

Figure 2. Effect of surfactant on gravity-capillary waves at low Bond number for $ak=0.3,\ \textit{Re}=4\times 10^4,\ \textit{Bo}=10$ and (a) $\beta =0.005$ , (b) $\beta =0.3$ and (c) $\beta =1$ . The top images show the instantaneous vorticity field at $t/T=1$ . The time evolution (from bottom to top, starting at $t/T=0$ and then at intervals $t/T=0.16$ ) of the interface is shown ,together with a colour map of the upper surface representing the distribution of normalised surfactant concentration, and the colour map of the lower surface represents the local Marangoni number $\textit{Ma}_{\textit{loc}}$ .

For a small value $\beta = 0.005$ , the influence of surfactants is weak, and the initial wave profile corresponds to a clean gravity-capillary wave. As the wave progresses, parasitic capillary waves gradually emerge on the forward face of the main wave, consistent with prior observations (Deike et al. Reference Deike, Popinet and Melville2015) (figure 2 a). Positive vorticity is located at the crests of both the main wave and the parasitic waves, while negative vorticity appears near the troughs. These parasitic capillary waves result from the localised high curvature at the wave crest, where surface tension effects are dominant, and are finally damped by viscosity (Fedorov & Melville Reference Fedorov and Melville1998). These parasitic capillaries are also found to significantly enhance dissipation (Fedorov & Melville Reference Fedorov and Melville1998; Melville & Fedorov Reference Melville and Fedorov2015).

The colour fields along the interface in the figure represent the normalised surfactant concentration and the normalised local surface tension gradient to visualise the influence of Marangoni forcing. The surfactant concentration illustrates the redistribution of surfactants from an initially uniform state. As time evolves, surfactants accumulate at the front face of the wave crest where parasitic waves originate. From the surfactant concentration and the EOS, the surface tension gradient can be calculated and indicates both the direction and strength of the Marangoni force, which is weak due to the low value of the parameter $\beta$ (effectively a very low surfactant concentration).

As $\beta$ increases (figure 2(b), $\beta =0.3$ ), surfactant effects become more pronounced. The formation of parasitic waves is nearly fully suppressed. This damping effect is consistent with previous experimental findings that show suppression of parasitic ripples with surfactant (Xu & Perlin Reference Xu and Perlin2023). The surfactant continues to accumulate at the front face of the crest, where it decreases the local surface tension. The non-uniform distribution of surfactant also induces stronger surface tension gradients compared to lower $\beta$ , which is shown in the colour field of surface tension gradients. It acts opposite to the direction of surfactant accumulation, and tends to stretch and smooth the wave crest. Correspondingly, near-surface vorticity intensifies compared to the $\beta = 0.005$ case, due to enhanced Marangoni flow, with the magnitude of the local Marangoni number being much higher.

For even larger values of $\beta$ (figure 2(c), $\beta =1$ ), parasitic capillary waves recur, but with reduced amplitude compared to the low- $\beta$ case ( $\beta = 0.005$ ). The surfactant again accumulates at both the main wave crest and the crests of the parasitic waves. However, the surface tension gradient is weaker than in the intermediate case ( $\beta =0.3$ ), indicating a weaker Marangoni effect. Similarly, the vorticity field is stronger than in the low- $\beta$ case, but weaker than in the intermediate- $\beta$ case, where parasitic suppression is most pronounced.

Because the surfactant concentration reaches local maxima at the crests of both the main wave and the parasitic ripples, and gradually decays away from these regions, the resulting surface tension distribution varies accordingly. The associated surface tension gradients drive opposing Marangoni flows on either side of each ripple, giving rise to the observed spatial variations in $\textit{Ma}_{\textit{loc}}$ . This variation appears as alternating blue–red spikes at the parasitic ripples in figure 2(a–c).

This non-monotonic behaviour in the surfactant parameter $\beta$ suggests a complex interaction between surfactant transport, Marangoni forcing and wave steepening, which we now analyse. Figure 3 shows the evolution of the wave interface at a larger Bond number, $\textit{Bo} = 200$ , and $ak = 0.3$ for varying values of $\beta$ . In this high- $\textit{Bo}$ regime, gravity dominates the wave dynamics. Similar to the $\textit{Bo} = 10$ cases, surfactant accumulates near the front of the wave crest, and the colour fields reveal a non-monotonic dependence of the surface tension gradient magnitude on $\beta$ . The wave surface at $\beta = 1$ remains smooth during time evolution, characteristic of a pure gravity wave, while for lower values of $\beta$ , small ripples can still be observed near the crest. The addition of surfactant introduces localised surface tension gradients that generate Marangoni shear along the crest. This shear can locally destabilise the interface and enhance wave steepening, particularly at moderately high $\textit{Bo}$ , where Marangoni stresses remain dynamically relevant despite the dominance of gravity. Nevertheless, the overall wave evolution is only weakly influenced by the presence of surfactant at high $\textit{Bo}$ . It was expected that the Marangoni effect is less pronounced in this high- $\textit{Bo}$ regime, which is highlighted by the relatively low local Marangoni numbers observed.

Figure 3. The instantaneous vorticity field and the time evolution of the interface for $ak=0.3,\ \textit{Re}=4\times 10^4,\ \textit{Bo}=200$ and (a) $\beta =0.005$ , (b) $\beta =0.3$ and (c) $\beta =1$ . The top images show the instantaneous vorticity field at $t/T=1.5$ . The time evolution (from bottom to top, starting at $t/T=0$ and then at intervals $t/T=0.16$ ) of the interface is shown together with a colour map of the upper surface representing the distribution of normalised surfactant concentration, and the colour map of the lower surface represents the local Marangoni number $\textit{Ma}_{\textit{loc}}$ .

3.2. Non-monotonic suppression on spilling breakers at higher $ak$

We now consider the effect of $\beta$ for larger $ak = 0.35$ , which can lead to wave breaking (Deike et al. Reference Deike, Popinet and Melville2015). The different wave breaking regimes are identified using slope-based criteria following Deike et al. (Reference Deike, Popinet and Melville2015): a spilling breaker is defined by the appearance of a vertical segment in the interface profile (i.e. a $-90^\circ$ slope), while a plunging breaker characterised by wave overturning is identified by the presence of both vertical and horizontal segments (i.e. a $-180^\circ$ configuration), as illustrated in figure 4.

Figure 5(a–c) show the time evolution of the wave interface for $ak = 0.35$ and $\textit{Bo} = 40$ at different values of $\beta$ . These snapshots illustrate a transition from spilling breakers to parasitic capillary waves, and then back to spilling breakers as $\beta$ increases. This transition can be seen from the time evolution of interface shapes shown in figure 5(a–c) and the slopes indicated by dashed lines at approximately one wave period (the 7th wave profile): at low and high $\beta$ , the leading wave crest is taller and sharper, while at intermediate $\beta$ , it becomes lower and broader, indicating the suppression of wave breaking. In figure 5(b,c), we see a transition from parasitic capillary wave to spilling breaker (vertical interface exceeds $90^\circ$ , with the spilling associated with more mixing and higher wave slope). Higher $\textit{Bo}$ will lead to plunging breakers following the criteria from figure 4(d). Note that at a later stage, the single wave crest splits into more crests, which is a representative spilling breaker as the bulge near the crest releases and transitions into a train of surface ripples (Liu & Duncan Reference Liu and Duncan2006). Surfactant accumulation on the front face could also enhance this deformation by reducing the local surface tension.

Figure 4. Illustration of typical wave shapes and breaking criteria: (a) gravity wave; (b) parasitic capillary wave; (c) spilling breaker with a vertical segment ( $-90^\circ$ ); (d) plunging breaker with both vertical and horizontal segments ( $-180^\circ$ ).

Figure 5. The instantaneous vorticity field and time evolution of the interface for $ak=0.35,\ \textit{Bo}=40 $ and (a) $\beta =0.005$ , (b) $\beta =0.3$ , (c) $\beta =1$ . Here, (a) and (c) are spilling breakers, and (b) is parasitic capillary waves following the crest shape criteria from Deike et al. (Reference Deike, Popinet and Melville2015). The top images show the instantaneous vorticity field at $t/T=1$ . The time evolution (from bottom to top, starting at $t/T=0$ and then at intervals $t/T=0.16$ ) of the interface is shown together with a colour map of the upper surface representing the distribution of normalised surfactant concentration, and the colour map of the lower surface represents local Marangoni number $\textit{Ma}_{\textit{loc}}$ .

The corresponding surfactant concentration distributions and local surface tension gradients for these cases are also shown in the colour field. These results are qualitatively similar to those observed for $ak = 0.3$ , except that the surfactant concentration is more focused compared to the non-breaking case due to stronger surface compression.

3.3. Mean square slope evolution and local Marangoni number

We have observed a modification of the dynamics at intermediate values of $\beta$ , leading to the suppression of parasitic capillary waves or breaking event due to strong Marangoni stresses, with such effects modulated by the relative strength of gravity and surface tension (Bond number) and wave slope.

To quantify the influence of surfactants on wave dynamics, we examine the time evolution of the mean square slope ( $s^2$ ) of the wave interface. Figure 6(a) shows the evolution of $s^2$ for various values of $\beta$ at low $\textit{Bo}$ ( $\textit{Bo}=10$ ). At low $\beta$ , $s^2$ initially increases during the first wave period due to the emergence of parasitic capillary waves, then gradually decays as the viscosity dissipates energy. In contrast, for intermediate values of $\beta$ , $s^2$ decreases nearly monotonically over time, reflecting the suppression of parasitic capillaries under strong Marangoni forcing. At even higher values of $\beta$ , $s^2(t)$ again exhibits an initial increase followed by decay – indicating the reappearance of parasitic capillaries. For the high- $\textit{Bo}$ case ( $\textit{Bo}=200$ ), as shown in figure 6(b), the effect of $\beta$ is negligible, consistent with the interface evolution in figure 3.

Figure 6. (a,b) The evolution $s^2$ for $ak=0.3$ and (a) $\textit{Bo}=10$ , (b) $\textit{Bo}=200$ , with different values of $\beta$ colour-coded. The points represent the global maxima. (c,d) The evolution of the spatially averaged local Marangoni number $\langle |\textit{Ma}_{\textit{loc}}|\rangle$ for $ak=0.3$ and (c) $\textit{Bo}=10$ , (d) $\textit{Bo}=200$ , with different $\beta$ .

Furthermore, we examine the time evolution of the spatially averaged local Marangoni number $\langle |\textit{Ma}_{\textit{loc}}|\rangle$ , defined in (1.2), at $\textit{Bo}=10$ and $200$ , as shown in figure 6(c,d). The results show a non-monotonic but opposite evolution trend with $\beta$ compared to the trend of $s^2$ . As $\beta$ increases, $\langle |\textit{Ma}_{\textit{loc}}|\rangle$ first increases and then decreases, consistent with our observation in figures 2 and 3. The trends are similar for different $\textit{Bo}$ but with different magnitudes – a much lower magnitude at higher $\textit{Bo}$ .

Then we compute the time-averaged $\overline {s^2}$ and time-averaged local Marangoni number $\overline {\langle \textit{Ma}_{{loc}}\rangle }$ , where the bar represents time averaging over the interval $t/T_0 \in [0.8, 1.2]$ , corresponding to the period when parasitic capillary waves or breaking waves begin to develop, which is the primary focus of our analysis, while significant wave damping has already occurred after that. Figure 7(a) present $\overline {s^2}$ versus $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ for all different $\textit{Bo}$ at $ak=0.3$ , with the circles representing each simulation case. They reveal the dependence of surface deformation on Marangoni forcing, with a strong modulation by the Bond number, and a clear collapse at high $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ , which shows that the surface roughness quantified by the mean square slope is dominated by Marangoni forcing. Figure 7(b,c) show both $\overline {s^2}$ and $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ as functions of $\beta$ with different $\textit{Bo}$ for $ak=0.3$ . They show non-monotonic trends of $\overline {s^2}$ and $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ with $\beta$ , with again a strong modulation by the Bond number: when capillary effects are stronger, the variations in Marangoni effects and roughness are stronger with the surfactant parameter $\beta$ .

Figure 7. (a) Time-averaged $\overline {s^2}$ as a function of time-averaged $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ with the colour map for different $\textit{Bo}$ at $ak=0.3$ . Plots of (b) $\overline {s^2}$ and (c) $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ as functions of $\beta$ with different $\textit{Bo}$ at $ak=0.3$ .

4.1. From gravity to parasitic capillary waves

Figure 8(a) presents the phase diagram with the symbols colour-coded with $\overline {\langle |\textit{Ma}_{\textit{loc}}| \rangle }$ , which shows a region of enhanced Marangoni forcing at low $\textit{Bo}$ and intermediate $\beta$ . We observe three wave regimes following the wave patterns described in figure 4: regime I of gravity waves (square markers), at mostly high Bond number; regime II of parasitic capillary waves (circle markers); and regime III at low Bond number and intermediate $\beta$ where parasitic capillary waves are suppressed by Marangoni effects (also circle markers, but with high values of Marangoni number). The classification is based on the presence of surface slope change, determined by a threshold on the local interface gradient.

Figure 8. The phase diagrams of wave patterns colour-coded by (a) $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ , (b) $\overline {s^2}$ , and (c) normalised $\overline {s^2}/\overline {s^2}_{\beta _{min}}$ , all for $ak=0.3$ . The squares represents the pure gravity waves, and the circles represents the parasitic capillary waves (PCW). In (a), the grey solid line represents the contour $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }=25$ from simulations. In (b), the grey solid line represents the contour $\overline {s^2}=0.06$ . In (c), the grey solid line represents the contour $\overline {s^2}/\overline {s^2}_{\beta _{min}}=0.9$ . The dashed lines represent the curve $\textit{Bo}_c(\beta )$ from (4.3), with $\textit{Bo}_c(0)=350$ . The snapshots I, II and III represent the typical wave shapes in each regime, which can be briefly described as smooth gravity wave, PCW, and suppressed PCW.

Regime I occurs at high $\textit{Bo}$ , and the transition from regime II to regime I agrees with the clean-water case when $\beta$ is small (Deike et al. Reference Deike, Popinet and Melville2015) – see more details in Appendix B – but shifts to lower Bond numbers as $\beta$ increases. We plot a contour line for $\overline {\langle |\textit{Ma}_{\textit{loc}}| \rangle }=25$ , which indicates the transition to regime III. Since the transition of surfactant effects is a gradual and continuous process rather than a sharp boundary, the threshold value $\overline {\langle |\textit{Ma}_{\textit{loc}}| \rangle }=25$ was chosen as a representative indicator to visualise the onset of strong Marangoni effects in the phase diagram. The contours serve as visual indicators to identify when the Marangoni effect becomes dominant.

To quantify the wave shape for different control parameters, we plot the phase diagram with $\overline {s^2}$ colour-coded in $(\textit{Bo},\beta )$ space at $ak=0.3$ , shown in figure 8(b). The results also reveal that the influence of $\beta$ on $\overline {s^2}$ is most pronounced at lower Bond numbers, corresponding to the region where $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ is high. At low $\textit{Bo}$ , as $\beta$ increases, $\overline {s^2}$ initially decreases, indicating a strong suppression of parasitic capillaries. Beyond a critical value ( $\beta \approx 0.3$ ), $\overline {s^2}$ begins to increase again with recurrence of parasitic waves. As $\textit{Bo}$ increases, the sensitivity of $\overline {s^2}$ to variations of $\beta$ becomes weaker, reflecting the reduced role of surface tension relative to gravity. We also plot a contour line for $\overline {s^2}=0.06$ , which behaves similarly and has a shape comparable to the contour of $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ in figure 6(a).

To focus on the Marangoni effect and eliminate the effect of $\textit{Bo}$ , we further show normalised $\overline {s^2}$ by $\overline {s^2}_{\beta _{min}}$ at the same $\textit{Bo}$ but smallest $\beta =0.005$ , as shown in figure 8(c). This reveals a significant suppression region, as in figure 8(a,b), down to as low as 40 % of the baseline value. We also include a contour line for $\overline {s^2}/\overline {s^2}_{\beta _{min}}=0.9$ , which closely resembles the contour of $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ and will be used in the following analysis.

More phase diagrams for wave patterns in the $(\textit{Bo}, \beta )$ space are shown in figure 9(a–c), with normalised $\overline {s^2}$ colour-coded, exploring different $ak$ and $\textit{Re}$ . A lower value of the initial wave slope $ak$ is shown in figure 9(a): in comparison to the phase diagram at $ak=0.3,\ \textit{Re}=4\times 10^4$ , the Marangoni effect is weaker in the low- $\textit{Bo}$ regime, leading to weaker suppression of wave deformation. In contrast, the results in the higher- $\textit{Re}$ case (figure 9 b) exhibit stronger Marangoni forcing, with the wave shape significantly suppressed across a broader range of $\textit{Bo}$ .

Figure 9. Wave pattern diagrams in the ( $\textit{Bo},\beta$ ) space with normalised $\overline {s^2}/\overline {s^2}_{\beta _{min}}$ colour-coded, for (a) $ak=0.25,\ \textit{Re}=4\times 10^4$ , (b) $ak=0.3,\ \textit{Re}=10^5$ , and (c) $ak=0.35$ , $\textit{Re}=4\times 10^4$ , where $\overline {s^2}_{\beta _{min}}$ is the case at the same $\textit{Bo}$ but smallest $\beta =0.005$ . The grey solid lines represent the contour $\overline {s^2}/\overline {s^2}_{\beta _{min}}=0.9$ . In (a,b,c), circles are parasitic capillary waves, squares are gravity waves, and the dashed lines represent the curve of $\textit{Bo}_c(\beta )$ from (4.3) for the transition from parasitic capillary wave to gravity wave, with $\textit{Bo}_c(0)=95$ in (a) and $\textit{Bo}_c(0)=350$ in (b). In (c), triangles are spilling breakers, and stars are plunging breakers. The black dashed and dotted lines represent the curves $\textit{Bo}_c(\beta )$ for different regime transitions, with $\textit{Bo}_c(0)=37$ for the dashed line, and $\textit{Bo}_c(0)=460$ for the dotted line, which is in good agreement with the transitions in previous study of wave breaking in clean water (Deike et al. Reference Deike, Popinet and Melville2015); see Appendix B.

For $ak=0.35$ , figure 9(c) exhibits a more complex behaviour, including cases of wave breaking, denoted as triangle symbols corresponding to spilling breakers, and star symbols for plunging breakers (see figure 4 for geometric definitions). The $\beta$ -dependent transition between regimes is again observed. At low $\textit{Bo}$ , strong Marangoni forcing suppresses the formation of parasitic capillary waves. At $\textit{Bo} = 40$ , the Marangoni forcing even inhibits wave breaking, revealing a non-trivial ‘spilling breaker–capillary wave–spilling breaker’ transition as $\beta$ increases, which is also shown in figure 5. This transition resembles that reported by Erinin et al. (Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023) between plunging breakers and spilling breakers when increasing the surfactant concentration. While the specific choice of EOS certainly influences the critical values separating the different regimes, another choice of EOS that would retain the ones discussed experimentally in Erinin et al. (Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023) (with similar effects on reduction of surface tension, almost linear decay, and then saturation of $\sigma$ as a function of $\varGamma$ ) would lead to a similar non-monotonic behaviour.

4.2. Scaling the boundaries

In this subsection, we provide a physical explanation for the mechanisms underlying the effects of surfactants on wave dynamics, which includes two main contributions: (i) the reduction of local surface tension due to the presence of surfactants; and (ii) the generation of Marangoni flows driven by surfactant concentration gradients, which act as additional forcing terms in the system. Consequently, two main outcomes are observed as $\beta$ varies: (1) a shift in the wave regime from capillary-dominated to gravity-dominated behaviour at high $\textit{Bo}$ , and (2) a non-monotonic variation in normalised $\rm \overline {s^2}$ , particularly at low $\textit{Bo}$ .

We first focus on the regime transition. At high $\textit{Bo}$ , the effect of Marangoni forcing becomes less significant, given the inverse dependence of $Ma$ on $\textit{Bo}$ . In this regime, the dominant mechanism is the reduction of surface tension by surfactants. While the input $\textit{Bo}$ is defined using the clean surface tension $\sigma _c$ , the actual surface tension $\sigma (\varGamma _0)$ depends on the surfactant distribution even for a uniform concentration $\varGamma = \varGamma _0$ , as specified by (2.5). Therefore, we introduce an effective Bond number $\textit{Bo}_{\textit{eff}}$ that accounts for the surfactant-modified surface tension, which is defined as

(4.1) \begin{equation} \textit{Bo}_{\textit{eff}}(\beta )=\frac {\Delta \rho\, g}{\sigma (\varGamma _0)k^2}=\textit{Bo}\frac {\sigma _c}{\sigma (\varGamma _0)}=\frac {\textit{Bo}}{1-\Delta \sigma _\infty \tanh (\beta /\Delta \sigma _\infty )}. \end{equation}

To characterise the transition between capillary and gravity wave regimes, we define a critical Bond number $\textit{Bo}_c(\beta )$ as a function of $\beta$ . When no surfactant is present ( $\beta = 0$ ), the critical value reduces to a constant $\textit{Bo}_c(0)$ . The regime transition occurs when the effective Bond number equals this critical value, $\textit{Bo}_{\textit{eff}}(\beta )=\textit{Bo}_c(0)$ , which gives

(4.2) \begin{equation} \frac {\textit{Bo}_c(\beta )}{1-\Delta \sigma _\infty \tanh (\beta /\Delta \sigma _\infty )}=\textit{Bo}_c(0). \end{equation}

Solving for $\textit{Bo}_c(\beta )$ in terms of $\beta$ yields

(4.3) \begin{equation} \textit{Bo}_c(\beta )= \textit{Bo}_c(0)\,(1-\Delta \sigma _\infty \tanh (\beta /\Delta \sigma _\infty )). \end{equation}

The resulting curve for $\textit{Bo}_c(\beta )$ is plotted as a dashed line in the phase diagrams in figures 8 and 9, and it agrees well with the observed transition boundaries in the phase diagrams. It shows that the observed regime transition is governed primarily by the change in effective Bond number due to surfactant-induced surface tension reduction, rather than by Marangoni stresses.

Then we focus on the non-monotonic dependence of $\overline {s^2}$ on $\beta$ observed at low $\textit{Bo}$ , which can be attributed to Marangoni flow induced from non-uniform surfactant distributions. As illustrated in figure 10(a,b), surfactants tend to accumulate at the front face of the wave crest, driven by surface compression associated with wave motion. We note that a dynamic gas pressure gradient exists on the water surface. This variation mainly reflects the dynamic pressure distribution consistent with the orbital motion, rather than the driver of it. The accumulation region coincides with the typical location where parasitic capillary waves would form in the absence of surfactants. However, the local surfactant accumulation reduces the local surface tension and increases the surface tension gradient. This gradient drives Marangoni flow directed opposite to the wave-induced surface motion, thereby reducing the curvature amplification that would otherwise trigger capillary ripples. In other words, instead of directly flattening the interface through curvature change, the Marangoni stresses damp the local surface velocity gradients responsible for increasing curvature. As a result, the formation of parasitic waves is suppressed.

To quantitatively show the non-monotonic dependence of Marangoni forcing on $\beta$ , we analytically show the dependence of the normalised mean surface tension gradient $(\Delta \sigma /\Delta \varGamma) \boldsymbol{\cdot }(\varGamma _0/\sigma _c)$ on $\beta$ based on the EOS (2.5) with different ranges of $\varGamma$ , as shown in figure 10(c). Although its magnitude and peak location vary depending on the averaging range of $\varGamma$ , the non-monotonic trend with respect to $\beta$ is robust. From our simulations, the maximum local $\varGamma /\varGamma _0$ reaches approximately $2.5$ for parasitic capillary waves, and $1.5$ for gravity waves. We choose the range $\varGamma /\varGamma _0 \in [1,2]$ to compute the mean surface tension gradient, which gives the peak of surface tension gradient at $\beta \approx 0.3$ . This peak coincides with the strongest suppression of $\overline {s^2}$ and is independent of $ak$ , $\textit{Bo}$ and $\textit{Re}$ .

Figure 10. (a) Illustration of the wave-driven surface flow pattern without surfactant. (b) Illustration of the redistribution of an initially uniform surfactant field due to wave-driven flow and the corresponding surfactant-driven flow. The surfactant concentration is focused at the crest front, which reduces the local surface tension, which further induces Marangoni flow as the dashed arrows. (c) The normalised characteristic surface tension gradient as function of $\beta$ for different compression ranges $\varGamma /\varGamma _0$ , calculated from the EOS.

Finally, we focus on the effect of Reynolds number $\textit{Re}$ , which scales positively with $\textit{Ma}_{local}$ as in (1.2), and implies that Marangoni effects become stronger at higher Reynolds numbers. We start with the force balance between viscous stress and Marangoni stress:

(4.4) \begin{equation} \mu _w\frac {\partial u}{\partial z}\approx \frac {\partial \sigma }{\partial x}. \end{equation}

Approximating the shear rate by the wave speed $c_0$ over a surface boundary layer thickness $\delta$ , and the surface tension gradient by ${\Delta \sigma }/{\Delta s}$ , this becomes

(4.5) \begin{equation} \mu _w\frac {c_0}{\delta }\approx \frac {{\Delta \sigma }}{{\Delta l}}. \end{equation}

By rearranging, we obtain a dimensionless form

(4.6) \begin{equation} \frac {\sigma _c}{\mu _w c_0}\frac {\Delta \sigma\, \lambda }{\Delta l\,\sigma _c}\frac {\delta }{\lambda }\approx \text{const}. \end{equation}

We further use the scaling $\delta / \lambda \sim \textit{Re}^{-1/2}$ for boundary layer thickness, which leads to an alternative expression

(4.7) \begin{equation} \frac {{\Delta \sigma }\,\lambda }{{\Delta l}\,\sigma _c}\, \textit{Bo}^{-1}\,\textit{Re}^{1/2}\approx \text{const}. \end{equation}

Furthermore, the surface tension gradient can be written as

(4.8) \begin{equation} \frac {{\Delta \sigma }}{{\Delta l}}=\frac {{\Delta \sigma }}{{\Delta \varGamma }}\frac {{\Delta \varGamma }}{{\Delta l}}, \end{equation}

where ${{\Delta \sigma }}/{{\Delta \varGamma }}$ is solely determined by $\beta$ through the EOS, shown in figure 10, while the surfactant concentration gradient ${{\Delta \varGamma }}/{{\Delta l}}$ is driven by surface advection. Figure 11 shows different values of ${{\Delta \sigma }}/{{\Delta \varGamma }}=(\sigma (\varGamma _1)-\sigma (\varGamma _2))/(\varGamma _1-\varGamma _2)$ for typical compression values $\Delta \varGamma /\varGamma _0=(\varGamma _1-\varGamma _2)/\varGamma _0$ . This advection can be approximately quantified by the transport by the surface velocity, which is proportional to $ak$ for a linear wave. Therefore, the final expression is

(4.9) \begin{equation} \frac {{\Delta \sigma }\,\varGamma _0}{{\Delta \varGamma }\,\sigma _c}(\beta )\cdot ak \cdot \textit{Bo}^{-1}\,\textit{Re}^{1/2}\approx \text{const}. \end{equation}

To verify the scaling in (4.9), we extract the contour plot of normalised $\overline {s^2}=0.9$ for all the different $ak$ and $\textit{Re}$ in figure 11(a), since $\overline {s^2}$ is a directly measurable indicator of the surfactant influence on surface roughness, and we have shown that the contours of normalised $\overline {s^2}$ are closely related to the Marangoni forcing. The contour shifts towards higher $\textit{Bo}$ as $\textit{Re}$ and $ak$ increase, indicating a broader region influenced by Marangoni stresses at higher $\textit{Re}$ and $ak$ , consistent with (4.9).

Figure 11. (a) The extracted contour lines of $\overline {s^2}/\overline {s^2}_{\beta _{min}}=0.9$ for different $\textit{Re}$ and $ak$ . (b) The rescaled contour lines by $\textit{Bo}\,\textit{Re}^{-1/2}\,(ak)^{-1}$ , where the dashed line represents the prediction from (4.9). The green circles represent the experimental data from Xu & Perlin (Reference Xu and Perlin2023), where suppression of parasitic capillary wave is observed with surfactant. The blue squares represent the experimental data from Erinin et al. (Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023), where plunging breakers occur for all $\beta$ . The red squares represent the experimental data from Liu & Duncan (Reference Liu and Duncan2006), where spilling breakers occur for all $\beta$ . The triangles show the cases without surfactant ( $\beta =0$ ).

We rescale the horizontal axis of the phase diagram using $\textit{Bo}\,\textit{Re}^{-1/2}\,(ak)^{-1}$ based on the scaling analysis from (4.9). As shown in figure 11(b), this rescaling collapses the threshold contours across different $\textit{Re}$ and $ak$ , and is well predicted by (4.9) represented by the dashed line. At small $\beta$ ( ${\lt } 0.1$ ), the EOS reduces to an approximately linear relation, and the dashed line exhibits the same behaviour. At large $\beta$ , the dashed line is governed by the full nonlinear EOS. The experimental results of Xu & Perlin (Reference Xu and Perlin2023) also show suppression of parasitic capillary waves in the presence of surfactant, indicated by the green circles in figure 11(b). In contrast, in the high- $\textit{Bo}$ regime where Marangoni forcing is weaker, the suppression becomes less pronounced, and plunging breaking persists across different $\beta$ . However, the detailed shape of the plunging breaker still changes due to local surfactant accumulation (Erinin et al. Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023), as shown by the blue squares in figure 11(b).

We further present an expected three-dimensional phase diagram in figure 12 with rescaled $\log (\textit{Bo}\,\textit{Re}^{-1/2})$ , $\beta$ and $ak$ . The $\log (\textit{Bo}\,\textit{Re}^{-1/2}){-}ak$ plane at $\beta =0$ represents the clean-water cases with gravity wave, parasitic capillary wave, spilling breaker and plunging breaker regimes, and corresponding transitions, following Deike et al. (Reference Deike, Popinet and Melville2015). As $\beta$ increases, these transitions shift depending on $\beta$ . The green and red surfaces corresponds to (4.3), where $\textit{Bo}_c(0)$ is obtained from $ak=K(1+\textit{Bo})^{1/3}$ with $K=1.12$ for no-breaking to spilling breaker (green surface), and $K=1.75$ for spilling breaker to plunging breaker (red surface) for the clean-water results from Deike et al. (Reference Deike, Popinet and Melville2015). The horizontal blue surface represents the breaking threshold boundary $ak=0.32$ , also from Deike et al. (Reference Deike, Popinet and Melville2015). The yellow surface represents (4.9); to its left lies the Marangoni-dominated regime where parasitic capillary waves, spilling breakers or plunging breakers are suppressed depending on the specific location within the diagram. We also include our simulation data and experimental data from Liu & Duncan (Reference Liu and Duncan2006), Xu & Perlin (Reference Xu and Perlin2023) and Erinin et al. (Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023), as in figure 11(b), to indicate their locations in this phase diagram (with estimation of $\beta$ based on the surface tension isotherms reported by Erinin et al. Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023).

Figure 12. A three-dimensional phase diagram with rescaled $\log (\textit{Bo}\,\textit{Re}^{-1/2})$ , $\beta$ and $ak$ . The surfaces with different colours represent the expected transitions among different regimes from (4.3) and (4.9). The different symbols represent our simulation data at different regimes and also experimental data from Liu & Duncan (Reference Liu and Duncan2006), Xu & Perlin (Reference Xu and Perlin2023) and Erinin et al. (Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023), as in figure 11(b).

Note that in real ocean environments or experiments, $\textit{Re}$ can be more than two orders of magnitude higher than those considered here (Erinin et al. Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023); the boundary layer becomes fully turbulent. For certain aspects of the effect of surfactant on waves in the ocean, three-dimensional effects would have to be taken into account. In particular, any study on air entrainment, bubble and droplet formation and statistics would have to account for surfactant effects on three-dimensional dynamics, including spanwise Marangoni effects. Spanwise Marangoni effects can also affect the dynamics of three-dimensional wave patterns. Numerically, as an example, assuming the required resolution scales with the boundary layer thickness as $1/2^L \sim \textit{Re}^{-1/2}$ , reaching $\textit{Re} \sim 4\times 10^6$ would demand a minimum grid refinement level $L \approx 13$ and also smaller time steps. This implies that in natural conditions, Marangoni effects might be important for the breaking during the first wave period and extend into higher- $\textit{Bo}$ regions compared to our simulation results, as illustrated by the experimental work of Erinin et al. (Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023) and Liu et al. (Reference Liu, Erinin, Liu and Duncan2025). Further exploration of three-dimensional effects would be a natural avenue for future work, targeting only specific ranges of the parameter space given the expected high computational cost.

Longuet-Higgins (Reference Longuet-Higgins1963, Reference Longuet-Higgins1995) provides a theoretical relationship between the shape of the interface and the formation of parasitic capillary waves (due to capillary pressure resulting in steep waves) and the link with underwater vorticity. In our study, the observed non-monotonic suppression of parasitic capillary waves with increasing $\beta$ is closely related to the same physical balance, but with the addition of variable surface tension and Marangoni stresses. Incorporating the effects of surfactant-dependent surface tension directly into the Longuet-Higgins analytical model would require reformulating the curvature–pressure balance and the surface boundary conditions to include spatially varying $\sigma (\varGamma )$ and its gradients. The non-monotonic trend observed in our results arises from two competing effects: (i) a reduction in mean surface tension that weakens the capillary restoring force; and (ii) enhanced Marangoni stresses at intermediate $\beta$ , which oppose surface compression and suppress the parasitic capillary waves. Developing a fully extended analytical model is non-trivial and beyond the present scope, but could be explored in future work.

5.1. Energy dissipation

We now discuss the effect of surfactant on energy dissipation in the water phase for the various wave regimes. To quantify the total energy dissipation, we consider the total energy $E=E_{\!p}+E_k+E_s$ , including potential energy $E_{\!p}$ , kinetic energy $E_k$ , and surface tension energy $E_s$ , which are calculated as

(5.1) \begin{equation} E_k=\frac {1}{2}\int \rho _w \boldsymbol{u}^2\,{\rm d}x\,{\rm d}y,\quad E_{\!p}=\int \rho _w gy\,{\rm d}x\,{\rm d}y,\quad E_s=\sigma (l-l_0), \end{equation}

where $l_0$ is the surface projection length. We quantify the total energy dissipation by the total energy change from the beginning to the end of the simulation, defined as $\Delta E$ .

We investigate the effect of surfactants on the total energy dissipation $\Delta E/E_0$ from the beginning to $t=2T_0$ . Figure 13(a–c) show the time evolution of the normalised total energy over two wave periods for $ak = 0.3$ at varying $\beta$ and $\textit{Bo}$ . In all cases, the total wave energy decreases in time due to dissipative processes.

Figure 13. The normalised total energy $E/E_0$ as a function of time for $ak=0.3$ and (a) $\textit{Bo}=10$ , (b) $\textit{Bo}=100$ and (c) $\textit{Bo}=1000$ . The amount of energy decay $\Delta E/E_0$ within $t=2T$ as a function of $\beta$ and $\textit{Bo}$ for (d) $ak=0.25$ , (e) $ak=0.3$ , (f) $ak=0.35$ .

At low $\textit{Bo} = 10$ (figure 13 a), significant energy dissipation is already visible even at low $\beta$ due to the presence of parasitic capillary waves (see Deike et al. Reference Deike, Popinet and Melville2015). The effect of $\beta$ on energy dissipation is then observed to be non-monotonic. For both low and high $\beta$ values, the energy decays more slowly than at intermediate $\beta$ values, particularly $\beta \approx 0.3$ . This trend contrasts with earlier observations of surface roughness, where rougher interfaces (associated with more prominent parasitic activity at low and high $\beta$ ) are typically expected to exhibit higher dissipation. However, the enhanced energy loss at intermediate $\beta$ is instead attributed to strong Marangoni flow, which results from larger surface tension gradients. This flow generates increased near-surface shear and vorticity, which is shown in the vorticity field snapshot in figure 2(b) where $\beta = 0.3$ produces more intense vorticity than other cases.

As $\textit{Bo}$ increases, shown in figure 13(b,c) for $\textit{Bo}=100$ and $\textit{Bo}=1000$ , the energy dissipation decreases due to the transition to purely gravity wave with weaker parasitic capillaries. The influence of $\beta$ on dissipation also becomes weaker, consistent with the weakening of Marangoni effects as $\textit{Bo}$ increases.

Figure 13(d–f) show the total energy loss over two wave periods as functions of both $\beta$ and $\textit{Bo}$ for different $ak$ . In all different $ak$ cases, the dissipation exhibits a non-monotonic dependence on $\beta$ , with the maximum at $\beta \approx 0.3$ , corresponding to the analytically strongest Marangoni forcing from EOS. This again supports the conclusion that surfactant-induced surface tension gradients are a key driver of enhanced dissipation at intermediate $\beta$ values.

Figure 14 shows the total energy dissipation $\Delta E/E_0$ as a function of $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ at $ak=0.3$ with different $\textit{Bo}$ . A positive correlation is observed between $\Delta E/E_0$ and $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ . As $\textit{Bo}$ increases, $\Delta E/E_0$ and $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ both decrease. Furthermore, for sufficiently large $\overline {\langle |\textit{Ma}_{\textit{loc}}|\rangle }$ , the data collapse onto a single curve, indicating that in this regime, the energy dissipation is primarily governed by Marangoni forcing. The same plot for other $ak$ and $\textit{Re}$ shows similar trends. This analysis demonstrates that as the Marangoni forcing increases, the energy dissipation also always increases for given values of other control parameters ( $\textit{Bo}$ , $\beta$ or $ak$ ), with the initial magnitude that is controlled by the baseline dissipation given by the wave slope and Bond number.

Figure 14. Plots of $\Delta E/E_0$ as functions of $\overline {\langle \textit{Ma}_{\textit{loc}}\rangle }$ with different $\textit{Bo}$ at $ak=0.3$ .

5.2. Mean surface vorticity

As mentioned above, Longuet-Higgins (Reference Longuet-Higgins1992) provides a link between the underwater generation of vorticity by nonlinear gravity-capillary waves (parasitic capillary waves) in the clean case, which was used by Deike et al. (Reference Deike, Popinet and Melville2015) to quantify the transition between non-breaking and breaking waves:

(5.2) \begin{equation} \varOmega =-2(ak)^2\omega , \end{equation}

where $\varOmega$ is the surface vorticity, and $\omega =\sqrt {gk+\sigma /\rho k^3}$ is the angular frequency from the linear dispersion relation. We examine how surfactant modifies the mean vorticity magnitude $\overline {|\varOmega |}$ .

We measure the mean absolute vorticity $\overline {|\varOmega |}$ averaged within a viscous boundary layer right below the interface (of thickness $\sim \lambda /\sqrt {\textit{Re}}$ ; see Fedorov & Melville Reference Fedorov and Melville1998; Deike et al. Reference Deike, Popinet and Melville2015). In figure 15, we plot $\overline {|\varOmega |}/\omega$ as a function of $ak$ for different values of $\beta$ . The exact value of $\overline {|\varOmega |}/\omega$ depends on the selected averaging region near the surface, but in the clean-water and small- $\beta$ cases, $\overline {|\varOmega |}/\omega$ follows the theoretical $(ak)^2$ scaling predicted by (5.2). As $\beta$ increases, $\overline {|\varOmega |}/\omega$ increases due to the enhanced Marangoni flow, consistent with the enhanced energy dissipation, and it shows a weaker and nearly zero dependence on $ak$ , as Marangoni flow dominates the vorticity at intermediate $\beta$ , which is independent of $ak$ .

Figure 15. Plots of $\overline {|\varOmega |}/\omega$ as functions of $ak$ with different $\beta$ at $ak=0.3$ and $\textit{Bo}=10$ . The solid line shows the $(ak)^2$ scaling (5.2), and the dashed horizontal line shows a relation independent of $ak$ .