4.1. Numerical procedure

To precondition the eigenvalue problem equations (3.7) for a more efficient numerical solution, we substitute the expression of the perturbed longitudinal velocity component $u$ from (3.7a ) into the rest of the equations of the eigenvalue problem (3.7). We also extract the expression for the interfacial deformation $\zeta$ from (3.7j ) and substitute it into the rest of the boundary conditions at $z=1$ . Hence, our eigenvalue problem is formulated in terms of $w, \theta, \varphi$ and it is treated numerically in this format. We note here that in the case of the non-deformable interface considered in § 4.4, the pressure will be eliminated without substituting its gradients into the normal stress boundary condition, since the latter is redundant, and $\zeta =0$ is imposed.

The numerical solution of the linear eigenvalue problem given by (3.7), or any of its simplified versions discussed in what follows, is obtained using the open-source code developed by Pearce et al. (Reference Pearce, Heil, Jensen, Jones and Prokop2018). It is based on the shooting method and numerical construction of the complex-valued Evans function (Evans & Feroe Reference Evans and Feroe1977) whose roots represent the eigenvalues of the given problem. The solution is initiated by specifying the range where the roots of the Evans function are searched and adjusted if needed.

Alternatively, it is possible to locate the eigenvalues graphically by sampling the Evans function on the complex plane of $\lambda$ , and by doing so, one obtains contour plots of the Evans function for a specified range of $(\lambda _r,\lambda _i)$ , where $\lambda _r \equiv \Re (\lambda )$ and $\lambda _i \equiv \Im (\lambda )$ . Then, the intersection points of the contour lines with the reference point of the plane or with the axis $\lambda _i=0$ are traced. These intersection points represent the eigenvalues of the problem among which one finds the eigenvalue relevant for the onset of either monotonic or oscillatory instability of the system, respectively. However, this approach is computationally costly, and the available routine to locate roots of the Evans function is preferable and is therefore used in our investigation.

The shooting method is implemented by solving the differential equations of the eigenvalue problem (3.7) or any of their simplified versions discussed below as a boundary-value problem using the compound matrix method (Ng & Reid Reference Ng and Reid1979) to construct the Evans function by satisfying the boundary conditions separately at each of the two endpoints $z=0$ and $z=1$ , and then matching the solutions at the interim point either in the middle of the domain or elsewhere.

To verify the convergence of the computation procedure, we check the accuracy of the results by changing the number of decimal digits used in the computation and by varying the matching point where the Evans function is created. We find that the values obtained for the critical values of the Marangoni number for various parameter sets of the problem based on three different matching points within the layer domain $ 0 \leqslant z \leqslant 1$ , namely $z=1/3,\,z=1/2$ and $z=2/3$ , converge up to the third digit.

In what follows, we concentrate on physical systems in which the solvent is water. Therefore, the corresponding value of the Prandtl number is taken in our investigation as $P=7$ . The value of the Lewis number used in most cases discussed below is $L=10^{-3}$ , which is quite representative for a moderately dense nanosuspension mixture. For reference, in the case of nanoparticles of diameter $d_p^*\sim 10$ nm, the Lewis number estimated using the Stokes–Einstein formula (2.4), is $L = 4.22\times 10^{-4}$ when the solvent temperature is $T^*=300\,$ K. Furthermore, we consider the range of the values for the Soret coefficient $\eta \in (0.31,3.1)$ corresponding to the temperature drop of $\Delta T^* \in (1,10)$ K across the layer. In addition, we determine the dimensionless parameters related to nanoparticle density via

\begin{align} \rho _{np} = \frac {\rho _{np}^*}{\rho _{bf}^*}\approx 4\quad \text{and}\quad (\rho c)_{np} = \frac {\rho _{np}^*c_{np}^*}{\rho _{bf}^*c_{bf}^*}\approx 0.68, \end{align}

which corresponds to alumina nanoparticles in water, also taking $a=7.47$ in this investigation, except for § 4.4 where the parameter $a$ is changing.

Before we present the results of the linear stability analysis, we note that the equilibrium surface tension $\sigma ^\ast _r$ ranges for aqueous solutions from ${\sim} 10^{-3} \ \textrm {to}\ 30 \times 10^{-3}$ Nm⁻¹ (De Wit, Gallez & Christov Reference De Wit, Gallez and Christov1994); therefore, the dimensionless surface tension for a layer of density $\rho ^\ast _{bf}=10^{3}\,\text{kg m}^{-3}$ , viscosity $\mu ^\ast _{bf}=10^{-3}\,\text{kg (m} \cdot \text{s)}^{-1}$ , thermal diffusivity $\kappa ^\ast _{bf}=0.146\times 10^{-6}\,\text{m}^2 \, \text{s}^{-1}$ and thickness $h_0^\ast =10^{-4}\,\text{m}$ is approximately $\Sigma _0=2\times 10^{4}$ . Similarly, the modified Galileo number for the terrestrial gravity of $g^\ast =9.81\,\text{m s}^{-2}$ yields $G=67.2$ for $h^\ast _0=10^{-4}$ m and $G \sim 10^{-4}$ for $h^\ast _0=10^{-6}$ m.  A list of typical values of the thermophysical properties of the nanofluid and the dimensionless numbers is given in table 1.

Table 1. Parameter nomenclature and their typical values used in this investigation.

4.2. Cooling at the substrate $(\mathcal{Q} = -1)$

In this section, we study the onset of the thermosolutocapillary instability when the nanofluid layer is cooled at the substrate. Studies of this configuration for pure fluids are quite rare in the literature (Deissler & Oron Reference Deissler and Oron1992; Oron & Rosenau Reference Oron and Rosenau1992; Burgess et al. Reference Burgess, Juel, McCormick, Swift and Swinney2001; Alexeev & Oron Reference Alexeev and Oron2007). Recently, Sarma & Mondal (Reference Sarma and Mondal2021b ) studied thermosolutal Marangoni instability in a layer of a viscoelastic binary fluid in the presence of the Soret effect. They also presented some results for the limiting case of a Newtonian fluid associated with a zero Deborah number. However, due to the fact they used a linearised version for the Soret component of the mass flux along with constant thermophysical properties of the fluid in contradistiction to the particle concentration-dependent Soret coefficient and other thermophysical properties of the nanofluid used in this paper, any quantitative comparison between the results is impossible. However, a qualitative comparison is possible and it will be made at the end of this section.

First, it must be clear that in the configuration at hand, thermocapillarity is stabilising for $M_T\gt 0$ , because any small elevation (depression) at the layer interface will have a higher (lower) temperature than in its vicinity along the interface, and thermocapillary interfacial tractions will lead to flattening of the interface. Figure 2(a) shows a stably stratified concentration profile and this configuration suggests a possibility of the emergence of the solutocapillary instability. Indeed, we find that solutocapillarity destabilises the system when the layer is cooled at the substrate and figure 4(a) displays the neutral curves for the pure monotonic solutocapillary instability for varying Biot numbers $B$ . We note that the solutal Marangoni number $M_S$ exhibits non-monotonic variation with the wavenumber $k$ . We also infer that the variation of the Biot number does not significantly affect the critical value of the solutal Marangoni number. However, we note that the neutral curves display the emergence of two local minima, one of them is associated with a low value of the wavenumber $k$ , whereas the second one belongs to the finite range of $k$ . Both of them are clearly seen in the neutral curve for $B=0.01$ in figure 4(a).

Figure 4. Monotonic solutocapillary instability in the case of cooling at the substrate $\mathcal{Q}=-1$ at $\Phi =0.01, a=7.47, \eta =0.31, L=10^{-3}, \Sigma _0\approx 2\times 10^{4}$ and $G=6.71$ . (a) Neutral curves $M_S(k)$ for the pure solutocapillarity instability, $M_T = 0$ , versus the wavenumber $k$ for various Biot numbers $B$ . (b) Neutral curves $M_S(k)$ for various values of $M_T$ with $B=0.01$ . The rise of the neutral curves with an increase in $M_T$ illustrates a stabilising effect of thermocapillarity on the monotonic solutocapillary instability. In both panels, the symbols $U$ and $S$ denote the unstable and stable domains of the system, respectively.

In addition, we note that when the layer is cooled at the substrate, the temperature base state increases with height $z$ , as shown in figure 2(b). This implies that at the interfacial hump, the temperature is higher than that at its depression. Thus, the emerging thermocapillary shear stress will lead to fluid flow from the hump into the depression, causing the tendency of interface flattening and therefore indicating a tendency to stabilisation (Deissler & Oron Reference Deissler and Oron1992; Alexeev & Oron Reference Alexeev and Oron2007). Figure 4(b) illustrates a significant stabilisation effect imparted by thermocapillarity on the solutocapillary instability threshold. In the range of $M_T$ presented in this figure, the critical value of $M_S$ slightly increases with $M_T$ ; however, the finite-wave modes are strongly stabilised with an increase in $M_T$ .

Figure 5(a) shows that the threshold value of the solutal Marangoni number $M_S$ in the pure solutocapillary instability increases linearly with the modified Galileo number $G$ . However, as we already mentioned, the neutral curves $M_S(k)$ exhibit two minima: one of them, $M_S=m_f$ , is in the finite-wave domain, whereas the second one, $M_S=m_l$ , belongs to the long-wave domain. The relative importance of these two minima with respect to the onset of the purely solutocapillary instability is now under scrutiny. The inset of figure 5(a) shows that in the domain of intermediate values of the modified Galileo number $14.8\leqslant G\leqslant 66.6$ , the long-wave minimum $M_S=m_l$ corresponding to the deformational mode prevails over the finite-wave minimum $M_S=m_f$ , namely $m_l\lt m_f$ , and corresponds to the instability onset. In addition, we find that in the rest of the $G$ range, i.e. for relatively low and sufficiently high values of $G$ , $m_l\gt m_f$ , and the finite-wave solutocapillary instability emerges. VanHook et al. (Reference VanHook, Schatz, Swift, McCormick and Swinney1997) and Golovin, Nepomnyashchy & Pismen (Reference Golovin, Nepomnyashchy and Pismen1997) observed an analogous effect of the competition between the long-wave deformational and finite-wave pattern forming modes for the thermocapillary instability in a layer with a deformable interface. Figure 5(b) displays the variation of the critical wavenumber $k_c$ with $G$ . We note that the critical wavenumber $k_c$ for the finite-wave minimum $M_T=m_f$ exhibits two distinct domains of variation with $G$ . For moderate values of the modified Galileo number $G$ , $k_c$ varies slowly; however, for sufficiently large values of $G$ , the critical wave number $k_c$ sharply increases. The inset of figure 5(b) shows a decrease in the critical wavenumber $k_c$ of the long-wave instability corresponding to the minimum $M_S=m_l$ with $G$ .

Figure 5. Onset of a purely solutocapillary monotonic instability in the case of cooling at the substrate $\mathcal{Q}=-1$ at $\Phi =0.01,\,a=7.47,\,\eta =0.31,\,L=10^{-3},\,M_T = 0,\,\Sigma _0\approx 2\times 10^{4}$ and $B = 0.01$ . (a) Variation of the critical value of the solutal Marangoni number $M_S$ versus the modified Galileo number $G$ that shows both the finite-wave range minimum $M_S=m_f$ and the long-wave minimum $M_S=m_l$ ; the inset presents the variation of the difference $m_l-m_f$ versus $G$ . FW and LW stand for domains of finite-wave and long-wave instability, respectively. (b) Variation of the critical wavenumber $k_c$ with $G$ ; the inset shows a decrease in the critical wavenumber $k_c$ corresponding to the long-wave minimum $M_S=m_l$ with $G$ . The symbols $U$ and $S$ denote the unstable and stable domains of the system, respectively.

Figure 6(a) shows that in the case of a purely solutocapillary instability, the critical solutal Marangoni number $M_S$ monotonically increases with the dimensionless surface tension number $\Sigma _0$ when the rest of the parameters are fixed and saturates at a very high value of $\Sigma _0$ that corresponds to the non-deformable interface. Slightly above it, a dashed curve corresponding to the secondary minimum of the corresponding neutral curve is located. The critical value of $M_S$ is attained at the local minimum of the neutral curves corresponding to $M_S=m_f$ and varies weakly with $\Sigma _0$ . The variation of the critical wavenumber $k_c$ is non-monotonic with $\Sigma _0$ in the finite but relatively low wavenumber range, as shown in figure 6(b), so the instability is finite-wave. The wavenumber corresponding to the secondary local minimum of the neutral curves $M_S(k)=m_l$ remains almost unchanging with $\Sigma _0$ and small, $k \approx 2.0 \times 10^{-3}$ ; therefore, this secondary mode is long-wave.

Figure 6. Variation of the onset of the pure solutocapillary monotonic instability with the dimensionless surface tension number $\Sigma _0$ for the case of cooling at the substrate $\mathcal{Q}=-1$ at $\Phi =0.01,\,a=7.47,\,\eta =0.31,\,L=10^{-3},\,M_T = 0,\,G=6.71$ and $B = 0.01$ . The $m_l$ and $m_f$ points correspond to the long-wave and finite-wave modes, respectively. The symbols $U$ and $S$ denote the unstable and stable domains of the system, respectively. (a) Variation of the values of the solutal Marangoni number $M_S$ corresponding to the two local minima $M_S=m_l$ and $M_S=m_f$ of the neutral curves $M_S(k)$ . The instability threshold $M_S$ increases with $\Sigma _0$ . The right vertical scale corresponds to $M_S$ related to the long-wave minimum   $M_S=m_l$ . (b) Variation of the wavenumbers corresponding to the two local minima of the neutral curves with $\Sigma _0$ . The critical wavenumber $k_c$ corresponding to the onset of the instability lies in the finite-wave range and varies non-monotonically with the dimensionless surface tension number $\Sigma _0$ . The wavenumber $k$ corresponding to the long-wave competing mode is almost constant with $\Sigma _0$ .

We next study the variation of the onset of the pure solutocapillary instability with the averaged nanoparticle bulk concentration $\Phi$ . Figure 7(a) displays a significant monotonic decrease in the threshold value of the solutal Marangoni number $M_S$ (for both modes $m_f$ and $m_l$ ) with an increase in $\Phi$ . The reason for this effect is mainly due to the Soret effect whose strength is proportional to the local concentration $\phi$ . In the case of cooling at the substrate, the temperature gradient is directed upward; therefore, the Soret effect with $\eta \gt 0$ drives the thermophoretic component of the mass flux downward and becomes stronger with an increase in the average particle concentration $\Phi$ . This leads to a stronger decrease in particle concentration in the neighbourhood of the interface and to an enhanced decrease in local viscosity there with an increase in $\Phi$ . Thus, the critical solutal Marangoni number $M_S$ driving the instability decreases with $\Phi$ .

It is natural to observe such a decrease of the solutocapillary instability threshold for a nanofluid suspension or a mixture with a surfactant where the surface tension linearly decreases with the interfacial concentration (Morozov et al. Reference Morozov, Oron and Nepomnyashchy2014; Shklyaev & Nepomnyashchy Reference Shklyaev and Nepomnyashchy2017). However, for an inorganic salt where the surface tension increases with the interfacial concentration, the opposite behaviour could be observed (Oron & Nepomnyashchy Reference Oron and Nepomnyashchy2004). Further, we note that the critical wavenumber $k_c$ for the finite-wave minimum gradually increases with $\Phi$ , whereas the wavenumber for the long-wave minimum remains almost constant around $k_c\approx 8\times 10^{-4}$ . Since the finite-wave minimum $M_S=m_f$ remains always slightly below the long-wave one $M_S=m_l$ , the instability remains finite-wave, with the wavenumber variation shown in figure 7(b). We also find that the structure of the neutral curves remains unchanged with the variation of $\Phi$ and a typical neutral curve is shown in the inset of figure 7(a).

Figure 7. Variation of the onset of the pure solutocapillary monotonic instability with the averaged bulk nanoparticle concentration $\Phi$ for the case of cooling at the substrate $\mathcal{Q}=-1$ and $a=7.47,\,\eta =0.31, L=10^{-3},\,M_T = 0,\,G=6.71,\,\Sigma _0=2\times 10^{4}$ and $B = 0.01$ . (a) Variation of the two minimal values $M_S=m_f$ and $M_S=m_l$ of the neutral curves $M_S(k)$ . The inset shows the neutral curve for $\Phi =0.03$ . (b) Variation of the wavenumbers corresponding to the two local minima of the neutral curves with $\Phi$ . The upper curve corresponding to $M_S=m_f$ represents the critical wavenumber $k_c$ . The symbols $U$ and $S$ denote the unstable and stable domains of the system, respectively.

Figure 8. (a) Variation of the critical value of the solutal Marangoni number $M_S$ for a pure solutocapillary monotonic instability with the Soret coefficient $\eta$ in the case of cooling at the substrate $\mathcal{Q}=-1$ with $a=7.47,\,\Phi =0.01,\,L=10^{-3},\,M_T = 0,\,G=6.71,\,\Sigma _0=2\times 10^{4}$ and $B = 0.01$ ; the inset shows the difference $m_f-m_l$ between the two minimal values of $M_S$ of the two solutocapillary modes versus $\eta$ . (b) Variation of the critical wavenumber $k_c$ corresponding to the short-wave minimum $M_S=m_f$ with $\eta$ ; the inset shows the variation of the wavenumber $k$ corresponding to the long-wave mode $M_S=m_l$ with $\eta$ . The dashed lines display various ranges of the data fit for the critical wavenumber as a function of $\eta$ . The symbols $U$ and $S$ denote the unstable and stable domains of the system, respectively. FW and LW stand for domains of finite-wave and long-wave instability, respectively.

Figure 8(a) illustrates the non-monotonic variation of the critical solutal Marangoni number $M_S$ with the Soret coefficient $\eta$ for both competing modes. We find that the temperature difference across the layer of the base state given by $\displaystyle \Delta T^*=q^*h_0^*/K_{bf}^*$ influences the value of the Soret coefficient $\eta$ defined in (2.17) and affects the competition between the long-wave deformational mode $M_S=m_l$ and the finite-wave pattern forming mode $M_S=m_f$ . The inset of figure 8(a) displays the domains of instability versus the parameter $\eta$ . We find that in the narrow range of the Soret coefficient $0.62 \leqslant \eta \leqslant 1.24$ corresponding to the temperature difference $2$ K $\leqslant \Delta T^\ast \leqslant 4$ K, the deformational mode $M_S=m_l$ dominates over the finite-wave mode $M_S=m_f$ in terms of the instability onset, namely $m_l\lt m_f$ . Outside of this range of $\eta$ , with a sufficiently low value of the modified Galileo number $G$ , the finite-wave mode $M_S=m_f$ prevails over the long-wave mode $M_S=m_l$ and determines the instability onset, namely $m_f\lt m_l$ . Figure 8(b) displays the variation of the critical wavenumber $k_c$ with the Soret coefficient $\eta$ . We observe that the critical wavenumber increases with $\eta$ and exhibits three different functional variations with $\eta$ . Notably, we observe that the critical wavenumber $k_c$ follows the scaling $k_c \sim \eta ^{1/2}$ up to the value of $\eta$ corresponding to the minimum of $k_c$ for $M_S=m_l$ , which is shown in the inset of figure 8(a). This is followed by two different scalings $k_c \sim \eta ^{4/3}$ and $k_c \sim \eta ^{2/3}$ for higher values of $\eta$ . The inset of figure 8(b) also shows that the critical wavenumber for the long-wave mode $M_S=m_l$ exhibits a non-monotonic variation with $\eta$ in its relevant subdomain.

Finally, we conclude that in the case of the layer cooled at the substrate, the pure solutocapillary instability is monotonic. It also remains monotonic when stabilising thermocapillarity is present along with gravity, the Soret effect and interfacial deformation. The type of this instability is finite-wave for both sufficiently low and sufficiently high values of the modified Galileo number $G$ , whereas in the intermediate range of $G$ , it is long-wave. A similar feature takes place with a variation of the Soret coefficient $\eta$ . Last, but not least, the fact that only the monotonic solutocapillary instability emerges in the case of cooling at the substrate is similar to what Sarma & Mondal (Reference Sarma and Mondal2021b ) found in this case making a statement that the monotonic mode is independent of the degree of viscoelasticity of the binary fluid. However, in contrast with the results obtained by Sarma & Mondal (Reference Sarma and Mondal2021b ) pointing to the emergence of the long-wave instability, our results show that the emerging solutocapillary instability may be either long-wave or finite-wave depending on the system parameters. This difference in the results between the current paper and Sarma & Mondal (Reference Sarma and Mondal2021b ) is possibly due to the presence of concentration-dependent thermophysical properties of the nanofluid.

4.3. Heating at the substrate $(\mathcal{Q} = 1)$

In this section, we investigate the onset of instability in the case of a heated substrate and reveal that it may be either monotonic or oscillatory. Joo (Reference Joo1995) and Oron & Nepomnyashchy (Reference Oron and Nepomnyashchy2004) found that in the case of a layer of a dilute binary mixture with either deformable or non-deformable interface, the system exhibits competition between destabilising thermocapillarity and stabilising solutocapillarity, and as a result, the oscillatory instability emerges.

In the case of a nanofluid layer heated from below, the base state displays a profile of the particle concentration increasing with height. Therefore, if the particle density is higher than that of the liquid, then, in the presence of gravity, there exists an ingredient of Rayleigh–Taylor instability (Chandrasekhar Reference Chandrasekhar1961), which, in what follows, will be referred to as a solutal buoyancy due to its qualitative similarity to the thermal buoyancy.

Figure 9. Eigenspectrum map of the two leading eigenvalues $\lambda$ in the complex plane spanned by their growth rate (the real part $\Re (\lambda ) \equiv \lambda _r$ the vertical axis on the left) and their frequency (the imaginary part $\Im (\lambda ) \equiv \lambda _i$ – the vertical axis on the right) versus $M_T$ in the case of heating at the substrate $\mathcal{Q}=1$ with $\Phi =0.01,$ $G=1,$ $\Sigma _0=10,$ $\eta =0.31,$ $L=10^{-3},$ $a=7.47,$ $M_S=0,$ $P=7,$ $B=0.01$ and $k=0.18$ . The domains $I{-}V$ are the domains with two negative real eigenvalues ( $2\Re _-$ ), two complex conjugate eigenvalues with a negative real part ( $2\mathbb{C}_-$ ), two complex conjugate eigenvalues with a positive real part ( $2\mathbb{C}_+$ ), two positive real eigenvalues ( $2\Re _+$ ), and two real eigenvalues (one positive ( $\Re _+$ ) and one negative ( $\Re _-$ ) ), respectively. The circle $(\circ )$ and star $(\star )$ points correspond to the real and imaginary parts of the complex growth rate $\lambda$ , respectively. The vertical dashed lines represent the transition from the subdomain $II$ to the subdomain $III$ , from $III$ to $IV$ and from $IV$ to $V$ . The horizontal dashed line represents $\lambda _r=0$ .

Figure 9 presents a typical example of the structure of the eigenspectrum of the problem given by (3.7) by following the two leading eigenvalue branches by increasing the thermal Marangoni number $M_T$ for $G=1$ near the critical wavenumber $k_c \approx 0.18$ for a parameter set specified in the caption. The entire variation range of $M_T$ is naturally subdivided into five different subdomains. With an increase in $M_T$ , the first subdomain (subdomain $I$ ) is where the two tracked eigenvalues are real and negative. The two branches of subdomain $I$ merge with an increase in $M_T$ and split into two branches, along which the eigenvalues are complex conjugate with the negative real part $\lambda _r$ ; this is subdomain $II$ . A combination of subdomains $I$ and $II$ belongs to the stability domain of the system. Following the two branches in the subdomain $II$ in figure 9, we observe that with an increase in $M_T$ , the frequency $\lambda _i$ of the two leading eigenvalues increases in its absolute value, whereas the growth rate $\lambda _r$ increases and crosses zero. We emphasise that the oscillatory instability emerges here in the case of $M_S=0$ ; hence, solutocapillarity is irrelevant. This is where the subdomain $II$ ends and the subdomain $III$ begins. In the subdomain $III$ , the two leading eigenvalues are still complex conjugates, but their real part $\lambda _r$ is positive and, therefore, the system undergoes oscillatory instability. In the subdomain $III$ , the growth rate $\lambda _i$ increases with $M_T$ , whereas the absolute value of the frequency $|\lambda _i|$ decreases until it reaches its zero value, and the two leading eigenvalues become real and positive. At this point, the subdomain $III$ ends at $M_T=M_2$ and the subdomain $IV$ begins. In the subdomain $IV$ , there are two real positive eigenvalues, one of them increases with an increase in $M_T$ , whereas the other one decreases until it reaches zero, where the subdomain $IV$ ends at $M_T=M_1$ , and the subdomain $V$ begins. In the subdomain $V$ , there is only one positive eigenvalue and the second leading eigenvalue is negative.

As an illustration for figure 9, figure 10 presents the neutral curves $M_T(k)$ in the case of a moderately low Galileo number $G,\,G=1$ , and the parameter set of figure 9, where the pure thermocapillary oscillatory instability sets in at a finite wavenumber $k \approx 0.18$ . Along with the neutral curve $M_T(k)$ , we present in figure 10 two more boundary curves that illustrate important features of the spectrum of the eigenvalue problem (3.7). The boundary located just above the threshold of the oscillatory instability in figure 10 marked as $M_2$ separates the domain where two complex conjugate eigenvalues with a positive real part exist from the domain of two real leading positive eigenvalues. The space between the $M_O$ and $M_2$ lines constitutes the subdomain $III$ . Slightly further above the boundary $M_2$ lies one more boundary marked as $M_1$ along which one of the leading positive eigenvalues turns to be zero, whereas the other one remains positive. The emergence of several boundaries in a close proximity illustrates a competition between thermocapillarity, gravity and other factors such as the Soret effect, interfacial deformability and possibly concentration-dependence of the thermal conductivity of the suspension which will be discussed below.

Figure 10. Pure thermocapillary oscillatory instability for the case of heating at the substrate $\mathcal{Q}=1$ and the parameter set of figure 9. The lowest curve marked as $M_O$ represents the onset of instability $M_T(k)$ showing the thermal Marangoni number versus the wavenumber $k$ . The subdomains $I {-} V$ correspond to the respective subdomains in figure 9. $M_1$ , $M_2$ and $M_O$ are transition values of $M_T$ between domains $IV$ to $V$ , from $III$ to $IV$ and from $II$ to $III$ , respectively.

The branch of the larger real eigenvalue in the subdomains $IV$ and $V$ corresponds to the thermocapillary mode which increases almost linearly with an increase in the Marangoni number $M_T$ . The branch of the smaller real eigenvalue, which remains positive but very small in the domain $IV$ and crossing zero at $M=M_1$ , is driven by particle diffusion, which is coupled with the unstable thermal mode. To check this observation, we have increased the Lewis number $L$ ten-fold, thereby effectively leading to a ten-fold increase in the smaller eigenvalue and, as a result of this, the positive decreasing branch of the second leading eigenvalue has decayed much faster and, therefore, $M_1$ has decreased. It turns out that the same decreasing branch weakly depends on the surface tension number $\Sigma _0$ as well.

Figure 11. Pure thermocapillary instability in the case of heating at the substrate $\mathcal{Q}=1$ with $\Phi =0.01,$ $\Sigma _0=10,$ $\eta =0.31,$ $L=10^{-3},$ $B=0.01,$ $M_S=0,$ $P=7.$ (a) Variation of the instability threshold $M_T$ with the modified Galileo number $G$ . The circle $(\circ )$ and star $(\star )$ points correspond to the monotonic and oscillatory modes, respectively; inset shows the upper line $M_T(G)$ above which the long-wave mode becomes unstable. Note that these values of $M_T$ are much larger than the critical values of $M_T$ for the oscillatory instability for $G \lt 55.6$ , so the values of the latter appear as zero. However, for $G \geqslant 55.6$ , the solutal buoyancy becomes dominant and drives the long-wave solutal buoyancy monotonic instability with the critical Marangoni number $M_T=0$ and the critical wavenumber $k_c=0$ . (b) Variation of the growth rate $\lambda$ versus wavenumber $k$ for $M_T=0$ for several values of $G$ in the domain of solutal buoyancy instability, $G \gt 55.6$ . The dashed lines illustrate data fit proportional to $k^2$ with the proportionality coefficient increasing with $G$ . The symbols $U,\,U^\dagger, $ and $S$ represent the domains where the two leading eigenvalues are real and positive, complex conjugates with a positive real part and complex conjugates with a negative real part, respectively.

Based on our earlier discussion, it is essential to understand the role of gravity in the onset of instability of the system at hand. Figure 11(a) displays the variation of the critical value of the thermal Marangoni number $M_T$ with the modified Galileo number $G$ in the case of a pure thermocapillary instability. Here, we note that for sufficiently low values of $G$ , the monotonic instability sets in. Note the difference with the case of $G=1$ presented in figures 9 and 10, where the oscillatory instability takes place. Following the monotonic branch of the $M_T{-}G$ plane, we observe that at $G \approx 10^{-1}$ , the monotonic branch becomes the $M_1$ - branch presented in figures 9 and 10 for $G=1$ , and the system undergoes the oscillatory instability. The left vertical dashed line in figure 11(a) represents the asymptotic line for the monotonic branch $\lambda =0$ . Note the $M_2$ branch is not presented here to not overburden the figure. Also, there exists a branch above which the long-wave mode becomes unstable (shown in the inset). However, since the values of $M_T$ along this branch are much higher than those along the monotonic and oscillatory branches discussed above, it is irrelevant for the onset of instability.

As mentioned above, in the case of heating at the substrate considered now, the density stratification in the layer is unstable for particles heavier than the carrier fluid, e.g., $\rho _{np}=4$ , and, in addition to the Marangoni effect, there exists at least one more instability ingredient referred to as solutal buoyancy. Also, gravity plays both stabilising and destabilising roles in the long-wavelength range, i.e. an increase in $G$ may promote the solutal buoyancy instability, but at the same time, the deformable interface tends to flatten out and to relieve a tendency to instability. We note that for relatively low values of $G$ , the solutal buoyancy mechanism is weak and the monotonic thermocapillary instability dominates. However, we find that with an increase in the value of $G$ , the solutal buoyancy mechanism becomes stronger, and, at a sufficiently large value of $G$ , $G \approx 55.6$ (marked by the right vertical dashed line shown also in the inset) for the given parameter set, it dominates over both the oscillatory thermocapillary instability and a long-wave solutal buoyancy instability emerging at $M_T=0$ marked by a horizontal thick line for $G \geqslant 55.6$ with the critical wavenumber $k_c=0$ . Figure 11(b) illustrates the growth rate $\lambda$ variation for different values of $G,\, G\gt 55.6$ with the data fit proportional to $k^2$ for $M_T=0$ . We also find that the proportionality constant increases with the value of $G$ within the unstable domain. As expected, the range of instability widens with an increase in $G$ .

Figure 12. Pure thermocapillary stability boundaries versus wavenumber $k$ in the domain of solutal buoyancy instability for the case of heating at the substrate $\mathcal{Q}=1$ with $\Phi =0.01,$ $\Sigma _0=10,$ $\eta =0.31,$ $L=10^{-3},$ $M_S=0,$ $a=7.47,$ $P=7$ and $B=0.01$ . The two solid (blue) curves on the left and the solid (red) curve on the right represent the monotonic and oscillatory boundaries with the vertical axes on the left and right vertical axes, respectively, for the two values of $G, G=55.8$ and $G=56$ . Note that, unlike the curves for the monotonic instability, the curves for the oscillatory instability almost overlap. Straight dashed lines represent the data fit $M_T \sim \alpha k^{-\beta }$ with $(\alpha =5.77\times 10^{-5}, \beta =2.2)$ for $G=55.8$ and $(\alpha =1.91\times 10^{-4}, \beta =2.12)$ for $G=56$ . The symbols $U,\,U^\dagger, $ and $S$ represent the domains of monotonic solutal buoyancy instability, of two complex conjugate leading eigenvalues with a positive real part, and that of stability, respectively.

Figure 12 presents more details related to the structure of the eigenspectrum of the system for the parameter set used in figure 11 in the range of $G$ adjacent to the onset of the solutal buoyancy instability. We observe that the long-wave monotonic branch $\lambda =0$ of the thermocapillary mode crosses zero at a finite but small wavenumber $k_c$ shown in figure 11(b). These curves represent boundaries of stabilisation of the thermocapillarity where the domain of instability and stability of the system are located beneath and above the corresponding curve. Therefore, higher values of $M_T$ stabilise the system. The instability range of the system in the solutal buoyancy regime expands with an increase in $G$ in terms of $M_T$ and $k_c$ . In the same range of $G$ , we also find that the oscillatory branch shown in the inset of figure 11(a) continues further into the solutal buoyancy domain, as also shown in figure 12 with the corresponding scale on the right-hand side of the figure. It is noted that the values of $M_T$ along the oscillatory finite-wave branch are much higher than those along the monotonic branches. Hence, in the range of finite wavenumbers $k$ , the oscillatory mode emerges in the solutal buoyancy instability domain and, within the range of the wavenumbers corresponding to this instability, the long-wave instability does not set in.

The variation of the critical wavenumber $k_c$ along the instability curves $M_T(G)$ shown in figure 11(a), are displayed in figure 13(a). We note that the critical wavenumber $k_c$ increases monotonically with $G$ along both the monotonic and oscillatory instability thresholds. Also, along the stretch of the monotonic mode when it does not represent the instability onset, the minimal wavenumber is lower than the critical wavenumber $k_c$ . In addition, we find that the critical wavenumber $k_c$ along the monotonic stabilisation boundary tends to zero with an increase in $G$ . Figure 13(b) displays a non-monotonic variation of the growth rate $\lambda$ with $k$ for several values of $G$ near the onset of the solutal buoyancy instability. Figure 13(c) shows the monotonic increase in the critical frequency along the oscillatory branch presented in figure 13(a) with the modified Galileo number $G$ .

Figure 13. Pure thermocapillary instability in the case of heating at the substrate $\mathcal{Q}=1$ with $\Phi =0.01,$ $\eta =0.31,$ $L=10^{-3},$ $a=7.47,$ $\Sigma _0=10,$ $M_S=0,$ $B=0.01$ . (a) Variation of the critical wavenumber $k_c$ versus the modified Galileo number $G$ . The $\circ$ , $\star$ and $\times$ points in the inset represent the critical wavenumber $k_c$ for the monotonic, oscillatory and solutal buoyancy instability, respectively. Inset shows the variation of the minimal wavenumber $k$ in the long-wave mode versus $G$ in the transition domain, $55 \leqslant G \leqslant 56.5$ . (b) Variation of the growth rate $\lambda$ versus wavenumber $k$ in the transition domain for various $G$ . (c) Variation of the frequency $\lambda _i$ versus the modified Galileo number $G$ along the oscillatory instability threshold of panel (a).

Figure 14(a) shows the neutral curves in the case of a pure thermocapillary monotonic instability for a small Galileo number $G=0.01$ for different Biot numbers $B$ . The results for the instability threshold shown in the $M_T{-}G$ plane in figure 11(a) suggest that the parameter set here for $B=0.01$ corresponds to the monotonic instability. The instability remains monotonic for the two other values of the Biot number $B$ . As in the case of the Marangoni instability in a simple fluid (Pearson Reference Pearson1958), we find that the stability domain located underneath the neutral curves expands with an increase in the Biot number accompanied with an increase of the critical Marangoni number $M_T$ at the minimum of the corresponding neutral curve. This is due to the fact that a higher Biot number implies that the system loses some of the heat to the atmosphere, and therefore it requires a higher temperature drop or the heat flux transferred into the system, equivalently, a higher $M_T$ to drive the Marangoni instability (Davis Reference Davis1987). Further, figure 14(b) presents the variation of the critical wavenumber $k_c$ with the Biot number $B$ . It follows the relationship $k_c \sim B^{1/4}$ which is a known asymptotic limit for the Marangoni instability in a simple Newtonian fluid (Sivashinsky Reference Sivashinsky1982) and is also one of the two known asymptotic limits for the long-wave instability in dilute binary fluids (Podolny et al. Reference Podolny, Oron and Nepomnyashchy2005) in the limit of small Biot numbers.

Figure 14. Pure monotonic thermocapillary instability in the case of heating at the substrate $\mathcal{Q}=1$ at $\Phi =0.01,$ $\eta =0.31,$ $L=10^{-3},$ $a=7.47,$ $M_S = 0,$ $\Sigma _0=10,$ $G=0.01$ . (a) Neutral curves $M_T(k)$ for various Biot numbers $B$ . The symbols $U$ and $S$ denote the unstable and stable domains, respectively. (b) Variation of the critical wavenumber $k_c$ with the Biot number $B$ . The dashed line represents the data fit $k_c \sim B^{1/4}$ with the factor of $0.178$ . The filled circle points correspond to the numerical solution of the EVP, (3.7).

Figure 15. Pure monotonic thermocapillary instability in the case of heating at the substrate $\mathcal{Q}=1$ at $B=0.01,$ $\eta =0.31,$ $L=10^{-3},$ $a=7.47,$ $M_S = 0$ and $G=0.01$ for several values of the inverse capillary number $\Sigma _0$ . (a) Variation of the critical thermal Marangoni number $M_T$ with the averaged bulk nanoparticle concentration $\Phi$ , where $ M_T/2$ and $M_T/10$ are presented for $\Sigma _0=10$ and $\Sigma _0=100$ , respectively. The symbols $U$ and $S$ denote the unstable and stable domains of the system, respectively. (b) Variation of the critical wavenumber $k_c$ with $\Phi$ , where $2 k_c$ and $3 k_c$ are presented for $\Sigma _0=10$ and $\Sigma _0=100$ , respectively.

Figure 15(a) presents the variation of the critical thermal Marangoni number $M_T$ with the mean particle concentration $\Phi$ . We note that the latter directly influences the viscosity and the heat conductivity of the fluid as given by (2.2f) and (2.2c ), respectively. We find that the enhanced nanoparticle concentration near the interface promoted by both the base state concentration profile and the Soret effect significantly affects the fluid viscosity near the layer interface, and thereby, notably leads to the delay of the pure thermocapillary instability and an increase in the critical value of the Marangoni number $M_T$ . We also find that the concentration-dependent thermal conductivity also delays the thermocapillary instability, but weakly, whereas the destabilising impact of the unstable stratification is also weak for the given small value of $G$ . Hence, the overall impact of an increase in $\Phi$ is stabilisation of the system, i.e. an ensuing increase in the critical value of $M_T$ . It is also revealed that an increase in the inverse capillary number $\Sigma _0$ leads to an increase in the rigidity of the interface and a significant increase in the critical value of the thermal Marangoni number. We also observe that the critical wavenumber $k_c$ of the purely thermocapillary instability gradually decreases with an increase in $\Phi$ and also with an increase in $\Sigma _0$ , as presented in figure 15(b).

We recall that in the base state when the system is heated at the substrate, the Soret effect with $\eta \gt 0$ causes an increase in the concentration of particles near the layer interface, as shown in figure 3. A change in $\eta$ may be caused by a change in the heat flux $q^\ast$ supplied to the system. Interestingly, an increase in the Soret coefficient $\eta$ leads to an increase in fluid viscosity in general and in the vicinity of the layer interface. As a result of this, the system stabilises, which is expressed by an increase in the critical Marangoni number $M_T$ , as illustrated in figure 16. The critical wavenumber $k_c$ is found to decrease with $\eta$ , similarly to the variation with $\Phi$ , as presented in the inset to figure 16. We note in passing that the solutal buoyancy instability is highly sensitive to the nanoparticle concentration at and near the interface, and indirectly to the strength of the Soret effect; nevertheless, we find that for low values of the modified Galileo number $G$ , the solutal buoyancy mechanism remains weak.

Figure 16. Variation of the critical values of the thermal Marangoni number $M_T$ in the case of a pure monotonic thermocapillary instability with the Soret coefficient $\eta$ in the case of heating at the substrate $\mathcal{Q}=1$ at $B=0.01,$ $\Phi =0.01,$ $L=10^{-3},$ $a=7.47,$ $M_S = 0,$ $\Sigma _0=10,$ $G=0.01$ ; inset shows variation of the critical wavenumber $k_c$ with $\eta$ . The symbols $U$ and $S$ denote the unstable and stable domains of the system, respectively.

Figure 17. Threshold of the combined soluto-thermocapillary instability shown in the plane $M_T{-} M_S$ in the case of heating at the substrate $\mathcal{Q}=1$ with $\Phi =0.01,\,\eta =0.31,\,L=10^{-3},\,\Sigma _0=10,\,B=0.01$ . Panels (a) and (b) correspond to $G=0.01$ and $G=1$ , respectively. The symbols $U,\,U^\dagger $ and $S$ represent the domains of one real positive eigenvalue, of two complex conjugate leading eigenvalues with a positive real part, and that of stability, respectively.

We now consider the onset of the combined soluto-thermocapillary instability. We find that the threshold value of the thermal Marangoni number $M_T$ increases with an increase in the value of the solutal Marangoni number $M_S$ , which implies that solutocapillarity has a stabilising effect, as shown in figure 17. As a reference, in the pure thermocapillary case ( $M_S=0$ ), the instability in the case of the given parameter set with $G=0.01$ is monotonic, as seen in figure 11. With an increase of the solutal Marangoni number, the instability turns from monotonic to oscillatory. With a very small value of the modified Galileo number $G$ , e.g., $G=0.01$ , and with a sufficiently small solutal Marangoni number $M_S$ , $M_S \lesssim 2$ , the emerging instability is monotonic suggesting that thermocapillarity dominates. However, the system undergoes transition from monotonic instability to oscillatory instability with an increase in the value of $M_S$ at approximately $M_S = 2$ , as shown in figure 17(a). The oscillatory instability is apparently induced by the disparity between the characteristic times of the thermal diffusivity and the nanoparticle diffusivity, i.e. $L \ll 1$ . For higher values of $G$ , e.g., $G=1$ , the instability is oscillatory, $M_S=0$ in the pure thermocapillary case, as shown in figure 11. In the combined soluto-thermocapillary case considered here, the instability is oscillatory, as the monotonic boundary remains above the oscillatory one for all $M_S\gt 0$ , as presented in figure 17(b). The critical wavenumber $k_c$ is presented in figure 18(a), which shows that $k_c$ weakly depends on $M_S$ in the case of the monotonic instability, whereas it significantly increases with $M_S$ for the oscillatory instability. It is also interesting to note that the value of $k_c$ exhibits a discontinuity at $M_S$ corresponding to the change in the instability type. This suggests that the transition from the monotonic to oscillatory instability takes place through the codimension-two point. The critical frequency $\lambda _i$ is presented in figure 18(b). It is found to increase monotonically and almost linearly with $M_S$ and collaterally with $M_T$ along the critical curve.

Figure 19 presents the neutral curves $M_T(k)$ for three values of $M_S$ with both monotonic and oscillatory branches shown. In the case of $M_S=0$ displayed in figure 19(a), the oscillatory branch bifurcates off the monotonic one at $M_T$ higher than the minimum of the neutral curve, and being a monotonically increasing function of $k$ , it remains above the critical value of $M_T$ corresponding to the onset of the monotonic instability of the system, as also shown in figure 17. The case of $M_S=2.1$ presented in figure 19(b), illustrates the emergence of a codimension-two point, so the minimal values of the monotonic and oscillatory branches are attained at the same $M_T$ , see the borderline between the two domains in figure 17(a). The critical wavenumber $k_c$ of the oscillatory instability is larger than that of the monotonic instability. This case is in particular interesting due to the dual mechanism of the selection of the type of the emerging instability. Near the threshold of the thermocapillary instability, the emerging instability may be monotonic or oscillatory or of both types, and thus exhibits a competition between the monotonic and oscillatory modes (Brand, Hohenberg & Steinberg Reference Brand, Hohenberg and Steinberg1984). With an increase in $M_S$ , the oscillatory branch bifurcates off the monotonic branch near the minimum of the monotonic branch and descends below it; hence, the emerging instability is oscillatory. The critical wavenumber $k_c$ of the oscillatory instability is larger than that for $M_S=2.1$ , which is consistent with an increase of $k_c$ with $M_S$ shown in figure 18(a) for the oscillatory instability.

Figure 18. The case of heating at the substrate $\mathcal{Q}=1$ at $\Phi =0.01,$ $\eta =0.31,$ $L=10^{-3},$ $a=7.47,$ $\Sigma _0=10,$ $G=0.01,$ $B=0.01$ . (a) Variation of the critical wavenumber $k_c$ with the solutal Marangoni number $M_S$ along the critical curve shown in figure 17(a). Note that the values of $k$ along the monotonic curve protruding into the domain of the oscillatory instability which are shown by the hollow circles are not critical. (b) Variation of the critical frequency $\lambda _i$ with the solutal Marangoni number $M_S$ along the critical curve presented in figure 17(a) in the domain of the oscillatory instability.

Figure 19. Neutral curves $M_T(k)$ in the case of heating at the substrate $\mathcal{Q}=1$ with $\Phi =0.01,$ $\eta =0.31,$ $L=10^{-3},$ $a=7.47,$ $\Sigma _0=10,$ $G=0.01$ and $B=0.01$ . Panels (a), (b) and (c) display the structure of the neutral curves for the monotonic $(\circ )$ and oscillatory branches $(\star )$ for $M_S=0,\,2.1$ and $20$ , respectively. Panel (b) shows the emergence of the co-dimension two point at $M_S=2.1$ . The symbols $U,\,U^\dagger $ and $S$ represent the domains of one real positive eigenvalue, two leading complex conjugate eigenvalues with a positive real part, and that of the system stability, respectively.

Figure 20. Variation of the critical thermal Marangoni number $M_T$ with the thermal conductivity stratification parameter $a$ in the case of heating at the substrate $\mathcal{Q}=1$ with $\Phi =0.01,$ $\eta =0.31,$ $L=10^{-3},$ $M_S = 0,$ $\zeta =0,$ $G=0$ and $B=0.01$ . The symbols $U,\,U^\dagger $ and $S$ denote the domains with one real positive eigenvalue, with two leading complex conjugate eigenvalues with a positive real part, and that of stability, respectively. The $\circ$ and $\star$ symbols denote, respectively, the onset of the monotonic and oscillatory instabilities in the case of a constant Brownian diffusivity $D_B$ , whereas $\times$ and $\diamond$ symbols represent, respectively, the onset of monotonic and oscillatory instability with nanoparticle concentration-dependent Brownian diffusion $D(\phi )$ given by (2.5b ).

4.4. Influence of the thermal conductivity stratification of the nanofluid

In this subsection, we study the effect of thermal conductivity stratification of the nanofluid on the type of Marangoni instability. To simplify the analysis, we assume here the absence of gravity, $G=0$ , the non-deformability of the interface, $\zeta =0$ , and the absence of solutocapillarity, $M_S=0$ . Figure 20 displays the variation of the critical thermal Marangoni number $M_T$ with the thermal conductivity stratification parameter $a$ . We find that the monotonic instability sets in for sufficiently small values of $a\lt a_c\sim 6\times 10^{-4}$ . However, the oscillatory instability emerges when the thermal conductivity stratification further increases. Therefore, fluids with nanoparticles with a low thermal conductivity undergo the monotonic instability, whereas those made of metals, i.e. with a higher thermal conductivity, for instance, alumina, copper, etc. (Buongiorno Reference Buongiorno2006; Coccia, Tomassetti & Di Nicola et al. Reference Coccia, Tomassetti and Di Nicola2021), are expected to exhibit oscillatory instability.

We recall that throughout this paper, we used a constant value for the Brownian diffusion coefficient. As follows from (2.19), with low values of the thermal conductivity parameter $a$ , the coefficient of the term linear in $\phi$ for $D(\phi )$ in equation (2.19) is larger than $a$ ; therefore, it may be inferred that the variation of the Brownian diffusivity with particle concentration must be accounted for. Figure 20 also demonstrates the threshold of the thermocapillary instability for the case in which the nanoparticle concentration-dependent Brownian diffusion coefficient $D(\phi )$ given by (2.19) is taken into account. We find that the critical value of $M_T$ based on the use of the concentration-dependent Brownian diffusion coefficient given by Batchelor (Reference Batchelor1976) differs from $M_T$ determined using a constant Brownian diffusion coefficient $D_B$ by less than $2\,\%$ . We also find that the instability type, either monotonic or oscillatory, and the transition value of $a=a_c$ from one to another are not affected by whether a constant or particle concentration-dependent Brownian diffusion coefficient is used. Moreover, we note that the variation of the threshold of the oscillatory thermocapillary instability up to the upper bound of the Hashin–Shtrikman interval (Hashin & Shtrikman Reference Hashin and Shtrikman1962; Keblinski et al. Reference Keblinski, Prasher and Eapen2008) remains smooth, and an increase in $a$ leads to a substantial increase in $M_T$ (not shown).

Figure 21(a) demonstrates that the critical wavenumber $k_c$ for the onset of Marangoni instability lies in the short-wave domain for both monotonic and oscillatory instabilities. The critical wavenumber slowly increases with $a$ continuously changing through the transition value of $a=a_c$ from the monotonic to oscillatory domains. Figure 21(b) shows an increase in the critical frequency $\lambda _i$ with an increase in $a$ . Again, the difference between the results based on a constant $D_B$ and nanoparticle concentration-dependent forms for the Brownian diffusivity $D(\phi )$ is minor.

Figure 21. The case of heating at the substrate $\mathcal{Q}=1$ at $\Phi =0.01,$ $\eta =0.31,$ $L=10^{-3},$ $M_S = 0,$ $\zeta =0,$ $G=0$ and $B=0.01$ . (a) Variation of the critical wavenumber $k_c$ with the thermal conductivity stratification parameter $a$ . The $\circ$ and $\star$ symbols show $k_c$ of the monotonic and oscillatory instabilities, respectively, when a constant Brownian diffusivity $D_B$ is used, whereas the $\times$ and $\diamond$ symbols denote $k_c$ of the monotonic and oscillatory instability, respectively, when the nanoparticle concentration-dependent Brownian diffusion $D(\phi )$ is used. The circles in the oscillatory domain denote the values of the wavenumber corresponding to the monotonic mode and they do not represent critical values. The $U,\,U^\dagger, $ and $S$ symbols represent the domains of monotonic instability, oscillatory instability and stability of the system, respectively. (b) Variation of the critical frequency $\lambda _i$ with the thermal conductivity stratification parameter $a$ in the domain where the oscillatory instability sets in. The $\star$ and $\diamond$ symbols represent the critical frequency $\lambda _i$ corresponding to a constant and nanoparticle concentration-dependent forms for the Brownian diffusivity, $D_B$ and $D(\phi )$ , respectively.

Figure 22. Variation of the critical thermal conductivity stratification parameter $a_c$ with the Lewis number $L$ in the case of heating at the substrate $\mathcal{Q}=1$ with $\Phi =0.01,$ $\eta =0.31,$ $M_S = 0,$ $\zeta =0,$ $G=0$ and $B=0.01$ . The thermal conductivity stratification parameter $a$ varies in the domain $a\in (0,7.47)$ . The dashed line represents the data fit $a_c\propto L^2$ with the factor of $8.80\times 10^2$ .

The variation of the critical thermal stratification parameter $a_c$ for the onset of oscillatory instability with the Lewis number $L$ is shown in figure 22 for the case of a layer with the non-deformable interface. It is found that in the limit of small $L$ , the value of $a_c$ is proportional to $L^2$ . We note that having a large value for the proportionality factor, the nanofluid layer with the particle concentration diffusion time scale comparable to the thermal diffusion time scale, e.g., with the Lewis number $L\sim O(10^{-1})$ , requires a significantly higher value of the thermal stratification parameter $a=a_c$ for the system to switch from the monotonic to oscillatory instability. However, with a sufficiently low value of the Lewis number $L\ll 1$ , the oscillatory instability occurs already at a very low thermal conductivity stratification parameter $a_c \ll 1$ . We also note that the critical value $a_c$ in figure 22 will be different for a layer when the interfacial deformation is present. The effect of the interfacial deformation via the dimensionless surface tension number $\Sigma _0$ on the variation of $a_c$ will be discussed elsewhere.

4.5. Eigenfunctions and physical mechanism

We examine the underlying physical mechanism driving the instabilities of our system by presenting sets of eigenfunctions near the criticality, i.e. the critical wavenumber $k_c$ , where the growth rate is the highest and the value of the control parameter, here the solutal or the thermal Marangoni number, reaches the minimal value along the neutral curve.

The purely solutocapillary mechanism of instability in the case of a nanofluid layer cooled at the substrate works in the following way. The concentration profile in this configuration exhibits a gravity stable stratification, e.g., a higher concentration of heavier nanoparticles near the substrate compared with their lower concentration near the interface. Now, suppose that due to an infinitesimal disturbance, a fluid packet with a higher concentration and lower temperature from underneath the interface is displaced towards it, where the surface tension decreases with the concentration. This creates a cooler spot with a higher concentration at the layer interface. If the thermocapillary effect is absent, $\sigma^\ast_{T^\ast}=0$ , a higher concentration spot represents a spot with a lower surface tension inducing flow along the interface emanating from it towards the domain with a higher surface tension which is that of a lower concentration. The flow is then fed by mass conservation bringing fresh fluid packets with even higher concentration from the bulk.

Figure 23. Normalised eigenfunctions of the EVP (3.7) in the case of cooling at the substrate $\mathcal{Q} = -1$ for the critical wavenumber $k_c = 0.16$ with $L=10^{-3},$ $\Phi =0.01,$ $\eta =0.31,$ $a=7.47,$ $B=0.01,$ $\Sigma _0\approx 2\times 10^{4},$ $G=6.71,$ $M_S = 23,\,M_T=0$ and $\lambda =1.1707\times 10^{-7}.$ The eigenfunctions $\bar {\phi }(x,z)$ and $\bar {T}(x,z)$ superimposed with the velocity vector field $\bar {\textbf{u}}(x,z)$ are shown in panels (a) and (b), respectively. The velocity vector field $\bar {\textbf{u}}(x,z)$ shows the convective flow driven from the low surface tension spot (high nanoparticle concentration) towards that of the high surface tension (low nanoparticle concentration).

Figure 23 displays a typical set of the eigenfunctions for the monotonic instability emerging in the case of the system cooled at the substrate. Panels (a) and (b) present the eigenfunctions for the concentration disturbances $\bar {\phi } (x,z)$ and those for the temperature $\bar {T}(x,z)$ , respectively, both superimposed with the flow velocity field presented by the velocity vectors. In figure 23(a) of nanoparticle concentration $\bar {\phi }(x,z)$ , we observe a rising flow in the high-concentration domain which is driven by the solutocapillary shear stress towards the low-concentration domain forming a descending flow there. Gravity, thermal diffusion and viscosity of the fluid all act against the destabilising solutocapillary effect and saturate the instability, so that the latter is monotonic. It is interesting to note that both fields of concentration and temperature disturbances $\bar {\phi }$ and $\bar {T}$ , respectively, form stripes. This is explained by the fact that the functions $\phi$ and $\theta$ found from the numerical solution of the eigenvalue problem (3.7) weakly vary with the height $z$ when the former retains its sign, whereas the latter changes its sign.

Figure 24. Normalised eigenfunctions of the EVP (3.7) in the case of heating at the substrate $\mathcal{Q} = 1$ for the critical wavenumber $k_c = 0.18$ with $L=10^{-3},$ $\Phi =0.01,$ $\eta =0.31,$ $a=7.47,$ $B=0.01,$ $\Sigma _0=10,$ $G=1,$ $M_S =0,\,M_T=1.22\ \text{and}\ \lambda =2.0773\times 10^{-6}+3.6965\times 10^{-5}i.$ The eigenfunctions for the concentration $\bar {\phi }(x,z)$ and temperature $\bar {T}(x,z)$ disturbances superimposed with the velocity vector field $|\bar {\textbf{u}}(x,z)|$ , are shown in panels (a) and (b), respectively.

We next present a set of eigenfunctions in the case of a pure oscillatory thermocapillary instability in a layer heated at the substrate in figure 24. Both $\varphi (z)$ and $\theta (z)$ are found to weakly depend on $z$ , so both eigenfunctions $\bar {\phi }(x,z)$ and $\bar {T}(x,z)$ display rolls contained between the extrema of the interfacial deformation. Since the depression of the interface is closer to a hot substrate and the interfacial elevation is farther away from it, their temperatures correspond to the maximum and the minimum of the interfacial temperature, respectively. This creates thermocapillary shear stresses directed away from the depression where the surface tension is the lowest to the elevation where the surface tension is the highest, thereby driving a flow within the entire bulk by means of fluid viscosity. In figure 24(b) of $\bar {T}(x,z)$ , we note the formation of the upwelling flow in the higher-temperature domain and the descending flow in the lower-temperature domain. As in the case shown in figure 23, the patterns shown in the profiles of the concentration and temperature disturbances eigenfunctions are stripes, since their amplitude functions $\phi$ and $\theta$ , respectively, depend weakly on $z$ . However, there is a phase shift between these stripes. This phase shift is due to the fact that in one of them, the real part dominates its imaginary part, whereas the opposite takes place in the other.

4.6. Incompressibility simplification

The fluid density, as described in § 2.1, varies with space, time and depends on the particle concentration which itself evolves in time and space. Therefore, generally speaking, the fluid in the problem at hand is ‘compressible’. It is possible to rewrite the non-dimensional form of the nanofluid continuity equation (2.13a ) as

(4.1) \begin{equation} \nabla \cdot \textbf{u} = -\left (\rho _{np}-1\right )\phi _m\left [P\phi _t + \nabla \cdot \left (\phi \textbf{u}\right )\right ]. \end{equation}

Based on (4.1) and on the nanoparticle mass flux balance equation (2.13d ), the divergence of the velocity vector is obtained in the form

(4.2) \begin{equation} \nabla \cdot \textbf{u} = -L\left (\rho _{np}-1\right ) \phi _m\left [\nabla ^2\phi + \nabla \cdot \left (\eta \phi \nabla T\right )\right ]. \end{equation}

The right-hand side of (4.2) represents the deviation of the full ‘compressible’ system from an incompressible one. Assuming that for a moderately dense nanofluid in the case of a small Lewis number $L$ satisfying the condition $L\in (10^{-4},10^{-2} )$ , the parameter $L(\rho _{np}-1)$ is $O (10^{-4}-10^{-2} )$ . Thus, (4.2) suggests that the nanofluid system at hand can be, by neglecting its right-hand side, simplified by imposing its ‘incompressibility’. The question is now about the implications of this simplification.

We now compare the results obtained for the problem governed in its full formulation by the EVP (3.7) with a non-uniform density depending on the space-time dependent particle concentration, and in this case, the problem is referred to as ‘compressible’, with those obtained for a simplified problem where the right-hand side of (4.2) is neglected and the rest of the governing equations and boundary conditions remain with no change. This simplified formulation will be referred to as an ‘incompressible’ one. In some sense, the latter is akin to the Boussinesq approximation where the density is assumed to be constant except for allowing for the buoyancy force arising from the variation in the fluid density.

Figure 25(a) illustrates the comparison between the neutral curve $M_S(k)$ for a pure solutocapillary instability in the case of cooling at the substrate for compressible and incompressible formulations of the problem. We observe that the neutral curve of the simplified incompressible problem displays a close match with the neutral curve obtained for the compressible formulation. Note the difference between the two which does not exceed 2 %. Further, figure 25(b) illustrates an excellent agreement between the variation of the critical thermal Marangoni number $M_T$ with the thermal conductivity stratification parameter $a$ for both monotonic and oscillatory instabilities in the case of heating at the substrate. Note that the curves for both the compressible and incompressible formulations almost fully overlap with the maximal difference of $1.66\,\%$ . Therefore, we infer that the incompressible simplification can be safely used for a nanofluid with stratification of thermophysical properties. We emphasise that all of the results presented here were obtained without using the incompressibility simplification, but note that this simplification may significantly reduce the numerical effort needed for treatment of the problem.

Figure 25. (a) Neutral curves for a pure solutocapillary monotonic instability $M_S(k)$ for the case of cooling at the substrate $\mathcal{Q}=-1$ with $B = 0.01,$ $\Phi =0.01,$ $\eta =0.31,$ $a=7.47,$ $L=10^{-3},$ $M_T = 0,$ $\Sigma _0\approx 2\times 10^{4}$ and $G=6.71$ . The $\circ$ and $\diamond$ points correspond to the cases of the full EVP (3.7) and a simplified incompressible formulation, respectively. (b) Variation of the critical thermal Marangoni number with the thermal conductivity stratification parameter $a$ for a pure thermocapillary instability in the case of heating at the substrate $\mathcal{Q}=1$ with $\Phi =0.01,$ $\eta =0.31,$ $L=10^{-3},$ $M_S = 0,$ $\zeta =0,$ $G=0,$ and $B=0.01$ . The $\circ$ , $\star$ , $\diamond$ and $\times$ points represent the values obtained for the monotonic and oscillatory instabilities based on the full EVP (3.7), respectively, and the monotonic and oscillatory instabilities obtained for a simplified incompressible formulation, respectively. The $U,\,U^\dagger, $ and $S$ symbols represent the domains of monotonic instability, oscillatory instability, and stability of the system, respectively.