Marangoni instabilities of drops of different viscosities in stratified liquids

Abstract

For an immiscible oil drop immersed in a stably stratified ethanol–water mixture, a downwards solutal Marangoni flow is generated on the surface of the drop, owing to the concentration gradient, and the resulting propulsion competes against the downwards gravitational acceleration of the heavy drop. In prior work of Li et al. (Phys. Rev. Lett., vol. 126, issue 12, 2021, 124502), we found that for drops of low viscosity, an oscillatory instability of the Marangoni flow is triggered once the Marangoni advection is too strong for diffusion to restore the stratified concentration field around the drop. Here we experimentally explore the parameter space of the concentration gradient and drop radius for high oil viscosities and find a different and new mechanism for triggering the oscillatory instability in which diffusion is no longer the limiting factor. For such drops of higher viscosities, the instability is triggered when the gravitational effect is too strong so that the viscous stress cannot maintain a stable Marangoni flow. This leads to a critical drop radius above which the equilibrium is always unstable. Subsequently, a unifying scaling theory that includes both the mechanisms for low and for high viscosities of the oil drops is developed. The transition between the two mechanisms is found to be controlled by two length scales: the drop radius $R$ and the boundary layer thickness $\delta$ of the Marangoni flow around the drop. The instability is dominated by diffusion for $\delta < R$ and by viscosity for $R<\delta$. The experimental results for various drops of different viscosities can well be described with this unifying scaling theory. Our theoretical description thus provides a unifying view of physicochemical hydrodynamic problems in which the Marangoni stress is competing with a stable stratification.

1. Introduction

Fluids in nature and technology are often multicomponent, with gradients in concentration when out of equilibrium (Lohse & Zhang 2020). The concentration gradients give rise to both density gradients and also solutal Marangoni stresses on an interface, thus leading to the competition between Marangoni convection and gravitational convection driven by density differences. Flows with such a competition are frequently encountered in technological applications. For example, in ink-jet printing, where the ink droplets are often multicomponent (Hoath 2016; Lohse 2022); in freezing emulsions (Ghosh & Coupland 2008; Degner et al. 2014; Deville 2017; Dedovets, Monteux & Deville 2018), where next to concentration gradients, thermal gradients also come into play; and in crystal growth (Chang & Wilcox 1976; Schwabe et al. 1978; Chang, Wilcox & Lefever 1979; Chun & Wuest 1979; Schwabe & Scharmann 1979; Schwabe, Scharmann & Preisser 1982; Preisser, Schwabe & Scharmann 1983; Kamotani, Ostrach & Vargas 1984), where the liquid column is subjected to a temperature gradient. The effect of gravity in these situations has often been ignored based on the perception that on small scales, gravitational effects are negligible owing to the small Bond number, which demonstrates the competition between gravity and capillary force (Nepomnyashchy, Legros & Simanovskii 2012). However, Edwards et al. (2018), Li et al. (2019a) and Diddens, Li & Lohse (2021) have recently discovered that gravitational effects can be important in evaporating binary droplets, in spite of the small scale reflected in small Bond numbers. It is therefore important to thoroughly understand the competing effect of gravity and Marangoni forces as some other flows where such a competition occurs may also need reconsideration.

The most basic example in which Marangoni forces compete with gravity is presumably an immiscible drop immersed in a stably stratified liquid, here made of an ethanol–water mixture. Such a simple droplet system allows for a systematic investigation of the parameter space of the concentration gradient, the drop size and the drop viscosity. Without gravitational effects (i.e. without a density stratification), the concentration gradient only leads to simple motion of the immersed drop (Lagzi et al. 2010; Maass et al. 2016; Jin, Krüger & Maass 2017; Dedovets et al. 2018). However, for a stable stratification, owing to density differences and gravity, intriguing and counter-intuitive phenomena can emerge, such as the oscillatory Marangoni flow around a drop/bubble (Zuev & Kostarev 2006; Viviani, Kostarev & Zuev 2008; Schwarzenberger et al. 2015; Bratsun et al. 2018) or the continuous bouncing of a drop (Blanchette & Shapiro 2012; Li et al. 2019b). For low viscosity (5 cSt) drops, it has been found that the continuous bouncing is actually caused by an oscillatory instability (Li et al. 2021) of the Marangoni flow, and the onset of the instability is triggered when the Marangoni advection is too strong for diffusion to restore the concentration field around the drop.

In this paper, we study the influence of the drop viscosity on the onset of the oscillatory instability in such a system by thoroughly exploring the three-dimensional parameter space spanned by the linearly stratified density gradient, the drop size and the drop viscosity. We start with very viscous silicone oil drops, for which we discover a new mechanism for the onset of the instability. We then develop a unifying scaling theory that includes both the newly discovered mechanism for high viscosity drops and the mechanism for low viscosity drops suggested by Li et al. (2021) and compare the results of this theory with the experimental results, which reveals good agreement. Finally, the instability thresholds for different viscosities are shown in one single phase diagram spanned by the Marangoni and Rayleigh numbers, thus clearly displaying the influence of the drop viscosity on the onset of the instability in the stratified flow.

2. Experimental procedure and methods

Linearly stratified ethanol–water mixtures were prepared in a cubic glass container (Hellma, 704.001-OG, Germany) with inner length $L=30\,\textrm {mm}$. Because the required liquid volume was very small (${<}27\,\textrm {ml}$), we used a slightly modified double-bucket method (Oster 1965) to generate the linear stratification, see figure 1(a). Syringe A was filled with a lighter liquid which was rich in ethanol (with ethanol weight fraction $w_{t}$), while syringe B was filled with a heavier liquid which was rich in water. The two syringes were filled to the same level and mounted on one syringe pump (Harvard Apparatus, PHD 2000, USA) to make sure the two liquids were injected at the same rate. The injection rate was no larger than 1 ml min$^{-1}$. During injection, the two liquids were constantly mixed in syringe A by a magnetic stirrer. The rod-like magnet in syringe A had a diameter of 4 mm. The resulting mixture flowed out of syringe A through a tube which led to the bottom of the glass container. A linear gradient was formed during the process, see figure 1(c). Strictly speaking, only an equal-mass flux rate could lead to a perfectly linear gradient in weight fraction. The equal-volume flux rate in our case and the volume reduction during the mixing of ethanol and water both led to an imperfect linear gradient. However, these imperfections were found to be negligible, see figure 2(e). Especially, these imperfections were negligible when compared to the small sizes of the drops, thus we still call it linear in the following text. Notice that by this method, the linear gradient could be varied continuously. However, in the experiments, we only chose several discrete values of the concentration gradient to span the parameter space.

Figure 1. Sketch of the experimental set-up. (a) The modified double-bucket method to generate the linearly stratified ethanol–water mixture in the glass container with inner length $L=30\,\textrm {mm}$. During injection, the two liquids are mixed by a magnetic stirrer. The resulting mixture flows out of syringe A through a tube which leads to the bottom of the container. A lid which holds the tube in position is used to cover the container to prevent evaporation. A linear gradient is then formed. A liquid layer of uniform ethanol concentration $w_{t}$ is first filled in the container before injection. (b) After the linear stratification has been formed, a layer of uniform ethanol concentration $w_{b}$ is injected at the bottom. Red indicates the mixture is rich in ethanol and blue rich in water. Silicone oil drops of different radii $R$ and different viscosities $\nu ^\prime$ are released in the stratified mixture. After releasing the drops, the container is again covered by a lid to prevent evaporation. Gravity is pointing downwards. (c) The ethanol weight fraction of the mixture $w_{e}$ linearly increases from $w_{b}$ at the bottom to $w_{t}$ at the top. (d) The density of the mixture $\rho$ linearly decreases from the bottom to the top. The density of the silicone oil is $\rho ^\prime$. The position $y=0$ marks the density-matched position where $\rho (w_{e}^\prime )=\rho ^\prime$. The height of the drop is $h$, counting from $y=0$.

Figure 2. (a) Successive snapshots of two 100 cSt silicone oil drops inside a linearly stratified ethanol–water mixture with $\mathrm {d}w_{e}/\mathrm {d}y\approx 40\,\textrm {m}^{-1}$, at times after they are released. The smaller drop ($R=54\pm 2\,\mathrm {\mu }\textrm {m}$) is levitating at $h\approx 4.8\,\textrm {mm}$, while the larger drop ($R=155\pm 2\,\mathrm {\mu }\textrm {m}$) is bouncing at $h\approx 3\,\textrm {mm}$. ($b$) Stream lines around a levitating drop ($R=143\pm 2\,\mathrm {\mu }\textrm {m}$) as revealed by the particle trajectories by superimposing 2300 images (covering 76.7 s). The red arrows indicate the flow directions. The surrounding mixture has a smaller concentration gradient $\mathrm {d}w_{e}/\mathrm {d}y\approx 10\,\textrm {m}^{-1}$. The scale bar is 0.5 mm. (c) Interfacial tension $\sigma (w_{e})$ between silicone oils of different viscosities and the ethanol–water mixture. Each point is an average of six measurements and the error bar is the standard deviation. The solid lines are polynomial fits to the data points. (d) Trajectories of 100 cSt silicone oil drops of different radii in a linearly stratified ethanol–water mixture. The trajectories all start from $t = 0\,\textrm {s}$. Here, $h = 0$ is the density-matched position, i.e. the position where $\rho =\rho ^\prime$. The filled circles represent the relative sizes of the drops. Drops with radius $R$ smaller than $66\,\mathrm {\mu }\textrm {m}$ are levitating and drops with radius $R$ larger than 89$\,\mathrm {\mu }$m are all bouncing. (e) The ethanol weight fraction $w_{e}$ of the mixture at corresponding height $h$ is measured by laser deflection. (f) The density $\rho$ of the mixture at corresponding height is calculated from the ethanol weight fraction $w_{e}$.

Before injection, the container was first filled with a liquid layer of uniform ethanol concentration $w_{t}$. After injection, a bottom layer of uniform ethanol concentration $w_{b}$ was injected with a third syringe, so that the ethanol weight fraction of the mixture $w_{e}$ increased continuously from $w_{b}$ at the bottom to $w_{t}$ at the top, see figure 1(b,c). The values of $w_{b}$ and $w_{t}$, and the depth of the linearly stratified layer were varied depending on the degree of stratification, which will be represented by $\mathrm {d}w_{e}/\mathrm {d}y$ in the following text. To avoid bubble formation during mixing, both ethanol (Boom B.V., 100 % (v/v), technical grade, the Netherlands) and Milli-Q water were degassed in a desiccator at ${\sim }$2000 Pa for 20 min before making the mixture. To reduce the disturbance to the linear gradient caused by the preferential evaporation of ethanol and subsequent Marangoni flows on the surface, the container was always covered by a lid except when releasing the drops. During injection, a lid which held the tube in place was used. Otherwise, a glass lid was used to cover the container.

The ethanol weight fraction $w_{e}(y)$ of the mixture at corresponding height $y$ was measured by laser deflection (Lin et al. 2013; Li et al. 2019b) immediately after the stratified mixture was made. The two uniform layers $w_{b}$ and $w_{t}$ were used to increase the accuracy of the laser deflection method. The density of the mixture $\rho (y)$ was calculated from $w_{e}(y)$ using an empirical equation (Khattab et al. 2012), and $y=0$ was set to the position where the density of the mixture $\rho (0)$ equalled that of the silicone oil $\rho ^\prime$, see figure 1(d). The ethanol weight fraction at this density-matched position $y=0$ was noted as $w_{e}^\prime$. As reported previously (Li et al. 2021), the maximum extent of the flow field induced by the drop became larger for weaker density stratifications, so that the container might not be large enough for weaker density stratifications. However, because we are not interested in the effect of the finite size of the container, experiments for weak density stratifications, i.e. smaller concentration gradients $\mathrm {d}w_{e}/\mathrm {d}y$, were repeated in a larger container (Hellma, 704.003-OG, Germany) with inner length $L=50\,\textrm {mm}$. Results that were different for the two different container sizes were discarded, thus all the results presented in this paper were not influenced by the finite size of the container.

Silicone oils (Sigma-Aldrich, Germany) of different viscosities $\nu ^\prime$ were used to generate drops of different radii $R$. They were injected from a $1\,\mathrm {\mu }$l syringe (Hamilton, KH7001) through a thin needle, whose outer-diameter was 0.515 mm. The properties of the silicone oils are listed in table 1. The drops were released from the top layer of the mixture, right after the laser deflection measurement. A collimated light-emitting diode (LED; Thorlabs, MWWHL4) was used to illuminate the drops in the mixture, and a camera (Nikon D850) connected to a long working distance lens system (Thorlabs, MVL12X12Z plus 0.25X lens attachment) was used to capture the side views of the drop. The radii and the trajectories of the drops were extracted from the sideview recordings. All images were recorded at 30 frames per second.

Table 1. Properties of the silicone oils used in the experiments.

For each concentration gradient, a large drop was first released. After it had bounced three times, the drop was carefully taken out with another thin needle with outer-diameter 0.515 mm, so as to minimize the disturbance to the stratified liquid. Then a second drop slightly smaller in size was released. In this way, we made sure that only one drop existed in the container at one time. This process was repeated until the smallest drop we could generate was released. One stratified mixture was used for no longer than 40 min. Otherwise, a new stratified liquid with the same concentration gradient was generated again to continue the process. The data were analysed after the first sweep in the drop radius $R$. If more data points were needed, the experiment was repeated at the same concentration gradient with the desired drop radii.

The flow field around the drop was revealed by adding tracer particles to the ethanol–water mixture. PSP (polyamide seeding particles, Dantec Dynamics, Denmark) particles of diameter 20$\,\mathrm {\mu }$m were added in both ethanol and water at 4 mg ml$^{-1}$ before making the linearly stratified mixture. These particles followed the flow faithfully, see Supplemental Material available at https://doi.org/10.1017/jfm.2021.983 for details. The interfacial tensions $\sigma$ between the silicone oils and the ethanol–water mixtures were measured by the pendant droplet method on a goniometer (OCA 15Pro, DataPhysics, Germany). Each measurement was repeated 6 times, and all the interfacial tension data are shown in figure 2(c). The solid lines are polynomial fits to the data points.

3. Results for 100 cSt silicone oil drops

For 100 cSt silicone oil drops, two typical states are observed after the initial sinking phase, that is, levitating or bouncing. In figure 2(a), we show the successive snapshots of two 100 cSt silicone oil drops in a mixture with $\mathrm {d} w_{e}/\mathrm {d}y\approx 40\,\textrm {m}^{-1}$ (see also Supplementary Movie 1). While the smaller drop ($R=54\pm 2\,\mathrm {\mu }\textrm {m}$) is levitating at $h\approx 4.8\,\textrm {mm}$, the larger drop ($R=155\pm 2\,\mathrm {\mu }\textrm {m}$) is bouncing at $h\approx 3\,\textrm {mm}$. The value $h=0$ marks the position where the density of the oil ($\rho ^\prime = 966\,\textrm {kg}\,\textrm {m}^{-3}$) equals that of the mixture (at $w_{e}^\prime =21.0\,\%$). The position of the drop $h$ as a function of time $t$ in the same stratified liquid, the ethanol weight fraction $w_{e}$ and the density of the mixture $\rho$ at the corresponding height are respectively shown in figure 2(d–f). The smallest drop ($R_1\approx 38\,\mathrm {\mu }\textrm {m}$) is levitating at $h\approx 5.9\,\textrm {mm}$. As the drop size increases, it levitates at a lower position until, above a critical radius $R_{cr}$, it starts to bounce instead of levitate. If its size is further increased, the drop bounces at a lower position (but still with $h>0$).

The smaller drop is able to levitate above the density-matched position ($h=0$) because of a stable Marangoni flow counteracting gravity, as revealed by the flow field around a levitating drop in a weaker gradient $\mathrm {d} w_{e}/\mathrm {d}y\approx 10\,\textrm {m}^{-1}$, see the particle trajectories shown in figure 2(b) (and also Supplementary Movie 2). The interfacial tension of the drop $\sigma$ decreases with increasing ethanol concentration of the mixture $w_{e}$, as shown in figure 2(c). The ethanol-concentration dependence of the interfacial tension of 20 and 50 cSt silicone oils, which will be used later in the paper, are also shown in this figure. This interfacial tension gradient at the surface of the drop pulls liquid downwards, generating a viscous force acting against gravity, which levitates the drop. The lighter liquids that are pulled down then flow sideways because of buoyancy. When the drop becomes large enough, however, the equilibrium becomes oscillatory and the drop starts to bounce between two different levels (see Supplementary Movie 3). Notice that the onset of this oscillatory instability is essentially the instability of the flow field around a drop which is not moving. Once the drop starts to move (bounce in this case), the flow field becomes drastically different because now the translational motion of the drop must also be considered. Consequently, the analysis of the bouncing behaviour is completely different from that of the instability. In this study, we are only interested in the onset of the instability, and the transition from a levitating drop to a bouncing one serves as a good indication.

Using this indication, we now explore the onset of this instability by independently varying $R$ and $\mathrm {d}w_{e}/\mathrm {d}y$. Because we cannot vary the two parameters continuously in the experiments, we use the following criterion to determine whether a drop is bouncing: if the bouncing amplitude $h_{A}$ of the drop is larger than its radius $R$, then the drop is considered to be bouncing (see Supplemental Material for more details). The results are shown in figure 3. The parameter space for $\mathrm {d}w_{e}/\mathrm {d}y\leq 40\,\textrm {m}^{-1}$ could not be explored because it was already contaminated by the finite size of the container. Different from the previously reported results of 5 cSt silicone oil drops (Li et al. 2021), where there is a critical gradient $(\mathrm {d} w_{e}/\mathrm {d}y)_{cr}$ above which the Marangoni flow is always unstable, no such critical gradient is observed for 100 cSt silicone oil drops up to $\mathrm {d}w_{e}/\mathrm {d}y\approx 150\,\textrm {m}^{-1}$. On the contrary, for all the tested gradients, there is always a critical radius $R_{cr}$ below which the flow is stable. This implies a different and new mechanism to trigger the Marangoni instability of 100 cSt drops compared with that of the 5 cSt drops. In the next section, we will develop a unifying scaling theory to account for both instability mechanisms, i.e. for low and high drop viscosity $\nu ^\prime$.

Figure 3. Phase diagram of the 100 cSt drops with a drop radius $R$ versus concentration gradient $\mathrm {d}w_{e}/\mathrm {d}y$ parameter space. Black triangles stand for bouncing drops, red circles for levitating ones. When experimentally determining the phase diagram, whether the drop bounces or not is decided within half an hour after its release to avoid any untraceable change of the gradient arising from long-term mixing induced by the drop motion itself.

4. A unifying scaling theory for the onset of the instability

For an immiscible drop immersed in the stably stratified ethanol–water mixture, a downwards Marangoni flow is generated on the surface of the drop, which also distorts the concentration field, see figure 4 for a sketch of the flow field and the concentration field around the drop. The system can be described in axisymmetrical spherical coordinates $(r, \theta )$ with its origin at the centre of the drop.

Figure 4. A sketch of the levitating drop (of radius $R$, density $\rho ^\prime$ and viscosity $\mu ^\prime$) and the flow field and ethanol concentration around it. Deeper red means higher ethanol concentration. The shaded ring inside the dashed circle represents the kinematic boundary layer with thickness $\delta$ set by the Marangoni flow, as represented by the solid arrows. The ethanol concentration inside this layer is enhanced by Maragnoni advection bringing down the ethanol rich liquid. The Marangoni flow velocity at the equator of the drop is $V_{M}$. The spherical coordinate $(r,\theta )$ has its origin at the centre of the drop. Here $\rho$ and $\mu$ are the density and viscosity of the liquid inside this layer, and $\rho ^*$ is the undisturbed density in the far field.

For the case of infinitely large solute diffusivity and zero density gradient, the analytical solution of the polar velocity $u_\theta$ outside the drop is given by (Young, Goldstein & Block 1959)

(4.1)\begin{equation} u_\theta={-}\frac{1}{4}\frac{1}{\mu+\mu^\prime}\frac{\mathrm{d} \sigma}{\mathrm{d} w_{e}} \frac{\mathrm{d} w_{e}}{\mathrm{d}y}R\left(\frac{R}{r}+\frac{R^3}{r^3}\right)\sin\theta, \end{equation}

where $\mathrm {d}\sigma /\mathrm {d}w_{e}$ is a material property (see figure 2c), $\mathrm {d}w_{e}/\mathrm {d}y$ the undisturbed ethanol gradient of the mixture, $\mu ^\prime$ the viscosity of the drop and $\mu$ the viscosity of the mixture. The Marangoni velocity at the equator of the drop is

(4.2)\begin{equation} u_\theta\vert_{r=R, \theta=90^{{\circ}}}={-}\frac{1}{2}\frac{\mathrm{d}\sigma}{\mathrm{d}w_{e}} \frac{\mathrm{d}w_{e}}{\mathrm{d}y}\frac{R}{\mu+\mu^\prime}. \end{equation}

In our case, the ethanol diffusivity $D$ and the density gradient both have a finite value, so the influence of advection and buoyancy cannot be neglected, see figure 2(b). Especially, the Marangoni advection tends to homogenize the concentration field around the drop (Li et al. 2021), effectively decreasing the concentration gradient $\mathrm {d}w_{e}/\mathrm {d}y$ around the drop, thus the Marangoni flow velocity is smaller than that predicted by (4.2). Let $V_{M}$ denote the Marangoni velocity at the equator of the drop in our case, i.e. $V_{M}=u_\theta \vert _{r=R, \theta =90^{\circ }}$, it is reasonable to drop the prefactor and write (also see Appendix A)

(4.3)\begin{equation} V_{M}\sim{-}\frac{\mathrm{d}\sigma}{\mathrm{d}w_{e}}\frac{\mathrm{d}w_{e}}{\mathrm{d}y} \frac{R}{\mu+\mu^\prime}. \end{equation}

The Marangoni flow imposes a viscous stress $\mu \nabla ^2\boldsymbol {u}\vert _{r=R}$ on the liquid close to the drop, which, in the axisymmetrical spherical coordinate, is expanded as

(4.4)\begin{equation} \mu\nabla^2\boldsymbol{u}\vert_{r=R}=\left.\mu\left[\frac{1}{r^2}\frac{\partial}{\partial r} \left(r^2\frac{\partial u_\theta}{\partial r}\right)+\frac{1}{r^2\sin\theta} \frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial u_\theta}{\partial\theta}\right)+ \frac{2}{r^2}\frac{\partial u_r}{\partial\theta}-\frac{u_\theta}{r^2\sin^2\theta}\right]\right\vert_{r=R}. \end{equation}

For the liquid close to the equator of the drop, i.e. at $\theta =90^{\circ }$, $u_\theta \vert _{r=R, \theta =90^{\circ }}=V_{M}$. Also notice that $u_r\vert _{r=R}=0$ and $\partial u_\theta /\partial \theta \vert _{\theta =90^{\circ }}=0$. The Marangoni viscous stress then becomes

(4.5)\begin{equation} \mu\nabla^2\boldsymbol{u}\vert_{r=R, \theta=90^{{\circ}}}\sim\mu \left(\frac{\partial^2V_{M}}{\partial r^2}+\frac{2}{R}\frac{\partial V_{M}}{\partial r}- \frac{V_{M}}{R^2}\right). \end{equation}

Let $\delta$ denote the boundary layer thickness of the Marangoni flow outside the drop, see figure 4. Because $\partial V_{M}/\partial r=\partial u_\theta /\partial r\vert _{r=R,\theta =90^{\circ }}\sim V_{M}/\delta$ and $\partial ^2 V_{M}/\partial r^2=\partial u_\theta /\partial r\vert _{r=R,\theta =90^{\circ }}\sim -V_{M}/\delta ^2$, we obtain

(4.6)\begin{equation} \mu\nabla^2\boldsymbol{u}\vert_{r=R, \theta=90^{{\circ}}}\sim\mu\left(-\frac{V_{M}}{\delta^2}+ \frac{2V_{M}}{\delta R}-\frac{V_{M}}{R^2}\right). \end{equation}

Two length scales are present in (4.6), $\delta$ and $R$, and they are both in the denominator. The next question is which of them is the more relevant length scale or, in other words, which of them is smaller. This can be inferred as follows: according to (4.3), an increase in the drop viscosity $\mu ^\prime$ would lead to a decrease in the Marangoni velocity $V_{M}$. Consequently, the Reynolds number $Re=\rho V_{M}R/\mu$ becomes smaller. Because within Prandtl–Blasius–Pohlhausen boundary layer theory (Schlichting & Gersten 2016), in general, the boundary layer thickness $\delta$ is inversely proportional to the square root of the Reynolds number, then $\delta$ is larger. That is to say, an increase in $\mu ^\prime$ would lead to an increase in $\delta$. Consequently, when $\mu ^\prime$ is very small, it is possible that $\delta < R$ so that $\delta$ is the more relevant length scale. However, when $\mu ^\prime$ is very large, it is possible that $\delta >R$ so that $R$ becomes the more relevant length scale. When $\mu ^\prime$ is taking some intermediate value, it is possible that $\delta \sim R$ so that the situation is more complicated and we are entering a transitional case. The two limiting cases and the transitional case are considered subsequently in the following subsections.

4.1. The limiting case when $\mu ^\prime$ is very small

When $\mu ^\prime$ is very small, $\delta < R$ so that (4.6) reduces to

(4.7)\begin{equation} -\mu\nabla^2\boldsymbol{u}\vert_{r=R, \theta=90^{{\circ}}}\sim\mu\frac{V_{M}}{\delta^2}. \end{equation}

This is the downward viscous stress acting on the liquid close to the equator of the drop. This liquid is lighter than its surroundings as it is entrained from the top. The buoyancy force acting on it is $g{\rm \Delta} \rho =g(\rho ^*-\rho )$, where $\rho$ is its density and $\rho ^*$ is the density of the undisturbed liquid in the far field, see figure 4. This liquid is brought down by the Marangoni flow along the surface of the drop, so ${\rm \Delta} \rho \sim -R\,\mathrm {d}\rho /\mathrm {d}y$. Now let us look at the boundary layer which is thin, $\delta < R$, i.e. all the liquid in the boundary layer is close to the drop. The boundary layer experiences a downward viscous force $\mu V_{M}/\delta ^2$ and an upward buoyancy force $-gR\,\mathrm {d}\rho /\mathrm {d}y$. For a levitating drop, the two forces balance, thus

(4.8)\begin{equation} \mu\frac{V_{M}}{\delta^2}\sim{-}gR\frac{\mathrm{d}\rho}{\mathrm{d}y}. \end{equation}

A very small $\mu ^\prime$ also leads to a very large $V_{M}$. This Marangoni advection acts to homogenize the concentration field inside the boundary layer while diffusion acts to restore it, see figure 4. When advection is too strong, the concentration field close to the drop will be smoothened, thus weakening the Marangoni flow itself, then the flow becomes unstable. The advection time scale is $\tau _{a}\sim R/V_{M}$ and the diffusion time scale is $\tau _{d}\sim \delta ^2/D$. The flow will become unstable when advection is faster than diffusion, i.e.

(4.9)\begin{equation} \tau_{a}<\tau_{d}. \end{equation}

Because diffusion is the limiting factor that causes this instability, we call it diffusion limited. Substituting the definitions of the two time scales into (4.9), we obtain

(4.10)\begin{equation} \frac{V_{M}R}{D}>\frac{R^2}{\delta^2}. \end{equation}

The left-hand side of (4.10) has the form of a Péclet number, which is referred to as the Marangoni number,

(4.11)\begin{equation} Ma=\frac{V_{M}R}{D}={-}\frac{\mathrm{d}\sigma}{\mathrm{d}w_{e}}\frac{\mathrm{d}w_{e}}{\mathrm{d}y}R^2 \frac{1}{(\mu+\mu^\prime)D}, \end{equation}

where we have used (4.3) with an equal sign. The instability criterion thus is

(4.12)\begin{equation} Ma>\frac{R^2}{\delta^2}. \end{equation}

Cancelling $\delta$ from (4.8) and (4.12), we obtain the diffusion-limited instability criterion for low viscosity $\mu ^\prime$ cases:

(4.13)\begin{equation} Ma/Ra^{1/2}>s, \end{equation}

where

(4.14)\begin{equation} Ra={-}\frac{\mathrm{d}\rho}{\mathrm{d}y}\, \frac{gR^4}{\mu D} \end{equation}

is the Rayleigh number for characteristic length $R$ and $s$ is a constant to be determined.

Equation (4.13) can also be expressed in dimensional quantities by substituting the definitions of $Ma$ and $Ra$:

(4.15)\begin{equation} \frac{\mathrm{d}w_{e}}{\mathrm{d}y}>\left(\frac{\mathrm{d}w_{e}}{\mathrm{d}y}\right)_{cr}=s^2 \left(\mu+\mu^\prime\right)^2\frac{gD}{\mu}\frac{\mathrm{d}\rho}{\mathrm{d}\sigma} \frac{\mathrm{d}w_{e}}{\mathrm{d}\sigma}. \end{equation}

Equation (4.15) actually predicts a critical concentration gradient above which the flow is always unstable.

Equations (4.13) and (4.15) were confirmed by results of 5 cSt silicone oil drops (Li et al. 2021) and the instability threshold is measured to be $s=275\pm 10$.

4.2. The limiting case when $\mu ^\prime$ is very large

When $\mu ^\prime$ is very large, $R<\delta$ so that (4.6) reduces to

(4.16)\begin{equation} -\mu\nabla^2\boldsymbol{u}\vert_{r=R, \theta=90^\circ}\sim\mu\frac{V_{M}}{R^2}. \end{equation}

Now the boundary layer is thick ($\delta >R$) and the force balance applied to the liquid close to the drop cannot be applied to the whole boundary layer. However, the liquid closest to the drop, or in other words the most inner layer of the boundary layer, still experiences the very same two forces: the downward viscous force, as shown by (4.16), and the upward buoyancy force $-gR\,\mathrm {d}\rho /\mathrm {d}y$. However, it also experiences a third force: an upward viscous force exerted by the outer layer. This is because the Marangoni velocity decreases radially outwards. Thus, for a stable Marangoni flow, we should have $\mu V_{M}/R^2>-gR\,\mathrm {d}\rho /\mathrm {d}y$. Otherwise, the Marangoni flow will become unstable because the viscous force on the inner layer cannot overcome its buoyancy. The instability criterion thus is

(4.17)\begin{equation} \mu\frac{V_{M}}{R^2}<{-}gR\frac{\mathrm{d}\rho}{\mathrm{d}y}. \end{equation}

Because viscosity is the limiting factor that causes this instability, we call it viscosity limited. Substituting (4.3) into (4.17) and reorganizing the left-hand side into the form of a Marangoni number, we obtain the viscosity-limited instability criterion

(4.18)\begin{equation} Ra/Ma>c, \end{equation}

where $c$ is a constant to be determined.

Equation (4.18) can also be expressed in dimensional quantities by substituting the definitions of $Ma$ and $Ra$:

(4.19)\begin{equation} R>R_{cr}=\sqrt{c\frac{\mu}{\mu+\mu^\prime}\frac{\mathrm{d}\sigma}{\mathrm{d}\rho}\frac{1}{g}}. \end{equation}

Equation (4.19) actually predicts a critical drop radius $R_{cr}$ above which the flow is always unstable. The fact that the diffusivity $D$ does not enter into this equation further indicates that this instability criterion is not dominated by diffusion. Note that the concentration gradient $\mathrm {d}w_{e}/\mathrm {d}y$ also does not enter into this equation, and the fluid properties $\mu$ and $\mathrm {d}\sigma /\mathrm {d}\rho$ depend on $w_{e}$ – the ethanol weight fraction at the levitation height – thus $R_{cr}$ is a function of $w_{e}$.

Before moving on, we want to stress that (4.17) is a force imbalance that only applies to the liquid closest to the drop and it does not determine the levitation height of the drop. For the levitation height of the drop, please see Appendix A. It is also found that the instability threshold $c$ should increase with $\mu ^\prime$, see Appendix B.

4.3. The transitional case when $\mu ^\prime$ is taking some intermediate value

In between the two limiting cases, we may also have a transitional case where $\delta \sim R$ so that the two mechanisms are intertwined. However, $R$ is an independent variable and an increase in $R$ also leads to an increase in the Reynolds number, thus decreasing $\delta$. Consequently, upon increasing the drop radius $R$ from very small to very large, the onset of the instability would change from the viscosity-limited regime to the diffusion-limited regime.

In the next section, we will compare the experimental results of drops of different viscosities ($\nu ^\prime =100$, 50 and 20 cSt) with the different cases mentioned above and measure the corresponding instability threshold $c$. Because it is found that $c$ is different for oils of different viscosities, we use $c_{100}$, $c_{50}$ and $c_{20}$ in the following text to represent the respective instability thresholds for 100, 50 and 20 cSt silicone oil drops. For general purposes, we will still use $c$, without specifying the drop viscosity.

To calculate the Marangoni and Rayleigh numbers, ethanol weight fractions at the positions where the drops levitate $w_{e}(h)$ are used to obtain the density $\rho$, viscosity $\mu$, diffusivity $D$ and the interfacial tension $\sigma$ (see Supplemental Material for the concentration dependence of $\rho$, $\mu$ and $D$). In the following, for bouncing drops, we use values corresponding to their highest positions.

5. Comparison with experimental results of different viscosities

The results of 100 cSt silicone oil drops, as shown in figure 3, are replotted in the $Ma$ versus $Ra$ parameter space in figure 5. As a comparison, the instability threshold of the diffusion-limited case (for 5 cSt silicone oil drops) $Ma/Ra^{1/2}\approx 275$ is shown as the red line. In the diffusion-limited case, no drop should be bouncing below the red line. However, while all the data points of 100 cSt drops are located below the red line, some of them still bounce. This indicates that diffusion is not the limiting factor to trigger the instability for the 100 cSt drops. Moreover, the blue dashed line which separates the bouncing drops from the levitating ones has a slope close to one (see (4.18)), which indicates that the instability is in the viscosity-limited regime. To further confirm this, we plot the results of 100 cSt drops in the $Ra/Ma$ versus $Ra$ parameter space in figure 6(a). As can be seen, there is indeed a critical value $(Ra/Ma)_{cr}$ above which the Marangoni flow is unstable, and the instability threshold is measured to be $c_{100}\approx 0.0177$ in the range $1\lesssim Ra\lesssim 5$. The critical radius $R_{cr}$ as a function of $w_{e}$ with $c_{100}=0.0177$ is calculated from (4.19) and shown as the blue curve in figure 6(b). The data shown in figure 6(a) are also replotted in figure 6(b) with ethanol weight fraction $w_{e}$ measured at the levitating height as the $x$-axis. As can be seen, the blue curve nicely separates the levitating drops and the bouncing ones. The dashed blue line in the range $w_{e}>77\,\textrm {wt}\%$ ($w_{e}<36\,\textrm {wt}\%$) corresponds to $Ra<1$ ($Ra>5$) and is not confirmed by experiments. We cannot confirm the region for $Ra>5$ because this would require much larger levitating drops for which container wall effects could affect the data. We also cannot confirm the range for $Ra<1$, i.e. $w_{e}>77\,\textrm {wt}\%$, because for concentration gradient $\mathrm {d}w_{e}/\mathrm {d}y$ larger than $150\,\textrm {m}^{-1}$, the drops tend to levitate close to the top layer ($w_{e}>77\,\textrm {wt}\%$), and the concentration gradient in this region is influenced by diffusion because it is in contact with the uniform top layer (see figure 1b). For $\mathrm {d}w_{e}/\mathrm {d}y=150\,\textrm {m}^{-1}$, the thickness of this region is ${\rm \Delta} h\approx 1.5\,\textrm {mm}$. The diffusion time scale to change the concentration in this region is $\tau \sim {\rm \Delta} h^2/D\approx 37\,\textrm {min}$. However, it normally takes ${\sim }20\,\textrm {min}$ to make sure that a small drop is indeed levitating. During this time, the concentration in this region would have changed owing to diffusion, which makes the measurement not accurate. For lower concentration gradients, the drops levitate in the middle of the linear region that is farther from the top layer. The concentration field in this region is not influenced by diffusion. Though in figure 6(a), only half a decade in the range of $Ra$ is confirmed by the experimental results, the confirmed range in figure 6(b) covers more than half of the available range: $21\,\textrm {wt}\%< w_{e}<100\,\textrm {wt}\%$ (because the drop cannot levitate below the density-matched position $w_{e}^\prime =21.0\,\textrm {wt}\%$).

Figure 5. Phase diagram of the 100 cSt drops replotted in the $Ma$ versus $Ra$ parameter space. Black triangles stand for bouncing drops, red circles for levitating ones. The red line is the diffusion-limited instability criterion $Ma/Ra^{1/2}=275$ for 5 cSt silicone oil drops, below which the drops of low viscosity are levitating. The red solid line in the range $6\times 10^{-3}\lesssim Ra\lesssim 3$ has been confirmed by the experiments of Li et al. (2021) and the red dashed lines are not. The blue line separates the bouncing drops from the levitating ones for the high viscosity drops of this present paper.

Figure 6. Levitating/bouncing phase diagrams for silicone oils of different viscosities: 100, 50 and 20 cSt. (a,c,e) $Ra/Ma$ versus $Ra$ phase diagrams for 100, 50 and 20 cSt silicone oils, respectively. As shown by the blue lines, the instability thresholds calculated by $(Ra/Ma)_{cr}$ (see (4.18)) are found to be $c_{100}\approx 0.0177$, $c_{50}\approx 0.0095$ and $c_{20}\approx 0.003$, respectively. Black triangles stand for bouncing drops and red circles for levitating ones. (b,d,f) Corresponding $R$ versus $w_{e}$ phase diagrams for 100, 50 and 20 cSt silicone oils, respectively. Black triangles stand for bouncing drops and red circles for levitating ones. The blue lines are calculated $R_{cr}$ with the corresponding instability thresholds from the left column of the figure. The blue solid lines in the range $36\,\textrm {wt}\%< w_{e}<77\,\textrm {wt}\%$, $39\,\textrm {wt}\%< w_{e}<86\,\textrm {wt}\%$ and $77\,\textrm {wt}\%< w_{e}<92\,\textrm {wt}\%$ are those confirmed by experiments in each case.

The parameter spaces for 50 and 20 cSt silicone oil drops are explored experimentally in the same way. Their interfacial tensions with the ethanol–water mixture are also measured in the same way, see figure 2(c). Their results are plotted in figure 6(c,d) for 50 cSt oil drops and figure 6(e,f) for 20 cSt oil drops. The instability threshold for 50 cSt silicone oil is measured to be $c_{50}\approx 0.0095$, and the blue line predicted from (4.19) with $c_{50}=0.0095$ again nicely separates the levitating drops from the bouncing ones. The dashed blue line in the range $w_{e}>86\,\textrm {wt}\%$ ($w_{e}<39\,\textrm {wt}\%$) corresponds to $Ra<1$ ($Ra>2.5$) and is not confirmed by experiments owing to the same reason as above. Again, while only a small range of $Ra$ in figure 6(c) is confirmed by the experimental results, the confirmed range in figure 6(d) covers more than half of the available range: $25\,\textrm {wt}\%< w_{e}<100\,\textrm {wt}\%$.

For 20 cSt silicone oil drops, no constant instability threshold is found, see figure 6(e,f), and we are probably entering the transitional case. Still, we cannot carry out experiments for $Ra<0.35$ and $Ra>16$ because of the same reasons as above. At $Ra\approx 0.35$ ($77\,\textrm {wt}\%< w_{e}<92\,\textrm {wt}\%$), the instability threshold is measured to be $c_{20}\approx 0.003$, as shown by the black solid line in figure 6(e) and the blue solid line in figure 6(f). In the range $Ra>0.35$ ($w_{e}<77\,\textrm {wt}\%$), the instability threshold deviates from $c_{20}\approx 0.003$.

To get a better understanding of this case, results of 20 cSt silicone oil drops are re-plotted in the $Ma/Ra^{1/2}$ versus $Ra$ parameter space in figure 7. The data points confirm that this is indeed the transitional case. The instability starts from the viscosity-limited regime $Ra/Ma=0.003$ when the drops are small ($Ra<0.35$), as shown by the blue dashed line. It slowly curves towards $Ma/Ra^{1/2}=275$ as the drop radius is increased ($0.35< Ra<3$) while still maintaining the viscosity-limited feature, as shown by the blue solid line. It finally turns into the diffusion-limited regime $Ma/Ra^{1/2}=275$ for much larger drops ($Ra>3$), as indicated by the red solid line.

Figure 7. Levitating/bouncing phase diagram for 20 cSt silicone oil drops. The instability threshold starts from the viscosity-limited regime $Ra/Ma=0.003$ for $Ra<0.35$ (blue dashed line) and it slowly changes to $Ma/Ra^{1/2}=275$ for $Ra>3$ (red solid line). The instability threshold changes from the viscosity-limited regime to the diffusion-limited regime in between (blue solid line). The solid lines are confirmed by experiments and the dashed lines are inferred from the trend of the experimental results.

We also performed experiments for 10 cSt silicone oil drops, where the onset of the instability becomes diffusion limited and the instability threshold is measured to be $s=265\pm 15$. This is indistinguishable from that of the 5 cSt drops $s=275\pm 10$ with the current experimental accuracy. We speculate that the instability threshold $s$ for the diffusion-limited case does not change with $\mu ^\prime$, because the balance between advection and diffusion happens in the bulk and the boundary layer thickness $\delta$ in (4.13) has cancelled out. However, the instability threshold $c$ for the viscosity-limited case is indeed found to increase with $\mu ^\prime$ (from $c_{20}\approx 0.003$ to $c_{100}\approx 0.0177$).

6. Instability thresholds for different viscosities

Now with all the results presented, we can have an overview of the instability thresholds for oils of different viscosities. They are plotted in figure 8. When the oil viscosity is very low, the instability is diffusion limited, which yields $Ma/Ra^{1/2}=s$. This is the case for 5 and 10 cSt silicone oil drops, and their instability thresholds are both $s\approx 275$. When the oil viscosity is very high, the instability is viscosity limited, which yields $Ra/Ma=c$. This is the case for 100 and 50 cSt silicone oil drops, and their instability thresholds are measured to be $c_{100}\approx 0.0177$ and $c_{50}\approx 0.0095$. When the oil viscosity is taking an intermediate value, the instability changes from the viscosity-limited regime to the diffusion-limited regime as the drop radius increases. This is the case for the 20 cSt silicone oil drops, whose instability thresholds are measured to be $c_{20}\approx 0.003$ and $s\approx 275$.

Figure 8. Phase diagram of the instability thresholds for different drop viscosities. The red line is the diffusion-limited instability threshold $Ma/Ra^{1/2}=275$ for 5 cSt silicone oil drops. Drops above the red line will be bouncing. The blue lines are the viscosity-limited instability thresholds, below which the drops will be bouncing. The instability threshold for 20 cSt silicone oil changes from $Ma/Ra^{1/2}\approx 275$ to $Ra/Ma\approx 0.003$. The instability thresholds for 50 and 100 cSt silicone oils are $Ra/Ma\approx 0.0095$ and $Ra/Ma\approx 0.0177$, respectively. Solid lines are confirmed by experiments, while the dashed lines are not.

7. Conclusions and outlook

In summary, the oscillatory instability of an immiscible oil drop immersed in a stably stratified ethanol–water mixture has been explored for oils of different viscosities. A unifying scaling theory has been developed, which predicts two different instability mechanisms depending on the two length scales in the system: the kinematic boundary layer thickness $\delta$ induced by the Marangoni flow and the drop radius $R$. (i) When $\delta < R$, the instability is triggered when advection is too strong so that diffusion cannot restore the concentration field around the drop. This is the diffusion-limited regime which yields $Ma/Ra^{1/2}=s$. (ii) When $R<\delta$, the instability is triggered when the gravitational effect is too strong so that the viscous stress cannot maintain a stable Marangoni flow. This is the viscosity-limited regime which yields $Ra/Ma=c$. An increase in the drop viscosity leads to an increase in $\delta$, which eventually moves the instability criterion from the pure diffusion-limited regime to the pure viscosity-limited regime, through a transitional case which connects both regimes. Because $\delta$ cancels out in the diffusion-limited regime, the instability threshold $s$ does not change with an increase in the drop viscosity $\mu ^\prime$, but the instability threshold $c$ in the viscosity-limited regime increases with increasing $\mu ^\prime$. The scaling theory is well supported by the experimental results.

Owing to the limitations of the scaling theory, we do not know how $\delta$ quantitatively depends on the drop viscosity $\mu ^\prime$ and drop radius $R$. Thus, we cannot predict beforehand in which regime will the instability be triggered. A quantitative understanding of $\delta$ may help to solve this problem. More detailed theoretical analysis may also be able to predict the instability thresholds beforehand. Nevertheless, our findings provide not only a deeper understanding to the simple droplet system under investigation, but also insight to other systems which allow for the direct competition between Marangoni and Rayleigh convection. For example, $Ra/Ma$ is called the dynamic Bond number in the context of Marangoni convections (Nepomnyashchy et al. 2012) $Bo_{d}=Ra/Ma=({(\mu +\mu ^\prime )}/{\mu })({\mathrm {d}\rho }/{\mathrm {d}\sigma })gR^2$. This number expresses the direct competition between Marangoni stress and gravity, as in the viscosity-limited regime. While the static Bond number $Bo_{s}={\rm \Delta} \rho gR^2/\sigma$ determines whether a drop/bubble maintains a spherical shape, the dynamic Bond number $Bo_{d}$ is more relevant for flows where the Marangoni flow is competing with a density gradient, such as in an evaporating binary sessile droplet (Edwards et al. 2018; Li et al. 2019a; Diddens et al. 2021).