Thermosolutal instabilities in a moderately dense nanoparticle suspension

J. Fluid Mech. (2025), vol. 1011, A52, doi:10.1017/jfm.2025.374 Thermosolutal instabilities in a moderately dense nanoparticle suspension Raj Gandhi1, Alexander Nepomnyashchy1and Alexander Oron2 1Department of Mathematics, Technion - Israel Institute of Technology, Haifa 3200003, Israel 2Department of Mechanical Engineering, Technion - Israel Institute of Technology, Haifa 3200003, Israel Corresponding author: Raj Gandhi, raj.g[email protected]hnion.ac.il (Received 2 August 2024; revised 27 March 2025; accepted 31 March 2025) We investigate the onset of thermosolutal instabilities in a moderately dense nanoparticle suspension layer with a deformable interface. The suspension is deposited on a solid substrate subjected to a specified constant heat flux. The Soret effect and the action of gravity are taken into account. A mathematical model for the system considered with nanoparticle concentration-dependent density, viscosity, thermal conductivity and the Soret coefficient is presented in dimensional and non-dimensional forms. Linear stability analysis of the obtained base state is carried out using disturbances in the normal mode, and the corresponding eigenvalue problem is derived and numerically investigated. The onset of various instabilities is investigated for cases of both heating and cooling at the substrate. The monotonic solutocapillary instability is found in the case of cooling at the substrate, which exhibits two competing mechanisms that belong to two different disturbance wavelength domains. We identify the occurrence of both monotonic and oscillatory thermocapillary instabilities when the system is heated at the substrate. Furthermore, we show the emergence of the solutal buoyancy instability due to density variation which is promoted by the Soret effect adding nanoparticles heavier than the carrier fluid in the proximity of the layer interface. Transitions from the monotonic to oscillatory thermocapillary instability are found with variation in the gravityand solutocapillarity-related parameters. Notably, we identify a previously unknown transition from monotonic to the oscillatory thermocapillary instability due to the variation in the strength of the thermal-conductivity stratification coupled with the Soret effect. Key words: Marangoni convection, coupled diffusion and flow © The Author(s), 2025. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/ licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited. 1011 A52-1 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron 1. Introduction Surface-tension-driven phenomena in pure simple-liquid systems are ubiquitous in nature, as well as in industrial applications. For example, tears of wine (Scriven & Sternling 1960), the footprint of a whale caused by the sweeping of the biomaterial to the surface (Levy et al. 2011), crystal growth techniques (Schwabe et al. 1978) and in experiments with heated fluid layers in microgravity (Smith & Davis 1983a,b; Schatz & Neitzel 2001)are some of these, just to mention a few. The inhomogeneity of scalar fields at the layer interface such as temperature that affects the local surface tension may induce, under certain conditions, a convective motion in a liquid layer, known in the literature as the Marangoni or thermocapillary instability (Davis 1987). Pearson (1958) was the first to theoretically investigate the onset of the monotonic Marangoni instability, which may be either long-wave or short-wave, by considering a layer of a pure Newtonian liquid with the non-deformable interface which is heated at its solid support, provided that the temperature drop across the layer is sufficiently large to overcome the dissipative properties of the liquid such as viscosity and thermal diffusivity. Scriven & Sternling (1964) later found that interfacial deflection significantly alters the stability boundary and may lead to the onset of the oscillatory Marangoni instability. However, their result was further corrected by Smith (1966) who showed the stabilisation of the long-wave gravitational waves predicted by Scriven & Sternling (1964). In a similar way, in an isothermal fluid layer, inhomogeneities of the local solute concentration at the layer interface also affect the local surface tension and may lead to the emergence of the solutocapillary instability (Davies & Rideal 1963). In a pure liquid layer, heating at the layer interface or, equivalently, cooling at the substrate, provides a stabilising mechanism (Deissler & Oron 1992;Oron&Rosenau1992; Alexeev & Oron 2007), which may lead to a full stabilisation or saturation at the nonlinear stage. The suppression of the Rayleigh–Taylor instability in an inverted air–oil system by the thermocapillarity was demonstrated experimentally by Burgess et al. (2001). In liquid mixtures such as binary mixtures, the Soret effect (Cross & Hohenberg 1993; Skarda, Jacqmin & McCaughan 1998) introduces an additional component to the expression of Fick’s law for the mass flux normally related to the gradient of the bulk solute concentration. It is associated with the temperature gradient in the mixture and is found to be important. By coupling between heat and mass transfer effects, the Soret effect by itself may affect and modify thermal instabilities in a layer of a binary mixture (Oron & Nepomnyashchy 2004). Various aspects of longand finite-wave thermosolutocapillary instabilities in dilute binary mixtures heated at the substrate have already been investigated (Oron & Nepomnyashchy 2004; Podolny, Oron & Nepomnyashchy 2005; Shklyaev et al. 2007,2009; Bestehorn & Borcia 2010; Podolny, Nepomnyashchy & Oron 2010; Morozov, Oron & Nepomnyashchy 2014). Sarma & Mondal (2021a) studied thermosolutocapillary instabilities in a viscoelastic binary fluid with a particular case of a Newtonian binary fluid at zero Deborah number. In all of these papers, the emergence of both monotonic and oscillatory instabilities was reported. The instability mechanism for such a system depends on the direction of heating or cooling. Joo (1995) revealed the onset of the monotonic solutocapillary instability in a layer of a binary mixture with heating imposed at the free interface. Furthermore, it was noted that oscillatory instability takes place due to the competition between the stabilising solutocapillary and destabilising thermocapillary instability in the case of heating at the substrate (Joo 1995). Sarma & Mondal (2021b) investigated thermosolutal Marangoni instability in a layer of a viscoelastic binary fluid heated at the interface. They showed the emergence of a long-wave monotonic instability in the cases of both a deformable and non-deformable layer interface, and demonstrated that this instability is driven by solutocapillarity. 1011 A52-2 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics Colloidal dispersions made of a mixture of a base fluid and nanoparticles of diameter d∗ pbetween 1 and 100 nm are known as nanofluids. The material properties of a nanofluid such as density, kinematic viscosity, thermal diffusivity and Soret diffusion coefficient naturally depend on the local bulk concentration of the particles (Buongiorno 2006), as well as the Brownian diffusion coefficient (Batchelor 1976). This fact introduces obvious challenges affecting theoretical work and practical applications such as inkjet printing (Lohse 2022), paint coating and microgravity experiments (Vailati et al. 2023). In heattransfer related technological applications, metallic particles, for instance, alumina, copper oxide, silica, titania, etc., are purposely employed to enhance the thermal conductivity of the nanofluid relative to that of the base fluid (Choi & Eastman 1995). There are also promising technological applications such as nanofluid fuel (Abramzon & Sirignano 1989; Basu & Miglani 2016), where thermosolutal Marangoni stresses are of importance (Vang &Shaw2020;Shaw2022). The thermophysical stratification in such systems introduces even more complex features into mathematical modelling and investigations due to that. It is important to note that instabilities may also be triggered in a simple fluid layer featuring a non-uniformity in one or more of the physical properties of the fluid. For instance, if the fluid density in a horizontal layer in the gravity field increases with height, a situation where a heavier fluid is above the lighter one, the system may become unstable (Rayleigh 1882; Chandrasekhar 1961). An analogue to this instability may arise in a nanofluid layer due to the presence of nanoparticles heavier than the carrier liquid in its upper stratum. We refer to this instability as to the solutal buoyancy instability since the local density of the fluid depends on the local particle concentration. The presence of the Soret effect with a positive thermodiffusivity coefficient promotes the formation of an unstable nanoparticle concentration stratification across the layer, and thereby enhances the possibility of the onset of the solutal buoyancy instability. Although off the scope of the current paper, in liquid metal batteries (Herreman et al. 2020), the onset of the solutal buoyancy instability during the charge phase was found due to the emergence of the unstable stratification of lithium. Interestingly, Herreman et al. (2020) found that the onset of solutal buoyancy convection actually helps to homogenise the alloy layer, and the same physical effect introduces complexities during the discharge phase by creating the stable stratification (Herreman et al. 2021). Solutal buoyancy instability was also found to create convective flow by dissolution from a soluble solid into a fluid (Berhanu et al. 2021). Furthermore, the viscosity of a nanofluid increases with the local particle concentration and can be approximated via an empirical model (Maron & Pierce 1956). We note that in shear-induced flows, additional contributions to the normal stress in a nanofluid may be important in the presence of a base flow (Phillips et al. 1992; Dhas & Roy 2022; Lavrenteva, Smagin & Nir 2024), but in the case considered in the current paper, the base state is quiescent, and hence, these effects may be safely omitted. Experimental data show that the thermal conductivity of a nanofluid varies with the local particle concentration. Maxwell (1873)andJeffrey(1973) derived analytical expressions for the thermal conductivity of a suspension of spherical particles in a fluid. Buongiorno (2006) proposed analytical expressions based on experimental data for this feature for two cases of nanofluids, namely those of alumina particles in water and titanium particles in water. It is interesting to note that thermal conductivity stratification in a two-layer system of Newtonian fluids may significantly influence its stability. For instance, Welander (1964) and Gershuni & Zhukhovitskii (1980) investigated the onset of an oscillatory instability induced by the thermal conductivity stratification in a stably stratified two-layer liquid system. Such instability may be driven by a disparity between the characteristic diffusion time scales of the two liquids. 1011 A52-3 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron The purpose of this paper is to investigate the combined thermo-solutocapillary instability in a layer of a moderately dense nanoparticle suspension subjected to the Soret effect, bounded by a deformable liquid–gas interface and supported by a horizontal solid substrate subjected to a prescribed heat flux in the gravity field. The cases of both heating and cooling at the substrate are considered. To simplify the analysis, we assume the carrier fluid to be a simple Newtonian liquid. In contradistinction with the Rayleigh–Bénard instability in a nanofluid layer considered by Chang & Ruo (2022), the thermosolutocapillary instability in nanofluids has not been investigated to date. The main challenge and the novelty of the current investigation is taking into account the dependence of all thermophysical properties of the nanofluid on the local particle concentration, which was not considered before in this context, e.g., by Chang & Ruo (2022). In our approach, because of the dependence of the fluid density on the local particle concentration which varies within the bulk, we adopt a ‘compressible’ approach to describe the dynamics of a nanofluid and show that the contribution of ‘compressibility’ is minor. As a result of accounting for the dependence of the thermophysical properties on the local particle concentration, we find that in the case of a layer heated at the substrate, such variation of thermal conductivity of the fluid leads to the change of the instability type from monotonic for weak variations to oscillatory for stronger variations. In the case of the layer cooled at the substrate, the emerging instability is monotonic, similar to what Sarma & Mondal (2021b) found for a Newtonian dilute binary fluid. Most of the instabilities found in this investigation are finite-wave, although narrow windows of long-wave instability also emerge in both cases of the heating direction. The plan of the paper is as follows. Section 2offers the problem formulation, and presents a set of governing equations and boundary conditions accounting for local particle concentration-dependent thermophysical properties of the system. The quiescent base state of the system is presented in § 2.4 with the details of its derivation given in Appendix A.1. Section 3is dedicated to the linear stability analysis of the determined base state and presents the eigenvalue problem which is numerically investigated. Section 4, subdivided into six subsections, presents the results of the investigation: §4.1 outlines the numerical procedure used for solving the linear eigenvalue problem; § 4.2 presents the results for the case of cooling at the substrate; whereas §4.3 delivers the results for the case of heating at the substrate. Further, §§ 4.4,4.5 and 4.6 explore the effect of the thermal conductivity stratification of the suspension on the system instability, presents typical eigenfunctions corresponding to the observed instabilities and discusses a possibility of using an incompressibility simplification, which is akin, in some sense, to the Boussinesq approximation for a heated layer of a simple liquid, respectively. A use of a simplified ‘incompressible’ formulation could significantly reduce the numerical effort associated with the solution of the linear eigenvalue problem in its full formulation. Finally, § 5 summarises the findings of the paper. 2. Problem formulation and governing equations We consider a nanofluid layer of a mean thickness h∗ 0, density ρ∗ nf, dynamic viscosity μ∗ nf, kinematic viscosity ν∗ nf =μ∗ nf/ρ∗ nf, thermal conductivity K∗ nf, heat capacity c∗ nf, and thermal diffusivity κ∗ nf =K∗ nf/ρ∗ nfc∗ nf. Here, the subscripts nf refer to the nanofluid. In what follows, the presence of nanoparticles in a nanofluid will be accounted for via the local concentration of particles. We further denote the physical properties of the base fluid and the nanoparticles with subscripts bf and np, respectively. The nanofluid layer is assumed to be of local thickness 1011 A52-4 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics z y tg x h0 Gas 0 Solid Constant heat flux Nanoflui d λ n Figure 1. Nanofluid (d∗ p=10−100 nm)layer on the solid substrate subjected to a constant heat flux at the substrate and exposed to the gas phase at its deformable interface. h∗and rests on a solid planar horizontal substrate being exposed at its deformable liquid–gas interface to the quiescent gas environment held at constant pressure p∗ ∞and temperature T∗ ∞in the gravity field g∗. The frame of reference is chosen so that the x∗ and y∗axes are located in the substrate, whereas the axis z∗is normal to the substrate and directed into the fluid layer opposite to the direction of gravity; hence, the substrate and the deformable liquid–gas interface are located at z∗=0andz∗=h∗, respectively (figure 1). The fluid is assumed to be a moderately dense nanofluid, i.e. a mixture of a Newtonian base fluid with nanoparticles of d∗ p≈10 nm whose volumetric concentration φ∗= φ∗(x∗,y∗,z∗,t∗)is varying in time t∗and space. The entire system is subjected to the prescribed heat flux Qq∗with Q=±1andq∗>0 at its solid bottom in the direction normal to the latter. Note that the value of Qis related to the direction of the heat flux, so that the cases of Q=1andQ=−1 correspond to heating and cooling at the solid–liquid interface, respectively. An imposed heat flux leads to the emergence of the temperature field T∗=T∗(x∗,y∗,z∗,t∗)varying with time and space within the layer. The surface tension at the liquid–gas interface is assumed to depend on both the interfacial temperature and nanoparticle concentration, and for small variations of temperature and particle concentration at the interface, is adequately approximated by a linear function σ∗(T∗,φ ∗)=σ∗ r−σ∗ T∗(T∗−T∗ r)−σ∗ φ∗(φ∗−φ∗ r), (2.1a) where σ∗ T∗=−∂σ∗/∂T∗>0andσ∗ φ∗=−∂σ∗/∂φ∗>0.(2.1b) Here, σ∗ ris the equilibrium reference value of the surface tension at the reference values of T∗ rand φ∗ r,soσ∗ r=σ∗(T∗ r,φ ∗ r), and both σ∗ T∗and σ∗ φ∗are positive, thus, the surface tension linearly decreases with both the temperature and particle concentration at the interface. We use, as an example, an alumina nanoparticle dispersion in distilled water which exhibits a vanishing value of the adsorption/desorption coefficient ratio; hence, it is possible to neglect interfacial kinetics mechanisms in this situation. However, we note that nanoparticle dispersions in different liquids such as non-stabilised water, n-decane, ndodecane, n-hexadecane, etc., exhibit a significant contribution of the interfacial kinetics. In those cases, consideration of the interfacial kinetics becomes necessary (Machrafi 2022), and instead of the interfacial value of the bulk particle concentration φ∗, a surface particle concentration Γ∗needstobeusedin(2.1a)and(2.1b). 1011 A52-5 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron 2.1. Thermophysical properties of a nanofluid In the case of non-dilute mixtures, the thermophysical properties of a nanofluid depend on the local particle concentration (Buongiorno 2006). The density and the specific heat capacity of a nanofluid are determined as the weighted sum of the respective properties of the base fluid and the nanoparticles, ρ∗ nf(φ∗)=(1−φ∗)ρ∗ bf +φ∗ρ∗ np,(2.2a) c∗ nf(φ∗)=(1−φ∗)(ρ∗c∗)bf +φ∗(ρ∗c∗)np ρ∗ nf .(2.2b) The thermal conductivity of nanofluids is also known to depend on the local particle concentration. In what follows, we use the empirical model K∗ nf(φ∗)=K∗ bf (1+aφ∗), (2.2c) where ais a constant describing the degree of the thermal conductivity variation with the local particle concentration. The value a=7.47 (Buongiorno 2006) is valid for an alumina Al2O3nanoparticles suspension in water. This value was extracted from the experimental data of Pak & Cho (1998) who measured thermal conductivity of an alumina–water nanofluid for various particle concentrations. The relationship between the thermal conductivity of the nanofluid and the particle concentration depends on both particles and solvent material. For instance, in the case of titania nanoparticles in water, the thermal conductivity varies with the particle concentration as K∗ nf(φ∗)=K∗ bf (1+2.92φ∗−11.99φ∗2)(2.2d) (Buongiorno 2006). We note that in the absence of an empirical relationship between the thermal conductivity K∗ nf of a nanosuspension and its particle concentration φ∗, one can employ the analytical models developed by Maxwell (1873)andJeffrey(1973), which estimated the value of K∗ nf for a dilute (φ∗/φm1withφmbeing the maximal packing volume fraction) suspension providing terms proportional to φ∗, whose coefficient could have yielded the value of a,and(φ∗)2, respectively. Our estimate for the value of a, based on Maxwell theory and the values of thermal conductivity of Al2O3and water, yields a3, which is quite far from the experimental value of a=7.47; therefore, in most of our results presented below, a=7.47 is adopted. Additionally, we mention that the limits of the effective thermal conductivity (K∗ nf/K∗ bf )of a monodisperse nanosuspension can be estimated by the Hashin– Shtrikman bounds (Hashin & Shtrikman 1962; Keblinski, Prasher & Eapen 2008) ⎡ ⎣1+ 3φ∗K∗ np −K∗ bf  3K∗ bf +(1−φ∗)K∗ np −K∗ bf ⎤ ⎦⩽K∗ nf K∗ bf ⩽⎡ ⎣1− 3(1−φ∗)K∗ np −K∗ bf  3K∗ np −φ∗K∗ np −K∗ bf ⎤ ⎦K∗ np K∗ bf . (2.2e) We find that in a low-concentration limit φ∗1, the effective thermal conductivity parameter afor Al2O3nanoparticles in water ranges therefore in the interval 2.84 ⩽a⩽38. Also, the effect of varying aon the properties of the instability will be briefly assessed in § 4.4. As for the nanofluid viscosity, we employ the empirical correlation model proposed by Maron & Pierce (1956) and de Kruif et al. (1985) for a concentrated suspension, μ∗ nf(φ∗)=μ∗ bf 1−φ∗ φm−2 ,(2.2f) 1011 A52-6 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics where φmis the maximal packing fraction. The value of φmvaries from 0.524 to 0.71. We use the value of φm=0.65 related to the random close packing (RCP) volume fraction of nanoparticles. The viscosity model given by (2.2f) illustrates the power-law variation of the viscosity with the nanoparticle concentration growing indefinitely when the nanoparticle concentration φ∗→φm. 2.2. Physical mechanisms relevant for our model Buongiorno (2006) described seven different slip mechanisms, i.e. mechanisms causing deviation of the particle velocity field from that of the carrier fluid, for the nanoparticle motion in a nanofluid. Out of these seven, we account for the two dominant slip mechanisms, namely, the Brownian and the Soret diffusion (thermophoresis). The Brownian diffusion coefficient is given by the generalised Stokes–Einstein formula valid for arbitrary mass fraction of nanoparticle concentration (Russel, Saville & Schowalter 1989; Bird, Stewart & Lightfoot 2002; Espín & Kumar 2014), D(φ∗)=DBK(φ∗)dφ∗Z(φ∗) dφ∗,(2.3) where DBis the diffusivity coefficient given by the Stokes–Einstein formula as DB=KBT∗ 3πμ∗ bf d∗ p ,(2.4) and KBis the Boltzmann constant, KB=1.380649 ×10−23 JK−1. The generalised Stokes–Einstein formula (2.3) exhibits a strong dependence of the diffusion coefficient on the local nanoparticle concentration via hydrodynamic and thermodynamic interaction given by the sedimentation coefficient K(φ∗)and the compressibility contribution Z(φ∗), respectively. Combining the Carnahan & Starling (1969) equation for the compressibility effect Z(φ∗)=1+φ∗+φ∗2−φ∗3 (1−φ∗)3and the semi-empirical expression for the sedimentation coefficient K(φ∗)=(1−φ∗)6.55, Russel et al. (1989) derived the following form for the generalised Stokes–Einstein formula: D∗(φ∗)=DB1−φ∗2.55 1+4φ∗+4φ∗2−4φ∗3+φ∗4,(2.5a) which reduces for low particle concentrations φ∗to D∗(φ∗)=DB1+1.45φ∗(2.5b) previously derived by Batchelor (1976). Despite this, in our investigation below, we will use a constant value for the Brownian diffusivity coefficient D∗(φ∗)=DBwith the justification given towards the end of § 2.3 andin§4.4. The second slip mechanism is due to the fact that the nanofluid layer is non-isothermal. The temperature gradient induces a flux of nanoparticles, and this phenomenon is known as the thermophoresis or the Soret effect. The Soret or thermal diffusion coefficient in a nanofluid is proportional to the particle concentration φ∗and is expressed by Whitmore & Meisen (1977) and Morozov (2002)as DT=0.26K∗ bf 2K∗ bf +K∗ np μ∗ bf ρ∗ bf φ∗.(2.6) 1011 A52-7 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron In what follows, we consider the total nanoparticle mass flux as a superposition of the Brownian and Soret diffusion processes via jp=−ρ∗ np CBT∗∇∗φ∗+CTφ∗∇∗T∗ T∗,(2.7) where CB=DB/T∗,CT=DT/φ∗and ∇∗=(∂x∗,∂ y∗,∂ z∗)with subscripts x∗,y∗,z∗ denoting partial differentiation with respect to the corresponding variable. Buzzaccaro et al. (2008) noted that the expression in the parentheses in (2.6) is valid for relatively large nanoparticles of diameter 1 μm in water. In fact, it was found that the Soret coefficient CT depends on the nanoparticle size (Braibanti, Vigolo & Piazza 2008; Michaelides 2015). It is also important to note that the total particle mass flux given by (2.7) ensures a consistently positive distribution of nanoparticle concentration across the nanofluid layer even when subjected to a strong thermophoresis. However, we also note that a use of a constant Soret coefficient in front of the ∇∗T∗term in (2.7) is constrained to a range bounded from above for this coefficient, since beyond this range, spurious unphysical negative values for the particle concentration φ∗emerge. This fact was also emphasised by Dastvareh & Azaiez (2018). The other five slip mechanisms mentioned by Buongiorno (2006) are inertia, diffusiophoresis, Magnus effect, fluid drainage and gravitational settling. Inertia has a negligible effect due to the homogeneous motion of nanoparticles with the surrounding continuum media. We neglect the diffusiophoresis effect for a one-component nanofluid; however, it may be important when the base fluid is subjected to an additional solute species gradient (Ruckenstein 1981; Anderson 1989; Morozov 2002). The Magnus effect arises due to a force perpendicular to the main flow direction, induced by the relative axial velocity between the nanoparticle and fluid flow. We neglect the Magnus effect as well considering an homogeneous motion of the nanoparticles with the surrounding fluid. The fluid drainage contribution is important for the distance between the wall and particle of the order of nanoparticle diameter and can be neglected for nanoparticles with a small diameter d∗ p. The relative strength of the particle flux due to gravitational settling (Mason & Weaver 1924; Shliomis & Smorodin 2005; Buzzaccaro et al. 2008; Cherepanov & Smorodin 2019) versus the flux due to Brownian diffusion may be estimated by the ratio Sg= h∗ 0(ρ∗ np −ρ∗ bf )g∗(πd∗3 p/6) KBT∗ r .(2.8) For nanoparticles of d∗ p≈10 nm and ρ∗ np ∼4gcm −3in a nanofluid layer of 0.1mm thickness at room temperature, and the terrestrial gravity is Sg≈3.7×10−4and, therefore, the mass flux induced by gravitational settling may be neglected. However, gravitational settling contributes significantly for nanoparticles with a larger diameter, say of d∗ p≈100 nm with Sg≈0.37. Therefore, in the latter case, the contribution of gravitational settling could not be disregarded, see also Chang & Ruo (2022)andtheir modification for the mass flux, their (2.3). For more details, see also Appendix A.2. Finally, the heat flux in a nanofluid is given by Fourier’s law of heat conduction, jT=−K∗ nf∇∗T∗.(2.9) The Dufour effect is exceedingly weak in liquid mixtures (Cross & Hohenberg 1993; Oron & Nepomnyashchy 2004; Morozov et al. 2014) and can be safely neglected in the case at hand. We also note that following Buongiorno (2006), the expression for the heat flux jTin (2.9) of Chang & Ruo (2022) had four more terms, with one of them arising from gravitational settling. All these terms are found to be negligible in our case in comparison 1011 A52-8 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics with the heat conduction term there, and will be omitted in what follows. A justification for this simplification will be presented in § 2.3 in the context of non-dimensional equations (2.13). 2.3. Governing equations The set of governing equations comprises the continuity, Navier–Stokes, advection– conduction heat transfer and nanoparticle mass transfer equations (Rohsenow, Hartnett &Cho1998; Batchelor 2000; Colinet, Legros & Velarde 2001;Birdet al. 2002), respectively (ρ∗ nf)t∗+∇∗·(ρ∗ nfu∗)=0,(2.10a) (ρ∗ nfu∗)t∗+∇∗·(ρ∗ nfu∗u∗)=−∇∗p∗+∇∗·τ∗−ρ∗ nf g∗ez∗,(2.10b) (ρ∗c∗)nf (T∗)t∗+(ρ∗c∗)nf(u∗·∇∗)T∗=∇∗·(K∗ nf∇∗T∗), (2.10c) (φ∗)t∗+∇∗·(u∗φ∗)=∇∗·CBT∗∇∗φ∗+CTφ∗∇∗T∗ T∗,(2.10d) where g∗is the gravity acceleration, ez∗is the unit vector in the z∗direction, u∗= (u∗,v ∗,w ∗)is the flow field vector and τ∗is the viscous part of the stress tensor τ∗=μ∗ nf(∇∗u∗+(∇∗u∗))−2 3μ∗ nf −K∗(∇∗·u∗)I, where the superscript denotes the transpose of the corresponding tensor, Iis the unity tensor, the subscript t∗stands for a partial derivative with respect to time t∗and K∗is the dilatational viscosity of the fluid. We note that the last term in (2.10b) 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.10). We impose three boundary conditions at the solid–liquid interface z∗=0. At the substrate, the fluid velocity exhibits no-slip and no-penetration, zero total mass flux implying the impermeability of the substrate, and a constant prescribed heat flux Qq∗. We note that it is quite natural to prescribe the heat flux at the substrate to better fit experimental settings in the case where the substrate is not made of material with a high thermal conductivity (Rohsenow et al. 1998): z∗=0:u∗=0,−K∗ nf T∗ z∗=Qq∗,CBT∗φ∗ z∗+CTφ∗T∗ z∗ T∗=0.(2.11a) At the deformable interface z∗=h∗(x∗,y∗,t∗), we impose the kinematic boundary condition, the continuity of the normal and tangential stresses, the continuity of the heat flux and impermeability for mass transfer, respectively h∗ t∗+u∗h∗ x∗+v∗h∗ y∗=w∗,(2.11b) n∗·T∗·n∗=2H∗σ∗,n∗·T∗·x∗=x∗·∇∗ sσ∗,n∗·T∗·y∗=y∗·∇∗ sσ∗,(2.11c) −K∗ nf∇∗T∗·n∗=q(T∗−T∗ ∞), (2.11d) CBT∗∇∗φ∗+CTφ∗∇∗T∗ T∗·n∗=0,(2.11e) where ∇∗ s=(I−n∗n∗)·∇∗is the surface gradient operator, n∗= −h∗ x∗ex∗−h∗ y∗ey∗+ez∗ 1+h∗2 x∗+h∗2 y∗ (2.11f) 1011 A52-9 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron where ε=⎛ ⎜ ⎝ ¯uz+¯wx ¯vz+¯wy 2¯wz−2 3∇·¯ u⎞ ⎟ ⎠(3.3e) is a vector containing the two off-diagonal and one diagonal components of the strain rate associated with the z-direction. The linearised boundary conditions at z=0are ¯ u=0,K0∂¯ T ∂z+¯ KdT0 dz=0,ηφ 0∂¯ T ∂z+ηdT0 dz ¯ φ+∂¯ φ ∂z=0.(3.4) Note that (3.3a), (3.3b), (3.3c), and (3.3d) represent the linearised versions of the continuity, three-dimensional momentum conservation, energy and mass diffusion equations, respectively. The linearised boundary conditions at z=1are P∂¯ ζ ∂t−¯w=0,(3.5a) ∂p0 ∂z ¯ ζ+¯p+2 3M0∇·¯ u−3∂¯w ∂z+Σ0∇2 ⊥¯ ζ=0,(3.5b) MS∇⊥¯ φ+dφ0 dz∇⊥¯ ζ+MT∇⊥¯ T+dT0 dz∇⊥¯ ζ+M0ε⊥=0,(3.5c) ¯ KdT0 dz+K0∂¯ T ∂z+dT0 dzB+dK0 dz+K0d2T0 dz2¯ ζ+B¯ T=0,(3.5d) ∂¯ φ ∂z+ηdT0 dz dφ0 dz+ηφ0d2T0 dz2+d2φ0 dz2¯ ζ+ηdT0 dz ¯ φ+ηφ0∂¯ T ∂z=0,(3.5e) where ∇⊥≡∂ ∂x,∂ ∂y and ε⊥=¯uz+¯wx ¯vz+¯wy.(3.5f) Equations (3.5a), (3.5b), (3.5c), (3.5d), and (3.5e) represent the linearised versions of the kinematic, balance of normal stresses, two-dimensional balance of tangential stresses, heat and mass fluxes boundary conditions, respectively. In what follows, we constrain our study to the two-dimensional (2-D) case in the x−zplane. By applying the ∇× operator to the momentum conservation equations and differentiating the normal stress balance boundary condition along the interface followed by substituting there the pressure gradient components obtained from the momentum conservation equations, we completely eliminate the pressure from the momentum conservation equations. Normal perturbation modes are introduced in the form ¯u,¯w, ¯ T,¯ φ, ¯ ζ;¯ X(x,z,t)=(u(z), w(z), θ(z), ϕ(z), ζ ;X)exp(ikx +λt), (3.6) where u,w,θ, ϕ, ζ,andX≡(R(z), M(z), K(z), H(z)) are the amplitudes of the horizontal and vertical fluid velocity components, temperature, concentration, interfacial deformation and material properties (all expressed via ϕ) perturbations, respectively, with kand λrepresenting the wavenumber and the growth rate of the disturbances, respectively. 1011 A52-16 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics In this representation, positive and negative values of the real part of λcorrespond to instability and stability regimes of the system, respectively. We note that Xrepresents the material properties of the nanofluid which are concentration-dependent and does not contain independent variables following (3.2). The resulting eigenvalue problem which determines the growth rate λof the disturbance as a function of its wavenumber kand the rest of the problem parameters reads R0w+iku+λPR+wR 0=0,(3.7a) −iGkR+k2M0u+2k2uM 0−ikwk2M0+λR0+M 0+ikM0w −uM0 −2uM 0+λR0u−uM 0+λuR 0=0,(3.7b) −T 0K+θk2K0+λPH0+T 0wH0−K−K0θ −K 0θ=0,(3.7c) ϕk2L−ηLT 0+λP+ηk2Lθφ0+φ0iku −ηLθ +w+φ 0w−ηLθ −ηLT 0ϕ−Lϕ =0.(3.7d) The boundary conditions at z=0are u=(aφmφ0+1)θ+aφmT 0ϕ=ηφ0θ+ηT 0ϕ+ϕ=0.(3.7e) The boundary conditions at the deformable interface z=hprojected onto z=1are ζλP−w=0,(3.7f) T 0aφmζφ 0+ϕ+Bζ+(aφmφ0+1)θ+ζT 0(aφmφ0+1)+Bθ=0,(3.7g) ζφ 0+ηφ0θ+ζT 0+ηT 0ζφ 0+ϕ+ϕ=0,(3.7h) kMSζφ 0+ϕ+kMTζT 0+θ+M0kw−iu=0,(3.7i) iζk3Σ0+u2k2M0+λR0+i−ζkp 0+kM0w−kM 0w+iM0u +iM 0u=0. (3.7j) Here and on, primes denote differentiation with respect to z. We reiterate that R0,M0,K0,H0represent the base state values of the thermophysical properties of the system given by (3.2), whereas R,M,K,Hare the amplitudes of their respective disturbances depending on z. Finally, we note that normal perturbations imposed on the bulk concentration φautomatically preserve the average bulk concentration due to periodicity of the x-component exp(ikx)of the perturbation ¯ φ. 4. Results 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 ufrom (3.7a) into the rest of the equations of the eigenvalue problem (3.7). We also extract the expression for the interfacial deformation ζ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, θ, ϕ 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 ζ=0isimposed. 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 1011 A52-17 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron developed by Pearce et al. (2018). It is based on the shooting method and numerical construction of the complex-valued Evans function (Evans & Feroe 1977) 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 λ, and by doing so, one obtains contour plots of the Evans function for a specified range of (λr,λi), where λr≡(λ)and λi≡(λ). Then, the intersection points of the contour lines with the reference point of the plane or with the axis λ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 1979)to construct the Evans function by satisfying the boundary conditions separately at each of the two endpoints z=0andz=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 ⩽z⩽1, namely z=1/3,z=1/2andz=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∼10 nm, the Lewis number estimated using the Stokes–Einstein formula (2.4),isL=4.22 ×10−4when the solvent temperature is T∗=300 K. Furthermore, we consider the range of the values for the Soret coefficient η∈(0.31,3.1)corresponding to the temperature drop of T∗∈(1,10) K across the layer. In addition, we determine the dimensionless parameters related to nanoparticle density via ρnp =ρ∗ np ρ∗ bf ≈4and(ρc)np =ρ∗ npc∗ np ρ∗ bf c∗ bf ≈0.68,(3.7j) which corresponds to alumina nanoparticles in water, also taking a=7.47 in this investigation, except for § 4.4 where the parameter ais changing. Before we present the results of the linear stability analysis, we note that the equilibrium surface tension σ∗ rranges for aqueous solutions from ∼10−3to 30 ×10−3Nm−1(De Wit, Gallez & Christov 1994); therefore, the dimensionless surface tension for a layer of density ρ∗ bf =103kg m−3, viscosity μ∗ bf =10−3kg (m ·s)−1, thermal diffusivity κ∗ bf =0.146 × 10−6m2s−1and thickness h∗ 0=10−4m is approximately Σ0=2×104. Similarly, the modified Galileo number for the terrestrial gravity of g∗=9.81 m s−2yields G=67.2for h∗ 0=10−4mandG∼10−4for h∗ 0=10−6m. A list of typical values of the thermophysical properties of the nanofluid and the dimensionless numbers is given in table 1. 1011 A52-18 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics Parameter Symbol Values Base fluid viscosity μ∗ bf 10−3kg (m ·s) −1 Base fluid density ρ∗ bf 1000 kg m−3 Base fluid kinematic viscosity ν∗ bf 10−6m2s−1 Base fluid thermal conductivity K∗ bf 0.61W(m·K) −1 Base fluid heat capacity c∗ bf ≈4.18 kJ (kg ·K) −1 Base fluid thermal diffusivity κ∗ bf 0.146 ×10−6m2s−1 Ambient gas temperature T∗ ∞≈300 K Equilibrium surface tension σ∗ r∼10−3−30 ×10−3Nm−1 Surface tension coefficient σ∗ T∗0.069 ×10−3N(mK)−1 Gravity g∗∼0.981−9.81 m s−2 Temperature difference across nanofluid layer T∗∼0.1−10 K Prandtl number P7 Modified Galileo number G∼10−4−67 Dimensionless surface tension Σ0∼10–2 ×104 Biot number B∼10−3−10−1 Random close-packing volume fraction (RCP) φm≈0.65 Nanoparticle bulk concentration Φ∼0.01−0.027 Nanoparticle diameter d∗ p≈10 nm Boltzmann constant KB1.380649 ×10−23 JK−1 Nanoparticle density ρ∗ np ∼4000 kg m−3 Nanoparticle heat capacity c∗ np ≈0.718 kJ (kg ·K) −1 Lewis number L∼10−3−4.22 ×10−4 Soret coefficient η∼0.31−3.1 Equilibrium nanofluid thickness h∗ 0∼10−6−10−4m Table 1. Parameter nomenclature and their typical values used in this investigation. 4.2. Cooling at the substrate (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 1992;Oron&Rosenau1992;Burgesset al. 2001;Alexeev&Oron2007). Recently, Sarma & Mondal (2021b) 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. 1011 A52-19 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron 22.70 90 80 70 60 50 40 30 20 22.65 22.60 22.55 22.50 MS (a)(b) MS 22.45 22.40 22.35 22.30 10–3 10–2 U S kk S U 10–1 10–2 10–1 B = 0.001 B = 0.01 B = 0.1 MT = 0 MT = 10 MT = 25 Figure 4. Monotonic solutocapillary instability in the case of cooling at the substrate Q=−1atΦ=0.01, a=7.47,η=0.31,L=10−3,Σ 0≈2×104and G=6.71. (a) Neutral curves MS(k)for the pure solutocapillarity instability, MT=0, versus the wavenumber kfor various Biot numbers B.(b) Neutral curves MS(k)for various values of MTwith B=0.01. The rise of the neutral curves with an increase in MTillustrates a stabilising effect of thermocapillarity on the monotonic solutocapillary instability. In both panels, the symbols Uand Sdenote the unstable and stable domains of the system, respectively. First, it must be clear that in the configuration at hand, thermocapillarity is stabilising for MT>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 MSexhibits 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.Bothofthemare clearly seen in the neutral curve for B=0.01 in figure 4(a). 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 1992;Alexeev&Oron2007). Figure 4(b) illustrates a significant stabilisation effect imparted by thermocapillarity on the solutocapillary instability threshold. In the range of MTpresented in this figure, the critical value of MS slightly increases with MT; however, the finite-wave modes are strongly stabilised with an increase in MT. Figure 5(a) shows that the threshold value of the solutal Marangoni number MSin the pure solutocapillary instability increases linearly with the modified Galileo number G. However, as we already mentioned, the neutral curves MS(k)exhibit two minima: one of them, MS=mf, is in the finite-wave domain, whereas the second one, MS=ml, 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 1011 A52-20 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics 10 20 30 40 50 60 70 G 22.5 25.0 27.5 30.0 32.5 35.0 37.5 40.0 42.5 U S 20 40 60 0 0.05 0.10 ml (LW ) mf (FW ) 10 20 30 40 50 60 70 G 0.2 0.3 0.4 0.5 0.6 0.7 0.8 kc kc mf ml 20 40 60 G 3 4 5 6 7 8 ×10 –4 MS (a)(b) mf–m l mf ml mf (FW ) Figure 5. Onset of a purely solutocapillary monotonic instability in the case of cooling at the substrate Q=−1atΦ=0.01,a=7.47,η=0.31,L=10−3,MT=0,Σ 0≈2×104and B=0.01. (a) Variation of the critical value of the solutal Marangoni number MSversus the modified Galileo number Gthat shows both the finite-wave range minimum MS=mfand the long-wave minimum MS=ml; the inset presents the variation of the difference ml−mfversus G. FW and LW stand for domains of finite-wave and long-wave instability, respectively. (b) Variation of the critical wavenumber kcwith G; the inset shows a decrease in the critical wavenumber kccorresponding to the long-wave minimum MS=mlwith G. The symbols Uand S denote the unstable and stable domains of the system, respectively. of figure 5(a) shows that in the domain of intermediate values of the modified Galileo number 14.8⩽G⩽66.6, the long-wave minimum MS=mlcorresponding to the deformational mode prevails over the finite-wave minimum MS=mf, namely ml<mf, and corresponds to the instability onset. In addition, we find that in the rest of the Grange, i.e. for relatively low and sufficiently high values of G,ml>mf, and the finite-wave solutocapillary instability emerges. VanHook et al. (1997) and Golovin, Nepomnyashchy & Pismen (1997) 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 kcwith G. We note that the critical wavenumber kcfor the finite-wave minimum MT=mfexhibits two distinct domains of variation with G. For moderate values of the modified Galileo number G,kcvaries slowly; however, for sufficiently large values of G, the critical wave number kcsharply increases. The inset of figure 5(b)shows a decrease in the critical wavenumber kcof the long-wave instability corresponding to the minimum MS=mlwith G. Figure 6(a) shows that in the case of a purely solutocapillary instability, the critical solutal Marangoni number MSmonotonically increases with the dimensionless surface tension number Σ0when the rest of the parameters are fixed and saturates at a very high value of Σ0that 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 MSis attained at the local minimum of the neutral curves corresponding to MS=mfand varies weakly with Σ0. The variation of the critical wavenumber kc is non-monotonic with Σ0in 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 MS(k)=mlremains almost unchanging with Σ0and small, k≈2.0×10−3; therefore, this secondary mode is long-wave. We next study the variation of the onset of the pure solutocapillary instability with the averaged nanoparticle bulk concentration Φ.Figure 7(a) displays a significant monotonic 1011 A52-21 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron 10 1 10 2 10 3 10 4 Σ0 10 1 10 2 10 3 10 4 Σ0 22.15 22.20 22.25 22.30 22.35 22.40 22.45 22.50 U S 22.41 22.42 22.43 22.44 22.45 22.46 22.47 m f m l 0 0.05 0.10 0.15 0.20 k = 2 × 10 –3 MS (a)(b) m f m l kc Figure 6. Variation of the onset of the pure solutocapillary monotonic instability with the dimensionless surface tension number Σ0for the case of cooling at the substrate Q=−1atΦ=0.01,a=7.47,η= 0.31,L=10−3,MT=0,G=6.71 and B=0.01. The mland mfpoints correspond to the long-wave and finite-wave modes, respectively. The symbols Uand Sdenote the unstable and stable domains of the system, respectively. (a) Variation of the values of the solutal Marangoni number MScorresponding to the two local minima MS=mland MS=mfof the neutral curves MS(k). The instability threshold MSincreases with Σ0. The right vertical scale corresponds to MSrelated to the long-wave minimum MS=ml.(b) Variation of the wavenumbers corresponding to the two local minima of the neutral curves with Σ0. The critical wavenumber kccorresponding to the onset of the instability lies in the finite-wave range and varies non-monotonically with the dimensionless surface tension number Σ0. The wavenumber kcorresponding to the long-wave competing mode is almost constant with Σ0. decrease in the threshold value of the solutal Marangoni number MS(for both modes mfand ml) with an increase in Φ. The reason for this effect is mainly due to the Soret effect whose strength is proportional to the local concentration φ. In the case of cooling at the substrate, the temperature gradient is directed upward; therefore, the Soret effect with η>0 drives the thermophoretic component of the mass flux downward and becomes stronger with an increase in the average particle concentration Φ.Thisleadstoastronger decrease in particle concentration in the neighbourhood of the interface and to an enhanced decrease in local viscosity there with an increase in Φ. Thus, the critical solutal Marangoni number MSdriving the instability decreases with Φ. 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. 2014; Shklyaev & Nepomnyashchy 2017). However, for an inorganic salt where the surface tension increases with the interfacial concentration, the opposite behaviour could be observed (Oron & Nepomnyashchy 2004). Further, we note that the critical wavenumber kcfor the finitewave minimum gradually increases with Φ, whereas the wavenumber for the long-wave minimum remains almost constant around kc≈8×10−4. Since the finite-wave minimum MS=mfremains always slightly below the long-wave one MS=ml, 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 Φand a typical neutral curve is shown in the inset of figure 7(a). Figure 8(a) illustrates the non-monotonic variation of the critical solutal Marangoni number MSwith the Soret coefficient ηfor both competing modes. We find that the temperature difference across the layer of the base state given by T∗=q∗h∗ 0/K∗ bf influences the value of the Soret coefficient ηdefined in (2.17) and affects the competition 1011 A52-22 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics 0.02 0.04 0.06 0.08 0.10 0.02 0.04 0.06 0.08 0.10 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 MS (a)(b) U S 10–2 10–1 k 9.94 9.96 9.98 10.00 MS U S Φ = 0.03 ΦΦ 0 0.05 0.10 0.15 0.20 0.25 m f m l kc ≈ 8 × 10 –4 m f m l kc Figure 7. Variation of the onset of the pure solutocapillary monotonic instability with the averaged bulk nanoparticle concentration Φfor the case of cooling at the substrate Q=−1anda=7.47,η=0.31, L=10−3,MT=0,G=6.71,Σ 0=2×104and B=0.01. (a) Variation of the two minimal values MS=mf and MS=mlof the neutral curves MS(k). The inset shows the neutral curve for Φ=0.03. (b) Variation of the wavenumbers corresponding to the two local minima of the neutral curves with Φ. The upper curve corresponding to MS=mfrepresents the critical wavenumber kc. The symbols Uand Sdenote the unstable and stable domains of the system, respectively. 0.5 1.0 1.5 2.0 2.5 3.0 12 14 16 18 20 22 24 26 MS (a)(b) US 12 η η ηη 3 –3 –2 –1 0 mf – ml ml (LW ) mf (FW ) 100 0.25 0.50 0.75 1.00 1.25 1.50 1.75 kc ∼ η4/3 3×10 13×10 0 0.5 1.0 1.5 2.0 0.50 0.75 1.00 1.25 1.50 ×10 –3 m f m l mf (FW ) m f m l kc ∼ η1/2 kc ∼ η2/3 kc kc Figure 8. (a) Variation of the critical value of the solutal Marangoni number MSfor a pure solutocapillary monotonic instability with the Soret coefficient ηin the case of cooling at the substrate Q=−1with a=7.47,Φ=0.01,L=10−3,MT=0,G=6.71,Σ 0=2×104and B=0.01;theinsetshowsthe difference mf−mlbetween the two minimal values of MSof the two solutocapillary modes versus η.(b) Variation of the critical wavenumber kccorresponding to the short-wave minimum MS=mfwith η;theinset shows the variation of the wavenumber kcorresponding to the long-wave mode MS=mlwith η. The dashed lines display various ranges of the data fit for the critical wavenumber as a function of η. The symbols U and Sdenote the unstable and stable domains of the system, respectively. FW and LW stand for domains of finite-wave and long-wave instability, respectively. between the long-wave deformational mode MS=mland the finite-wave pattern forming mode MS=mf. The inset of figure 8(a) displays the domains of instability versus the parameter η. We find that in the narrow range of the Soret coefficient 0.62 ⩽η⩽1.24 corresponding to the temperature difference 2K ⩽T∗⩽4K, the deformational mode MS=mldominates over the finite-wave mode MS=mfin terms of the instability onset, namely ml<mf. Outside of this range of η, with a sufficiently low value of the modified 1011 A52-23 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron Galileo number G, the finite-wave mode MS=mfprevails over the long-wave mode MS=mland determines the instability onset, namely mf<ml.Figure 8(b) displays the variation of the critical wavenumber kcwith the Soret coefficient η. We observe that the critical wavenumber increases with ηand exhibits three different functional variations with η. Notably, we observe that the critical wavenumber kcfollows the scaling kc∼η1/2 up to the value of ηcorresponding to the minimum of kcfor MS=ml, which is shown in the inset of figure 8(a). This is followed by two different scalings kc∼η4/3and kc∼η2/3 for higher values of η. The inset of figure 8(b) also shows that the critical wavenumber for the long-wave mode MS=mlexhibits a non-monotonic variation with η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 η. 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 (2021b) 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 (2021b) 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 (2021b) is possibly due to the presence of concentration-dependent thermophysical properties of the nanofluid. 4.3. Heating at the substrate (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 (1995) and Oron & Nepomnyashchy (2004) 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 1961), which, in what follows, will be referred to as a solutal buoyancy due to its qualitative similarity to the thermal buoyancy. 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 MTfor G=1 near the critical wavenumber kc≈0.18 for a parameter set specified in the caption. The entire variation range of MTis naturally subdivided into five different subdomains. With an increase in MT, the first subdomain (subdomain I)is where the two tracked eigenvalues are real and negative. The two branches of subdomain Imerge with an increase in MTand split into two branches, along which the eigenvalues are complex conjugate with the negative real part λr; this is subdomain II. A combination of subdomains Iand 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 MT,the frequency λiof the two leading eigenvalues increases in its absolute value, whereas the 1011 A52-24 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 MT –1 0 1 2 3 4 λrλi ×10 –4 I : (2–) II : (2C–) III : (2C+) IV : (2+) V : (+, –) M2M1 MO III III IV V –4 –3 –2 –1 0 1 2 3 4 ×10 –5 Figure 9. Eigenspectrum map of the two leading eigenvalues λin the complex plane spanned by their growth rate (the real part (λ)≡λrthe vertical axis on the left) and their frequency (the imaginary part (λ)≡λi– the vertical axis on the right) versus MTin the case of heating at the substrate Q=1withΦ=0.01,G=1, Σ0=10,η=0.31,L=10−3,a=7.47,MS=0,P=7,B=0.01 and k=0.18. The domains I−Vare the domains with two negative real eigenvalues (2−), two complex conjugate eigenvalues with a negative real part (2C−), two complex conjugate eigenvalues with a positive real part (2C+), two positive real eigenvalues (2+), and two real eigenvalues (one positive (+) and one negative (−) ), respectively. The circle (◦)and star () points correspond to the real and imaginary parts of the complex growth rate λ, respectively. The vertical dashed lines represent the transition from the subdomain II to the subdomain III,fromIII to IV and from IV to V. The horizontal dashed line represents λr=0. growth rate λrincreases and crosses zero. We emphasise that the oscillatory instability emerges here in the case of MS=0; hence, solutocapillarity is irrelevant. This is where the subdomain II ends and the subdomain III begins. In the subdomain III,thetwo leading eigenvalues are still complex conjugates, but their real part λris positive and, therefore, the system undergoes oscillatory instability. In the subdomain III,thegrowth rate λiincreases with MT, whereas the absolute value of the frequency |λ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 MT=M2and the subdomain IV begins. In the subdomain IV, there are two real positive eigenvalues, one of them increases with an increase in MT, whereas the other one decreases until it reaches zero, where the subdomain IV ends at MT=M1, and the subdomain Vbegins. 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 MT(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≈0.18. Along with the neutral curve MT(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 M2separates 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 MOand M2lines constitutes the subdomain III. Slightly further above the boundary M2lies one more boundary marked as M1along which one of the leading positive eigenvalues turns to be zero, whereas the other one remains positive. The emergence of 1011 A52-25 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron 0.35 0.30 U U†U† Monotonic Monotonic Oscillatory Monotonic Oscillatory Oscillatory S S 0.25 0.20 0.15 10–1 100101101 100 MSMS MT 3.5 3.0 2.5 2.0 1.5 MT (a)(b) Figure 17. Threshold of the combined soluto-thermocapillary instability shown in the plane MT−MSin the case of heating at the substrate Q=1withΦ=0.01,η=0.31,L=10−3,Σ 0=10,B=0.01. Panels (a)and (b) correspond to G=0.01 and G=1, respectively. The symbols U,U†and Srepresent the domains of one real positive eigenvalue, of two complex conjugate leading eigenvalues with a positive real part, and that of stability, respectively. instability, whereas it significantly increases with MSfor the oscillatory instability. It is also interesting to note that the value of kcexhibits a discontinuity at MScorresponding 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 λiis presented in figure 18(b). It is found to increase monotonically and almost linearly with MSand collaterally with MTalong the critical curve. Figure 19 presents the neutral curves MT(k)for three values of MSwith both monotonic and oscillatory branches shown. In the case of MS=0 displayed in figure 19(a), the oscillatory branch bifurcates off the monotonic one at MThigher than the minimum of the neutral curve, and being a monotonically increasing function of k, it remains above the critical value of MTcorresponding to the onset of the monotonic instability of the system, as also shown in figure 17. The case of MS=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 MT, see the borderline between the two domains in figure 17(a). The critical wavenumber kcof 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 1984). With an increase in MS, 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 kcof the oscillatory instability is larger than that for MS=2.1, which is consistent with an increase of kcwith MSshown in figure 18(a) for the oscillatory instability. 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, ζ=0, and the absence of solutocapillarity, MS=0. Figure 20 displays the variation of the critical thermal 1011 A52-32 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics 9.0 ×10 –2 ×10 –6 8 7 6 5 4 3 2 1 λi 8.5 Monotonic Monotonic Oscillatory Oscillatory 8.0 7.5 7.0 6.5 6.0 5.5 0510 MS kc MS 15 20 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 (a)(b) Figure 18. The case of heating at the substrate Q=1atΦ=0.01,η=0.31,L=10−3,a=7.47,Σ 0=10, G=0.01,B=0.01. (a) Variation of the critical wavenumber kcwith the solutal Marangoni number MSalong the critical curve shown in figure 17(a). Note that the values of kalong 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 λiwith the solutal Marangoni number MSalong the critical curve presented in figure 17(a) in the domain of the oscillatory instability. Marangoni number MTwith the thermal conductivity stratification parameter a.Wefind that the monotonic instability sets in for sufficiently small values of a<ac∼6×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 2006; Coccia, Tomassetti & Di Nicola et al. 2021), 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 φfor D(φ) in equation (2.19)islarger 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 concentrationdependent Brownian diffusion coefficient D(φ) given by (2.19) is taken into account. We find that the critical value of MTbased on the use of the concentration-dependent Brownian diffusion coefficient given by Batchelor (1976) differs from MTdetermined using a constant Brownian diffusion coefficient DBby less than 2 %. We also find that the instability type, either monotonic or oscillatory, and the transition value of a=ac from one to another are not affected by whether a constant or particle concentrationdependent 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 1962; Keblinski et al. 2008) remains smooth, and an increase in aleads to a substantial increase in MT(not shown). Figure 21(a) demonstrates that the critical wavenumber kcfor the onset of Marangoni instability lies in the short-wave domain for both monotonic and oscillatory instabilities. The critical wavenumber slowly increases with acontinuously changing through the transition value of a=acfrom the monotonic to oscillatory domains. Figure 21(b)shows 1011 A52-33 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron 0.5 0.19 0.18 0.17 0.16 0.15 0.14 0.13 0.4 0.3 0.2 0.1 2.0 1.5 1.0 0.5 0 0.050 0.10 0.15 0.20 0.25 0.30 0 0.05 0.10 kk k MT (a)(b) (c) MS = 0 MS = 2.1 MS = 20 UU U S S S U† U† U† MT 0.15 0.20 0.02 0.04 0.06 0.08 0.10 MT Figure 19. Neutral curves MT(k)in the case of heating at the substrate Q=1withΦ=0.01,η=0.31, L=10−3,a=7.47,Σ 0=10,G=0.01 and B=0.01. Panels (a), (b)and(c) display the structure of the neutral curves for the monotonic (◦)and oscillatory branches () for MS=0,2.1 and 20, respectively. Panel (b) shows the emergence of the co-dimension two point at MS=2.1. The symbols U,U†and Srepresent the domains of one real positive eigenvalue, two leading complex conjugate eigenvalues with a positive real part, and that of the system stability, respectively. an increase in the critical frequency λiwith an increase in a. Again, the difference between the results based on a constant DBand nanoparticle concentration-dependent forms for the Brownian diffusivity D(φ) is minor. The variation of the critical thermal stratification parameter acfor the onset of oscillatory instability with the Lewis number Lis 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 acis proportional to L2. 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∼O(10−1),requires a significantly higher value of the thermal stratification parameter a=acfor the system to switch from the monotonic to oscillatory instability. However, with a sufficiently low value of the Lewis number L1, the oscillatory instability occurs already at a very low thermal conductivity stratification parameter ac1. We also note that the critical value acin 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 Σ0on the variation of acwill be discussed elsewhere. 1011 A52-34 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics 70.0 67.5 65.0 62.5 60.0 MTU S Monotonic Monotonic Oscillatory Oscillatory (D(φ)) Monotonic (D(φ)) Oscillatory 57.5 55.0 52.5 50.0 10–4 10–3 10–2 a 10–1 100 U† Figure 20. Variation of the critical thermal Marangoni number MTwith the thermal conductivity stratification parameter ain the case of heating at the substrate Q=1withΦ=0.01,η=0.31,L=10−3,MS=0,ζ=0, G=0andB=0.01. The symbols U,U†and Sdenote the domains with one real positive eigenvalue, with two leading complex conjugate eigenvalues with a positive real part, and that of stability, respectively. The ◦ and symbols denote, respectively, the onset of the monotonic and oscillatory instabilities in the case of a constant Brownian diffusivity DB, whereas ×and symbols represent, respectively, the onset of monotonic and oscillatory instability with nanoparticle concentration-dependent Brownian diffusion D(φ) given by (2.5b). 1.2 5 ×10 –3 4 3 2 1 0 1.1 1.0 0.9 kc Monotonic Oscillatory 0.8 0.7 0.6 10–4 10–3 10–2 10–1 aa 10001234567 Monotonic Oscillatory Oscillatory (D(φ)) Monotonic (D(φ)) D(φ) DB λi (a)(b) Figure 21. The case of heating at the substrate Q=1atΦ=0.01,η=0.31,L=10−3,MS=0,ζ=0, G=0andB=0.01. (a) Variation of the critical wavenumber kcwith the thermal conductivity stratification parameter a.The◦and symbols show kcof the monotonic and oscillatory instabilities, respectively, when a constant Brownian diffusivity DBis used, whereas the ×and symbols denote kcof the monotonic and oscillatory instability, respectively, when the nanoparticle concentration-dependent Brownian diffusion D(φ) 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†,and Ssymbols represent the domains of monotonic instability, oscillatory instability and stability of the system, respectively. (b) Variation of the critical frequency λiwith the thermal conductivity stratification parameter ain the domain where the oscillatory instability sets in. The and symbols represent the critical frequency λicorresponding to a constant and nanoparticle concentration-dependent forms for the Brownian diffusivity, DBand D(φ), respectively. 1011 A52-35 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron 101 100 10–1 10–2 10–3 10–3 10–2 L ac ac∝L2 a∈(0, 7.47) Monotonic Oscillatory 10–1 Figure 22. Variation of the critical thermal conductivity stratification parameter acwith the Lewis number L in the case of heating at the substrate Q=1withΦ=0.01,η=0.31,MS=0,ζ=0,G=0andB=0.01. The thermal conductivity stratification parameter avaries in the domain a∈(0,7.47). The dashed line represents the data fit ac∝L2with the factor of 8.80 ×102. 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 kc, 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, σ∗ T∗=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 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)presentthe eigenfunctions for the concentration disturbances ¯ φ(x,z)and those for the temperature ¯ T(x,z), respectively, both superimposed with the flow velocity field presented by the velocity vectors. In figure 23(a) of nanoparticle concentration ¯ φ(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 ¯ φand ¯ T, respectively, form stripes. This is explained by the fact that the functions φand θfound from the numerical solution 1011 A52-36 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics 1.0 0.8 0.6 0.4 0.2 0 1.0 0.8 0.6 0.4 0.2 10 20 x z (a)(b) z x 30 40 1.0 2 T ¯(x, z)φ ¯(x, z)|u¯(x, z)||u¯(x, z)| 0.6 0.4 0.2 0 1 0 –1 –2 0.6 0.4 0.2 0 0.5 0 –0.5 –1.0 010203040 Figure 23. Normalised eigenfunctions of the EVP (3.7) in the case of cooling at the substrate Q=−1for the critical wavenumber kc=0.16 with L=10−3,Φ=0.01,η=0.31,a=7.47,B=0.01,Σ 0≈2×104, G=6.71,MS=23,MT=0andλ=1.1707 ×10−7.The eigenfunctions ¯ φ(x,z)and ¯ T(x,z)superimposed with the velocity vector field ¯ u(x,z)are shown in panels (a)and(b), respectively. The velocity vector field ¯ 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). 1.0 0.8 0.6 0.4 0.2 0 1.0 0.8 0.6 0.4 0.2 10 20155 x z (a)(b) z 3025 35 010 20155 x 3025 35 0.2 T ¯(x, z)φ ¯(x, z)|u¯(x, z)||u¯(x, z)| 0 0.2 0.4 0.6 0.8 1.01.0 0.8 0.6 0.4 0.2 0 0.1 0 –0.1 –0.2 –4 –2 0 2 4 Figure 24. Normalised eigenfunctions of the EVP (3.7) in the case of heating at the substrate Q=1for the critical wavenumber kc=0.18 with L=10−3,Φ=0.01,η=0.31,a=7.47,B=0.01,Σ 0=10,G=1, MS=0,MT=1.22 and λ=2.0773 ×10−6+3.6965 ×10−5i.The eigenfunctions for the concentration ¯ φ(x,z)and temperature ¯ T(x,z)disturbances superimposed with the velocity vector field |¯ u(x,z)|,areshown in panels (a)and(b), respectively. of the eigenvalue problem (3.7) weakly vary with the height zwhen the former retains its sign, whereas the latter changes its sign. 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ϕ(z)and θ(z)are found to weakly depend on z, so both eigenfunctions ¯ φ(x,z)and ¯ 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 ¯ T(x,z), we note the formation of the upwelling flow in the highertemperature 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 φand θ, 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. 1011 A52-37 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron 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 ∇·u=−ρnp −1φm[Pφt+∇·(φu)].(4.1) 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 ∇·u=−Lρnp −1φm!∇2φ+∇·(ηφ∇T)".(4.2) 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 Lsatisfying the condition L∈(10−4,10−2), the parameter L(ρ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 MS(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 MTwith the thermal conductivity stratification parameter afor 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. 5. Summary and conclusions In this paper, we present a set of model governing equations and boundary conditions describing the dynamics of a moderately dense heated nanofluid layer with a deformable gas–liquid interface when the carrier fluid is Newtonian. This set of equations is based on continuity, momentum conservation, energy conservation and particle mass conservation equations, and represents an extension of the governing equations valid for dilute binary mixtures. Since, in the case at hand, the mixture is moderately dense, it is unrealistic to assume that its thermophysical properties, e.g., density, dynamic viscosity, thermal 1011 A52-38 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics 22.7 Compressible Simplification 70.0 67.5 65.0 62.5 60.0 57.5 55.0 52.5 50.0 22.6 22.5 22.4 22.3 MS (a)(b) MT 22.2 22.1 22.0 10–3 10–2 10–1 10–4 10–3 10–2 10–1 100 ka SS U UU† Monotonic Monotonic Oscillatory Oscillatory (simplification) Monotonic (simplification) Oscillatory Figure 25. (a) Neutral curves for a pure solutocapillary monotonic instability MS(k)for the case of cooling at the substrate Q=−1withB=0.01,Φ=0.01,η=0.31,a=7.47,L=10−3,MT=0,Σ 0≈2×104and G=6.71. The ◦and 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 afor a pure thermocapillary instability in the case of heating at the substrate Q=1with Φ=0.01,η=0.31,L=10−3,MS=0,ζ=0,G=0,and B=0.01. The ◦,,and ×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†,and Ssymbols represent the domains of monotonic instability, oscillatory instability, and stability of the system, respectively. conductivity and heat capacity, are constant. Therefore, here they are assumed to depend on the local particle concentration (Maron & Pierce 1956; Krieger & Dougherty 1959;de Kruif et al. 1985; Buongiorno 2006). We also assume that the Soret effect is present and the thermodiffusion coefficient is also local-particle-concentration dependent (Scriven & Sternling 1964; Platten & Legros 1984). We apply our equations to investigate the stability of a nanofluid layer open to the atmosphere at its deformable interface and subjected to a specified heat flux whether heating or cooling at its underneath support in the gravity field. We also assume that surface tension is both temperatureand particle-concentration dependent, so the thermocapillary and solutocapillary effects are present and accounted for. However, we assume that the considered layer is sufficiently thin, so buoyancy effects arising from an unstable temperature distribution with height is neglected. It is important to emphasise that since the fluid density in our model depends on the local particle concentration varying in both time and space, the system considered here is compressible, so the fluid velocity field is not solenoidal. We find a steady base state of the system which is written out analytically in terms of the Lambert W function. In the case of cooling at the substrate, the base state exhibits stable stratification in terms of both temperature and particle concentration, which affects directly the fluid density. In contrast, in the case of heating at the substrate, the temperature decreases with height, whereas the particle concentration increases with height; therefore, both conceive several instability mechanisms. We carry out the linear stability analysis of the base state of the system based on normal mode disturbances. We find that in the case of the system cooled at the substrate, the solutocapillary effect destabilises the system, whereas the thermocapillarity provides a stabilising effect. Interestingly, we note that the neutral curves vary weakly with the disturbance wavenumber and exhibit the emergence of two local minima, one of them 1011 A52-39 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press R. Gandhi, A. Nepomnyashchy and A. Oron located in the long-wave domain, whereas the other is in the finite-wave domain. These two minima compete with each other owing to the variation of the modified Galileo number G and the Soret coefficient η. We also find that only the finite-wave minimum of the neutral curves is sensitive to variation of the inverse capillary number which is related to the interfacial deformability. In the case of heating at the substrate, the linear stability properties of the system are by far more diverse than in the case of a cooled substrate. We find the emergence of both monotonic and oscillatory instabilities. The former is more typical for low values of the modified Galileo numbers Gequivalent to thinner nanofluid layers, mainly driven by thermocapillarity, whereas the latter is typical for moderate values of the Galileo number equivalent to thicker layers and displays a competition between thermocapillarity, solutocapillarity, and gravity. It is interesting to note that the low-Galileo number monotonic instability is long-wave with the critical wavenumber kcthat follows the wellknown scaling with the Biot number Bfor a pure fluid and is one of the two possible scalings found in the literature for a dilute binary mixture, namely kc∼B1/4.Further,we find that the long-wave thermocapillary instability is stabilised with an increase in the averaged bulk nanoparticle concentration Φand in the Soret coefficient η. We also reveal that the oscillatory instability for moderate values of Gis finite-wave. For sufficiently high values of G, the instability becomes again monotonic and long-wave driven solely by gravity, referred as to solutal buoyancy instability emerging from an unstable density stratification due to the number of particles which increases with height. We have also elucidated the details of the structure of a typical eigenspectrum by following the two leading eigenvalues for different values of the thermal Marangoni number MT. Among other details, we find that near the emergence of the monotonically growing mode driven by thermocapillarity, there exists another, slowly decreasing with the Marangoni number diffusional mode with a small growth rate which depends on the Lewis number and which is small. Near the inception point where the two modes are comparable, oscillations may emerge. When both thermocapillarity and solutocapillarity are active, the instability is monotonic for lower values of the solutal Marangoni numbers, whereas with an increase in the latter, the instability becomes oscillatory via the codimension-two point. Both of these instabilities occur at non-zero wavenumbers with the critical wavenumber for the oscillatory instability higher than that for the monotonic instability. As mentioned above, since the nanofluid density depends on the local nanoparticle concentration which is timeand space-dependent, the problem is essentially compressible, so the fluid velocity is not solenoidal. This fact poses additional difficulties in the numerical treatment of the linear eigenvalue problem solved to carry out the linear stability analysis. However, despite the fact that all of the results presented in this paper are obtained by solving the full formulation of the problem, we find that in many cases, the difference between the results obtained by solving the problem in its full version and, alternatively, in its simplified version in which the continuity equation is written in the form identical to that of an incompressible fluid is small. The presence of nanoparticles in a fluid causes a variation in the thermal conductivity of the nanofluid. As mentioned by Buongiorno (2006), thermal conductivity of the nanofluid in the case of alumina particles in water is linear with the local particle concentration with the parameter areferred to here as the thermal conductivity stratification parameter. Our results show that the value of aaffects the type of instability, namely, a pure thermocapillary instability is monotonic for small aand oscillatory when aexceed a certain critical value which is found to be proportional to the square of the Lewis number when the layer interface is non-deformable and in the absence of gravity. This oscillatory 1011 A52-40 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press Journal of Fluid Mechanics instability adds a new mechanism leading to the oscillatory instability induced by the thermal conductivity stratification in addition to other cases known in the literature (Joo 1995; Nepomnyashchy & Simanovskii 1995). Most of the instabilities found and discussed here are finite-wave instabilities. However, many of the instabilities, both monotonic and oscillatory dealt with in the similar setting but in the context of dilute binary mixtures (Oron & Nepomnyashchy 2004; Podolny, Oron & Nepomnyashchy 2006; Shklyaev et al. 2007,2009) with both constant Soret coefficient and the thermal conductivity of the mixture, were long-wave. To resolve this apparent mismatch, we emphasise that in our setting of a moderately dense mixture with concentration-dependent thermal conductivity and the Soret coefficient, the instabilities become also long-wave in the combined limit of mean particle concentration Φ, the Soret coefficient η, the thermal conductivity stratification parameter aand the Biot number Bbeing all very low. In this case, our present theory matches the critical values of the Marangoni number derived by Oron & Nepomnyashchy (2004) and Podolny et al. (2005) for the layers with either non-deformable or deformable interface. This issue will be further discussed in detail elsewhere. Finally, we emphasise that various analytical approximations have been derived over the years for the thermophysical properties of dilute monodisperse suspensions. At first order with respect to the low local particle concentration, they introduce factors which depend on the respective thermophysical properties of the base fluid and the suspended hard spherical particles. To mention the most prominent ones, these are the expressions for the viscosity (Einstein 1906), thermal conductivity (Maxwell 1873) and Brownian diffusivity (Batchelor 1976) of a suspension. In contrast, there are various empirical fits for these properties which are based on experiments conducted with different nanofluids. To bridge between the differences, we have employed the expressions arising from the experimental data, but also compared the results with those based on the theoretical expressions mentioned above. We have found that in terms of the thresholds of the thermosolutal instabilities investigated here, the differences appear to be very small within several percent and without any qualitative discrepancies. Acknowledgements. We are indebted to the anonymous referees whose valuable comments have contributed to the improvement of the quality of this paper. Funding. This work is supported by the grant from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska–Curie grant agreement number 955612 (NanoPaInt). Declaration of interests. The authors report no conflict of interest. Appendix A. Base state of the system A.1. Nanoparticle concentration-dependent Soret coefficient and thermal conductivity First, we note that since we seek a quiescent (u0=0) equilibrium state of the system, the fluid viscosity does not affect the result. To proceed with this task, we rewrite here for the reader’s convenience the system of equations and boundary conditions given by (2.21), (2.22), and (2.23) with prime denoting differentiation with respect to z: p 0=−G−Gφmρnp −1φ0,(A1a) (1+aφmφ0)T 0=0,(A1b) φ 0+ηφ0T 0=0,(A1c) z=0:(1+aφmφ0)T 0=−Q,φ  0+ηφ0T 0=0,(A2a) z=1:p0=0,(1+aφmφ0)T 0+BT0=0,φ  0+ηφ0T 0=0,(A2b) 1011 A52-41 https://doi.org/10.1017/jfm.2025.374 Published online by Cambridge University Press