2.1. Nonlinear equations

The fluid motion is governed by the continuity and Navier–Stokes equations while the temperature field of the system is governed by the energy equation, i.e.

(2.2)$$\begin{gather} \boldsymbol{\nabla^*} \boldsymbol{\cdot} \boldsymbol{v^*} =0, \end{gather}$$

(2.3)$$\begin{gather}\rho \left(\frac{\partial \boldsymbol{v^*}}{\partial t^*}+\boldsymbol{v^*}\boldsymbol{\cdot} \boldsymbol{\nabla^*} \boldsymbol{v^* }\right) =-\boldsymbol{\nabla^*} {p^*} + \mu \nabla^{*2} \boldsymbol{v^*} - \rho ( A^*\omega{^*}^2 \cos(\omega^* t^*))\boldsymbol{e_z} \end{gather}$$

and

(2.4)\begin{equation} \rho C_p \left(\frac{\partial T^*}{\partial t^*}+\boldsymbol{v^*}\boldsymbol{\cdot}\boldsymbol{\nabla^*} {T^*} \right) = \lambda\nabla^{*2} T^*. \end{equation}

In (2.2)–(2.4), $\boldsymbol {e_z}$ is the unit vector along the $z$ direction. Here, $A^*$ and $\omega ^*$ are the amplitude and the angular frequency of the parametric forcing.

At the bottom wall, $z^* = -d^*$, we impose no-slip and no-penetration conditions for velocity and a Dirichlet condition for temperature i.e.

(2.5a,b)\begin{equation} \boldsymbol{v^*}=\boldsymbol{0}\quad \textrm{and}\quad T^*=T^*_H. \end{equation}

At the impenetrable interface, $z^*=h^*({x^*},t^*)$, the jump mass balance yields

(2.6)\begin{equation} (\boldsymbol{v^*}- \boldsymbol{u^*})\boldsymbol{\cdot} \boldsymbol{n^*}=0, \end{equation}

for which the interface velocity vector, $\boldsymbol {u^*}$, the unit normal, $\boldsymbol {n^*}$, and unit tangent vectors, $\boldsymbol {t^*}$, are given by

(2.7a–c)\begin{align} \boldsymbol{u^*}\boldsymbol{\cdot} \boldsymbol{n^*}=\dfrac{\dfrac{\partial h^*}{\partial t^*}}{\left(1+\left(\dfrac{\partial h^*}{\partial x^*}\right)^{2} \right)^{1/2}},\quad \boldsymbol{n^*}= \dfrac{-\dfrac{\partial h^*}{\partial x^*} \boldsymbol{e_{x}} +\boldsymbol{e_z}}{\left(1+\left(\dfrac{\partial h^*}{\partial x^*}\right)^2 \right)^{1/2}}, \quad \boldsymbol{t^*}= \dfrac{\boldsymbol{e_x}+\dfrac{\partial h^*}{\partial x^* } \boldsymbol{e_{z}} }{\left(1+\left(\dfrac{\partial h^*}{\partial x^*}\right)^2 \right)^{1/2}}. \end{align}

The fluid loses heat from the free surface through a purely conducting gas layer via Newton's law of cooling, i.e.

(2.8)\begin{equation} \mathcal{H}(T^*-T^*_{a})=-\lambda \boldsymbol{\nabla^*} T^*\boldsymbol{\cdot}\boldsymbol{n^*}, \end{equation}

where the heat transfer coefficient is defined via $\mathcal {H}={\lambda _a}/{d^*_a}$, with $\lambda _a$ and $d^*_a$ being the thermal conductivity and height of the gas layer. The force balance at the interface is

(2.9)\begin{equation} \boldsymbol{n^*}\boldsymbol{\cdot}\boldsymbol{\mathcal{T}}^*+\gamma \boldsymbol{n^*}\boldsymbol{\nabla^*}\boldsymbol{\cdot}\boldsymbol{n^*}-\boldsymbol{\nabla^*_s}\boldsymbol{\gamma}=0, \end{equation}

where, the stress tensor $\boldsymbol{\mathcal{T}}^*$ is given by $\boldsymbol{\mathcal{T}}^*= -p^*\boldsymbol {I}+\mu ({\boldsymbol{\nabla^*} \boldsymbol {v^*}+\boldsymbol{\nabla^{*}} \boldsymbol{v^{*}}}^T)$. Here $\boldsymbol{\nabla^*_s}$ is the surface gradient operator and is given by $\boldsymbol{\nabla ^*}-\boldsymbol {n^*}(\boldsymbol {n^*}\boldsymbol {\cdot }\boldsymbol{\nabla ^*})$. We now turn our focus to obtaining the non-dimensional governing equations.

2.2. Non-dimensional equations

The following characteristic scales are used to obtain non-dimensionalized variables:

(2.10)\begin{equation} \left.\begin{array}{c} \displaystyle x=\dfrac{x^*}{W^*},\quad z=\dfrac{z^*}{d^*}, \quad u=\dfrac{u^*}{U}, \quad w=\epsilon\dfrac{w^*}{U}, \quad t= \dfrac{t^*}{\left(\dfrac{d^*}{\epsilon U}\right)}, \\ \displaystyle \omega=\omega^* \left(\dfrac{d^*}{\epsilon U}\right)\!, \quad p= \dfrac{p^*}{\left(\dfrac{\mu U}{d^*}\right)} \quad \mbox{and}\quad T= \dfrac{T^*-T^*_{i0}}{\Delta T^*}, \end{array}\right\} \end{equation}

where $U={\kappa }/{d^*}$ is a characteristic horizontal velocity based on the time scale for thermal diffusion, $\kappa ={\lambda }/{\rho C_p}$ is the thermal diffusivity and $\Delta T^*=T^*_{H}-T^*_{i0}$ is the temperature difference across the unperturbed fluid layer.

Observe that we render $x$ and $z$ dimensionless using $W^*$ and $d^*$, respectively. This leads to a length scale separation parameter, $\epsilon$, given by $\epsilon =d^*/W^*$, which will play a role when we consider the long-wave limit of the problem in § 4. Upon non-dimensionalizing equations (2.2)–(2.9) we get

(2.11)$$\begin{gather} \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0, \end{gather}$$

(2.12)$$\begin{gather}\epsilon Re \left(\frac{\partial u}{\partial t}+ u\frac{\partial u}{\partial x} +w\frac{\partial u}{\partial z}\right) =-\epsilon \frac{\partial p}{\partial x} + \epsilon^2\frac{\partial^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial z^2}, \end{gather}$$

(2.13)$$\begin{gather}\epsilon^3 Re \left(\frac{\partial w}{\partial t}+ u\frac{\partial w}{\partial x} +w\frac{\partial w}{\partial z}\right) =-\epsilon \frac{\partial p}{\partial z} + \epsilon^4\frac{\partial^2 w}{\partial x^2} +\epsilon^2\frac{\partial^2 w}{\partial z^2} - \epsilon \mathcal{A} \cos(\omega t) \end{gather}$$

and

(2.14)\begin{equation} \epsilon \left(\frac{\partial T}{\partial t}+ u\frac{\partial T}{\partial x} +w\frac{\partial T}{\partial z}\right) = \epsilon^2\frac{\partial^2 T}{\partial x^2} +\frac{\partial^2 T}{\partial z^2}, \end{equation}

where $Re={\rho U d^*}/{\mu }$ is the Reynolds number and $\mathcal {A}={\rho A^* \omega {^*}^{2} d{^*}^{2}}/{\mu U}$. The scaled boundary conditions at the bottom wall, $z=-1$, are

(2.15a–c)\begin{equation} u=0,\quad w=0 \quad \mbox{and}\quad T=1. \end{equation}

The non-dimensional form of the mass balance at the impenetrable interface, $z=h(x,t)$, along with the tangential and normal force balance conditions are

(2.16)$$\begin{gather} \epsilon w-\epsilon u \frac{\partial h}{\partial x} = \epsilon \frac{\partial h}{\partial t}, \end{gather}$$

(2.17)$$\begin{gather}\left(\frac{\partial u}{\partial z}+\epsilon^2 \frac{\partial w}{\partial x}\right)\left(1-\epsilon^2 \left(\frac{\partial h}{\partial x}\right)^2\right)+ 4 \epsilon^2 \frac{\partial w}{\partial z}\frac{\partial h}{\partial x}\notag\\ =-\epsilon Ma \left(\frac{\partial T}{\partial x}+\frac{\partial T}{\partial z}\frac{\partial h}{\partial x}\right)\left[1+\epsilon^2\left(\frac{\partial h}{\partial x}\right)^2\right]^{1/2} \end{gather}$$

and

(2.18)\begin{align} &-\!\epsilon p +2 \epsilon^2\left[ \frac{\partial w}{\partial z}\left(1-\epsilon^2 \left(\frac{\partial h}{\partial x}\right)^2\right)-\left(\frac{\partial u}{\partial z}+\epsilon^2\frac{\partial w}{\partial x}\right)\left(\frac{\partial h}{\partial x}\right)\right]\left[ 1+\epsilon^2\left(\frac{\partial h}{\partial x}\right)^{2}\right]^{-1}\nonumber\\ &\quad =\epsilon^3 \left(\frac{1}{Ca}-Ma T\right)\frac{\partial^2 h}{\partial x^2}\left[ 1+\epsilon^2\left(\frac{\partial h}{\partial x}\right)^2\right]^{-{3}/{2}}, \end{align}

where $Ca={\mu U}/{\gamma }$ and $Ma={\gamma _T \Delta T^*}/{\mu U}$ are the capillary and Marangoni numbers, respectively. Newton's law of cooling at the free surface becomes

(2.19)\begin{equation} Bi ( T-T_a) + \left(\dfrac{\partial T}{\partial z}-\epsilon^2\dfrac{\partial T}{\partial x}\dfrac{\partial h}{\partial x}\right)\left[ 1+\epsilon^2\left(\dfrac{\partial h}{\partial x}\right)^2\right]^{-{1}/{2}}=0, \end{equation}

where $T_a=({T^*_a-T^*_{io}})/({T^*_H-T^*_{io}})$ and where $Bi={\mathcal {H} d^*}/{\lambda }$ is the Biot number.

In what follows we shall consider the nonlinear (2.11)–(2.19) in different limits. In § 3, we discuss the stability analysis for arbitrary wavenumbers by choosing $d^*$ to be the common length scale for the horizontal and vertical directions. This amounts to setting $\epsilon =1$ without any loss of generality. By contrast, in § 4, we will retain the original scaling and consider the limit $\epsilon \ll 1$.

3.1. Pure Marangoni instability

The linear stability calculation for the pure Marangoni problem, done for the parameters given in tables 1 and 2, is depicted in figures 2(a) and 2(b) as plots of critical Marangoni number versus the wavenumber. The initial rise and fall in the curves are due to the competition between the stabilizing effect of surface tension and the destabilizing effect of surface tension gradients. This is followed by a final rise due to stabilization afforded by viscosity and thermal diffusion. Due to the finiteness of the horizontal extent of the container, the allowable wavenumbers, $k$, are given by $k=2m{\rm \pi} /W$, where $m\in \mathbb {N}$ and $W$ is the dimensionless width, i.e. $W^*/d^{*}$. For example, when $1/W=0.075$, the allowable wavenumbers are depicted by the vertical lines drawn in figure 2(a) i.e. at $k=2{\rm \pi} /W$, $4 {\rm \pi}/W$, $6 {\rm \pi}/W$, etc. Clearly, the minimum critical Marangoni number, denoted $Ma_c$, occurs at the wavenumber, $k=8 {\rm \pi}/W$, i.e. corresponding to four full waves within the container and denoted by the star in figure 2(a). However, if $1/W =0.003$ we arrive at figure 2(b) where now the minimum critical Marangoni number occurs at $k=2 {\rm \pi}/W$, i.e. corresponding to one full wave, and is lower than the short-wave critical Marangoni number depicted in figure 2(b), denoted by an asterisk.

Table 2. Values of dimensionless groups used in the calculations.

Figure 2. Critical $Ma$ versus $k$ curves, marking allowable wavenumbers $k=2m{\rm \pi} /W$ with vertical lines. The numbers between parentheses correspond to $m$. Here (a) ${1}/{W} =0.075$ and (b) ${1}/{W} =0.003$. The inset in (b) is a magnification of the low wavenumber range.

3.2. Pure Faraday instability

We now turn to the second of the two primary calculations, where the critical amplitude of shaking, $\mathcal {A}$, for the onset of Faraday instability is obtained in the absence of Marangoni instability. The Faraday instability principally arises when the parametric forcing frequency is commensurate with the system's natural frequency. The natural frequency, however, depends on $Ca$, $Re$ and on the wavenumber of the disturbance, $k$. It is determined by calculating $Im(\sigma )$ i.e. $\sigma _i$ via (3.3)–(3.5) and (3.7a–c)–(3.10) modified here by taking $Ma=0$, $\omega =0$ and $\mathcal {A}=0$. Now, depending on the frequency range of oscillation, the instability can give rise to different waveforms which are marked in the figures 3(a) and 3(b). The calculation is done by setting the frequency and then determining the minimum critical amplitude, denoted $A^*_c$, as we sweep through the allowable waveforms, say, $k=2{\rm \pi} /W$, $k=4 {\rm \pi}/W$, $k=6{\rm \pi} /W$, etc. The lowest of the minimum amplitudes and the corresponding waveforms constitute the critical diagram, figure 3, displayed in terms of dimensioned quantities to emphasize the practical attainability of the results. This figure, drawn again for the two cases, i.e. $1/W=0.075$ and $1/W=0.003$, shows the effect of sidewall proximity. From figure 3(a) we see that there are minimum points in the frequency range where resonance is strongest for each interfacial mode. These minimum points in each of these frequency ranges are near, but not precisely at, the natural frequency associated with the waveforms. The critical curve is concatenated from the individual waveforms. However, this feature becomes less discernible as the dimensionless width becomes very large. We see this behaviour in figures 3(a) and 3(b), with figure 3(b) being drawn for a very large value of the dimensionless width compared with figure 3(a).

Figure 3. Critical forcing amplitude $(A^*_c)$ versus frequency $(\,f^*)$ for the pure Faraday instability. The numbers between parentheses correspond to $m$. The vertical line corresponds to $f^{*}={1}$ Hz. Here (a) ${1}/{W} =0.075$ and (b) ${1}/{W} =0.003$.

Unlike the Marangoni instability, the Faraday instability is strongly dependent on the fluid's inertia i.e. the Reynolds number, $Re$. The inertial motion acts on the deflecting interface and just as gravity acts to stabilize the long-wave Marangoni instability (Vanhook et al. Reference Vanhook, Schatz, Swift, McCormick and Swinney1997) we might expect that inertial action on the fluid interface can stave off the Marangoni instability with parametric forcing. To see if this is so we turn to the next stability calculation which shows the effect of resonance on the Marangoni instability.

3.3. The effect of Faraday forcing on the long-wave Marangoni instability

We now consider the effect of parametric forcing on the Marangoni instability focusing first on the long-wave regime. In this setting, there are two leading instability modes. The first mode, with critical Marangoni number $Ma_s$, gives rise to predominantly monotonic temporal growth beyond criticality and produces a steady response at criticality, i.e. the dominant term in the Floquet expansion is $n=0$. The second mode, with threshold $Ma_o$, is predominantly oscillatory as a result of resonance. It is subharmonic in nature, i.e. the dominant term in the Floquet expansion is $n=1$.

We first focus on mode $s$, which connects to the classical Marangoni instability mode in the absence of forcing ($\mathcal {A}=0$). We perform stability calculations for different forcing amplitudes, $\mathcal {A}$, at fixed frequency $f^*={1}$ Hz and $1/W=0.003$. Figure 4 represents the results of these calculations. The solid curve corresponds to the pure Marangoni instability ($\mathcal {A}=0$) while the dashed line corresponds to $\mathcal {A}=0.5\mathcal {A}_{c}$ and the dot–dashed line corresponds to $\mathcal {A}=0.9\mathcal {A}_{c}$. Observe that the critical Marangoni number, $Ma_s$, increases with amplitude of shaking over a considerable range of wavenumbers that includes the maximum of the $Ma$ versus $k$ curve. The stabilization in the low wavenumber region is illustrated in figure 4(b), where the wavenumber $k=2{\rm \pi} /W$, corresponding to a single wave within the container width, is demarcated by a vertical solid line. Thus, the long-wave monotonic Marangoni instability mode can be strongly stabilized by parametric forcing.

Figure 4. Critical $Ma$ versus $k$ showing the stabilizing effect of Faraday forcing for a fixed frequency of ${1}$ Hz. Here the amplitude of shaking is increased from $\mathcal {A}=0$ (solid curve) to $\mathcal {A}=0.5\mathcal {A}_{c}$ (dashed curve) and $\mathcal {A}=0.9\mathcal {A}_{c}$ (dot–dashed curve). The corresponding critical amplitude is $A^*_c = {16.4}$ mm and is marked by the vertical line in figure 3(b) ($1/W=0.003$). Here $\boldsymbol {S}$ represents the stable region and $\boldsymbol {U}$ represents the unstable region.

To understand what happens when $\mathcal {A}>\mathcal {A}_c$ we turn to figure 5, which represents the corresponding stability diagram. While the $Ma_s$ are obtained for $k=2{\rm \pi} /W$, the most unstable possible wavenumber for the container width chosen, the $Ma_o$ are obtained for $k=40{\rm \pi} /W$, the lowest wavenumber at which oscillatory modes first appear, as seen from figure 3(b). The dark and light grey regions indicate stable and unstable regimes which are identified based on the real part of $\sigma$, for the two leading eigenvalues. Observe that as $\mathcal {A}$ increases, $Ma_{s}$ (shown by the solid dots in figure 5a) increases, which is in accordance with figure 4. The second instability mode, $Ma_o$, represented by the open circles in the magnified view shown in figure 5(b), appears in the proximity of $\mathcal {A}_c$. For $\mathcal {A}<\mathcal {A}_c$, $Ma_o$ rapidly approaches negative infinity as $\mathcal {A}$ is reduced. When $\mathcal {A}=\mathcal {A}_{c}$, $Ma_o$ reaches zero while $Ma_s$ continues to increase, but at a slower rate. As a result, when $\mathcal {A}>\mathcal {A}_c$, $Ma_o$ can overtake $Ma_s$. It is clear from figures 5(a) and 5(b) that when $Ma < Ma_o$ the flow is oscillatory but if $Ma > Ma_s$ the flow is monotonic. Clearly, if $Ma_o < Ma_s$, then the system is stable for $Ma_o < Ma < Ma_s$ and if $Ma_o > Ma_s$, then the system is always unstable.

Figure 5. Stability diagrams depicting the emergence of a second instability mode for a fixed frequency, $f^*={1}$ Hz. Critical Marangoni numbers versus forcing amplitude, $\mathcal {A}$, are shown along with the sign of the real part of $\sigma$, $\mathcal {R}(\sigma$), for the long-wave regime. The symbol $\bullet$ represents $Ma_{s}$ and the symbol $\circ$ represents $Ma_{o}$. Panel (b) is the magnification of (a) in the vicinity of $\mathcal {A}_c$. The dark region represents $\mathcal {R}(\sigma )<0$ and the light region represents $\mathcal {R}(\sigma )>0$.

We now turn to the task of analysing the effect of Faraday forcing on the short-wave regime.

3.4. Effect of Faraday forcing on the short-wave Marangoni instability

To access the short-wave Marangoni instability we consider, as an example, $1/W=0.075$ and refer to figure 2(a), where the star locates the onset of short-wave instability and corresponds to $k = 8{\rm \pi} /W$. Figure 4(a), although drawn for a wider container, suggests that the effect of forcing on the short-wave Marangoni instability is negligible. Let us verify whether this holds for the current container width, where the short-wave instability is dominant. As in the discussion of the long-wave case, figure 6 represents the corresponding stability diagram. To make a consistent comparison with the long-wave case, illustrations are shown for the same frequency, i.e. $f^* = {1}$ Hz. Unlike the long-wave case, the $Ma_s$ are obtained for $k = 8{\rm \pi} /W$ (marked by a star in figure 2a) and the $Ma_o$ are obtained for $k = 2{\rm \pi} /W$ (marked by the vertical line in figure 3a), i.e. the most unstable mode for each instability type.

Figure 6. Illustration of $Ma$ versus $\mathcal {A}$ depicting the sign of the real part of $\sigma$, $\mathcal {R}(\sigma )$ for the short-wave regime. The symbol $\bullet$ represents $Ma_{s}$ and the symbol $\circ$ represents $Ma_{o}$. The dark region represents $\mathcal {R}(\sigma )<0$ and the light region represents $\mathcal {R}(\sigma )>0$.

Observe from figures 6(a) and 6(b) that $Ma_s$ is not affected by parametric forcing, unlike the long-wave case, while $Ma_o$ displays the same behaviour as seen in figure 5. We thus focus on the effect of parametric forcing on the long-wave instability mode in the remainder of the manuscript.

3.5. A practical perspective: the effect of increasing forcing amplitude at fixed Marangoni number where $Ma>Ma_c$

In § 3.3, we saw that raising the forcing amplitude, $\mathcal {A}$, increases the threshold Marangoni number, $Ma_{s}$, in the long-wave region. This implies that an increase in $\mathcal {A}$ at fixed Marangoni number, $Ma>Ma_c$, can render the flow stable. In the current subsection, we investigate whether this holds over a wide frequency range. To this end, we consider a system with physically realizable properties given by table 1, now setting the width to $1/W=0.003$. The Marangoni number is set to $Ma = 1.4 Ma_{c}$, where $Ma_{c}$ corresponds to a waveform with $k = 2{\rm \pi} /W$. This number, $Ma\sim 25$, can be read off the inset of figure 4(b). It may be seen from the inset of figure 2(b) that for $Ma=25$, only the waveform corresponding to one wave, i.e. $m=1$ is unstable, all other allowable waveforms being stable.

Having assigned $Ma$, we now seek the stability of the flow as we increase the dimensioned forcing amplitude, $A^*$, until $\sigma _{r}=0$, while incrementing the input frequency, $f^*$. Results are displayed in figure 7, where we once again use dimensioned quantities to emphasize the practical feasibility of the calculations. In this figure, the dashed curve corresponds to the threshold amplitude, $A^*_s$, required for the stabilization of the monotonic, long-wave Marangoni mode. On this curve, the assigned $Ma$ equals $Ma_s$. A connection between figures 4 and 7 can be made by considering a frequency of ${1}$ Hz, where the critical amplitude is approximately ${15}$ mm according to figure 7. This corresponds to $\mathcal {A} \sim 0.9 \mathcal {A}_c$, i.e. to the intersection of the horizontal and vertical lines in figure 4.

Figure 7. Frequency dependence of the stability bounds for the monotonic ($s$) and oscillatory ($o$) instability modes: $1/W=0.003$. The $A^*$ versus $f^*$ curves at fixed $Ma=1.4 Ma_c\vert _{m=1}$ (marked by circle in figure 2b). Integers between parentheses, ($m$), identify the most unstable wavenumber $k=2 m{\rm \pi} /W$. The upward arrow marks the cutoff frequency, where $A^*_o$ becomes equal to $A^*_s$.

Returning to figure 7, the solid curve represents a second threshold, $A^*_o$, where the interface deformation is predominantly subharmonic in nature and reflective of resonance. It is here that the assigned $Ma=1.4 Ma_c$ equals $Ma_o$, as denoted earlier by open circles in figure 5. The integer, $m$, between parentheses in figure 7 characterizes the most-unstable wavenumber, $k = m 2{\rm \pi} /W$, for a given frequency.

The plot, thus obtained, demarcates three regions of instability/stability, which are separated by the dashed and solid curves. Region $\mathrm {I}$ is unstable to a monotonic mode related to the Marangoni instability, region $\mathrm {II}$ is fully stable and region $\mathrm {III}$ is unstable to a subharmonic oscillatory mode related to Faraday instability.

It should be observed that figure 7 is similar to figure 3(b) even though now $Ma\neq 0$. In fact, the numerical values of the critical amplitude, $A^*_o$, given by the solid lines in the two figures are almost indistinguishable from one another. This implies that the Marangoni number has little influence on the onset of resonance, even though the parametric forcing has a strong influence in quenching the Marangoni instability.

Region II exists only in a certain frequency range, until the solid and dashed curves intersect. We designate the frequency, $f^{*}$, at the intersection point of the three regions, indicated by the vertical arrow in figure 7, as the cutoff frequency. Beyond this frequency, we move from an unstable Marangoni regime (region $\mathrm {I}$) directly to an unstable resonant regime (region $\mathrm {III}$) eluding any stable transition. Within the thin region between the solid and dashed curves, both the monotonic and oscillatory instability modes grow, but the oscillatory one is dominant. For example, at a forcing frequency of ${1.56}$ Hz, the oscillatory mode with $m=24$ ought to be observed in an experiment. This means that to the right of the cutoff frequency an increase in amplitude will elicit a drastic change in the modal response from one wave to more than $21$ waves. This sudden transition is interpreted from the perspective of $A^{*}$ versus $k$ curves, in Appendix A.

Let us now consider how the cutoff frequency changes when increasing the assigned Marangoni number. A higher Marangoni number will enhance surface tension gradient flows and to suppress these flows, a higher amplitude of oscillation is required at a given frequency. Hence, the dashed line, which represents the threshold amplitude to quench Marangoni instability, will shift upwards, increasing region $\mathrm {I}$, where the Marangoni instability is dominant and decreasing region $\mathrm {II}$, where the fluid system is stable. This will cause a shift in the cutoff frequency as shown by the graph in figure 8.

Figure 8. The $Ma/Ma_c$ versus cutoff frequency.

The above observations are important as it is known that the pure Marangoni problem in the long-wave region leads to dry spots when the Marangoni number exceeds its critical value (Vanhook et al. Reference Vanhook, Schatz, Swift, McCormick and Swinney1997). To stave off such dry spots with Faraday forcing our calculations would have us hypothesize that we may either force at low frequency leading to a stable flat surface or force at high frequency at the risk of producing large amplitude interfacial oscillatory waves. To check this hypothesis, we proceed to a nonlinear model that allows us to track the interface evolution. This is an involved task but we can gain much insight by appealing to a reduced-order model that uses separation of length scales in the case of thin films, where $\epsilon = 1/W$ is much less than unity. This is done in the next section.

4.1. Long-wave approximation

As we are concerned with the long-wave Marangoni instability mode, we invoke the long-wave approximation, $\epsilon \ll 1$, where $\epsilon = d^*/W^*$ denotes the ratio of normal to transverse length scales (cf. figure 1), in order to derive a low-dimensional WRIBL model, allowing efficient, high-fidelity nonlinear computations (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2011).

We start with the non-dimensional governing equations (2.11)–(2.19), which we truncate at $O(\epsilon ^2)$. Thus, both inertia, which plays an important role in parametric instabilities enters at $O(\epsilon )$, and transverse viscous and thermal diffusion, which enter at $O(\epsilon ^2)$, are accounted for. We thus obtain the truncated continuity equation

(4.1)\begin{equation} \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z} = 0, \end{equation}

the momentum equations

(4.2)$$\begin{gather} \epsilon Re \left(\frac{\partial u}{\partial t}+ u\frac{\partial u}{\partial x} +w\frac{\partial u}{\partial z}\right) =-\epsilon \frac{\partial p}{\partial x} + \epsilon^2\frac{\partial^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial z^2}, \end{gather}$$

(4.3)$$\begin{gather}- \epsilon \frac{\partial p}{\partial z}+\epsilon^2\frac{\partial^2 w}{\partial z^2}-\epsilon \mathcal{A} \cos(\omega t)=0 \end{gather}$$

and the energy equation

(4.4)\begin{equation} \epsilon \left(\frac{\partial T}{\partial t}+ u\frac{\partial T}{\partial x} +w\frac{\partial T}{\partial z}\right) = \epsilon^2\frac{\partial^2 T}{\partial x^2} +\frac{\partial^2 T}{\partial z^2}, \end{equation}

where we have assumed $Re$ and $\mathcal {A}$ to be of at least $O(\epsilon )$. Further, at $z=h(x,t)$, the interphase coupling conditions are

(4.5a)$$\begin{gather} \epsilon\frac{\partial h}{\partial t}=-\epsilon u\frac{\partial h}{\partial x}+\epsilon w, \end{gather}$$

(4.5b)$$\begin{gather}\frac{\partial u}{\partial z}\left(1-\epsilon^2 \left(\frac{\partial h}{\partial x}\right)^2\right)+\epsilon^2 \frac{\partial w}{\partial x}+ 4 \epsilon^2 \frac{\partial w}{\partial z}\frac{\partial h}{\partial x} =-\epsilon Ma \left(\frac{\partial T}{\partial x}+\frac{\partial T}{\partial z}\frac{\partial h}{\partial x}\right)\!, \end{gather}$$

(4.5c)$$\begin{gather}-\epsilon p+2 \epsilon^2 \left(\frac{\partial w}{\partial z}-\frac{\partial u}{\partial z}\frac{\partial h}{\partial x}\right)=\frac{\epsilon^{3}}{Ca}\frac{\partial^{2} h}{\partial x^{2}} \end{gather}$$

and

(4.5d)$$\begin{gather}Bi (T-T_a)\left[ 1+\frac{\epsilon^2}{2}\left(\dfrac{\partial h}{\partial x}\right)^2\right] + \frac{\partial T}{\partial z}-\epsilon^2\dfrac{\partial T}{\partial x}\dfrac{\partial h}{\partial x}=0, \end{gather}$$

where we have assumed the capillary number, $Ca$, to be at most of $O(\epsilon )$, and the Marangoni number, $Ma$, to be at least of $O(\epsilon )$. We designate the system of truncated equations given by (4.1)–(4.5) as a $1+\epsilon ^2$ model.

Next, we eliminate the pressure, $p$, from (4.2) via (4.3) and (4.5c). Integration of (4.3) from $z$ to $z = h(x,t)$, and substitution of (4.5c) for ${p}|_{h}$ yields the pressure profile in the liquid film, i.e.

(4.6)\begin{equation} \epsilon p\vert_{z}=-\frac{\epsilon^{3}}{Ca}\frac{\partial^{2} h}{\partial x^{2}}+\epsilon \mathcal{A} \cos(\omega t) (h-z)+\epsilon^2 \left[\left(\frac{\partial w}{\partial z}-2 \frac{\partial u}{\partial z}\frac{\partial h}{\partial x}\right)\left\vert_{h}+\frac{\partial w}{\partial z}\right\vert_{z}\right]\!, \end{equation}

which is then substituted into (4.2). As a result, we obtain the so-called boundary layer equation governing momentum transport within the liquid film in the long-wave limit, which replaces (4.2) and (4.3)

(4.7)\begin{align} &\epsilon Re \left(\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+w \frac{\partial u}{\partial z}\right) = \frac{\epsilon^{3}}{Ca}\frac{\partial^{3} h}{\partial x^{3}}-\epsilon \mathcal{A} \cos(\omega t)\frac{\partial h}{\partial x}+\frac{\partial^{2} u}{\partial z^{2}}+2\epsilon^2\frac{\partial^{2} u}{\partial x^{2}}\nonumber\\ &\quad +\epsilon^2\frac{\partial}{\partial x} \left(2\frac{\partial u}{\partial z}\frac{\partial h}{\partial x}-\frac{\partial w}{\partial z}\right)\Big\vert_{h}. \end{align}

4.2. The WRIBL model

We reduce the dimension of the system of governing equations (4.1), (4.4), (4.5) and (4.7) by eliminating the $z$-coordinate from the problem and describing the flow in terms of the local instantaneous film surface deflection $h(x,t)$, flow rate $q(x,t)$ and film surface temperature $\theta (x,t)$. This is achieved by integrating (4.4) and (4.7) in the $z$-direction across the film thickness, i.e. from $z=-1$ to $z=h(x,t)$, following the WRIBL approach (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2011).

First, the dependent variables, $u$ and $T$, are decomposed into a leading-order contribution (highlighted by a hat) and an $O(\epsilon )$ correction (highlighted by a tilde), i.e.

(4.8)\begin{equation} u(x,z,t)=\underbrace{\hat{u}(x,z,t)}_{O(1)}+ \underbrace{\tilde{u}(x,z,t)}_{O(\epsilon)},\end{equation}

and

(4.9)\begin{equation} T(x,z,t)=\underbrace{\hat{T}(x,z,t)}_{O(1)}+ \underbrace{\tilde{T}(x,z,t)}_{O(\epsilon)}, \end{equation}

and $w$ follows from (4.1). For this decomposition to hold, $\hat {u}$ and $\hat {T}$ should remain sufficiently close to the exact solution at all times. This can be achieved by requiring $\hat {u}$ and $\hat {T}$ to satisfy the governing equations (4.1), (4.4), (4.5) and (4.7) truncated at $O(1)$ and by introducing appropriate gauge conditions.

In particular, $\hat {u}$ is governed by the boundary value problem

(4.10a–c)\begin{equation} \frac{\partial^{2}\hat{u}}{\partial z^{2}}=K, \quad \hat{u}\vert_{z=-1}=0, \quad \mbox{and}\quad \left.\frac{\partial \hat{u}}{\partial z}\right\vert_{z=h(x,t)}=-\epsilon Ma \left(\frac{\partial \hat{T}}{\partial x}+\frac{\partial \hat{T}}{\partial z}\frac{\partial h}{\partial x}\right)\!, \end{equation}

where we have further restricted $Ma$ to be at least of $O(1/\epsilon )$. We have additionally restricted $Ca$ to be at most of $O(\epsilon ^{3})$, in order to retain the inhomogeneity, $K$. Further, we impose the gauge condition,

(4.11)\begin{equation} \int_{-1}^{h(x,t)} \hat{u} \,{\rm d} z=q(x,t), \end{equation}

which requires that the base velocity profile yield the local instantaneous flow rate $q(x,t)$. Solving (4.10a–c) for $\hat {u}$ subject to (4.11), yields

(4.12)\begin{equation} \hat{u}=-\frac{3 q (z+1-2 (h+1))(z+1)}{2 (h+1)^3}-\epsilon Ma \frac{\partial\theta}{\partial x} \frac{(z+1) (3(z+1)-2 (h+1))}{4 (h+1)}, \end{equation}

which implies for the velocity correction, $\tilde {u}$,

(4.13a–c)\begin{equation} \int_{-1}^{h(x,t)} \tilde{u} \,{\rm d} z=0,\quad \tilde{u}\vert_{z=-1}=0 \quad\text{and} \quad \left.\frac{\partial \tilde{u}}{\partial z}\right\vert_{z=h}=-\left(\epsilon^2 \frac{\partial \hat{w}}{\partial x}+ 4 \epsilon^2 \frac{\partial \hat{w}}{\partial z}\frac{\partial h}{\partial x}\right)\!. \end{equation}

Similarly, $\hat {T}$ is governed by the boundary value problem

(4.14a–c)\begin{equation} \frac{\partial^{2} \hat{T}}{\partial z^{2}}=0,\quad \hat{T}\vert_{z=-1}=1 \quad \mbox{and}\quad \hat{T}\vert_{z=h(x,t)}=\theta, \end{equation}

where $\theta (x,t)$ is the temperature of the film surface. Solving (4.14a–c) for $\hat {T}$, we obtain

(4.15)\begin{equation} \hat{T}=\frac{ (\theta-1)}{h+1}(z+1)+1, \end{equation}

which implies for the temperature correction, $\tilde {T}$,

(4.16a)\begin{equation} \tilde{T}\vert_{z=-1}=0, \quad \tilde{T}\vert_{z=h}=0 \end{equation}

and

(4.16b)\begin{equation} \frac{\partial \tilde{T}}{\partial z}\vert_{z=h}=-\frac{h+\theta }{h+1}-Bi \theta+\epsilon ^2 \left[\frac{\partial h}{\partial x} \frac{\partial \theta}{\partial x}-\left(\frac{\partial h}{\partial x}\right)^2\left(\frac{ Bi \theta +1}{2 }+\frac{ \theta -1}{h+1}\right)\right]\!. \end{equation}

Next, we introduce (4.12) and (4.15) in (4.4) and (4.7), truncate the resulting equations at $O(\epsilon ^2)$, multiply with the weight functions $F(x,z,t)$ and $\mathcal {W}(x,z,t)$, and integrate across the film height, $1+h(x,t)$. This yields the integral momentum equation

(4.17a)\begin{align} &\int_{-1}^{h} \epsilon Re \left(\frac{\partial \hat{u}}{\partial t}+\hat{u} \frac{\partial \hat{u}}{\partial x}+\hat{w} \frac{\partial \hat{u}}{\partial z}\right) F\,{\rm d} z= \boxed{ \int_{-1}^{h} \frac{\partial^{2} }{\partial z^{2}}(\hat u+\tilde u) F\,{\rm d} z}+\int_{-1}^{h}2\epsilon^2\frac{\partial^{2} \hat{u}}{\partial x^{2}} F\,{\rm d} z\nonumber\\ &\quad - \left(\epsilon \mathcal{A} \cos(\omega t)\frac{\partial h}{\partial x}-\frac{\epsilon^{3}}{Ca}\frac{\partial^{3} h}{\partial x^{3}} -\epsilon^2\frac{\partial}{\partial x} \left[ (2\frac{\partial \hat{u}}{\partial z}\frac{\partial h}{\partial x}-\frac{\partial \hat{w}}{\partial z})\vert_{h}\right]\right) \int_{-1}^{h} F\,{\rm d} z, \end{align}

and the integral energy equation

(4.17b)\begin{align} \int_{-1}^{h} \epsilon \left(\frac{\partial \hat{T}}{\partial t}+ \hat{u}\frac{\partial \hat{T}}{\partial x} +\hat{w}\frac{\partial \hat{T}}{\partial z}\right) \mathcal{W} \,{\rm d} z= \boxed{\int_{-1}^{h} \frac{\partial^2 }{\partial z^2}(\hat T+ \tilde T) \mathcal{W} \,{\rm d} z}+\int_{-1}^{h} \epsilon^2\frac{\partial^2 \hat T}{\partial x^2} \mathcal{W} \,{\rm d} z, \end{align}

where the unknown corrections, $\tilde {u}$ and $\tilde {T}$, only appear in the boxed terms, expressing wall-normal diffusion. At most other places, $\tilde {u}$ and $\tilde {T}$ drop out due to truncation at $O(\epsilon ^2)$. Only for the convective terms have we made additional assumptions by dropping the contributions of the corrections, even though $O(\epsilon \tilde {u}) = O(\epsilon \tilde {T}) = O(\epsilon ^2)$, implying weak (but finite) inertial and heat convection effects.

The weight functions, $F$ and $\mathcal {W}$, are chosen such as to cancel $\tilde {u}$ and $\tilde {T}$ from the boxed terms in (4.17a) and (4.17b). For $F$, we require

(4.18a–c)\begin{equation} \frac{\partial^{2}F}{\partial z^{2}}=1,\quad F\vert_{z=-1}=0,\quad \left.\frac{\partial F}{\partial z}\right\vert_{z=h(x,t)}=0, \end{equation}

yielding

(4.19)\begin{equation} F=\tfrac{1}{2} (z+1) (-2 h+z-1), \end{equation}

and for $\mathcal {W}$, we require

(4.20a–c)\begin{equation} \frac{\partial^{2} \mathcal{W}}{\partial z^{2}}=0, \quad \mathcal{W}\vert_{z=-1}=0 \quad \mbox{and} \quad \left. \frac{\partial \mathcal{W}}{\partial z}\right\vert_{z=h(x,t)}=\mathcal{C}, \end{equation}

yielding

(4.21)\begin{equation} \mathcal{W}=\mathcal{C}(z+1), \end{equation}

where $\mathcal {C}$ is an arbitrary constant. Because the weight functions, according to (4.19) and (4.21), correspond to the leading-order profiles $\hat {u}$ (4.12) and $\hat {T}$ (4.15), the weighted integration we have performed corresponds to a Galerkin projection.

As a result of these definitions, the boxed terms in (4.17a) and (4.17b) can be rewritten using integration by parts

(4.22)\begin{equation} \int_{-1}^{h} \frac{\partial^{2} }{\partial z^{2}}(\hat u+\tilde u) F\,{\rm d} z=\left[\frac{\partial (\hat u+\tilde u) }{\partial z} F-(\hat u+\tilde u) \frac{\partial F}{\partial z} \right]_{-1}^{h}+\int_{-1}^{h} (\hat{u}+\tilde{u}) \,{\rm d} z, \end{equation}

and

(4.23)\begin{equation} \int_{-1}^{h} \frac{\partial^2 }{\partial z^2}(\hat T+ \tilde T) \mathcal{W} \,{\rm d} z= \left[\frac{\partial(\hat T+\tilde T) }{\partial z} \mathcal{W}- (\hat T+\tilde T) \frac{\partial \mathcal{W}}{\partial z} \right]_{-1}^{h}. \end{equation}

Then, substituting (4.12), (4.13a–c), (4.15), (4.16), (4.19) and (4.21) into (4.17a) and (4.17b), we obtain the final integral momentum and energy equations of our WRIBL model

(4.24a)\begin{align} &\epsilon Re\left\{\frac{18 q^2}{35}\frac{\partial h}{\partial x}-\frac{2\hbar^2}{5} \frac{\partial q}{\partial t} -\frac{34\hbar q}{35} \frac{\partial q}{\partial x}\right\}+\epsilon^2 {Re} Ma\left\{\frac{\hbar^2 q}{56}\left( \frac{3}{2}\hbar\frac{\partial ^2\theta }{\partial x^2} +\frac{\partial h}{\partial x} \frac{\partial \theta }{\partial x} \right)\right.\nonumber\\ &\qquad +\left.\frac{\hbar^3}{120}\left(\hbar \frac{\partial ^2\theta }{\partial x\partial t} +\frac{19}{7}\frac{\partial \theta }{\partial x} \frac{\partial q}{\partial x} \right)\right\} -\epsilon^3{Re}Ma^2\frac{\hbar^4 }{336}\frac{\partial \theta}{\partial x}\left\{\frac{19}{5}\frac{\partial h}{\partial x}\frac{\partial \theta}{\partial x}+\frac{5}{2}\hbar\frac{\partial^2 \theta}{\partial x^2}\right\}\nonumber\\ &\quad =q+\frac{\hbar^3}{3}\left\{\epsilon A \frac{\partial h}{\partial x} \cos (\omega t)-\frac{\epsilon ^3}{Ca}\frac{\partial ^3h}{\partial x^3}\right\}+\epsilon Ma\frac{\hbar^2}{2}\frac{\partial \theta }{\partial x}\nonumber\\ &\qquad +\epsilon ^2\left\{\frac{9 \hbar}{5} \left(\frac{\partial h}{\partial x} \frac{\partial q}{\partial x} - \frac{\partial ^2q}{\partial x^2} \hbar \right)+\frac{4 q}{5}\left(3\hbar \frac{\partial ^2h}{\partial x^2} -2 \left(\frac{\partial h}{\partial x}\right)^2\right)\right\}\nonumber\\ &\qquad +\epsilon^3Ma \hbar^2\left\{\frac{\hbar}{10}(3 \frac{\partial ^2\theta }{\partial x^2} \frac{\partial h}{\partial x} + \frac{\partial ^3\theta }{\partial x^3} \hbar)+\frac{1}{15} \frac{\partial ^2h}{\partial x^2} \frac{\partial \theta }{\partial x} \hbar-\frac{4}{5}\frac{\partial \theta }{\partial x} \left(\frac{\partial h}{\partial x}\right)^2 \right\}, \end{align}

and

(4.24b)\begin{align} &\epsilon \left\{ \frac{1}{3} \hbar^2 \frac{\partial \theta}{\partial t}+\frac{7}{120} \hbar \frac{\partial q}{\partial x} (\theta-1)+\frac{9}{20} \hbar q \frac{\partial \theta}{\partial x} \right\}\nonumber\\ &\quad +\frac{\epsilon^2 Ma \hbar^2}{40}\left[(1-\theta)\left(\frac{1}{2} \frac{\partial \theta ^2}{\partial x^2} \hbar + \frac{\partial h}{\partial x} \frac{\partial \theta }{\partial x}\right) - \hbar\left(\frac{\partial \theta }{\partial x}\right)^2 \right]=1-\theta -\hbar(Bi \theta +1)\nonumber\\ &\quad -\epsilon^2\hbar\left\{\frac{(Bi \theta+1)}{2} \left(\frac{\partial h}{\partial x}\right)^2-\frac{1}{3} \frac{\partial ^2\theta }{\partial x^2} \hbar + \frac{(\theta-1)}{3}\left(\frac{\partial ^2h}{\partial x^2} +\frac{1}{\hbar}\left(\frac{\partial h}{\partial x}\right)^2\right)-\frac{1}{3} \frac{\partial h}{\partial x} \frac{\partial \theta }{\partial x} \right\}, \end{align}

where we have substituted $\hbar = 1+h$. These equations are completed by the integral continuity equation

(4.24c)\begin{equation} \frac{\partial h}{\partial t}+\frac{\partial q}{\partial x}=0, \end{equation}

obtained by integrating (4.1) from $z=-1$ to $z=h(x,t)$, and invoking the interphase mass balance (4.5a). Thus, our final model (4.24) consists of a nonlinear system of three evolution equations in terms of the three unknowns $h(x,t)$, $q(x,t)$ and $\theta (x,t)$.

Observe that (4.24a) contains a mixed derivative, $\partial _{xt}\theta$, at one place. The latter appears as a result of the base velocity profile $\hat {u}$ chosen in (4.12), and, in particular, as a result of retaining the leading-order Marangoni stress in the boundary condition (4.10a–c). We find that this choice greatly improves the linear stability predictions of our model. Neither in the linear stability calculations nor in the nonlinear computations (where we use a spectral method for spatial discretization) do the mixed derivative pose any numerical challenge.

4.3. Numerical procedure

Numerically, to solve the system of nonlinear equations, (4.24), we discretize variables in space along the horizontal direction using a Fourier spectral scheme (Hesthaven, Gottlieb & Gottlieb Reference Hesthaven, Gottlieb and Gottlieb2007). This automatically satisfies periodicity in the $x$-direction, which is the setting considered in the current manuscript. In all our nonlinear computations, the domain length is set to the container width, $W$. Time integration is realized via the NDSolve function in Mathematica (version 12.1). Each of our nonlinear computations is started from an initial perturbation constructed with the eigenfunctions ${h'}|_{t=0}$, ${q'}|_{t=0}$, and ${\theta '}|_{t=0}$, defined according to (3.2), and obtained from a preliminary linear stability calculation for $\mathcal {A}=0$ and $Ma>Ma_{c}$. The initial deflection amplitude was set to be $5\,\%$ of the fluid depth. The grid resolution for our simulations was determined via grid independence tests by a stepwise doubling of the number of collocation points until the change in the minimum interface deflection, $h_{min}(t)$, was less than 1 $\%$.

4.4. Nonlinear suppression of dry spots

To validate our long-wave model (4.24), we start by reproducing the neutral stability bounds from figure 7. For this, we linearly perturb the dependent variables $h(x,t)$, $q(x,t)$ and $\theta (x,t)$ according to (3.2), and solve the resulting eigenvalue problem to obtain the critical forcing amplitude, $\mathcal {A}$, at fixed $\omega$, $Ma$ and $\epsilon =0.003$, using the parameters in table 1. Results are represented in figure 9, where, as in figure 7, the dashed curve represents the monotonic instability mode ($\hat {h}_0 \gg \hat {h}_n$ for all $n \neq 0$) and the solid curve represents the oscillatory instability mode ($\hat {h}_n \gg \hat {h}_0$ for at least one $n \neq 0$). Agreement with the corresponding curves in figure 7 is excellent, due to the smallness of the chosen $\epsilon =0.003$, which warrants the long-wave approximation underlying our model.

Figure 9. Linear effect of parametric forcing on Marangoni instability, as predicted by the long-wave model (4.24). Critical forcing amplitude $A^*$ versus forcing frequency $f^*$. The same parameters as in figure 7: $d^*=3$ mm; $\epsilon = d^*/W^*=0.003$; $Ma=1.4 Ma_c|_{k=2{\rm \pi} }$. Dashed curve, stability bound of monotonic mode ($\hat {h}_0\gg \hat {h}_n$ for all $n \neq 0$); solid curve, stability bound for oscillatory modes ($\hat {h}_n\gg \hat {h}_0$ for at least one $n \neq 0$). Numbers ($m$) designate wavenumber $k = 2 m{\rm \pi}$ of critical spatial mode. For the chosen parameters, only $m=1$ yields an unstable Marangoni mode at $A^*=0$.

We next perform nonlinear computations based on (4.24) for parameters lying in the three characteristic regions of figure 9. In region I, the liquid film is always unstable to a Marangoni mode with $m=1$, i.e. $k = 2{\rm \pi}$. Recall that, in this section, the wavenumber, $k$, is scaled by the container width, $W^{*}$, i.e. $k = k^{*} W^{*}$. In region II, the film is rendered linearly stable via parametric forcing. In region III, the film is always unstable to Faraday waves.

Figure 10 represents the time evolution of the film thickness profile, $\hbar =1+h(x,t)$, in the limit $A^*=0$, where only the Marangoni instability is active. The value of the Marangoni number is $Ma=1.4 Ma_{c}\vert _{k=2{\rm \pi} }=25$. Here $k = 2{\rm \pi}$ is the only possible instability mode (cf. figure 2b), as a result of our choice for the container width, i.e. $\epsilon =0.003$. It is observed that the deflection of the film surface increases as time progresses, until the liquid–gas interface at the trough reaches the lower wall. The film surface buckles, forming two secondary troughs, which drastically slow down further draining of liquid from the trough region. As shown by Boos & Thess (Reference Boos and Thess1999) and Dietze, Picardo & Narayanan (Reference Dietze, Picardo and Narayanan2018), additional buckling events may follow (at higher $Ma$), but, eventually, the liquid film will always rupture due to spinodal dewetting (Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009). Thus, it is safe to say that the nonlinear evolution in figure 10, will lead to the formation of a dry spot, in accordance with the experimental observations of Vanhook et al. (Reference Vanhook, Schatz, Swift, McCormick and Swinney1997) for very thin liquid films. Supplementary movie 1 shows the evolution towards dry spot formation in action.

Figure 10. Nonlinear computation for parameters according to region I in figure 9: $A^*=0$; $Ma=25$. Interface profile of dry spot formation at different times $t$. Dashed curve, $t=0$; dashed-dotted curve, $t=7$; solid curve, $t=14$. See also Supplementary movie 1, available at https://doi.org/10.1017/jfm.2024.58.

We now consider how the nonlinear evolution in figure 10 is altered by applying additional parametric forcing with increasing forcing amplitude, $A^*$. For this, we consider two forcing frequencies, i.e. $f^*={1}$ Hz, which corresponds to the natural frequency at $m = 20$, and $f^*={1.56}$ Hz, which corresponds to the natural frequency at $m = 24$. We start by discussing our computations for $f^*={1}$ Hz which are represented in figure 11. At $A^* = 0.9 A^*_c$ (figure 11a), parameters lie in region II of figure 9, and our nonlinear computation indicates a full attenuation of the initial interface deflection. Thus, parametric forcing indeed allows us to obtain a perfectly flat film surface in this parameter range. Supplementary movie 2 shows the evolution towards flatness in action.

Figure 11. Nonlinear effect of parametric forcing on Marangoni instability. Numerical computations with WRIBL model (4.24), showing transition from regions II to III in figure 9: $f^*={1}$ Hz; $Ma=1.4 {Ma_c}|_{k=2{\rm \pi} }$. Film surface profiles at different times $t$: dashed curves, $t=0$. (a) Linear suppression of surface deformations: $A^*=0.9 A^*_c$ (region $\mathrm {II}$ in figure 9). Dash–dotted curve, $t=45$; solid curve, $t=90$. (b) Saturated Faraday waves: $A^*=1.1 A^*_{c}$ (region $\mathrm {III}$ of figure 9). Dashed-dotted curve, $t=0.001$; solid curve, $t=0.002$. See also Supplementary movies 2 and 3. Here $A^*_c$ corresponds to solid curve in figure 9.

When the forcing amplitude is increased to $A^* = 1.1 A^*_c$, which lies in region III of figure 9, Faraday waves appear. Here, the resonant action takes over and the interface develops $m=20$ waves per container width, as shown in figure 11(b). Although dry spot formation is still prevented under these conditions, we see that large interface deflections result from the parametric forcing. Supplementary movie 3 shows the formation of resonant waves in action.

A different picture emerges for the higher forcing frequency, $f^*=1.56$ Hz, as shown in figure 12. It is known from figure 9 that linear stabilization of the liquid film cannot be achieved under these conditions, as the solid curve for the onset of oscillatory instability lies below the dashed curve demarcating stabilization of the monotonic instability mode. Nonetheless, figure 12(a) shows evidence of a very strong nonlinear saturation of the interface deflection at $A^* = 0.9 A^*_c$, even though this point lies in region I of figure 9. Here, the interface remains almost flat, notwithstanding that it is subject to the long-wave Marangoni instability mode, responsible for forming dry spots in the absence of forcing. Thus, dry spot formation can be prevented and the film can be rendered virtually flat also at higher forcing frequencies, albeit due to a nonlinear mechanism (figure 12a), in contrast to the linear mechanism observed at low frequencies (figure 11a).

Figure 12. Nonlinear effect of parametric forcing on Marangoni instability. Numerical computations with WRIBL model (4.24), showing transition from regions I to III in figure 9: $f^*={1.56}$ Hz; $Ma=1.4{Ma_c}|_{k=2{\rm \pi} }$. Film surface profiles at different times $t$: dashed curves, $t=0$. (a) Nonlinear attenuation of surface deformations: $A^*=0.9 A^*_c$ (region $\mathrm {I}$ in figure 9). Dash–dotted curve, $t=2$; solid curve, $t=3.5$. (b) Saturated Faraday waves: $A^*=1.1 A^*_{c}$ (region $\mathrm {III}$ in figure 9). Dashed–dotted curve, $t=0.001$; solid curve, $t=0.002$. See also Supplementary movies 4 and 5.

This nonlinear saturation mechanism is associated with small standing waves or ripples forming on the main humps of the interface deformation. These are only slightly visible in figure 12(a), but their oscillatory motion can be distinguished clearly in the Supplementary movie 4. These standing waves or ripples form at the main humps or crests rather than the main troughs on account of the inertial action promoted by the height of fluid below the large deflecting crests. We believe that these standing waves are responsible for arresting further drainage of liquid from the deformation trough, as a result of a cellular flow pattern that develops underneath the deformation hump. This is seen in figure 13, which displays streamlines in the laboratory reference frame for the two cases from figures 10 and 12(a).

Figure 13. Change in flow structure within the liquid film. Streamlines in the laboratory reference frame for parameters corresponding to figure 10 and figure 12(a). (a) Dry spot formation (figure 10): $A^*=0.1 A^*_c$; $Ma=1.4 Ma_c|_{k=2{\rm \pi} }$; and $f={1.56}$ Hz. (b) Nonlinear saturation due to parametric forcing (figure 12a): $A^*=0.9 A^*_c$; $Ma=1.4 Ma_c|_{k=2{\rm \pi} }$; and $f^*={1.56}$ Hz. Only one half of the domain is represented in each panel, for convenience.

Finally, we turn to figure 12(b), which shows that the nonlinear saturation of the surface deformation is lost, when the forcing amplitude is increased to $A^* = 1.1 A^*_c$, where large-amplitude short Faraday waves form once again, as a result of oscillatory instability (see also Supplementary movie 5).

From an application point of view, one may ask at what forcing amplitude, $A^*$, does the nonlinear saturation of the interface deflection set in? We have quantified this for the current parameters in figure 14, which represents time traces of the minimum film thickness, $\hbar _{min}=1+h_{min}$, for different subcritical (figure 14a) and supercritical (figure 14b) forcing amplitudes $A^*$, all other parameters remaining the same as in figure 12. Based on figure 14(a), we may conclude that a significant nonlinear attenuation of the interface deflection is achieved for $A^*/A^*_c \ge 0.7$.

Figure 14. Effect of relative forcing amplitude $A' = A^*/A^*_c$ on long-time nonlinear evolution of the heated liquid film. Parameters according to figure 9: $f^*={1.56}$ Hz. Time traces of the minimum film thickness $\hbar _{min}=1+h_{min}$. (a) Nonlinear attenuation of surface deformation in region I of figure 9: $A' < 1$. From bottom to top: $A'=0$, 0.5, 0.6, 0.7, 0.8 and 0.9. (b) Saturated Faraday waves: $A' \ge 1$. From bottom to top: $A'=1.1$ and 1.