1. Introduction
The use of ‘Storm oil’ to calm ocean waves was widely practised by Roman sailors as reported by Franklin et al. (Reference Franklin, Brownrigg and Farish1774) and Lohse (Reference Lohse2023). Aitken (Reference Aitken1884) described how nearly insoluble oil spilled onto water can act as surfactant that calms stormy seas in Aberdeen Harbour. Inspired by such observations, this work focuses on surfactant-covered parametric gravity-capillary standing waves, known as Faraday waves.
Faraday (Reference Faraday1831) noticed that vertically vibrating a fluid layer produces surface waves oscillating at half the driving frequency. Crossing a threshold amplitude, these waves usually organise into patterns such as squares, hexagons and triangles (Westra et al. Reference Westra, Binks and van de Water2003). Complications arise from factors such as contact line dissipation and surface contamination. In this work, we focus on the effects of surface contamination on the Faraday-wave patterns.
Henderson (Reference Henderson1998) found that the measured damping rates with insoluble surfactants agreed with the damping constants of elastic interfacial films (Miles Reference Miles1967). Kumar & Matar (Reference Kumar and Matar2002) presented a linear stability theory for surfactant-covered Faraday waves in the lubrication approximation. Subsequent research (Kumar & Matar Reference Kumar and Matar2004) emphasised the role of the phase difference that influences the Marangoni stresses. Depending on the phase difference, the Marangoni stresses may oppose (in phase) or support (out of phase) the fluid flow. Ubal et al. (Reference Ubal, Giavedoni and Saita2005a,Reference Ubal, Giavedoni and Saitab) computed the two-dimensional numerical simulations of surfactant-covered Faraday waves. However, these studies are limited to linearised one- or two-dimensional models, with some being carried out using lubrication theory. Strickland et al. (Reference Strickland, Shearer and Daniels2015) conducted experiments to investigate the spatiotemporal evolution of the surfactant-covered Faraday waves in a cylindrical container of 
$10$ wavelengths in diameter. Due to the confined region of interest, they observed the combined effects of meniscus and Faraday waves. As a result, the formation of patterns such as squares and hexagons was suppressed because of the superposition of these waves. In a related study, Lau et al. (Reference Lau, Westerweel and van de Water2020) used surfactant-covered Faraday waves to measure ultralow interfacial tension between two immiscible fluids (water-dodecane system). They found that the predicted wavelength from linear stability analysis is two orders of magnitude higher than the experiments, which shows that the Marangoni stresses are crucial in surfactant-covered interfacial waves. Despite this, Lau et al. (Reference Lau, Westerweel and van de Water2020) obtained square patterns in a surfactant-covered system. No further studies have been conducted on pattern formation in surfactant-covered Faraday waves.
Périnet et al. (Reference Périnet, Juric and Tuckerman2009) were the first to perform full three-dimensional direct numerical simulations for the study of Faraday waves. Kahouadji et al. (Reference Kahouadji, Périnet, Tuckerman, Shin, Chergui and Juric2015) further exploited the highly parallelised front tracking code, BLUE (Shin et al. Reference Shin, Chergui and Juric2017), to find supersquare patterns. Ebo-Adou et al. (Reference Ebo-Adou, Tuckerman, Shin, Chergui and Juric2019) employed BLUE to study Faraday waves on a sphere. Recently, Panda et al. (Reference Panda, Kahouadji, Tuckerman, Shin, Chergui, Juric and Matar2023, Reference Panda, Kahouadji, Abdal, Tuckerman, Shin, Chergui and Matar2024) used the same code for studying surface waves on a water drop. Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018) further extended BLUE by including modules to solve surfactant dynamics on the interface as well as in the bulk medium.
In this work, we report the results of simulations of three-dimensional surfactant-covered Faraday waves; we focus on the influence of Marangoni effects on the surface-wave patterns. Our study reveals that the dominance of Marangoni flow leads to transitions away from the square patterns to asymmetric squares, weakly wavy stripes, and ridges and hills. These ridges and hills are new features that occurred on a highly elastic surface. Ridges are found to rise non-uniformly and fall by forming a hill. Our direct numerical simulations help to uncover the rich physics of the dynamics of these newly observed ridges and hills.
This paper is organised as follows. First, we briefly present the problem, scaling and the numerical method. We then present the numerical threshold acceleration which is validated by comparison with the two-dimensional simulations of Ubal et al. (Reference Ubal, Giavedoni and Saita2005b). After that, we present a phase diagram that highlights the influence of Marangoni stresses in the pattern transition of surfactant-covered Faraday waves. These patterns are analysed spectrally. Finally, we explain the newly observed ridges and hills in detail.
Figure 1. (a) Schematic representation of the computational domain: the height of the domain 
$\tilde H = 5.00\, \textrm {mm}$, and the lateral dimensions 
$\tilde {\lambda }_c \times \tilde {\lambda }_c$, where 
$\tilde {\lambda }_c$ is the critical wavelength. No-penetration and no-slip boundary conditions are applied at the bottom and top of the domain and periodic boundaries on the sides. Here, 
$\tilde F$ is the oscillatory driving acceleration, 
$\tilde A$ is the acceleration amplitude, and 
$\omega$ is the angular frequency. (b) Critical acceleration 
$F_c$ for a surfactant-free interface where the solid lines represent the neutral curves for the hydrodynamic parameters listed in Ubal et al. (Reference Ubal, Giavedoni and Saita2005b) and the present work, evaluated using the method of Kumar & Tuckerman (Reference Kumar and Tuckerman1994). Here ‘SH’ and ‘H’ refer to the subharmonic and harmonic tongues. (c,d) Temporal evolution of the total kinetic energy 
$E_k$ for a (c) surfactant-free and (d) surfactant-covered (
$\beta _s = 1.0,\,\Gamma _0 = 0.2$) interface at different acceleration amplitudes 
$F$. The wavelength in both cases is the critical wavelength 
$\tilde {\lambda }_c = 5.3023\,\textrm {mm}$ for the surfactant-free case.
2. Problem formulation, non-dimensionalisation and numerical method
Our computational domain is shown in figure 1(a), which contains a layer of heavy fluid overlaid by light fluid under the gravitational acceleration, 
$g$. We choose a simulation set-up and hydrodynamic parameters based on Ubal et al. (Reference Ubal, Giavedoni and Saita2005b), where the lower heavy fluid is a water–glycerine mixture of depth 
$\tilde h = 1\,\textrm {mm}$, density 
$\tilde \rho _w = 1000\,{\mathrm kg}\, {\mathrm m}^{-3}$ and viscosity 
$\tilde \mu _w = 0.025\, {\mathrm kg}\,{\mathrm m}^{-1}\, {\mathrm s}^{-1}$. Unlike Ubal et al. (Reference Ubal, Giavedoni and Saita2005b), we include an upper air layer of height 
$4\, \textrm {mm}$, density 
$\tilde \rho _a =1.206\,{\mathrm kg}\, {\mathrm m}^{-3}$ and viscosity 
$\tilde \mu _a = 1.82 \times 10^{-5}\,{\mathrm kg}\,{\mathrm m}^{-1}\, {\mathrm s}^{-1}$. The surface tension of the liquid–gas surfactant-free interface is 
$\tilde \sigma _0 = 70 \times 10^{-3}\,{\mathrm kg}\, {\mathrm s}^{-2}$. Due to the low density ratio (
$10^{-3}$) and capillary length 
$l_c = \sqrt{\tilde \sigma_0/((\tilde \rho_w - \tilde \rho_a)g)} = 2.67\ \textrm{mm}$ being smaller than the air layer height, the upper fluid minimally influences the Faraday instability, allowing comparison with Ubal et al. (Reference Ubal, Giavedoni and Saita2005b). The frequency of the external vibration is 
$100$ Hz (angular frequency 
$\omega =2\pi \: 100\ \rm rad\ s^{-1}$).
We consider an insoluble surfactant that is present only on the interface since we consider that the time scale of surfactant desorption from the interface into the bulk is larger than the vibratory time scale. The saturated surfactant concentration at the critical micelle concentration is 
$\tilde \Gamma _\infty \sim {O}(10^{-6})\ \rm mol\ m^{-2}$; the range of surfactant elasticity parameter 
$\beta _s$ (whose definition is discussed in the following section) is 
$0.1 \lt \beta _s \lt 0.9$. The diffusivity for the surfactant 
$\mathcal {D}$ is set to 
$2.5\times 10^{-9}\,{\mathrm m}^{2}\, {\mathrm s}^{-1}$ to align with the work of Ubal et al. (Reference Ubal, Giavedoni and Saita2005a,Reference Ubal, Giavedoni and Saitab). Unless otherwise specified, we set the initial surfactant coverage to 
$\tilde \Gamma _0 = 0.5\tilde \Gamma _\infty$.
We list the major time scales in the problem: (i) the capillary time scale 
$\Delta \tilde t_c = (\tilde \rho _w \tilde h^3/\tilde \sigma _0 )^{1/2}$ of natural capillary oscillations of the perturbed planar interface; (ii) the imposed vibrational time scale 
$\Delta \tilde t_i = 1/\omega$; and (iii) the Marangoni time scale 
$\Delta \tilde t_m = \tilde \mu _w \tilde h/(\tilde \sigma _0-\tilde \sigma (\tilde \Gamma _0))$, where 
$\tilde \sigma$ denotes the surface tension of a surfactant-laden interface, which quantifies the surfactant dynamics on the interface. Our choice of parameters leads to 
$\Delta \tilde t_c \sim {O}(10^{-3})\,\textrm {s}$, 
$\Delta \tilde t_i \sim {O}(10^{-3})\, \textrm {s}$ and 
$\Delta \tilde t_m \sim {O}(10^{-4}{-}10^{-3})\, \textrm {s}$. This choice ensures that we observe a competition between the vibrational, capillary and Marangoni effects.
We choose the height of the liquid 
$\tilde h$ as the length scale, the inverse angular frequency 
$1/\omega$ as the time scale, and 
$\tilde \rho _w \omega ^2 \tilde h^2$ as the pressure scale. Finally, the interfacial concentration 
$\tilde \Gamma$ is scaled by the saturated interfacial concentration 
$\tilde \Gamma _\infty$. The dimensionless hydrodynamic equations are then written as
\begin{equation} {\boldsymbol \nabla} \cdot {\mathbf u} = 0, \end{equation}
\begin{eqnarray} {\rho }\left (\frac {\partial {\mathbf {u}}}{\partial {t}}+ {\mathbf {u}}\cdot {\boldsymbol{\nabla}}{\mathbf {u}}\right ) &=& -{\boldsymbol{\nabla}}{p} - \frac {\rho }{Fr^2}\left (1-F\cos {{t}}\right )\mathbf {i}_z+ \frac {1}{Re}{\boldsymbol{\nabla}}\cdot \left [{\mu }({\boldsymbol{\nabla}}{\mathbf {u}}+{\boldsymbol{\nabla}}{\mathbf {u}}^T)\right ] \nonumber \\ && \mbox {}+ \frac {1}{We}\int _{{\mathcal {A}}({t})} \left (\sigma \kappa \mathbf {n} + {\boldsymbol{\nabla}}_s {\sigma }\right ) \delta ({\mathbf {x}}-{\mathbf {x}}_f) {\rm d}{\mathcal {A}}. \end{eqnarray}
Here, the dimensionless density and dynamic viscosity are given by
\begin{equation} \rho = \tilde{\rho}_a/\tilde \rho _w + (1 - \tilde{\rho} _a/\tilde{\rho}_w)\mathcal {H}({{\tilde{\mathbf {x}}}}, \tilde t), \,\mu = \tilde{\mu}_a/\tilde{\mu}_w + (1 - \tilde{\mu}_a/\tilde{\mu} _w)\mathcal {H}({\tilde{\mathbf{x}}},\tilde{t}), \end{equation}
where 
$\mathcal {H}({\tilde{\mathbf{x}}}, \tilde t)$ is the Heaviside function, which is set to 
$0$ for air (subscript 
$a$) and 
$1$ for water (subscript 
$w$). The last term on the right-hand side of (2.2) corresponds to the surface force at the interface 
${\mathbf {x}}={\mathbf {x}}_f$ and 
$\mathbf n$ is the vector normal to the interface. Inside the integral, the first and second terms account for forces arising from the normal and tangential stresses; the latter are the Marangoni stresses induced by the presence of surface tension gradients. Here 
$\mathcal {A}(t)$ denotes the dimensionless time-dependent interfacial area. The interfacial concentration 
$\Gamma$ evolves according to
\begin{equation} \frac {\partial \Gamma }{\partial t} + \boldsymbol{\nabla} _s \cdot ( \Gamma {\mathbf {u}}_s) = \frac {1}{Pe} \boldsymbol{\nabla} _s^2 {\Gamma }, \end{equation}
where 
$\mathbf u_s$ is the surface velocity, 
${\boldsymbol{\nabla}}_s$ (
$\boldsymbol{\nabla} _s \cdot$) is the gradient (divergence) in the plane locally tangent to the interface. As justified by Stone (Reference Stone1990) and used by Muradoglu & Tryggvason (Reference Muradoglu and Tryggvason2008) and Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018), in the Lagrangian frame of reference the velocity satisfies 
$\mathbf u \cdot \mathbf n = 0$, which gives rise to (2.4). The dimensionless parameters in (2.2) and (2.4) are the Reynolds, Weber, Péclet and Froude numbers, and the ratio of imposed acceleration 
$\tilde A$ to gravitational acceleration 
$g$:
\begin{equation} Re = \frac {\omega \tilde \rho _w \tilde h^2 }{\tilde \mu _w}, \,\,\, We = \frac {\omega ^2\tilde \rho _w \tilde h^3}{\tilde \sigma _0}, \,\,\, Pe = \frac {\omega \tilde h^2}{\mathcal {D}}, \,\,\, Fr = \omega \sqrt {\frac {\tilde h}{g}},\,\,\, F=\frac {\tilde A}{g}. \end{equation}
The surfactant dynamics is coupled with the hydrodynamics through the nonlinear Langmuir equation of state given by
\begin{equation} \sigma = \max \left [0.05, 1 + \beta _s \ln {\left (1 - {\Gamma }\right )}\right ], \,\,\,\,\, \beta _s \equiv \frac {\mathcal {R}\tilde T\tilde \Gamma _\infty }{\tilde \sigma _0}, \end{equation}
where 
$\beta _s$ is the surfactant elasticity number measuring the sensitivity of the surface tension to the surfactant concentration and where the lower limit of 
$\sigma$ has been set to 
$0.05$, below which the Langmuir equation of state may diverge. Here, 
$\tilde T$ is the temperature and 
$\mathcal R$ is the universal gas constant.
The interplay between the surfactant-dependent elastic interface, capillary forces and vibratory inertial forces are captured by 
$\beta _s$ and 
$We$. Scaling the stresses by vibratory inertial stress (
$\tilde \rho _w \omega ^2 \tilde h^2$) and surface tension by the surfactant-free surface tension (
$\tilde \sigma _0$), we obtain the Marangoni stress, 
$\tau$, along the tangent 
$\mathbf t$ to the interface:
\begin{equation} \tau \equiv \frac {1}{We}{\boldsymbol{\nabla}}_s {\sigma }\cdot \mathbf {t} = -\frac {Ma}{(1-\Gamma )}{\boldsymbol{\nabla}}_s {\Gamma }\cdot \mathbf {t}, \end{equation}
where 
$Ma = \beta _s/ We$ is the Marangoni number. Here 
$Ma$ describes the competition between the surfactant-dependent surface elastic forces and the vibratory inertial forces. Rescaling the surface tension by the difference between the surfactant-free and surfactant-covered surface tension (
$\tilde \sigma (\tilde \Gamma _0) + \sigma (\tilde \sigma _0 - \tilde \sigma (\tilde \Gamma _0))$) and stresses by the viscous stresses (
$\tilde \mu _w \omega$), we obtain a dimensionless parameter 
$B$, given by,
\begin{equation} B \equiv \frac {\tilde \sigma _0 - \tilde \sigma (\tilde \Gamma _0)}{\omega \tilde \mu _w \tilde h} = -\frac {\tilde \sigma _0 \beta _s \ln {(1-\tilde \Gamma _0/\tilde \Gamma _\infty )}}{\omega \tilde \mu _w \tilde h}. \end{equation}
Kumar & Matar (Reference Kumar and Matar2004) refer to 
$B$ as the Marangoni number. Because 
$B$, unlike 
$Ma$, includes the surface tension gradient forces, as well as the vibratory inertial and viscous forces, we will use 
$B$ as a control parameter in the remainder of the paper.
We refer to Shin et al. (Reference Shin, Chergui and Juric2017, Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018) for an exhaustive description of the numerical implementation, parallelisation and validation of the numerical framework which we briefly outline here. The spatial derivatives on the Eulerian grid are calculated using a standard cell-centred scheme, except for the nonlinear convective term for which we implemented an essentially non-oscillatory (ENO) procedure on a staggered grid. Peskin’s immersed boundary method is used to couple the Eulerian and the Lagrangian grids. The advection of the Lagrangian field 
$\mathbf x_f(t+\Delta t) = \int _t^{t+\Delta t} \mathbf u_f(t) {\rm d}t$, where 
$\mathbf u_f(t)$ is the interpolated velocity at the interface at time 
$t$, is accomplished by second-order Runge–Kutta numerical integration. A resolution of 
$\Delta x = \Delta y = \lambda _c/44$ was found to be necessary to capture the Faraday-wave dynamics in Périnet et al. (Reference Périnet, Juric and Tuckerman2009) and Kahouadji et al. (Reference Kahouadji, Périnet, Tuckerman, Shin, Chergui and Juric2015). We choose a finer resolution of 
$\lambda _c/128$ to capture the coupling with the surfactant dynamics. Adaptive time stepping is used, where the time step is restricted by the surfactant diffusion, advection, capillary, viscous and vibratory numerical time-step criteria. More details can be found in Kahouadji et al. (Reference Kahouadji, Périnet, Tuckerman, Shin, Chergui and Juric2015) and Muradoglu & Tryggvason (Reference Muradoglu and Tryggvason2008).
3. Results and discussion
We begin by computing the Faraday-wave threshold on the surfactant-free (clean) surface using the method for linear stability analysis detailed in Kumar & Tuckerman (Reference Kumar and Tuckerman1994). We determined that the critical acceleration amplitude 
$F_c$ and wavelength 
$\lambda _c$ are 
$12.34$ and 
$5.3023$, respectively (see figure 1b). We can also compute a threshold from our nonlinear numerical simulations by computing the initial growth rates of the total kinetic energy 
$E_k$ for several values of 
$F$ near 
$F_c$. Since the growth rate varies linearly with the acceleration near the threshold, we can compute the threshold 
$F_c^N$ by linear interpolation. Furthermore, since our focus is on evaluating the growth rate, we selected a lateral surface area of 
$\lambda _c \times \lambda _c$. For a surfactant-free interface, we considered three acceleration amplitudes 
$F = (0.9,1,1.1)F_c$, as shown in figure 1(b). Interpolation to zero growth rate yields 
$F_c^N=12.32$, which differs by only 
$0.16\,\%$ from the theoretical 
$F_c$, as shown in the first line of table 1.
Table 1. Numerical threshold acceleration 
$F_c^N(B)$ for surfactant-free and surfactant-covered interfaces for varying initial surfactant coverage 
$\Gamma _0$ and elasticity number 
$\beta _s$ and a fixed wavelength 
$\lambda _c = 5.3023$. The surfactant-free critical acceleration 
$F_c = 12.34$ is obtained by using the linear stability method of Kumar & Tuckerman (Reference Kumar and Tuckerman1994). The table demonstrates the agreement of our thresholds with those of Ubal et al. (Reference Ubal, Giavedoni and Saita2005b) via 
$\delta ^{\textrm {Ubal}}(B) \equiv |F_c^N(B)-F_c^{\rm Ubal}(B)|/F_c^N(B)$. The last column presents the increase in the Faraday threshold due to surfactant coverage via 
$\Delta \equiv (F^N_c(B)-F_c)/F_c$.
A theoretical linear stability analysis such as that of Kumar & Tuckerman (Reference Kumar and Tuckerman1994) for a surfactant-covered interface would require linearising the Langmuir equation of state (Kumar & Matar 2002, 2004), a task that has not yet been carried out. However, we can compute the acceleration of the numerical threshold 
$F_c^N(B)$ using the procedure described earlier. We compute growth rates from numerical simulations with surfactant-covered 
$\beta _s = 1$ interfaces for different initial surfactant coverage 
$\Gamma _0$ (and corresponding values of 
$B$). Although the critical wavelength varies with the elasticity number (Kumar & Matar Reference Kumar and Matar2004), we approximate it by its surfactant-free value. The resulting thresholds 
$F_c^N$ are displayed in the next three rows of table 1. The same computations were carried out by Ubal et al. (2005a) and Ubal et al. (Reference Ubal, Giavedoni and Saita2005b) using a two-dimensional finite-element technique. Their values are displayed as 
$F_c^{\textrm {Ubal}}(B)$ in table 1. The relative errors 
$\delta ^{\textrm {Ubal}}\equiv |F_c^N(B)-F_c^{\textrm {Ubal}}(B)|/F_c^N$ between our results and those of Ubal et al. (Reference Ubal, Giavedoni and Saita2005b) are less than 
$0.7\,\%$. The last column of table 1 shows the strong dependence of the Faraday threshold on the surfactant coverage via the relative increase 
$\Delta \equiv |F_c^N(B)-F_c|/F_c$. Our results show that increasing 
$B$ stabilises the interface, as observed in previous studies (Henderson Reference Henderson1998; Ubal et al. Reference Ubal, Giavedoni and Saita2005a).
Table 2 shows the increase in the Faraday threshold for many other values of elasticity number 
$\beta _s$ and surfactant coverage 
$\Gamma _0$. The damping rate increases with either of these parameters, leading to an increase in the threshold of Faraday waves. We note that the threshold depends almost entirely on their combination, 
$B$; that is, when 
$\beta _s$ and 
$\Gamma$ are varied so as to produce the same value of 
$B$, then 
$F_c^N$ is unchanged. See, for example, the parameter pairs 
$(\beta _s=0.45, \Gamma _0=0.40)$, which yield 
$B=1.02$, 
$F_c^N=43.0$ and 
$(\beta _s=0.80,\Gamma _0=0.25)$, which yield 
$B=1.03$, 
$F_c^N=43.6$. Other pairs of 
$(\beta ,\Gamma _0)$ values that yield very close values of 
$B$ and 
$F_c^N$ can also be seen in table 2.
Table 2. Numerical threshold acceleration 
$F_c^N$ for wavelength 
$\lambda _c$ and varying 
$\beta _s$, 
$\Gamma _0$ and 
$B$, and its relative increase 
$\Delta \equiv (F_c^N-F_c)/F_c$ from the surfactant-free case. The bold data are used in figure 2.
After a transient phase, Faraday waves appear, which correspond to subharmonic waves whose amplitude is steady and whose response period 
$T$ is twice that of the forcing period. We set 
$t=0$ to be an instant within the steady-amplitude Faraday-wave regime. The computations for assessing the influence of 
$B$ on interfacial dynamics in the nonlinear regime are conducted for 
$F = 1.1F_c^N$ on a larger lateral surface area of 
$3\lambda _c \times 3\lambda _c$. In the surfactant-free case, we observe square patterns.
As shown in figure 2(a), for 
$B \lt 1$ (dark blue dots, purple region), the interface exhibits square symmetry. In a narrow band of 
$1\leqslant B\leqslant 1.23$ (light blue dots), the vertical and horizontal directions differ slightly; we refer to these patterns as asymmetric squares. Within 
$1.23\leqslant B \leqslant 1.46$ (orange dots), the asymmetric square pattern undergoes a transition to weakly wavy stripes. Ridges (ellipses whose major axes are in the 
$y$-direction) appear very faintly as dots for 
$B = 1.23$, 
$t = 3T/4$, and more prominently on the wavy stripes for 
$B=1.51$. For 
$B=1.83$, 
$t=0$, one can also see circular hills between each set of ridges. The hills are the continuation of the ridges formed in the previous half-period. One such instance is shown at 
$t = 3T/4$, where the ridges have disappeared but the hills are present. We explore below the role of Marangoni stresses in the formation of these patterns.
Figure 2. (a) Phase diagram in the 
$\beta _s{-}\Gamma _0$ parameter plane showing the inertia-dominated (violet) and Marangoni-dominated (pink) regions. The solid, dotted, dot-dashed and dashed lines correspond to the 
$B = 1, 1.23, 1.46, 2$ contours, respectively. The four typical patterns are squares, asymmetric squares, weakly wavy stripes, and ridges and hills. The phase boundaries are accurate to within 
$\Delta B = \pm 0.1$. The corresponding values of 
$B$ and 
$F_c^N$ are reported in table 2. (b) Spatiotemporal evolution of the surface deflection 
$\zeta$ over one time period is shown from left to right; squares (
$B = 0.92$), asymmetric squares (
$B = 1.23$), weakly wavy stripes (
$B = 1.51$), and ridges and hills (
$B = 1.83$) are shown from top to bottom rows, respectively. Panels (c,d) show 
$\zeta _{mn}$ and 
$\Gamma _{mn}$, the maximal magnitudes over time of the 
$\zeta$ and 
$\Gamma$ fourier coefficients, respectively, as a function of 
$B$.
Figure 3. (a–d) Three-dimensional visualisation of the surface. (a) Rise of ridges and necking process at 
$t = T/4$ and (b) maximum rise of the ridge at 
$t = 3T/8$. (c) Prominent hill on the ridge at 
$t = T/2$. (d) Falling hill at 
$t = 3T/4$. (e–h) Two-dimensional projections on 
$x{-}z$ slice containing interface curve 
$s$ (indicated in (a)) for 
$t = T/4, 3T/8, T/2$ and 
$3T/4$, respectively. A half-wavelength (ridge to trough) is shown. Colour-coding of the plane indicates 
$y-$vorticity 
$\omega _y$, while streamlines show flow in the 
$x{-}z$ plane. The interface curve 
$s$ is coloured according to the surfactant concentration. Red dots indicate the point of maximum curvature. (i–l) Tangential (see (f)) Marangoni stress and velocity along 
$s$ at 
$t = T/4, 3T/8, T/2$ and 
$3T/4$ shown as black and red curves, respectively. When the sign of one of these quantities is positive (negative), its direction points rightwards (leftwards) from the apex (trough) through the neck to the trough (apex) of the ridge, as indicated at the top (bottom) of figure 3(i). The vertical dashed line indicates the necking region, shown as the red dot in the corresponding 
$x{-}z$ projection. The length of 
$s$ decreases from approximately 
$5$ at 
$t=T/4, 3T/8,$ 
$T/2$ to approximately 
$4$ at 
$t=3T/4$, as can be seen in the curves in (e–h).
To quantify the patterns, we evaluate the spatial Fourier spectra for the surface height, 
$\zeta$, and surfactant concentration, 
$\Gamma$, defining 
$\hat {\zeta }_{mn}(t)$ and 
$\hat {\Gamma }_{mn}(t)$ to be the Fourier coefficients associated with the 
$(x,y)$ wavevector 
$\boldsymbol{k}_{mn}$. We then set 
$\zeta _{mn} \equiv \max _{ [t,t+T ]}{|\hat \zeta _{mn}(t)|}$ and 
$\Gamma _{mn} \equiv \max _{ [t,t+T ]}{|\hat \Gamma _{mn}(t)|}$. Figures 2(c,d) present an overview of the spatial Fourier spectra of 
$\zeta$ and 
$\Gamma$ as a function of 
$B$ in the range 
$B \in [0, 1.51]$. At higher 
$B$, ridges and hills emerge, where steep spatial gradients and many higher spatial harmonics appear.
For 
$B \lt 1$, the square pattern is characterised by comparable amplitudes of 
$\zeta _{10}$ and 
$\zeta _{01}$, as shown in figure 2(c). For 
$B \lt 0.5$, where Marangoni effects are weak, the 
$\zeta _{mn}$ modes have magnitudes similar to those associated with the clean case corresponding to 
$B = 0$, consistent with previous findings (Constante-Amores et al. Reference Constante-Amores, Batchvarov, Kahouadji, Shin, Chergui and Matar2021). For 
$B\gt 1$, Marangoni-driven stresses dominate over inertial effects. The square symmetry is broken, and by 
$B = 1.23$, 
$\zeta _{10}$ surpasses 
$\zeta _{01}$, with an increase in higher-order modes, such as the 
$\zeta _{20}$ mode. As 
$B$ increases further, strong 
$x$-dependent modes emerge, leading to a transition from asymmetric squares to stripes (see figure 2c).
A parallel change occurs in the 
$\Gamma$-spectrum. For 
$B \lt 1$, the surfactant is advected without being significantly hindered by Marangoni stresses, aligning the 
$\Gamma$-spectrum with the 
$\zeta$-spectrum, where 
$\Gamma _{10}$ and 
$\Gamma _{01}$ dominate (see figure 2d). For 
$B\gt 1$, 
$\Gamma _{10}$ begins to surpass 
$\Gamma _{01}$. Thus, 
$B \approx 1$ is a pivotal point in the dynamics, at which there is an equilibrium between the opposing mechanisms of advection-driven surfactant inhomogeneity and Marangoni-driven homogeneity.
We now turn to the formation of hills and ridges on the interface. Figures 3(a–d) illustrate the evolution of a small portion of the interface, colour-coded by surfactant concentration. During the first half-cycle, the ridges rise, and the fluid and surfactant flow up from the troughs, advecting the surfactant to the apex of the ridge. Figures 3(e–h) show two-dimensional projections containing arc 
$s$, as indicated in figure 3(a). As the surfactant is advected towards the apex, a 
$\Gamma$-deficit (higher 
$\sigma$) is created at the trough.
The capillary force resulting from the 
$\Gamma$-deficit leads to the emergence of a bulb on the ridge, surrounded by a narrow region of negative curvature, which we call a neck, and which is highlighted by a red spot on the interface in figures 3(e–h); 
$\Gamma$ accumulates at the ends of the ridge as shown in figures 3(a,b). Marangoni stresses along 
$s$ counteract the 
$\Gamma$-inhomogeneity caused by the surface advection. This is shown in figure 3(i), where 
$\nabla_s \sigma \cdot \mathbf t = \tau > 0$ and 
$\mathbf u_s \cdot \mathbf t = u_t < 0$ along the arc 
$s$. We call this a barrier. This barrier rigidifies the surface during the first half-cycle, leading to 
$|\mathbf u_s| \approx 0$ at 
$t = 3T/8$, as shown in figure 3(b).
The negative vorticity along the surface in figure 3(f) indicates that 
$\tau$ opposes the surface advection. Due to this barrier, a backflow develops on the surface from the apex towards the neck, as indicated by 
$u_t \gt 0$ in the inset of figure 3(j). This drives surfactants from the apex towards the neck. Simultaneously, the accumulated surfactant at the ends of the ridge flows towards the neck due to a similar mechanism, as illustrated by the red arrows in figure 3(b). During this process, the midpoint of the ridge rises to form a bulb; see figure 3(c). By 
$t=T/2$, 
$\Gamma$ is maximal (so 
$\sigma$ is minimal) at the neck.
The accumulated surfactant causes Marangoni stresses, with distinct peaks of 
$\tau \gt 0$ and 
$\tau \lt 0$ across the neck (figure 3k). The barrier is now formed at the neck (shown as a white dotted region in figure 3c) where these stresses in the region between the apex and the neck begin to oppose the flow reversal at half-cycle. Meanwhile, surface tension decreases at the neck. As a result, the neck begins to reopen (see the streamlines in figure 3g) as is commonly observed in surfactant-laden neck reopening phenomena, discussed in detail in Constante-Amores et al. (Reference Constante-Amores, Batchvarov, Kahouadji, Shin, Chergui and Matar2021).
In the next half-cycle (
$t \geqslant T/2$), the ridge begins to fall. However, the opposing Marangoni stress between the neck and the apex (
$\tau \lt 0$ in figure 3k) slows the collapse of this region. This slower descent of 
$u_t$ (
${\rm d}u_t/{\rm d}s \approx 0$) leads to the formation of the hill on the ridge. Meanwhile, at 
$t = 3T/4$, the region between the neck and the trough continues to fall more quickly than the hill. This accelerated fall is driven by the surfactant gradients towards the trough (
$\tau \gt 0$ as shown in figure 3l) which, instead of opposing the bulk flow as before, now begin to support it due to 
$u_t \gt 0$. The presence of two high-vorticity regions (blue zones) along the interface in figure 3(h) is an effect of the two distinct roles of Marangoni stresses at the neck. As a consequence, new ridges develop while the hills of the previous cycle are still present, as seen in figure 2 at 
$t = 0$, 
$B=1.83$.
Figure 4. Marangoni-influenced ridge formation: (a) 
$x{-}z$ projection containing 
$s$, as defined in figure 3, at 
$t = T$; the colour-coding used here is that of figure 3. (b) Three-dimensional visualisation of the interface colour-coded by the magnitude of Marangoni stresses 
$|\boldsymbol{\nabla} _s \sigma |$, indicating the barriers around the rising ridge.
Figure 4 further elucidates the mechanism of ridge formation. The surfactant accumulates on the developing ridge due to a combination of Marangoni-driven surface flow from the neck to the trough (
$\tau \gt 0$ in figure 3h,l) as previously discussed, and advection through bulk flow in the second half-cycle leading to strong surface compression at the ridge. This accumulation (see 
$\Gamma$-surplus region highlighted in figure 4a) generates a Marangoni stress, directed from the newly developed ridge toward the falling hill (as highlighted by the arrow indicating the direction of 
$\tau$ in figure 4a). The magnitude of the Marangoni stress, 
$|\boldsymbol{\nabla} _s \sigma |$, is shown in figure 4(b). This high-stress region, which surrounds the developing ridge, highlights the strength of the barrier to ridge formation. Close inspection of this region reveals that the barrier is weaker at the midpoint of the ridge, allowing stronger inward-directed surface flow to this region (viz. the velocity glyphs in figure 4b). This, in turn, leads to a higher elevation at the midpoint of the ridge than at its ends, as shown in figure 4(b).
4. Conclusion
The study highlights the role of Marangoni stresses in Faraday-wave-pattern transitions. Numerical simulations were validated against previously reported two-dimensional simulation. Using the parameter 
$B$ to compare the Marangoni and inertial time scales, we found that the threshold acceleration increases with 
$B$. After we evaluated 
$B$, we increment the acceleration by 10 % of their respective threshold acceleration. Square patterns are observed for the surfactant-free interface. For the surfactant-covered interface, we found four different patterns as we increased 
$B$. We showed that at 
$B \approx 1$, square patterns transition to asymmetric squares. Increasing Marangoni strength further, asymmetric squares change to weakly wavy stripes. The novel finding highlighted here is the fact that at further higher 
$B$ values, ridges and hills appear. Due to strong Marangoni flow during a cycle of forcing, surfactant and surface flow compete (
$\tau \gt 0$ and 
$u_t \lt 0$), which we call a barrier. The barrier slows down a rising ridge which then reaches its maximum height, resembling a bulb, at 
$t = 3T/8$. While the bulb falls in the next half-cycle, a 
$\Gamma$-surplus region forms at the neck of the ridge. This creates a bi-directional Marangoni stress, where the flow is opposed (supported) between the apex (neck) and the neck (trough). This led to a faster collapse of the ridge between the neck and the trough. However, the bulb falls at a slower rate resembling a hill structure on the ridge. In the next cycle, 
$\Gamma$ accumulates at the newly forming crest. The barrier is weaker at the midpoint than at the sides of the rising crest. This creates a faster rise of the midpoint of the crest, resembling a ridge structure. The existence of such a barrier at the newly forming crest and at the neck are the cause of the formation of these interesting ridges and hills.
Acknowledgements.
D.P. and L.K. acknowledge HPC facilities provided by the Imperial College London Research Computing Service. D.P. acknowledges Imperial College President's scholarship for his doctoral studies. D.J. and J.C. acknowledge support through 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 Equipement National de Calcul Intensif) grant 2024 A0162B06721. The numerical simulations were performed with code BLUE (Shin et al. Reference Shin, Chergui and Juric2017, Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018) and the visualisations were generated using ParaView.
Funding.
This work was supported by the Engineering and Physical Sciences Research Council, UK, through the PREMIERE (EP/T000414/1) programme grant and the ANTENNA Prosperity Partnership (EP/V056891/1). O.K.M. acknowledges funding from PETRONAS and the Royal Academy of Engineering for a Research Chair in Multiphase Fluid Dynamics.
Declaration of interests.
The authors report no conflict of interest.
Open access
