Marangoni-driven patterns, ridges and hills in surfactant-covered parametric surface waves

Abstract

Parametric oscillations of an interface separating two fluid phases create nonlinear surface waves, called Faraday waves, which organise into simple patterns, such as squares and hexagons, as well as complex structures, such as double hexagonal and superlattice patterns. In this work, we study the influence of surfactant-induced Marangoni stresses on the formation and transition of Faraday-wave patterns. We use a control parameter, $B$, that assesses the relative importance of Marangoni stresses as compared with the surface-wave dynamics. Our results show that the threshold acceleration required to destabilise a surfactant-covered interface through vibration increases with increasing $B$. For a surfactant-free interface, a square-wave pattern is observed. As $B$ is incremented, we report transitions from squares to asymmetric squares, weakly wavy stripes and ultimately to ridges and hills. These hills are a consequence of the bidirectional Marangoni stresses at the neck of the ridges. The mechanisms underlying the pattern transitions and the formation of exotic ridges and hills are discussed.

1. Introduction

The use of ‘Storm oil’ to calm ocean waves was widely practised by Roman sailors as reported by Franklin et al. (1774) and Lohse (2023). Aitken (1884) 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 (1831) 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. 2003). 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 (1998) found that the measured damping rates with insoluble surfactants agreed with the damping constants of elastic interfacial films (Miles 1967). Kumar & Matar (2002) presented a linear stability theory for surfactant-covered Faraday waves in the lubrication approximation. Subsequent research (Kumar & Matar 2004) 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. (2005a,b) 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. (2015) 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. (2020) 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. (2020) 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. (2009) were the first to perform full three-dimensional direct numerical simulations for the study of Faraday waves. Kahouadji et al. (2015) further exploited the highly parallelised front tracking code, BLUE (Shin et al. 2017), to find supersquare patterns. Ebo-Adou et al. (2019) employed BLUE to study Faraday waves on a sphere. Recently, Panda et al. (2023, 2024) used the same code for studying surface waves on a water drop. Shin et al. (2018) 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. (2005b). 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. (2005b) and the present work, evaluated using the method of Kumar & Tuckerman (1994). 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. (2005b), 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. (2005b), 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. (2005b). 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. (2005a,b). 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

(2.1) \begin{equation} {\boldsymbol \nabla} \cdot {\mathbf u} = 0, \end{equation}

(2.2) \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

(2.3) \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

(2.4) \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 (1990) and used by Muradoglu & Tryggvason (2008) and Shin et al. (2018), 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$:

(2.5) \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

(2.6) \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:

(2.7) \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,

(2.8) \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 (2004) 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. (2017, 2018) 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. (2009) and Kahouadji et al. (2015). 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. (2015) and Muradoglu & Tryggvason (2008).

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 (1994). 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 (1994). The table demonstrates the agreement of our thresholds with those of Ubal et al. (2005b) 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 (1994) 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 2004), 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. (2005b) 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. (2005b) 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 1998; Ubal et al. 2005a).

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. 2021). 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. (2021).

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.