2.1. Thermophysical properties and physical mechanisms

The thermophysical properties of a nanofluid exhibit a strong dependence on the local nanoparticle concentration (Buongiorno Reference Buongiorno2006). We consider the density and specific heat of a nanofluid layer as a weighted sum of the respective properties, with the particle volume fraction as a weight

(2.1) \begin{equation} \rho ^\ast _{\textit{nf}}(\phi ^\ast ) =(1-\phi ^\ast )\rho ^\ast _{\textit{bf}}+\phi ^\ast \rho ^\ast _{\textit{np}},\quad c^\ast _{\textit{nf}}(\phi ^\ast ) = \frac {(1-\phi ^\ast )(\rho ^\ast c^\ast ) _{\textit{bf}}+ \phi ^\ast (\rho ^\ast c^\ast ) _{\textit{np}}}{\rho _{\textit{nf}}^\ast }. \end{equation}

Figure 1.

Nanofluid layer with particles of diameter $d_p^* \sim 50 \,\text{nm}$ deposited on the solid substrate, subjected to a constant heat flux at the substrate and exposed to the gas phase at its deformable interface with interfacial adsorption/desorption particle kinetics.

Pak & Cho (Reference Pak and Cho1998), Buongiorno et al. (Reference Buongiorno2009) and Keblinski, Prasher & Eapen (Reference Keblinski, Prasher and Eapen2008) showed experimentally a significant enhancement of thermal conductivity with the Al $_2$ O $_3$ nanoparticles in water as a base liquid. Buongiorno (Reference Buongiorno2006) suggested the empirical linear fit for the thermal conductivity variation with local nanoparticle concentration, which is valid for a moderately dense nanofluid. This empirical model for the thermal conductivity reads

(2.2) \begin{equation} K^\ast _{\textit{nf}}(\phi ^\ast ) = K^\ast _{\textit{bf}}(1 + a \phi ^\ast ), \end{equation}

where $K_{\textit{bf}}^*$ is the thermal conductivity of the base fluid and the parameter $a$ represents the strength of variation of the thermal conductivity with the local nanoparticle concentration  $\phi ^\ast$ . We note that the thermal conductivity model given by (2.2) is similar in its form to the classical model given by Maxwell (Reference Maxwell1873) in the limit of a dilute suspension of monodispersed hard spheres in a single-component liquid. Jeffrey (Reference Jeffrey1973) extended the model for the case of the thermal conductivity of a suspension of Maxwell (Reference Maxwell1873) valid up to $O(\phi ^*)$ to the higher order in local nanoparticle concentration $O(\phi ^{*2})$ . The empirical fit given by Buongiorno (Reference Buongiorno2006) for the experiments of Pak & Cho (Reference Pak and Cho1998) in the case of alumina nanoparticles in water leads to the value of the parameter $a = 7.47$ , on the other hand, the corresponding value of $a$ obtained from the Maxwell (Reference Maxwell1873) model yields $a \leqslant 3$ , which is quite far from the experimental value $a=7.47$ . This difference could be associated with the possibility of the formation of aggregates. In what follows, we use the value $a=7.47$ . In addition, it is possible to estimate the value of the thermal conductivity variation via the Hashin–Shtrikman bounds (Hashin & Shtrikman Reference Hashin and Shtrikman1962; Keblinski et al. Reference Keblinski, Prasher and Eapen2008)

(2.3) \begin{equation} \left [1+\frac {3\phi ^* \big (K_{\textit{np}}^*-K_{\textit{bf}}^*\big )}{3K_{\textit{bf}}^*+\left (1-\phi ^*\right )\big (K_{\textit{np}}^*-K_{\textit{bf}}^*\big )}\right ]\leqslant \frac {K_{\textit{nf}}^*}{K_{\textit{bf}}^*}\leqslant \left [1-\frac {3(1-\phi ^*)\big (K_{\textit{np}}^*-K_{\textit{bf}}^*\big )}{3K_{\textit{np}}^*-\phi ^*\big (K_{\textit{np}}^*-K_{\textit{bf}}^*\big )}\right ]\frac {K_{\textit{np}}^*}{K_{\textit{bf}}^*}. \end{equation}

For the low nanoparticle concentration limit, $\phi ^*\ll 1$ , the thermal conductivity parameter $a$ for the Al $_2$ O $_3$ nanoparticles in water ranges in the interval $2.84\le a\le 38$ , so the value $a=7.47$ belongs to this interval.

For the nanofluid viscosity, we use the model suggested by de Kruif et al. (Reference de Kruif, van Iersel, Vrij and Russel1985) for the hard-sphere nanoparticle suspension in an electrostatically stabilised base liquid, i.e.

(2.4) \begin{equation} \mu ^\ast _{\textit{nf}}(\phi ^\ast ) = \mu ^\ast _{\textit{bf}} \left (1-\frac {\phi ^\ast }{\phi _m}\right )^{-2}, \end{equation}

where $\phi _m$ stands for random close packing. The values of $\phi _m$ vary in the interval between $0.524-0.71$ . We use the value $\phi _m=0.65$ satisfying the experimental observations of de Kruif et al. (Reference de Kruif, van Iersel, Vrij and Russel1985) for the zero-shear limit.

The absolute velocity of a nanoparticle is described by the sum of the carrier fluid velocity and the relative (slip) velocity. Buongiorno (Reference Buongiorno2006) describes seven different slip mechanisms for a nanofluid suspension, which are inertia, Brownian diffusion, thermophoresis or the Soret effect, diffusiophoresis, Magnus effect, fluid drainage and gravitational settling. Out of these seven, we consider the two dominant slip mechanisms, i.e. Brownian diffusion and the Soret effect.

The Brownian diffusivity is determined from the Stokes–Einstein formula, i.e.

(2.5) \begin{equation} D_B = \frac {K_B^*T^\ast }{3\pi \mu ^\ast _{\textit{bf}}d_p^\ast }, \end{equation}

where $K_B^*$ is the Boltzmann constant, $K_B^*= 1.380649 \times 10^{-23}$ J/K. We note that Brownian diffusivity also depends on the nanoparticle concentration (Batchelor Reference Batchelor1976; Russel et al. Reference Russel, Saville and Schowalter1989; Bird, Stewart & Lightfoot Reference Bird, Stewart and Lightfoot2002); however, for a low particle concentration, the particle concentration-dependent Brownian diffusivity imparts negligible influence on the stability (Gandhi et al. Reference Gandhi, Nepomnyashchy and Oron2025c ); therefore, in this condition, the constant Brownian diffusivity provides a good approximation.

In a non-isothermal nanofluid layer, the temperature gradient induced between the solid substrate and the deformable interface can induce another slip mechanism known as thermophoresis that is similar to the Soret effect in the kinetic theory of gases. We assume the Soret effect for the single-component nanofluid mixture given by (Whitmore & Meisen Reference Whitmore and Meisen1977; Buongiorno Reference Buongiorno2006)

(2.6) \begin{equation} D_T = \left (\frac {0.26K^\ast _{\textit{bf}}}{2K^\ast _{\textit{bf}}+K^\ast _{\textit{np}}}\right )\frac {\mu ^\ast _{\textit{bf}}}{\rho ^\ast _{\textit{bf}}}\phi ^\ast . \end{equation}

Thence, the total nanoparticle mass flux is expressed by the mass flux induced by the Brownian diffusion and the Soret effect as

(2.7) \begin{equation} \boldsymbol{j}_p = -\rho ^\ast _{\textit{np}}\left ( C_BT^\ast \boldsymbol{\nabla} ^\ast \phi ^\ast + C_T\phi ^\ast \frac {\boldsymbol{\nabla} ^\ast T^\ast }{T^\ast }\right ), \end{equation}

where $C_B = D_B/T^\ast$ , $C_T = D_T/\phi ^\ast $ and $\displaystyle \boldsymbol{\nabla} ^\ast = (\partial _{x^\ast }, \partial _{y^\ast },\partial _{z^\ast } )$ , with subscripts ${x^\ast }, {y^\ast },{z^\ast }$ representing partial differentiation with respect to the corresponding variable. We note in passing that the total mass flux of nanoparticles given by (2.7) ensures a positively definite distribution of nanoparticle concentration, even subjected to strong thermophoresis. The other five slip mechanisms are negligible for a thin nanofluid layer with the diameter $d_p^*=50$ nm; however, the gravitational setting becomes important for the thick layer or in suspensions with larger nanoparticles (Chang & Ruo Reference Chang and Ruo2022; Gandhi et al. Reference Gandhi, Nepomnyashchy and Oron2025c ).

We assume that the heat flux in a nanofluid is given by the Fourier’s law of heat conduction, i.e.

(2.8) \begin{equation} \boldsymbol{j}_T = - K_{\textit{nf}}^\ast \boldsymbol{\nabla} ^\ast T^\ast , \end{equation}

where we neglect the heat transfer contribution due to the Dufour effect and the nanoparticle mass transfer. The Dufour effect is negligible for the liquid mixture (Cross & Hohenberg Reference Cross and Hohenberg1993; Oron & Nepomnyashchy Reference Oron and Nepomnyashchy2004); however, important for the heat transfer in gases. We note that other components of the mass flux of nanoparticles (Buongiorno Reference Buongiorno2006), such as sedimentation (Chang & Ruo Reference Chang and Ruo2022), can also contribute to the total heat flux. However, in the case of a thin nanofluid layer with small-sized nanoparticles, they are negligible (Gandhi et al. Reference Gandhi, Nepomnyashchy and Oron2025c ).

2.2. Governing equations and boundary conditions

The set of governing equations describing the dynamics of a nanofluid layer contains the continuity, Navier–Stokes, advection–conduction heat transfer and particle mass transfer equations, respectively (Rohsenow, Hartnett & Cho Reference Rohsenow, Hartnett and Cho1998; Batchelor Reference Batchelor2000; Colinet, Legros & Velarde Reference Colinet, Legros and Velarde2001; Bird et al. Reference Bird, Stewart and Lightfoot2002):

(2.9a) \begin{align} \big(\rho _{\textit{nf}}^*\big)_{t^*} + \boldsymbol{\nabla} ^\ast & \boldsymbol{\cdot }\big(\rho _{\textit{nf}}^*\boldsymbol{u}^*\big) = 0, \end{align}

(2.9b) \begin{align} \big(\rho _{\textit{nf}}^*\boldsymbol{u}^*\big)_{t^*} + \boldsymbol{\nabla} ^\ast \boldsymbol{\cdot }\big(\rho _{\textit{nf}}^*\boldsymbol{u}^*\boldsymbol{u}^*\big) &= -\boldsymbol{\nabla} ^\ast p^* + \boldsymbol{\nabla} ^\ast \boldsymbol{\cdot }\boldsymbol{\tau }^* -\rho _{\textit{nf}}^*g^*\boldsymbol{e}_{z^*}, \end{align}

(2.9c) \begin{align} \big(\rho ^* c^*\big)_{\textit{nf}}T^*_{t^*} + \big(\rho ^* c^*\big)_{\textit{nf}} \big(\boldsymbol{u}^\ast \boldsymbol{\cdot }\boldsymbol{\nabla} ^\ast \big) T^\ast &= \boldsymbol{\nabla} ^\ast \boldsymbol{\cdot }\big(K_{\textit{nf}}^*\boldsymbol{\nabla} ^\ast T^*\big), \end{align}

(2.9d) \begin{align} \phi ^*_{t^*} + \boldsymbol{\nabla} ^\ast \boldsymbol{\cdot }\big(\boldsymbol{u}^*\phi ^*\big) &= \boldsymbol{\nabla} ^\ast \boldsymbol{\cdot }\biggl (C_BT^*\boldsymbol{\nabla} ^\ast \phi ^* +C_T\phi ^*\frac {\boldsymbol{\nabla} ^\ast T^*}{T^*}\biggr ). \end{align}

Here, $g^*$ is the gravity acceleration, $\boldsymbol{e}_{z^*}$ is the unit vector in the $z^\ast$ direction, ${\boldsymbol u}^* = (u^*,v^*,w^* )$ is the flow velocity field vector and $\boldsymbol{\tau }^*$ is the viscous part of the stress tensor $\displaystyle \boldsymbol{\tau }^* =\mu _{\textit{nf}}^\ast [\boldsymbol{\nabla} ^\ast {\boldsymbol u^*}+ (\boldsymbol{\nabla} ^\ast {\boldsymbol u^*} )^\intercal ] - ({2}/{3}\mu _{\textit{nf}}^\ast - \mathcal{K}^\ast ) (\boldsymbol{\nabla} ^\ast \boldsymbol{\cdot }\boldsymbol{u}^\ast ) \boldsymbol{I}$ , where the superscript $^\intercal$ denotes the transpose of the corresponding tensor, $\boldsymbol{I}$ is the unity tensor, the subscript $t^\ast$ stands for a partial derivative with respect to time $t^\ast$ and $\mathcal{K}^\ast$ is the dilatational viscosity of the fluid. We note that the last term in (2.9b ) represents the buoyancy force, which is due to the fluid density varying within the layer with the particle concentration that is, in turn, coupled to the fluid temperature via the governing equations (2.9).

We impose three boundary conditions at the solid–liquid interface $z^*=0$ . The boundary conditions at the solid substrate are

(2.10) \begin{equation} z^\ast =0: \quad \boldsymbol{u}^* = 0, \quad -K_{\textit{nf}}^{*}T^*_{z^\ast } = -q^*, \quad C_BT^* \phi ^*_{z^\ast }+C_T\phi ^*\frac {T^*_{z^\ast }}{T^*} = 0. \end{equation}

The first condition represents the no-slip and no-penetration for the flow velocity at the substrate. The second condition describes the heat transfer to the solid substrate made of a material with a low thermal conductivity, which is represented by a specified heat flux $-q^*$ (Rohsenow et al. Reference Rohsenow, Hartnett and Cho1998). Finally, the third condition corresponds to the impermeability of the substrate in terms of the total mass flux of the nanoparticles.

In what follows, we consider the deformable liquid–air interface exposed to the gas phase in the gravity field. The motion of the interface located at $z^*=h^*(x^*,y^*,t^*)$ is governed by the kinematic boundary condition

(2.11) \begin{equation} h_{t^*}^*+u^*h_{x^*}^*+v^*h_{y^*}^*=w^*. \end{equation}

The next boundary condition at the deformable interface $z^* = h^*(x^*,y^*,t^*)$ is the continuity of the normal and tangential stresses

(2.12) \begin{equation} \boldsymbol{n}^\ast \boldsymbol{\cdot }\boldsymbol{\mathcal{T}}^*\boldsymbol{\cdot }\boldsymbol{n}^\ast = 2\boldsymbol{\mathcal{H}}^\ast \sigma ^\ast ,\quad \boldsymbol{n}^\ast \boldsymbol{\cdot }\boldsymbol{\mathcal{T}}^*\boldsymbol{\cdot }\boldsymbol{x}^\ast = \boldsymbol{x}^\ast \boldsymbol{\cdot }\boldsymbol{\nabla} ^\ast _s\sigma ^*,\quad \boldsymbol{n}^\ast \boldsymbol{\cdot }\boldsymbol{\mathcal{T}}^*\boldsymbol{\cdot }\boldsymbol{y}^\ast = \boldsymbol{y}^\ast \boldsymbol{\cdot }\boldsymbol{\nabla} ^\ast _s\sigma ^*, \end{equation}

where $\boldsymbol{\mathcal{T}}^*= -p^*\boldsymbol{I} + \boldsymbol{\tau }^\ast$ is the total stress tensor. Here, $\displaystyle \boldsymbol{\nabla} ^\ast _s = (\boldsymbol{I} - \boldsymbol{n}^\ast \boldsymbol{n}^\ast )\boldsymbol{\cdot }\boldsymbol{\nabla} ^\ast$ is the surface gradient operator,

(2.13) \begin{equation} \boldsymbol{n}^\ast = \frac {-h_{x^*}^*\boldsymbol{e}_{x^*}-h_{y^*}^*\boldsymbol{e}_{y^*}+\boldsymbol{e}_{z^*}}{\sqrt {1+h_{x^*}^{*2}+h_{y^*}^{*2}}} \end{equation}

is the unit vector normal to the interface, $\boldsymbol{e}_{x^\ast }$ and $\boldsymbol{e}_{y^\ast }$ are the orthonormal vectors in the $x^\ast$ and $y^\ast$ directions, respectively,

(2.14) \begin{equation} \boldsymbol{x}^\ast = \frac {\boldsymbol{e}_{x^*}+h_{x^*}^*\boldsymbol{e}_{z^*}}{\sqrt {1+h_{x^*}^{*2}}},\qquad \boldsymbol{y}^\ast =\frac {\boldsymbol{e}_{y^*}+h_{y^*}^*\boldsymbol{e}_{z^*}}{\sqrt {1+h_{y^*}^{*2}}}, \end{equation}

are unit vectors tangent to the interface and $2\boldsymbol{\mathcal{H}}^*$ is the mean curvature of the interface

(2.15) \begin{equation} 2\boldsymbol{\mathcal{H}}^* = -\boldsymbol{\nabla} ^\ast _s\boldsymbol{\cdot }\boldsymbol{n}^\ast . \end{equation}

Henceforth, following Stokes’ hypothesis (Buresti Reference Buresti2015), we assume a zero value of the dilatational viscosity $\mathcal{K}^\ast$ . This signifies that the absolute value of $\displaystyle \mathcal{K}^\ast \boldsymbol{\nabla} ^\ast \boldsymbol{\cdot }\boldsymbol{u}^\ast$ is negligible compared with the thermodynamic pressure $p^\ast$ , i.e. $\displaystyle |\mathcal{K}^\ast \boldsymbol{\nabla} ^\ast \boldsymbol{\cdot }\boldsymbol{u}^\ast |\ll p^\ast$ . We note that the value of the surface tension $\sigma ^\ast$ will be derived in (2.25).

The heat transfer boundary condition at the interface $z^*=h^*(x^*,y^*,t^*)$ is expressed by Newton’s law of cooling

(2.16) \begin{equation} -K_{\textit{nf}}^*\boldsymbol{\nabla} ^\ast T^*\boldsymbol{\cdot }\boldsymbol{n}^\ast =\widehat {q}(T^*-T_\infty ^*), \end{equation}

where $\widehat {q}$ is the rate of heat transfer from the liquid to gas phase by convection. This boundary condition effectively couples the temperature and local particle concentration at the interface since $K_{\textit{nf}}^\ast$ depends on the latter.

Additionally, at the liquid–air interface, the interfacial concentration $\varGamma ^\ast$ defined hereon in units of mol m $^{-2}$ is given by (Machrafi Reference Machrafi2022)

(2.17) \begin{equation} \varGamma ^* = \vartheta \,\varGamma ^*_\infty , \end{equation}

where $\vartheta$ and $\displaystyle {\varGamma ^*_\infty = \vartheta _\infty \mathcal{L}_{\textit{np}}^*({\rho _{\textit{np}}^*}/{\mathcal{M}_{\textit{np}}^*}})$ are the surface coverage and the maximum surface concentration in units of mol m $^{-2}$ , respectively. For spherical particles, the maximum geometrical coverage fraction $\vartheta _\infty = 0.547$ , the characteristic length $\mathcal{L}^*_{\textit{np}} = d_p^*/3$ and the factor

(2.18) \begin{equation} \frac {\rho _{\textit{np}}^*}{\mathcal{M}_{\textit{np}}^*} = \left [\frac {4\pi ^2}{9\sqrt {2}}\left (\frac {d_p^*}{2}\right )^3N_A^*\right ]^{-1} \end{equation}

yields an equivalent bulk nanoparticle concentration volume fraction to the interfacial molar mass, in units of mol m $^{-3}$ , where $N_A^* = 6.02\times 10^{23}$ mol $^{-1}$ is Avogadro’s number. Note that nanoparticles are neither molecules nor atoms, yet they are treated as a polymer made of many monomers, a collection of chemically connected atoms or molecules. In what follows, the molar mass of a nanoparticle is given by (Machrafi Reference Machrafi2022)

(2.19) \begin{equation} \mathcal{M}_{\textit{np}}^* = \mathcal{M}_{\textit{np}}^{*\prime }f_{\textit{np}}\frac {\mathcal{V}_{\textit{np}}^*}{\mathcal{V}_{\textit{np}}^{*\prime }}, \end{equation}

where $\mathcal{M}_{\textit{np}}^{*\prime }$ and $\displaystyle {f_{\textit{np}} = {\pi }/{3\sqrt {2}}}$ represent the molar mass of one molecule and maximum three-dimensional nanoparticle packing density, e.g. face-centred cubic of spheres, respectively. Furthermore, the volumes of a nanoparticle and a spherical molecule are defined, respectively, as

(2.20) \begin{equation} \mathcal{V}_{\textit{np}}^{*} = \frac {4\pi }{3}\left (\frac {d_p^*}{ 2}\right )^3\quad \text{and}\quad \mathcal{V}_{\textit{np}}^{*\prime } = \frac {4\pi }{3}\ell _{\textit{np}}^{*3};\quad \ell _{\textit{np}}^{*} = \left (\frac {3\mathcal{M}_{\textit{np}}^{*\prime }}{4\pi \rho _{\textit{np}}^*N_A^*}\right )^{1/3}. \end{equation}

Here, $\ell _{\textit{np}}^{*}$ depicts the radius of a molecule; in a similar way, we define in what follows the effective radius of a base fluid molecule $\mathcal{L}_{\textit{bf}}^*$ .

The spatiotemporal evolution of the interfacial nanoparticle concentration $\varGamma ^*(x^*,y^*,t^*)$ at the deformable interface $z^* = h^*(x^*,y^*,t^*)$ is given by the extension of the conservation equation derived for the case without interfacial kinetics $\mathcal{J}_a^* = \mathcal{J}_d^*=0$ (Scriven Reference Scriven1960; Edwards et al. Reference Edwards, Brenner and Wasan1991; Wong, Rumschitzki & Maldarelli Reference Wong, Rumschitzki and Maldarelli1996; Pereira & Kalliadasis Reference Pereira and Kalliadasis2008; Manikantan & Squires Reference Manikantan and Squires2020):

(2.21a) \begin{eqnarray} \frac {\partial \varGamma ^\ast }{\partial t^*} &&- \frac {\partial h^*}{\partial t^*}\left (\boldsymbol{e}_{z^*}\boldsymbol{\nabla} _s^*\varGamma ^\ast \right ) +\boldsymbol{\nabla} _s^*\boldsymbol{\cdot }\left (\boldsymbol{u}_s^*\varGamma ^*\right ) + \varGamma ^*\left (\boldsymbol{\nabla} _s^*\boldsymbol{\cdot }\boldsymbol{n}^*\right )\left (\boldsymbol{u}^*\boldsymbol{\cdot }\boldsymbol{n}^*\right ) \nonumber \\ && =-\boldsymbol{\nabla} _s^*\mathcal{J}_{Ds}^* + \varGamma _\infty ^*\mathcal{S}_a^*\left [\mathcal{J}_a^* - \mathcal{J}_d^*\right ]. \end{eqnarray}

Here $\boldsymbol{u}_s^* = (\boldsymbol{I} - \boldsymbol{n}^\ast \boldsymbol{n}^\ast )\boldsymbol{\cdot }\boldsymbol{u}^*$ is the surface velocity component and the specific surface area covered by the interfacial nanoparticles per their unit mole $\mathcal{S}_a^*$ is given by

(2.21b) \begin{equation} \mathcal{S}_a^* = \frac {\mathcal{M}_{\textit{np}}^{*}}{\rho _{\textit{np}}^*\mathcal{L}_{\textit{np}}^*}. \end{equation}

The values $\mathcal{J}_{Ds}^*, \mathcal{J}_{a}^*$ and $\mathcal{J}_{d}^*$ are, respectively, the fluxes due to the surface diffusion, adsorption and the desorption. We express the surface diffusion flux via Fick’s law of diffusion with a constant surface diffusivity $D_s$ :

(2.21c) \begin{equation} \mathcal{J}_{Ds}^* = -D_s\boldsymbol{\nabla} _s^*\varGamma ^*. \end{equation}

We note that, similarly to the bulk diffusion coefficient, the surface diffusivity also exhibits the dependence on the interfacial concentration, as noted by Valkovska & Danov (Reference Valkovska and Danov2000) and Manikantan & Squires (Reference Manikantan and Squires2020). However, the variation of surface diffusion for the regime of moderate nanoparticle concentration is quite small. Thus, in what follows, we proceed with the constant surface diffusivity.

The adsorption/desorption flux characterises the nanoparticle exchange between the bulk and the air–liquid interface. We adopt the Langmuir isotherm model for a moderately dense colloidal suspension. The adsorption flux $\mathcal{J}_a^*$ is thus expressed by the adsorption of the bulk molar concentration $C^*$ to the available empty interfacial space $(1-\varGamma ^*/\varGamma _\infty ^*)$ with the adsorption rate $k_a^*$ , thereby, the adsorption flux is given by

(2.21d) \begin{equation} \mathcal{J}_a^* = k_a^* C^*\left (1-\frac {\varGamma ^*}{\varGamma ^*_\infty }\right ), \end{equation}

where ${C^* =({\rho _{\textit{np}}^*}/{\mathcal{M}_{\textit{np}}^*})\phi ^*}$ is the bulk molar nanoparticle concentration. The desorption flux $\mathcal{J}_d^*$ of the adsorbed nanoparticles per specific area with a desorption rate $k_d^*$ is given by

(2.21e) \begin{equation} \mathcal{J}_d^* = \frac {k_d^*}{\mathcal{S}_a^*}\left (\frac {\varGamma ^*}{\varGamma _\infty ^*}\right ). \end{equation}

Hence, upon substitution of (2.21c )–(2.21e ) into (2.21a ), we obtain the generalised interfacial spatiotemporal evolution equation for the interfacial nanoparticle concentration.

At the deformable interface, the total molar mass flux balances the total interfacial kinetic flux:

(2.22) \begin{equation} -\left (C_BT^*C_{z^*}^* +C_TC^*\frac {T_{z^*}^*}{T^*}\right )\boldsymbol{\cdot }\boldsymbol{n}^* = \mathcal{J}_a^* - \mathcal{J}_d^*. \end{equation}

Before we formulate the generalised surface tension variation model with heat and mass transfer, we begin with the formulation of the surface tension variation with the interfacial and bulk nanoparticle concentrations. The surface tension variation given by the Gibbs adsorption isotherm (Manikantan & Squires Reference Manikantan and Squires2020; Machrafi Reference Machrafi2022) is

(2.23) \begin{equation} \text{d}\sigma ^*\big(\phi ^*,\varGamma ^*\big) = -\varGamma _{e}^*\text{d}\,\varXi _{\textit{np}}^* (\vartheta ), \end{equation}

where the excess interfacial concentration $\varGamma _e^*$ is obtained by subtracting the maximum interfacial concentration in the occupied space and the surface equivalent imaginary bulk concentration, thence

(2.24a) \begin{equation} \varGamma ^*_e = \vartheta \varGamma _\infty ^*-\phi ^*\varGamma _b^* = \left (\vartheta K_\varSigma -\phi ^*\right )\varGamma _b^*, \end{equation}

where $\varGamma _b^* =( { {\mathcal{L}_{\textit{np}}^{*2}\rho _{\textit{np}}^*}/{\mathcal{L}_{\textit{bf}}^{*}\mathcal{M}_{\textit{np}}^*}})$ and $\displaystyle {K_\varSigma = \varGamma _\infty ^\ast /\varGamma _b^\ast } =({\vartheta _\infty \mathcal{L}_{\textit{bf}}^{*}/}{\mathcal{L}_{\textit{np}}^*})$ accounts for the ability of a nanoparticle to adsorb to the interface or stay in the bulk. The effective radius of the base fluid molecule is given by

(2.24b) \begin{align} \mathcal{L}_{\textit{bf}}^\ast = \frac {1}{3} \left (\frac {3\mathcal{M}_{\textit{bf}}^\ast }{4\pi \rho _{\textit{bf}}^\ast N_A^\ast } \right )^{1/3}, \end{align}

where $\mathcal{M}_{\textit{bf}}^\ast$ is the molar mass of the base fluid. Moreover, the surface chemical potential of the interfacial nanoparticle concentration is expressed by (Machrafi Reference Machrafi2022)

(2.24c) \begin{equation} \varXi _{\textit{np}}^* = K_B^*T_\infty ^*N_A^*\left [\ln (\vartheta )- \omega _{\textit{np}}\ln \left (1-\vartheta \right )\right ]. \end{equation}

Here, the parameter

(2.24d) \begin{equation} \omega _{\textit{np}} = \frac {\big(\rho _{\textit{bf}}^*/\mathcal{M}_{\textit{bf}}^*\big)\mathcal{L}_{\textit{bf}}^*}{\big(\rho _{\textit{np}}^*/\mathcal{M}_{\textit{np}}^*\big)\mathcal{L}_{\textit{np}}^*} \end{equation}

represents the ratio between the specific interfacial area per mole occupied by the interfacial nanoparticle concentration and the specific interfacial area per mole covered by the base fluid molecules.

The variation of the surface tension $\sigma ^\ast$ with temperature $T^*$ is assumed to be linear when the particle concentration at the interface is kept constant, i.e. $\displaystyle \sigma ^\ast =\sigma ^\ast _r - \sigma ^\ast _ {T^\ast } (T^\ast - T^\ast _r ),$ with $\sigma ^\ast _{T^\ast } \equiv -\partial \sigma ^\ast /\partial T^\ast \gt 0$ , where the subscript $r$ denotes some reference values. Along with the Gibbs isotherm relation, we derive the equation of state for the surface tension of a nanofluid with interfacial heat and mass transfer, by integrating (2.23) with respect to $\vartheta$ keeping $\phi ^*$ constant to obtain a general expression for $\sigma ^\ast$

(2.25) \begin{align} \sigma ^*(T^*,\varGamma ^*,\phi ^*)& = \sigma _r^* - \sigma _{T^*}^* (T^*-T_\infty ^*) + T_\infty ^* N_A^* \varGamma _b^* K_B^* \Biggl \{-K_{\varSigma } \frac {\varGamma ^*}{\varGamma _\infty ^*} + \phi ^* \ln \left (\frac {\varGamma ^*}{\varGamma _\infty ^*}\right ) \nonumber \\ &\quad +\, \omega _{\textit{np}} \biggl [\left (K_{\varSigma }-\phi ^* \right ) \ln \left (1-\frac {\varGamma ^*}{\varGamma _\infty ^*}\right )+K_{\varSigma } \frac {\varGamma ^*}{\varGamma _\infty ^*}\biggr ] \Biggr \}. \end{align}

Here, $\sigma ^\ast _r$ is the equilibrium reference value of the surface tension of the base fluid with the gas, and since $\sigma _{T^*}^*$ is positive, the surface tension decreases linearly with temperature. Interestingly, we note that the equation of state varies with both the bulk $\phi ^*$ and interfacial nanoparticle concentrations $\varGamma ^*$ . The surface tension exhibits a non-monotonic variation with the nanoparticle concentration, since for the low-concentration adsorption of nanoparticles from the bulk, the surface tension initially decreases with the concentration. However, with an increase in the concentration and depending on the values of the adsorption/desorption kinetics, the advancing nanoparticles fill the available sites, hence tending to increase the surface tension further.

We comment on the variation of surface tension with the nanoparticle concentration at equilibrium between the adsorption and desorption of nanoparticles. In equilibrium,

(2.26) \begin{equation} \begin{aligned} \mathcal{J}_a^* ={}& \mathcal{J}_d^*,\\ k_a^* \frac {\rho _{\textit{np}}^*}{\mathcal{M}_{\textit{np}}^*}\phi ^*\left (1-\frac {\varGamma ^*}{\varGamma ^*_\infty }\right ) = {}& \frac {k_d^*}{\mathcal{S}_a^*}\left (\frac {\varGamma ^*}{\varGamma _\infty ^*}\right ). \end{aligned} \end{equation}

Therefore, the dimensionless surface coverage $\displaystyle {\vartheta = \varGamma ^\ast /\varGamma _\infty ^\ast }$ is given by

(2.27) \begin{equation} \vartheta = \frac {\phi ^*k_a^*\mathcal{S}_a^*}{\left (\frac {\mathcal{M}_{\textit{np}}^*}{\rho _{\textit{np}}^*}\right )k_d^* + \phi ^*k_a^*\mathcal{S}_a^*}. \end{equation}

In what follows, we substitute (2.27) into (2.25), and deduce the variation of the equation of surface tension with the nanoparticle concentration $\phi ^*$ , $\sigma ^* (\phi ^* )$ , whose variation with the nanoparticle concentration is illustrated in figure 2. As mentioned above, the surface tension displays a minor drop at low nanoparticle concentrations, while it increases steadily with concentration thereafter. In figure 2 the circles and triangles denote the experimental data of surface tension variation for the nanofluid containing Al $_2$ O $_3$ nanoparticles in water obtained from Tanvir & Qiao (Reference Tanvir and Qiao2012) and Harikrishnan et al. (Reference Harikrishnan, Dhar, Agnihotri, Gedupudi and Das2017). We note that the experimental results display proximity to the theoretical model. A ‘mild decrease’ in the surface tension variation is observed at low nanoparticle concentrations, which is due to the low surface energetics of alumina nanoparticles in water, i.e. a low value of the ratio between the adsorption and desorption rates. On the other hand, for different base liquids, such as ethanol, n-decane, etc., a non-monotonic variation of surface tension with the nanoparticle concentration expressed by the equation of state is more pronounced in the cases of high surface energetics (Tanvir & Qiao Reference Tanvir and Qiao2012; Harikrishnan et al. Reference Harikrishnan, Dhar, Agnihotri, Gedupudi and Das2017; Machrafi Reference Machrafi2022).

Figure 2.

Variation of the surface tension $\sigma ^*(\phi ^*)$ Nm⁻¹ with the nanoparticle concentration $\phi ^*$ . The thick line $(\boldsymbol{-})$ , $\circ$ and $\triangle$ points represent the surface tension values obtained theoretically using (2.25) and (2.27), and from the experiments $1$ (Tanvir & Qiao Reference Tanvir and Qiao2012) and $2$ (Harikrishnan et al. Reference Harikrishnan, Dhar, Agnihotri, Gedupudi and Das2017), respectively. The inset shows a close-up of the surface tension variation in the domain of a low nanoparticle concentration.

Finally, to obtain the closure of the problem, one more condition must be satisfied. The average molar mass of the particles $({\rho _{\textit{np}}^*}/{\mathcal{M}_{\textit{np}}^*})\varPhi ^*$ . with $\varPhi ^\ast$ being the average volume concentration of particles, is conserved, therefore,

(2.28a) \begin{equation} \begin{aligned} \left (\frac {\rho _{\textit{np}}^*}{\mathcal{M}_{\textit{np}}^*}\right )I(\mathcal{D}^*)h_0^* \varPhi ^* &= \left [\int _{0}^{h_0^*}\left (\iint _{\mathcal{D}^*}\frac {\rho _{\textit{np}}^*}{\mathcal{M}_{\textit{np}}^*}\phi ^*\big(x^*,\,y^*,\,z^*,t^*\big)\text{d}x^*\text{d}y^*\right )\text{d}z^* \right. \nonumber \\ &\quad \left. +\, \iint _{{\mathscr{L}}^*}\varGamma ^*\big(x^*,\,y^*,\,t^*\big)\text{d}\textsf{S}^* \right ], \end{aligned} \end{equation}

with

(2.28b) \begin{equation} \text{d}\textsf{S}^* = \sqrt {1+h_{x^*}^{*2}+h_{y^*}^{*2}}\,\text{d}x^*\text{d}y^*, \end{equation}

where $\mathcal{D}^*$ is a projection of the flow domain onto the $x^\ast -y^\ast$ plane and $I({\mathcal{D}^*})$ is the area of this projection, whereas ${\mathscr{L}}^*$ represents the projection of the gas–liquid interface onto the $x^*-y^*$ surface. Note that the average particle concentration $\varPhi ^*$ is multiplied by the volume, i.e. $I(\mathcal{D}^*)h_0^*$ , which expresses the total mass of nanoparticles across the volume of a nanofluid layer (Edmonstone et al. Reference Edmonstone, Craster and Matar2006). This implies that the total mass of nanoparticles in the bulk and at the interface is conserved; it also suggests an inherent inhomogeneity of nanoparticle distribution in the bulk and at the interface. A needed closure of the problem is obtained by imposing a constraint (2.28a ) with a prescribed and fixed value of $\varPhi ^*$ . We note that the condition of the constant nanoparticle averaged bulk concentration $\varPhi ^*$ is appropriate for the case of the quiescent base state. However, in the presence of the base state flow, the conservation of the nanoparticle mass flux or the nanoparticle concentration current should be employed (Krishnan, Beimfohr & Leighton Reference Krishnan, Beimfohr and Leighton1996; Frank et al. Reference Frank, Anderson, Weeks and Morris2003; Ramachandran & Leighton Reference Ramachandran and Leighton2008; Morris Reference Morris2020).

2.3. Dimensionless formulation

We non-dimensionalise the governing equation and the boundary conditions with the following normalisation:

(2.29) \begin{align} \big (x^*,y^*,z^* \big )=h_0^\ast \left (x,y,z\right ), \quad h^\ast =h_0^\ast h,\quad {\boldsymbol{u}^*}=\frac {\kappa ^\ast _{\textit{bf}}}{h_0^\ast } {\boldsymbol{u}}, \quad t^*=\frac {h_0^{*2}}{\nu ^\ast _{\textit{bf}}} t,\quad \phi ^* = \phi _m \phi , \nonumber \\ p^*=p^\ast _\infty +\frac {\mu ^\ast _{\textit{bf}}\kappa ^\ast _{\textit{bf}}}{h_0^{*2}} p, \quad T^*=T^\ast _\infty +\frac {q^\ast h_0^\ast }{K^\ast _{\textit{bf}}}T,\quad \varGamma ^* = \varGamma _\infty ^* \varGamma , \quad \varPhi ^*=\phi _m\varPhi , \nonumber \\ K^\ast _{\textit{nf}}=K^\ast _{\textit{bf}}K_{\textit{nf}},\quad \rho _{\textit{nf}}^* = \rho ^\ast _{\textit{bf}}\rho _{\textit{nf}}, \quad \mu _{\textit{nf}}^*=\mu ^\ast _{\textit{bf}} \mu _{\textit{nf}},\quad \big(\rho ^\ast c^\ast \big)_{\textit{nf}} =\big(\rho ^\ast c^\ast \big) _{\textit{bf}} (\rho c)_{\textit{nf}}.\quad \end{align}

From here and on, the quantities without an asterisk represent the variables in dimensionless form. A set of dimensionless governing equations is written as

(2.30a) \begin{align} P(\rho _{\textit{nf}})_t+\boldsymbol{\nabla }\boldsymbol{\cdot }(\rho _{\textit{nf}}\boldsymbol{u}) &= 0, \end{align}

(2.30b) \begin{align} (\rho _{\textit{nf}}\boldsymbol{u})_t + P^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }(\rho _{\textit{nf}}\boldsymbol{u}\boldsymbol{u}) &= -\boldsymbol{\nabla }p +\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau } - \rho _{\textit{nf}} G \boldsymbol{e}_z, \end{align}

(2.30c) \begin{align} (\rho c)_{\textit{nf}}\left [PT_t + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })T\right ] &= \boldsymbol{\nabla }\boldsymbol{\cdot }( K_{\textit{nf}}\boldsymbol{\nabla }T) , \end{align}

(2.30d) \begin{align} PL^{-1}\phi _t + L^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u}\phi ) &= \boldsymbol{\nabla} ^2\phi + \boldsymbol{\nabla }\boldsymbol{\cdot }(\eta \phi \boldsymbol{\nabla }T), \end{align}

where $\displaystyle \boldsymbol{\nabla }\equiv ({\partial }/{\partial x}, {\partial }/{\partial y},{\partial }/{\partial z} )$ and $\displaystyle \boldsymbol{\tau } =\mu _{\textit{nf}} [\boldsymbol{\nabla }{\boldsymbol u}+ \boldsymbol{\nabla }{\boldsymbol u}^\intercal -{2}/{3} (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} ) \boldsymbol{I} ]$ is the viscous part of the dimensionless stress tensor.

The non-dimensional form of the boundary conditions at the substrate $z=0$ reads

(2.31a) \begin{equation} \boldsymbol{u} = 0,\quad K_{\textit{nf}} T_z = 1, \quad \phi _z + \eta \phi T_z =0. \end{equation}

At the deformable interface $z=h(x,y,t)$ , the boundary conditions are

(2.31b) \begin{align} Ph_t+uh_x &+vh_y=w, \end{align}

(2.31c) \begin{align} \boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{\mathcal{T}}\boldsymbol{\cdot }\boldsymbol{n} = 2\boldsymbol{\mathcal{H}}\sigma ,\quad \boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{\mathcal{T}}\boldsymbol{\cdot }\boldsymbol{x} &= \boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{\nabla} _s \sigma ,\quad \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\mathcal{T}}\boldsymbol{\cdot }\boldsymbol{y} = \boldsymbol{y}\boldsymbol{\cdot }\boldsymbol{\nabla} _s \sigma , \end{align}

where the dimensionless total stress tensor is $\boldsymbol{\mathcal{T}} = - p\boldsymbol{I} + \boldsymbol{\tau }$ and the normal and tangential vectors to the interface are

(2.31d) \begin{equation} \boldsymbol{n} = \frac {-h_{x}\boldsymbol{e}_{x}-h_{y}\boldsymbol{e}_{y}+\boldsymbol{e}_{z}}{\sqrt {1+h_{x}^{2}+h_{y}^{2}}},\quad \boldsymbol{x} = \frac {\boldsymbol{e}_{x}+h_{x}\boldsymbol{e}_z}{\sqrt {1+h_{x}^{2}}},\quad \boldsymbol{y} =\frac {\boldsymbol{e}_{y}+h_{y}\boldsymbol{e}_z}{\sqrt {1+h_{y}^{2}}},\quad \text{with}\quad 2\mathcal{H} = -\boldsymbol{\nabla} _s\boldsymbol{\cdot }\boldsymbol{n}, \end{equation}

and the surface tension is given by

(2.31e) \begin{equation} \sigma = \varSigma _0 - M_T T + E \left [\omega _{\textit{np}} \left [\left (K_{\varSigma }-\phi \phi _m \right ) \ln \left (1-\varGamma \right )+K_{\varSigma } \varGamma \right ]-K_{\varSigma } \varGamma + \phi \phi _m \ln \left (\varGamma \right )\right ]. \end{equation}

The boundary conditions for heat and mass transfer at the deformable interface $z=h(x,y,t)$ are

(2.31f) \begin{align} K_{\textit{nf}}\boldsymbol{\nabla }T\boldsymbol{\cdot }\boldsymbol{n} & +BT=0, \end{align}

(2.31g) \begin{align} P\big [\varGamma _{t}- h_{t}(\boldsymbol{e}_z\boldsymbol{\cdot }\boldsymbol{\nabla} _s )\varGamma \big ]+\boldsymbol{\nabla} _s\boldsymbol{\cdot }(\boldsymbol{u}_s\varGamma ) &+\varGamma (\boldsymbol{\nabla} _s\boldsymbol{\cdot }\boldsymbol{n})(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}) = L_S\big (\boldsymbol{\nabla} _s^{2}\varGamma \big ) \nonumber\\ & + \mathcal{B}_A \big[K_{\textit{ad}}\phi (1-\varGamma ) - \varGamma \big ] , \end{align}

(2.31h) \begin{align} -\left (\boldsymbol{\nabla }\phi + \eta \phi \boldsymbol{\nabla }T\right )\boldsymbol{\cdot }\boldsymbol{n} &= \mathcal{B}_D\left [K_{\textit{ad}}\phi (1-\varGamma ) - \varGamma \right ]. \end{align}

The problem possesses a long list of parameters, and the dynamics of the system is described by the following dimensionless numbers:

(2.32a) \begin{equation} \begin{aligned} & P=\frac {\nu ^\ast _{\textit{bf}} \rho ^\ast _{\textit{bf}} c_{\textit{bf}}^\ast }{K^\ast _{\textit{bf}}}, \quad G=\frac {g^\ast h_0^{*3}}{\nu ^\ast _{\textit{bf}}\kappa ^\ast _{\textit{bf}}}, \quad L=\frac {C_B T^\ast _\infty \rho ^\ast _{\textit{bf}} c_{\textit{bf}}^\ast }{K^\ast _{\textit{bf}}}, \quad L_s=\frac {D_s \rho ^\ast _{\textit{bf}} c_{\textit{bf}}^\ast }{K^\ast _{\textit{bf}}}, \nonumber \\ & B=\frac {\widehat {q}h_0^\ast }{K^\ast _{\textit{bf}}},\qquad \varSigma _0 = \frac {\sigma ^\ast _rh_0^\ast }{\mu ^\ast _{\textit{bf}}\kappa ^\ast _{\textit{bf}}},\qquad \eta = \frac {C_T q^\ast h_0^\ast }{C_BT_\infty ^{\ast 2} K^\ast _{\textit{bf}}}. \end{aligned} \end{equation}

Here, $P,$ $G,$ $L,$ $L_s,$ $B,$ $\varSigma _0$ and $\eta ,$ are the Prandtl, modified Galileo, Lewis, surface Lewis, Biot, dimensionless surface tension numbers and the Soret coefficient, respectively. The convective heat transfer induces the interfacial surface tension gradient, which is quantified by the dimensionless thermal Marangoni number

(2.32b) \begin{equation} M_T = \frac {\sigma _{T^*}^*\Delta T^{*}h_0^*}{\mu _{\textit{bf}}^*\kappa _{\textit{bf}}^*}, \quad \text{where}\quad \Delta T^*\equiv \frac {q^\ast h_0^\ast }{K^\ast _{\textit{bf}}}. \end{equation}

The interfacial adsorption and desorption processes are controlled by

(2.32c) \begin{equation} K_{\textit{ad}} = \frac {k_a^* \rho _{\textit{np}}^*\phi _m\mathcal{S}_a^*}{\mathcal{M}_{\textit{np}}^*k_d^*},\quad \mathcal{B}_D = \frac {k_d^*h_0^*\mathcal{M}_{\textit{np}}^{*}}{C_BT_\infty ^*\phi _m\mathcal{S}_a^*\rho _{\textit{np}}^*},\quad \mathcal{B}_A = \frac {k_d^*h_0^{*2}}{\kappa _{\textit{bf}}^*}, \end{equation}

depicting the ratio of the adsorption to desorption rate, the surface diffusion and advection Biot numbers, respectively. The parameters $\mathcal{B}_D$ and $\mathcal{B}_A$ denote the ratio of the desorption flux to the diffusion flux and to the advection flux, respectively. We finally note that the ratio between the two surface Biot numbers yields the dimensionless Péclet number $ Pe=\mathcal{B}_D/\mathcal{B}_A$ .

In what follows, we define the dimensionless quantity, the Elasticity number

(2.33) \begin{equation} E = \frac {h_0^*K_B^*N_A^*T_\infty ^*\varGamma _b^*}{\mu _{\textit{bf}}^*\kappa _{\textit{bf}}^*}, \end{equation}

which describes the ratio of interfacial elasticity moduli at a fixed $\varGamma _b^*$ to viscous stress (Stebe et al. Reference Stebe, Lin and Maldarelli1991; Seiwert, Dollet & Cantat Reference Seiwert, Dollet and Cantat2014; Manikantan & Squires Reference Manikantan and Squires2020). Keeping in mind (2.25) for the surface tension implies that the Elasticity number intrinsically quantifies the variation of surface tension with interfacial nanoparticle concentration $ \partial \sigma /\partial \varGamma$ and with bulk nanoparticle concentration $\partial \sigma /\partial \phi$ . In fact, there exist two dynamic solutal Marangoni ‘numbers’ (more precisely, Marangoni functions), namely, the surface concentration Marangoni number

(2.34a) \begin{equation} M_\varGamma (\varGamma ,\phi ) \equiv \frac {\partial \sigma }{\partial \varGamma } = \frac { E (1 -\varGamma +\varGamma \omega _{\textit{np}} ) (\varGamma K_{\varSigma }-\phi \phi _m )}{(\varGamma -1) \varGamma } \end{equation}

and the bulk concentration Marangoni number

(2.34b) \begin{equation} M_\phi (\varGamma ) \equiv \frac {\partial \sigma }{\partial \phi } = E \phi _m \left [\ln (\varGamma )-\omega _{\textit{np}}\ln (1-\varGamma ) \right ]. \end{equation}

It is interesting to note that the surface concentration Marangoni number (2.34a ) exhibits dependence on both the interfacial nanoparticle concentration and the bulk nanoparticle concentration; however, the bulk concentration Marangoni number (2.34b ) depends only on the interfacial nanoparticle concentration. Remarkably, this signifies that the sign and amplitude of the surface tension gradient determine the mode of solutal Marangoni instability.

Figure 3.

Variation of solutal Marangoni numbers $M_\phi$ and $M_\varGamma$ with nanoparticle concentration $\phi$ obtained from (2.31e ) at equilibrium determined by (2.27) at $E\approx 2.53,$ $K_{\textit{ad}} = 3.84,$ $K_\varSigma =2.11\times 10^{-3}$ and $\omega _{\textit{np}}=6.23\times 10^3$ . The solid and dashed curves represent the variation of the surface concentration Marangoni number $M_\varGamma$ and the bulk concentration Marangoni number $M_\phi$ , respectively.

Figure 3 displays the variations of the surface tension surface and bulk concentration gradients with the nanoparticle concentration obtained from (2.31e ) at an equilibrium interfacial nanoparticle concentration

(2.35) \begin{equation} \varGamma = \phi K_{\textit{ad}}/\left (1+\phi K_{\textit{ad}}\right ), \end{equation}

following from (2.27). We observe that $M_\varGamma$ increases linearly with the bulk nanoparticle concentration $\phi$ . The positivity of $M_\varGamma$ for all values of the nanoparticle concentration suggests that the system inhibits a possible intrinsic solutocapillary instability due to the surface tension gradient of interfacial nanoparticle concentration $M_\varGamma$ . However, we note that $M_\phi$ is negative for low values of the bulk concentration $\phi$ and positive otherwise.

We further emphasise that the variation of surface tension due to the solutal effects is given by

(2.36) \begin{equation} \Delta \sigma = M_\varGamma \Delta \varGamma + M_\phi \Delta \phi = M_\varGamma \left ( \Delta \varGamma + \frac {M_\phi }{M_\varGamma } \Delta \phi \right ). \end{equation}

Note that, since $\Delta \varGamma$ and $\Delta \phi$ may be independently of either sign, it is impossible to claim whether $\Delta \sigma$ given by (2.36) leads to the emergence of the solutocapillary instability before obtaining the solution of the relevant EVP; see § 3.

It follows from figure 3 that there exists an intersection point $\tilde {M}$ corresponding to $\phi = \tilde {\phi }$ , so that $M_{\varGamma } \gt M_{\phi }$ for $\phi \lt \tilde {\phi }$ and $M_{\varGamma } \lt M_{\phi }$ for $\phi \gt \tilde {\phi }$ . The Elasticity number $E$ reflects the impact of the variation of the surface tension with the bulk and interfacial concentration of particles under non-equilibrium conditions, and its value as an eigenvalue of the EVP referred to above is associated with the possible onset of the solutocapillary instability.

The adsorption/desorption rates in the case of alumina nanoparticles in water can be derived using the interaction between a nanoparticle and the layer interface, as outlined by Ruckenstein & Prieve (Reference Ruckenstein and Prieve1976), Grow & Shaeiwitz (Reference Grow and Shaeiwitz1982) and Machrafi (Reference Machrafi2022) and detailed in Appendix A. We find the adsorption and desorption rates for Al $_2$ O $_3$ nanoparticles of diameter $d_p^*=50$ nm to be $k_a^* = 2.27\times 10^{-4}$ ms⁻¹ and $k_d^* = 2.3\times 10^3$ s $^{-1}$ , respectively. Thereby, we obtain the values of the surface Biot numbers for a thin nanofluid layer $h_0^*=10^{-6}$ m as $\mathcal{B}_D = 6.7236$ and $\mathcal{B}_A = 0.0157$ . The representative Péclet number obtained by the ratio of the surface Biot number yields $ Pe \approx 426$ , which implies that the advective flux dominates the diffusive flux; hence, we refer to this case as a high-Péclet-number flow $ Pe \gg 1$ . Moreover, for a thicker layer $h_0\approx 10^{-4}$ m or for a layer with a higher desorption rate $(k_d^* \approx 10^5)$ , the Péclet number reduces significantly to the case of $ Pe \sim O(1)$ . Besides, there also exists a regime of $ Pe \ll 1$ , where the diffusion is faster than the advection (Manikantan & Squires Reference Manikantan and Squires2020). The other parameter ranges for the investigation of the current paper are provided in table 1.

Table 1.

Parameter nomenclature and their typical values used in this investigation.

The dimensionless thermophysical properties appearing in the set of non-dimensional governing equations and boundary conditions (2.30) and (2.31) are defined as

(2.37) \begin{equation} \begin{aligned} \rho _{\textit{nf}}(\phi ) =(1-\phi \phi _m)+\phi \phi _m\rho _{\textit{np}},\quad (\rho c)_{\textit{nf}}(\phi ) = (1-\phi \phi _m)+ \phi \phi _m(\rho c)_{\textit{np}},\\ K_{\textit{nf}}(\phi ) = 1 + a \phi \phi _m,\quad \mu _{\textit{nf}}(\phi ) = (1-\phi )^{-2}. \end{aligned} \end{equation}

We finally express the dimensionless form of the conservation of the total particle volume across the nanofluid layer. Hence, following (2.28), the dimensionless conservation of the total nanoparticle volume across the layer is given by

(2.38a) \begin{equation} I(\mathcal{D})\varPhi = \int _{0}^{1}\left (\iint _{\mathcal{D}}\phi (x,\,y,\,z,\,t)\text{d}x\,\text{d}y\right )\text{d}z + \mathcal{A} \iint _{{\mathscr{L}}}\varGamma (x,\,y,\,t)\, \text{d}\textsf{S} , \end{equation}

with

(2.38b) \begin{equation} \text{d}\textsf{S} = \sqrt {1+h_{x}^{2}+h_{y}^{2}}\,\text{d}x\,\text{d}y, \end{equation}

where ${\mathcal{A} = ({\mathcal{M}_{\textit{np}}^*\varGamma _\infty ^*}/{\rho _{\textit{np}}^*h_0^*\phi _m}})$ represents the degree of particle adsorption factor at the liquid–gas interface.

2.4. Base state solution

We now obtain the quiescent $(\boldsymbol{u}=\boldsymbol{u}_0 = 0)$ steady state of the system with the flat interface $z=1$ in terms of the components of temperature $T_0(z)$ , bulk nanoparticle concentration $\phi _0(z)$ , pressure $p_0(z)$ and interfacial particle concentration $\varGamma _0$ , which depend only on the vertical coordinate $z$ . We note that for the base state, the variables are strongly coupled via the concentration-dependent thermophysical properties; therefore, the appropriately reduced (2.30) and (2.31) need to be solved simultaneously, i.e.

(2.39a) \begin{align} p_{0,z}= - G -G \phi _m (\rho _{\textit{np}}-1) \phi _0, \end{align}

(2.39b) \begin{align} \left [(1+a \phi _m \phi _0)T_{0,z}\right ]_z =0, \end{align}

(2.39c) \begin{align} \left (\phi _{0,z} + \eta \phi _0 T_{0,z}\right )_z=0, \end{align}

with the following boundary conditions:

(2.40a) \begin{align} \left (1+a\phi _m\phi _0\right )T_{0,z} = 1,\qquad \phi _{0,z} + \eta \phi _0 T_{0,z} =0\quad \text{at}\quad z=0, \end{align}

(2.40b) \begin{align} p_0 =0, \qquad (1+a \phi _m \phi _0)T_{0,z} +B T_0=0, \qquad \phi _{0,z} + \eta \phi _0 T_{0,z}=0\quad \text{at}\quad z=1. \end{align}

At the equilibrium, the adsorption and desorption kinetic fluxes balance each other; therefore, it follows from (2.31h ) that

(2.40c) \begin{equation} K_{\textit{ad}}\phi _0(1)(1-\varGamma _0)-\varGamma _0 = 0. \end{equation}

The base state of the system satisfying (2.39) with boundary conditions (2.40) reads

(2.41a) \begin{align} T_0 &= \frac {1}{\eta B}\biggl [B \ln \left (\frac {1}{a\phi _m}\mathcal{W}\left (\gamma (\varPhi ) a\phi _m \exp {\left (\gamma (\varPhi ) a\phi _m - \eta \right )}\right )\right )\nonumber\\ & -B \ln \left (\frac {1}{a\phi _m}\mathcal{W}\left (\gamma (\varPhi ) a\phi _m \exp {\left (\gamma (\varPhi ) a\phi _m - \eta z\right )}\right )\right ) - \eta \biggr ], \end{align}

(2.41b) \begin{align} \phi _0 &=\frac {\mathcal{W}\left (\gamma (\varPhi ) a\phi _m \exp {\left (\gamma (\varPhi ) a\phi _m - \eta z\right )}\right )}{a\phi _m}, \end{align}

(2.41c) \begin{align} \varGamma _0 &= \frac {K_{\textit{ad}} \mathcal{W}\left (\gamma (\varPhi ) a\phi _m \exp \left (\gamma (\varPhi ) a\phi _m - \eta \right )\right )}{K_{\textit{ad}} \mathcal{W}\left (\gamma (\varPhi ) a\phi _m \exp \left (\gamma (\varPhi ) a\phi _m - \eta \right )\right ) + a\phi _m}, \end{align}

(2.41d) \begin{align} p_0 &= G (1-z) - \frac {1}{2 \eta a}\biggl \{G (\rho _{\textit{np}}-1) \biggl [\mathcal{W}\left (\gamma (\varPhi ) a\phi _m \exp {\left (\gamma (\varPhi ) a\phi _m - \eta \right )}\right ) \nonumber \\ & -\mathcal{W}\left (\gamma (\varPhi ) a\phi _m \exp {\left (\gamma (\varPhi ) a\phi _m - \eta z\right )}\right )\biggr ]\times \biggl [\mathcal{W}\left (\gamma (\varPhi ) a\phi _m \exp {\left (\gamma (\varPhi ) a\phi _m - \eta \right )}\right ) \nonumber \\ & +\mathcal{W}\left (\gamma (\varPhi ) a\phi _m \exp {\left (\gamma (\varPhi ) a\phi _m - \eta z\right )}\right )+2\biggr ]\biggr \}, \end{align}

where $\gamma (\varPhi )$ is a constant of integration to be determined below and $\mathcal{W}(z)$ is the principal branch of the Lambert W function (Corless et al. Reference Corless, Gonnet, Hare, Jeffrey and Knuth1996; Vallis, Parker & Tobias Reference Vallis, Parker and Tobias2019; Gandhi et al. Reference Gandhi, Nepomnyashchy and Oron2025c ) given by

(2.42) \begin{equation} z={\mathcal{W}}(z) \exp {\mathcal{W}}(z). \end{equation}

Finally, following (2.38), the conservation of the total particle volume in the bulk and at the interface is imposed by specifying the constant value of the average concentration of the nanoparticles $\varPhi$ , i.e.

(2.43) \begin{equation} \varPhi = \int _{0}^{1}\phi _0(z)\text{d}z + \mathcal{A} \varGamma _0, \end{equation}

whereby, we deduce the value of $\gamma (\varPhi )$ by substituting (2.41b ) and (2.41c ) into (2.43). We note that the value of the constant $\gamma (\varPhi )$ depends on the system parameters $a$ , $\phi _m$ , $\eta ,$ $K_{\textit{ad}}$ , $\varPhi$ , $\mathcal{A}$ and $\rho _{\textit{np}}^*/\mathcal{M}_{\textit{np}}^*$ .

Figure 4.

Variation of the quiescent base state in terms of (a) the bulk concentration $\phi _0$ , and (b) temperature components presented as $T_0(z)-T(z=1)$ with the transverse coordinate $z$ , with $B=0.01,$ $\eta = 1.66,$ $a=7.47,$ $\varPhi = 0.01$ and different values of the temperature difference $\Delta T^*$ . The inset in panel (a) represents the variation of the interfacial nanoparticle concentration $\varGamma _0$ with the Soret effect $\eta$ .

Figure 4 illustrates the variation of the base state components in terms of the nanoparticle concentration in the bulk $\phi _0(z)$ in figure 4(a) and the temperature $T_0(z)$ in figure 4(b) with the coordinate $z$ at different values of the temperature difference $\Delta T^*$ in the layer.

We note that the Soret effect with a positive Soret coefficient drives the nanoparticle flux from the hot interface $z=h$ to the cold substrate $z=0$ , thus creating a stable density stratification profile. We find that particle concentration profiles of the base state exhibit a nonlinear variation with height $z$ , whereas for low values of $\Delta T^\ast$ , the temperature profiles of the base state in figure 4(b) display an almost linear variation with  $z$ . The nonlinear variation of $T_0(z)$ becomes visible for higher values of the temperature difference $\Delta T^*$ (not shown in figure 4 b). Also, in contrast to the concentration profiles $\phi _0(z)$ , the temperature $T_0(z)$ remains almost independent of the temperature difference $\Delta T^*$ . The negative values of the temperature difference $T_0(z)-T_0(z=1)$ express the fact that the values of the temperature in the bulk are lower than $T_\infty ^*$ . Furthermore, we note that the interfacial adsorption/desorption kinetics expressed by $K_{\textit{ad}}$ describes the relative nanoparticle adsorption to desorption ratio from/to the layer adjacent to the interface. The inset of figure 4(a) displays the variation of the interfacial concentration $\varGamma _0$ component of the base state with the Soret coefficient $\eta$ . We find that for the nanofluid layer cooled at the substrate, the Soret effect and the desorption kinetics complement each other; therefore, an increase in the Soret coefficient $\eta$ via an increase in the temperature difference $\Delta T^*$ substantially reduces the interfacial concentration $\varGamma _0$ .

We note that the mathematical framework presented in this paper formally addresses a moderately dense nanofluid, i.e. $\varPhi \leqslant 0.1$ and $\phi _0(z) \le 0.1$ for all $0 \le z \le 1$ . As observed in figure 4(a), due to the presence of the Soret effect, the particle concentration near a cool bottom is the highest, so the value of $\phi _0(0)$ is higher than the average particle concentration $\varPhi$ by more than three times for the set of parameters given. Therefore, to rule out the emergence of high nanoparticle concentrations near the bottom $z=0$ , we have to set a specific limitation, hereafter referred to as limitation (i), for our theory to be valid. Another limitation, hereafter referred to as limitation (ii), should be set for the upper bound of the temperature drop $\Delta T^*$ across the layer to exclude the temperature sensitivity of the thermophysical properties of the nanofluid. To first address limitation (ii), we assume that the temperature drop across the layer $\Delta T^*$ should not exceed $10$ K, thus limiting the heat flux $q^*$ from above accordingly. Now, with this upper bound at hand, we seek the conditions to satisfy limitation (i). We note that the upper bound for limitation (i) depends on the full set of parameters of the system. For example, for the fixed parameter set of $a=7.47,\, \phi _m=0.65,\, K_{\textit{ad}} = 3.84,\, \mathcal{A} = 0.014$ , $\rho _{\textit{np}}^*/\mathcal{M}_{\textit{np}}^* = 0.034$ along with $d_p^* = 50$ nm and the average particle concentration $\varPhi =0.01$ considered here, the temperature drop across the nanofluid layer should be $\Delta T^*\leq 7$ K, which corresponds in this case to $\phi _0(0)\leqslant 0.095$ . Similarly, in terms of $\varPhi$ , the upper bound due to limitation (i) should increase with a decrease in $\Delta T^*$ . For instance, for the same set of parameters, as specified above, with the temperature drop of $\Delta T^*=1$ K, the average particle concentration is limited to $\varPhi \leqslant 0.05$ . Accordingly, an upper bound for $\Delta T^*$ imposes an upper bound for the admissible thermal Marangoni number $M_T$ . For example, for a nanofluid layer considered here with thickness $h_0^*=10^{-6}$ m, the admissible values of the thermal Marangoni number must satisfy $M_T\leq 3.3$ .

3.1. Numerical procedure and the incompressibility condition

We further simplify the EVP (3.10) by substituting the expression of the perturbed longitudinal velocity component $u$ from (3.10a ) into the rest of the equations of the EVP (3.10). We also extract the expression of the perturbed deformation $\zeta$ from (3.10g ) and substitute it into the rest of the boundary conditions at $z=1$ . We note that the system consists of two concentration fields, the bulk concentration $\varphi$ and the interfacial concentration $\varGamma$ , to accommodate the numerical implementation of the EVP, we combine the boundary condition for the bulk mass flux balance equation (3.10j ) with the evolution boundary condition for the interfacial concentration equation (3.10k ), and obtain an explicit expression of the disturbances of the interfacial concentration $\varGamma$ . In what follows, we substitute the obtained expression of $\varGamma$ into the rest of the boundary conditions at $z=1$ . Hence, our simplified EVP reduces to contain only three dependent perturbed variables $w,$ $\theta $ and $\varphi$ .

The solution of the linear EVP gives us the threshold values of the instabilities. We solve the EVP and its simplifications using open-source code (Pearce et al. Reference Pearce, Heil, Jensen, Jones and Prokop2018). It is based on the shooting method and the numerical construction of the complex-valued Evans function (Evans & Feroe Reference Evans and Feroe1977) using the compound matrix method (Ng & Reid Reference Ng and Reid1979) whose roots represent the eigenvalues of the given problem. The solution is initiated by specifying the range in which the roots of the Evans function are searched and adjusted if needed. More details related to the numerical implementation are provided in Gandhi et al. (Reference Gandhi, Nepomnyashchy and Oron2025c ); Appendix C presents a validation of the numerical method by comparing the results of the fully numerical solution of the EVP with the analytical solutions of the toy problem of § 4.3.

We note that the density of the nanofluid system varies with space and time and also depends on the local nanoparticle concentration, which in turn varies spatiotemporally. Hence, under this consideration, the system is ‘compressible’. As in Gandhi et al. (Reference Gandhi, Nepomnyashchy and Oron2025c ), a further simplification of the problem considered here can be achieved. We rewrite the non-dimensional form of the nanofluid continuity equation (2.30a ) as

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

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

(3.12) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = -L\big (\rho _{\textit{np}}-1\big ) \phi _m\big [\boldsymbol{\nabla} ^2\phi + \boldsymbol{\nabla }\boldsymbol{\cdot } (\eta \phi \boldsymbol{\nabla }T )\big ]. \end{equation}

The right-hand side of (3.12) represents the deviation of the full ‘compressible’ system from an incompressible one. For nanosuspension with small concentration, e.g. $\varPhi \leqslant 0.05$ , and the small Lewis number $L\in (10^{-5},10^{-3} )$ , the right-hand side of (3.12) effectively vanishes and yields $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}\approx 0$ . This implies that under these conditions, the system at hand can be treated as ‘incompressible’. In our previous studies Gandhi et al. (Reference Gandhi, Nepomnyashchy and Oron2025c ), we validated the underlying incompressibility assumption via a direct comparison between the solution of the full EVP in its compressible and incompressible forms. Henceforth, we solve the EVP (3.10) under the assumption of an incompressible fluid system.