2.1. Experimental set-up

Experiments were performed on sessile droplets deposited on a superhydrophobic substrate, consisting of a micro-structured glass slide coated with a fluorinated monolayer to increase its hydrophobicity. This yielded static macroscopic contact angles larger than 150 $^{\circ }$ , and very low roll-off angles (low hysteresis). The flow inside the droplets was visualised by fluorescently labelled polystyrene particles (provided by microParticles GmbH) with diameter approximately 1 ${\unicode{x03BC} }{\textrm {m}}$ , which makes them suitable as flow tracers. Such colloidal particles are stabilised by sulphate groups at their surface, which conveys stability to the suspension with no need of surfactants, which would interfere with our measurements. The flow was visualised and quantified by particle image velocimetry (PIV), for which we used a very low particle concentration ( $5\times 10^{-5}$ w/v %). Illumination was performed by a thin laser sheet, aimed at the central cross-section of the droplet. Imaging was performed with a Ximea CCD camera, mounted on a long-distance microscope, at 0.75 frames per second, which was enough to capture the slow motion of the particles. The PIV images were processed using our own algorithm. We used a multi-pass process with decreasing window size, from $128\times128$ pixels down to $32\times32$ pixels, with a 50 % window overlap in each pass.

Due to the spherical-cap shape of the droplet, the light coming from the particles into the sensor experience refraction at the liquid–air interface, therefore the particles appear at different positions in the image compared to in the drop. After the images had been processed with the PIV algorithm, we applied an optical correction using ray-tracing to the centres of the interrogation windows and the measured displacements. Velocities that were measured close to the edges of the drop ( $r/R\gt 0.9$ ) were omitted, since there the optical distortion becomes very large, leading to unreliable results.

The droplet was carefully deposited over the superhydrophobic substrate inside a closed (not pressurised) chamber. The substrate was previously cleaned by rinsing it once in water and carefully drying it with a nitrogen gas flow. To reduce the impact of contamination from the air in the room, the chamber was cleaned with a wet cloth, and it then remained closed until the experiment was performed. The temperature inside the chamber was not controlled, but remained constant at 21 $^{\circ }$ C throughout the whole experiment. The relative humidity was passively controlled using reservoirs of pure water or salt-saturated water solutions in the chamber, which allowed us to vary the relative humidity from $35\,\%$ to $90\,\%$ . The relative humidity and temperature were measured inside the chamber using a sensor (HGC 30 from DataPhysics) near the droplet.

2.2. Experimental results

Figures 1(a–c) show three different experiments conducted at different relative humidity, namely 90 %, 74 % and 50 %, therefore having different evaporation rates. In these experiments, the flow starts as a single convection roll vortex within the droplet (upper row of figure 1), eventually transitioning to an axisymmetric flow pattern (lower row of figure 1). The initial single-roll flow pattern remains stable for several minutes (see supplementary movies 1–3), indicating that it is not merely a transient phenomenon.

Figure 1. Experimental measurements of the flow in an evaporating drop. Each column corresponds to a measurement at a different relative humidity ( $\textit{RH}$ ). The plane $y=0$ is shown. The evaporation speed increases from left to right. Snapshots at two different times are shown. Initially, the flow is asymmetric, with a single roll in a different direction for each experiment (upper row). After some time, the flow becomes axisymmetric in all experiments, with a flow from the apex towards the contact line, indicating dominant thermal buoyant forces (lower row).

The direction of the observed axisymmetric interfacial velocity field in the late phase, from the apex towards the contact line, and directed upwards at the symmetry axis, suggests that the flow is driven by a thermal buoyant circulation (from here on, axisymmetric Rayleigh flow), and not a thermally driven Marangoni flow (from here on, axisymmetric Marangoni flow).

However, such a flow pattern transition is not always observed, even though the experiments are run under identical experimental conditions. The flow transition illustrated in figure 1 occurs in 40 % of the experiments done (a total of 50 experiments), and an additional 43 % of the experiments show either only a single-roll flow pattern or only the axisymmetric Rayleigh flow, with no observable transition. A minority of them (4 %) showed an inverted transition, from axisymmetric Rayleigh flow to single roll. Finally, a small but not negligible fraction of the experiments (13 %) showed a transition from an axisymmetric Marangoni flow, characterised by a downward-directed flow in the symmetry axis, into an axisymmetric Rayleigh flow, characterised by upward-directed flow. The variety of observed flow patterns reveals the high sensitivity of the transition due to the presence of contaminants in the environment. Despite the diversity of flow patterns observed, the most frequently observed flow patterns are clearly dominated by thermally driven convective flows, such as those shown in figure 1.

Such experimental results contrast with the theoretical predictions (Hu & Larson Reference Hu and Larson2005; Sáenz et al. Reference Sáenz, Sefiane, Kim, Matar and Valluri2015; Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b ) that suggest that thermal Marangoni flows should be dominant, manifesting interfacial flows directed from the contact line towards the apex, and downward-directed flow in the symmetry axis.

This feature resembles the case of sessile droplets with lower contact angles, in which thermally driven Marangoni flows are expected to develop flow patterns even stronger than the capillary flows. However, in practice, such thermally driven flows are rather weak and only inconsistently observed (Hu & Larson Reference Hu and Larson2005; Gelderblom et al. Reference Gelderblom, Diddens and Marin2022). The most accepted explanation of this surprising feature is that the presence of contaminants at the surface hinders the thermally driven surface tension gradient at the interface.

Similarly, we hypothesise that the presence of contaminants at large contact angles is the cause of the dominance of convective flows, in the same way that capillary flows dominate the evaporation-driven flows in low contact angle sessile droplets (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten1997; Gelderblom et al. Reference Gelderblom, Diddens and Marin2022). The presence of contaminants also explains the large diversity of flow patterns observed in our experimental results. In the experiments of figure 1, the time of the transition from single roll to axisymmetric Rayleigh flow is significantly different. We theorise that this is due mainly to different amounts of contaminants being present. Experimentally, we have not observed a strong dependence of the flow pattern on the relative humidity.

In the remainder of the paper, we explore the influence of the presence of contaminants at the droplet’s interface by making use of numerical methods, as has been done in previous studies (Savino et al. Reference Savino, Paterna and Favaloro2002; Hu & Larson Reference Hu and Larson2005, Reference Hu and Larson2006; Xu & Luo Reference Xu and Luo2007; Girard et al. Reference Girard, Antoni and Sefiane2008; van Gaalen et al. Reference van Gaalen, Wijshoff, Kuerten and Diddens2022).

3.1. Transient model

We consider a simple model that captures the essential physics of the experiments. A similar model has been used by Diddens et al. (Reference Diddens, Li and Lohse2021) for the prediction of the flow inside an evaporating binary glycerol–water droplet. Here, however, the temperature field inside a pure water droplet is considered, rather than the compositional change of the species.

A droplet of initial volume $V_0$ , density $\rho$ , viscosity $\mu$ , thermal conductivity $k$ and specific heat capacity $c_p$ is deposited on a substrate, forming a contact angle $\theta$ ; see figure 2. The substrate is approximated to be perfectly thermally conductive, such that its temperature is set to the ambient temperature $T_0$ , which is far below the boiling point. Considering finite thermal conductibility would result in lower thermal gradients within the droplet, weakening the thermal Marangoni effects. However, our numerical tests with the experimental thermal conductibility of the substrate suggest that the interfacial velocity remains of the same order of magnitude as by considering a perfectly thermally conductive substrate, thereby justifying neglecting this effect. The droplet is surrounded by air at $T_0$ and relative humidity $RH = c_{\infty }/c_{{sat}}$ , where $c_\infty$ is the vapour concentration far from the droplet, and $c_{{sat}}$ is the vapour saturation concentration at atmospheric pressure and constant $T_0$ . All liquid properties are temperature-dependent.

Figure 2. Schematics of the model used in the numerical simulations.

We consider the gas surrounding the water to be in a quiescent state. The solutal and thermal Rayleigh numbers of the water vapour around the droplet are defined as (Dietrich et al. Reference Dietrich, Wildeman, Visser, Hofhuis, Kooij, Zandvliet and Lohse2016) ${Ra}_c^g = (\beta _c^g g V)/(\mu ^g \kappa _c^g)$ and ${Ra}_T^g = (\beta _T^g g V)/(\mu ^g \kappa _T^g)$ , respectively. Here, $\beta _c^g = (\partial \rho ^g / \partial c)/\rho ^g$ is the solutal expansion coefficient, $\rho ^g$ is the gas density, $c$ is the gas vapour concentration, $g$ is gravitational force, $\mu ^g$ is the dynamic viscosity of the gas, $\kappa _c^g$ is the diffusivity of the solutal vapour, $\beta _T^g = (\partial \rho ^g / \partial T)/\rho ^g$ is the thermal expansion coefficient, and $\kappa _c^g$ is the diffusivity of the thermal vapour. When calculating ${Ra}_c$ and ${Ra}_T$ for the gas and vapour around the droplet, we obtain ${Ra}_s = 1.2$ and ${Ra}_T = 0.5$ , both much smaller than the critical Rayleigh value 12 for convective effects to become relevant. Since there are no external heat sources, the water vapour concentration $c$ can be considered to quasi-stationarily diffuse far from the droplet (see (3.1)) and the evaporation rate $j$ is given by the diffusive flux at the liquid–gas interface (see (3.2)) (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten1997, Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000; Deegan Reference Deegan2000; Popov Reference Popov2005):

(3.1) \begin{align}&\qquad \nabla ^2 c = 0 , \end{align}

(3.2) \begin{align}& j = - D^{{vap}}\, \boldsymbol \nabla c \boldsymbol \cdot \boldsymbol{n} , \end{align}

where $D^{{vap}}$ is the vapour diffusion coefficient, and $\boldsymbol{n}$ is the normal to the interface. As the droplet evaporates, the droplet–gas interface moves according to the kinematic boundary condition

(3.3) \begin{equation} \rho (\boldsymbol{u} - \boldsymbol{u}_{{i}}) \boldsymbol \cdot \boldsymbol{n} = j. \end{equation}

Here, $\boldsymbol{u}$ is the velocity field, and $\boldsymbol{u}_{{i}}$ is the velocity of the interface. On account of the latent heat of evaporation $\Lambda$ , evaporative cooling occurs. Due to the presence of the isothermal substrate and the spatially non-uniform evaporation rate, temperature gradients along the interface are expected. In the liquid phase, if convective motion is initiated through e.g. thermal Marangoni forces, then it is expected to have significant effects on the droplet temperature profile due to the low thermal diffusivity $\kappa = k/(\rho c_p)$ of water (or any similar liquid), which is of order $\kappa \sim$ ${10^{-7}}\ {\textrm m^2\ \textrm s^{-1}}$ . Therefore, we consider a convection–diffusion equation for the temperature $T$ in the liquid phase,

(3.4) \begin{equation} \rho c_p ( \partial _t T + \boldsymbol{u} \boldsymbol\cdot \boldsymbol \nabla T ) = k\, \nabla ^2 T , \end{equation}

with a boundary condition for the evaporative cooling in a frame co-moving with the interface,

(3.5) \begin{equation} k\, \boldsymbol \nabla T \boldsymbol \cdot \boldsymbol{n} = - j \Lambda , \end{equation}

and $T=T_0$ at the substrate. Thermal transport in the gas phase is disregarded due to its small contribution in the droplet’s flow features.

We consider the $z$ -axis to be the vertical axis pointing to the apex of the droplet, and gravity $g$ acts in the negative $z$ -direction for a sessile droplet, a configuration to which this study is restricted. The velocity field and pressure $p$ are governed by the Navier–Stokes equation

(3.6) \begin{equation} \rho ( \partial _t \boldsymbol{u}+ \boldsymbol{u} \boldsymbol \cdot \boldsymbol \nabla \boldsymbol{u} ) = - \boldsymbol \nabla p + \mu\, \nabla ^2 \boldsymbol{u} + \rho \boldsymbol{g}, \end{equation}

where $\boldsymbol{g} = - g \boldsymbol{e}_z$ is the gravitational acceleration ( $\boldsymbol{e}_z$ is the unit vector in the $z$ -direction, positive upwards). The mass conservation is ensured by the continuity equation

(3.7) \begin{equation} \partial _t \rho + \boldsymbol \nabla \boldsymbol \cdot (\rho \boldsymbol{u})=0 . \end{equation}

Similarly to what is considered by Diddens et al. (Reference Diddens, Li and Lohse2021), we take a small slip length to simulate the assumed pinned contact line in Constant Contact Radius mode (Picknett & Bexon Reference Picknett and Bexon1977) to resolve the incompatibility of no-slip at the substrate and evaporation at the contact line. As long as the slip length is orders of magnitude smaller than the droplet size, the results are independent of the slip length (Diddens et al. Reference Diddens, Li and Lohse2021), which is also confirmed by our simulations.

The liquid–gas interface must account for the Laplace pressure and Marangoni stresses, which are given by

(3.8) \begin{equation} \boldsymbol{n} \boldsymbol \cdot ( - p \boldsymbol{I} + \mu (\boldsymbol \nabla \boldsymbol{u} + (\boldsymbol \nabla \boldsymbol{u})^{\textrm T}) ) \boldsymbol \cdot \boldsymbol{n} = \sigma C \boldsymbol{n} \end{equation}

and

(3.9) \begin{equation} \boldsymbol{n} \boldsymbol \cdot (\mu (\boldsymbol \nabla \boldsymbol{u} + (\boldsymbol \nabla \boldsymbol{u})^{\textrm T}) ) \boldsymbol \cdot \boldsymbol{t} = \boldsymbol \nabla _t \sigma \boldsymbol \cdot \boldsymbol{t} \, . \end{equation}

Here, $\boldsymbol{t}$ is the tangential vector to the interface, $\boldsymbol{I}$ is the identity matrix, $C$ is the curvature of the interface, and $\boldsymbol \nabla _t$ is the surface derivative at the interface.

3.2. Quasi-stationary model

In a droplet evaporating at ambient conditions, the velocity $\boldsymbol{u}$ associated with the circulation within the droplet is typically much faster than the velocity of the interface $\boldsymbol{u}_{{i}}$ (Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b ). Buoyancy and thermal Marangoni forces form a quasi-equilibrium, thereby justifying the quasi-stationary approximation. The kinematic boundary condition (3.3) is then reduced to $\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} = 0$ , which can be imposed numerically via a Lagrange multiplier field at the interface. Due to its small size, the temperature variations within the droplet are typically a few per cent of the absolute ambient temperature. It is then reasonable to consider an expansion of the temperature into $ T(\boldsymbol{x},t) = T_0 + \Delta T(\boldsymbol{x},t)$ , i.e. the ambient temperature summed with the small local variations $\Delta T(\boldsymbol{x},t)$ . Taking these approximations, it is logical to expand the density and surface tension into first-order Taylor series around the ambient temperature $T_0$ :

(3.10) \begin{align} \rho (T) &= \rho _0 + (T - T_0)\, \left . {\partial _T \rho } \right |_{T_0} , \end{align}

(3.11) \begin{align} \sigma (T) &= \sigma _0 + (T - T_0)\, \left . {\partial _T \sigma } \right |_{T_0} , \end{align}

where $\rho _0$ and $\sigma _0$ are the density and surface tension at $T_0$ , and $\left . {\partial _T \rho } \right |_{T_0}$ and $\left . {\partial _T \sigma } \right |_{T_0}$ are the first-order derivatives of the density and surface tension with respect to the temperature, evaluated at $T_0$ . We simplify the notation to $\partial _T \rho$ and $\partial _T \sigma$ in the following. All other properties are considered to be constant and are evaluated at $T_0$ , i.e. $\mu =\mu _0$ , $k=k_0$ , $c_p=c_{p0}$ , $D^{{vap}}=D^{{vap}}_0$ , $\Lambda = \Lambda _0$ , $c_{{sat}}=c_{{sat0}}$ – that is, we assume the Oberbeck–Boussinesq approximation for these properties.

We scale the vapour concentration such that $\tilde {c}=1$ at the droplet–gas interface, and $\tilde {c}=0$ far away. For the temperature differences $T-T_0$ , one should take into consideration that no temperature difference should be expected when no evaporation happens, i.e. when relative humidity is 100 % (or $c_{{sat}}=c_\infty$ ). We choose the factors in evaporative cooling (3.5) for the scaling of the temperature variation:

(3.12) \begin{equation} \begin{aligned} (c-c_\infty ) &= (c_{{sat}}-c_\infty )\tilde {c} , \quad (T-T_0) = \dfrac {(c_{{sat}}-c_\infty ) \Lambda _0 D^{{vap}}_0}{k_0} \tilde {T} . \end{aligned} \end{equation}

We take the cubic root of the volume to non-dimensionalise length scales. The thermal diffusivity $\kappa _0 = k_0/(\rho _0 c_{p0})$ is taken to non-dimensionalise time, since thermal effects dictate the flow. We then have

(3.13) \begin{equation} \begin{aligned} \boldsymbol{x} = V^{1/3} \tilde {\boldsymbol{x}}, \quad t = \frac {V^{2/3}}{\kappa _0} \tilde {t} , \quad \boldsymbol{u}= \frac {\kappa _0}{V^{1/3}} \tilde {\boldsymbol{u}}. \end{aligned} \end{equation}

A dimensionless mass transfer rate $\tilde {j}$ can then be defined as $\tilde {j} = V^{1/3}/(D^{vap}_0 (c_{sat}-c_\infty )) j = \boldsymbol \nabla \tilde {c} \boldsymbol \cdot \boldsymbol{n}$ . In non-dimensional form, (3.4) and (3.5) can be rewritten as

(3.14) \begin{align}& \partial _{\tilde {t}} \tilde {T} + \tilde {\boldsymbol{u}} \boldsymbol \cdot \tilde {\boldsymbol \nabla } \tilde {T} = \tilde {\nabla }^2 \tilde {T} , \end{align}

(3.15) \begin{align}&\qquad \tilde {\boldsymbol \nabla } \tilde {T} \boldsymbol \cdot \boldsymbol{n} = - \tilde {j} . \end{align}

The Oberbeck–Boussinesq approximation is taken for the Navier–Stokes equations since $(T-T_0)\, \partial _T \rho$ is small compared to $\rho _0$ . The density appears as in (3.10) only in the gravity term $\rho \boldsymbol{g}$ , while considered constant $\rho _0$ in the remaining terms. Due to the small Reynolds number ( $\sim 0.1$ at experimental conditions), Stokes flow is assumed to be sufficient. Equations (3.6) and (3.7) can then be reduced to

(3.16) \begin{align}&\qquad\qquad\quad \tilde {\boldsymbol \nabla } \boldsymbol \cdot \tilde {\boldsymbol{u}} = 0 , \end{align}

(3.17) \begin{align}& - \tilde {\boldsymbol \nabla } \tilde {p} + \tilde {\nabla }^2\boldsymbol{u} + {Ra}_T\, \tilde {T} \boldsymbol{e}_z = 0, \end{align}

where

(3.18) \begin{equation} \tilde {p} = \frac {V^{2/3}}{\kappa _0 \mu _0} (p - \rho _0 g z) \, \end{equation}

is the non-dimensional pressure, and

(3.19) \begin{equation} {Ra}_T = \frac {V g\, |\partial _{\tilde T} \rho |}{\kappa _0 \mu _0} \end{equation}

is the thermal Rayleigh number.

To remove the null space of the pressure field, we use a Lagrange multiplier, which ensures that the average pressure $\tilde {p}$ in the liquid domain is zero. Note that we derive $\rho$ with the non-dimensional temperature $\tilde T$ to ensure that $RaT=0$ in the absence of evaporation, e.g. when $c_{{sat}}=c_\infty$ . The non-dimensional capillary (3.8) and Marangoni (3.9) forces are given by

(3.20) \begin{align}& \tilde {p} + \boldsymbol{n} \boldsymbol \cdot (\tilde {\boldsymbol \nabla } \tilde {\boldsymbol{u}} + (\tilde {\boldsymbol \nabla } \tilde {\boldsymbol{u}})^{\textrm T}) \boldsymbol \cdot \boldsymbol{n} = \frac {1}{Ca} (\tilde {C} + Bo\, \tilde {z}) , \end{align}

(3.21) \begin{align}&\quad \boldsymbol{n} \boldsymbol \cdot (\tilde {\boldsymbol \nabla } \tilde {\boldsymbol{u}} + (\tilde {\boldsymbol \nabla } \tilde {\boldsymbol{u}})^{\textrm T}) \boldsymbol \cdot \boldsymbol{t} = - {Ma}_T\, \tilde {\boldsymbol \nabla }_t \tilde {T} \boldsymbol{\cdot} \boldsymbol{t} , \end{align}

where the capillary number $Ca$ , the Bond number $Bo$ and the thermal Marangoni number ${Ma}_T$ are given by

(3.22) \begin{equation} \begin{aligned} Ca = \frac {\kappa _0 \mu _0}{V^{1/3} \sigma _0}, \quad Bo = \frac {\rho _0 g V^{2/3}}{\sigma _0} , \quad {Ma}_T = \frac {V^{1/3} \,|\partial _{\tilde T} \sigma |}{\kappa _0 \mu _0} . \end{aligned} \end{equation}

In small droplets, the capillary number is small ( ${\sim}10^{-6}$ at the experimental conditions). Since the Bond number is also small ( ${\sim}10^{-1}$ at the experimental conditions), the last term of (3.20) is assumed to be irrelevant, leading to an equilibrium spherical-cap shape with a constant contact angle $\theta$ .

In our model, we have considered only the case of decreasing density and surface tension with increasing temperature, which is the case for water at ambient lab conditions. This allows us to consider only positive ${Ra}_T$ and ${Ma}_T$ values, following the sign convection in the definition of the system of equations.

3.3. Procedure

The system of equations of §§ 3.1 and 3.2 are solved using a finite element method on triangular elements with the package pyoomph (https://pyoomph.github.io/) (Diddens & Rocha Reference Diddens and Rocha2024), which is based on oomph-lib (Heil & Hazel Reference Heil and Hazel2006) and GiNaC (Bauer et al., Reference Bauer, Frink and Kreckel2002). The code was mutually validated by an analogous implementation in ngsolve (Schöberl Reference Schöberl1997). The mesh motion in the transient model of § 3.1 is modelled by using an arbitrary Lagrangian–Eulerian method. Since the droplet geometry is assumed to be axially symmetric, a two-dimensional mesh is used in all simulations, even when analysing the azimuthal symmetry-breaking phenomenon (§ 4.3). This can be done by assuming that the axial symmetry is broken in a bifurcation, due to instabilities that act on an azimuthal wavenumber $m$ . We employ the method described and validated by Diddens & Rocha (Reference Diddens and Rocha2024), a linear stability analysis method that allows investigation of symmetry-breaking instabilities. Together with bifurcation tracking and pseudo-arc-length continuation methods, this tool allows identification of the critical parameters for which the symmetry is broken.

The vapour diffusion equation (3.1) in the gas phase is discretised using first-order continuous Lagrangian elements. Its solution then provides the value of the evaporation rate $\tilde {j}$ at the interface. In order to spatially approximate the velocity–pressure pair in the Stokes equation, MINI-elements (Arnold, Brezzi & Fortin Reference Arnold, Brezzi and Fortin1984) are used. The temperature field is approximated using first-order continuous Lagrangian elements. The choices of the latter two spaces are made in order to avoid spurious oscillations in the tangential velocity at the droplet–gas interface. Further details concerning the choice of field spaces are provided in Appendix A.

In § 4.2, we calculate the percentage of the total volume of liquid associated with the stream function $\psi$ that flows in the direction from the apex towards the contact line. This percentage is denoted as $\psi ^+$ . The value of $\psi ^+$ provides insight into the dominant forces driving the flow. Specifically, if $\psi ^+ \gt 90\,\%$ , then the flow is considered to be dominated by buoyancy forces, and if $\psi ^+ \lt 10\,\%$ , then the flow is considered to be dominated by thermocapillary forces. By fixing ${Ma}_T$ , we can formally consider ${Ra}_T$ as a Lagrange multiplier to obtain the unique ${Ra}_T$ value that satisfies each $\psi ^+$ constraint. Subsequently, we employ pseudo-arc-length continuation to trace the curves that delineate the regime where there is competition between Marangoni and Rayleigh effects.

4.1. At experimental conditions: analysis of the axisymmetric flow direction

We take the initial experimental conditions of figure 1(c) as input for the numerical simulations to compare the results. More specifically, we consider here: $RH=50\,\%$ , $T_0=21\,^\circ$ C, $V_0={5.6}\ {\unicode{x03BC} }{\textrm l}$ and $\theta =150^\circ$ . The transient model of § 3.1 is used here to describe the time evolution of the evaporating droplet. The Young–Laplace equation is solved to obtain the slightly gravity-deformed droplet shape at deposition time, used as the initial geometry condition for our simulations.

Figure 3. Numerical results for the temperature and velocity magnitude fields at the initial experimental conditions of figure 1(c), i.e. $RH=50\,\%$ , $T_0=21^\circ$ C, $V_0={5.6}\ {\unicode{x03BC} }\text{l}$ and $\theta =150^\circ$ , 1500 s after deposition on the substrate. (a) Considering all effects generating the flow, the flow direction is axisymmetric, going from the rim to the apex, indicating dominant thermocapillary forces. (b) If the surface tension is considered constant, then the flow direction of the axisymmetric flow is the other way around, namely from the apex to the rim, indicating dominance of the buoyancy forces.

We analyse the flow direction obtained numerically in a time frame where the experiments present an axisymmetric flow from the apex towards the contact line, in particular 1500 s after deposition, as depicted in figure 1(cii). At this instant, the numerical calculations show a vortex in the droplet from the contact line towards the apex, indicating dominant thermocapillary forces; see figure 3(a). The velocity magnitude is of the order of $\sim {10}\ {\textrm {mm s}}^{-1}$ , which is significantly larger than the experimental measurements, $\sim {1}\ {\unicode{x03BC} }{\textrm {m}\ {\textrm s}}^{-1}$ . One can theoretically consider a constant surface tension to disregard any Marangoni flow; see figure 3(b). Here, despite the theoretical assumption, the flow direction agrees well with the experimental results, and the velocity magnitude is much closer to that of the experiments, $\sim {10}\ {\unicode{x03BC} }{\textrm{m s}}^{-1}$ . The discrepancy in flow direction and velocity magnitude between the numerical and experimental results is evident, indicating that there is an additional factor that significantly reduces the thermal Marangoni forces, resulting in a thermal-buoyancy-driven flow.

4.2. Phase diagram of the axisymmetric flow direction

Hereinafter, the quasi-stationary model of § 3.2 is used for the following calculations. We present a comprehensive analysis of the flow direction inside a droplet of $\theta =150^\circ$ , in a ${Ma}_T$ – ${Ra}_T$ phase diagram (figure 4), assuming axial symmetry. To ensure the accuracy of our findings, we first examine the stability of the stationary solution across the entire parameter space. By calculating the eigenvalues of the linearised system of equations, we determine whether the obtained stationary solution is stable. From figure 4, one can identify a small bistable region where multiple stable solutions coexist for the same ${Ma}_T$ and ${Ra}_T$ values. For simplicity, we exclude this region from our analysis, as indicated by the grey area in the phase diagram. Then three remaining distinct regions can be identified in the phase space: the blue region, where buoyancy forces predominantly drive the flow, corresponding to ${Ma}_T$ and ${Ra}_T$ values such that $\psi ^+ \gt 90\,\%$ ; the orange region, where thermocapillary forces primarily drive the flow, bounded by the curve $\psi ^+ \lt 10\,\%$ ; and the yellow region, where two competing vortices are observed, defined by $10\,\% \leqslant \psi ^+ \leqslant 90\,\%$ .

Figure 4. The ${Ma}_T$ – ${Ra}_T$ phase diagram for flow direction for droplet of $\theta =150^\circ$ . The grey area represents the bistable region, where multiple stable solutions coexist for the same ${Ma}_T$ and ${Ra}_T$ values. The blue region is dominated by thermal buoyancy, while the orange region is dominated by thermocapillary forces. In the yellow region, one observes rolls driven by both thermal Marangoni flow and thermal buoyancy. The black dashed line represents the $\psi ^+=50\,\%$ curve, where the competing vortices occupy equal volumes.

The black dashed line in figure 4 represents the $\psi ^+ = 50\,\%$ curve, indicating where the competing vortices occupy equal volumes. Note the slope of 1 for this curve in this log-log phase space, which results from the linear combination of the effects of ${Ma}_T$ and ${Ra}_T$ at the onset of flow as ${Ma}_T \rightarrow 0$ and ${Ra}_T \rightarrow 0$ .

The phase diagram (figure 4) highlights a significant shift in flow direction with varying ${Ma}_T$ , evident in the thin transitional yellow region. This shift can be attributed to the presence of warmer fluid in the bulk. When a disturbance locally reduces the surface tension at the interface, a Marangoni flow transports fluid away from the perturbed region. Incompressibility requires that this fluid is replenished from the bulk. Warmer fluid, with corresponding lower liquid–gas surface tension, is pushed towards the perturbed region, further increasing the Marangoni flow. As a result, a self-enhancing thermal Marangoni flow is generated, leading to a rapid transition into a vortex from the contact line towards the apex as ${Ma}_T$ increases. The phase diagram also reveals that within the examined range of ${Ma}_T$ and ${Ra}_T$ , Marangoni forces predominantly drive the flow for sufficiently large ${Ma}_T$ (larger than $\sim 10$ ). Typically, ${Ra}_T$ is significantly smaller than ${Ma}_T$ , rendering the blue region of the phase diagram unrealistic for water droplets. Thereby, the phase diagram also suggests that an additional factor that mitigates the influence of thermocapillary forces under the given experimental conditions must be considered to replicate the experimental observations.

4.3. Axisymmetrybreaking

In the previous subsection, we assumed axial symmetry to determine the flow direction. Here, we investigate the stability of these stationary solutions under azimuthal symmetry. To do so, we use the method described in § 3.3 to estimate the critical parameters at which azimuthal stability is broken, and at which, simultaneously, the corresponding eigenfunction resembles a single-roll convection. We begin by setting ${Ma}_T=1$ and iteratively adjusting ${Ra}_T$ to find a value close to the critical ${Ra}_T^c$ for axial symmetry breaking at azimuthal wavenumber $m=1$ . This is done by calculating the solutions of the eigenproblem at each iteration until the largest eigenvalue approaches zero, through a shift–invert Arnoldi method (Arnoldi Reference Arnoldi1951; Saad Reference Saad2011). Then we activate the solver described by Diddens & Rocha (Reference Diddens and Rocha2024) to accurately find ${Ra}_T^c$ at ${Ma}_T=1$ . Once ${Ra}_T^c$ is determined, we obtain, through pseudo-arc-length continuation, the critical ${Ra}_T^c$ – ${Ma}_T$ curve, represented by the black solid line of figure 5.

Figure 5. Phase diagram of the axial symmetry-breaking phenomenon for a droplet of 150 $^\circ$ contact angle. The black solid line represents the critical ${Ra}_T^c$ – ${Ma}_T$ curve for which the flow loses stability of the axial symmetry solution (blue region) and can exhibit single-roll convection (purple region). The black dashed lines represent chosen experimental conditions, including those shown in figure 1, initially (circular marker) and after $90\,\%$ of volume is evaporated (tringular marker). The experimental conditions (see supplementary movies 1–5) are coloured according to the observed flow at the corresponding time instance, i.e. axisymmetric (blue) and single-roll convection (purple). Clearly, in three of the five cases presented here, the experimentally found single-roll convection at the beginning of the evaporation process disagrees with the theoretical expectation of this plot, reflecting that it is insufficient to consider only thermal effects as in § 4. In § 5, considering also solutal Marangoni forces due to surfactants, this discrepancy is resolved. The right- and left-hand side images show the flow fields for axial symmetry and for single-roll convection, respectively.

The bifurcation curve of figure 5 represents the critical parameters for which the flow loses the axial symmetry (blue region) and can exhibit single-roll convection (purple region). Note that the method used here (see § 3.3) does not provide the full nonlinear dynamics after symmetry breaking, which would require full three-dimensional simulations. Bearing that in mind, the representations of possible flow fields in each region shown on the right- and left-hand sides of the figure are merely illustrative. These flow representations depict axial symmetry and single-roll convection, respectively, and they are obtained by revolving the two-dimensional flow field around the $z$ -axis and assigning the total value of the perturbed velocity at each point, i.e. of $\boldsymbol{u}(r,\phi ,z, t) = \boldsymbol{u}^0(r,z) + \epsilon\, \boldsymbol{u}^m(r,z)\,{\textrm e}^{{\textrm i} m \phi + \lambda _m t}+\text{c.c.}$ , where $m=1$ , $r$ and $\phi$ are the radial and azimuthal directions, $\epsilon$ is an attributed small amplitude of the perturbation, $\boldsymbol{u}^0$ is the axial symmetric part of the velocity field, $\boldsymbol{u}^m$ is the eigenfunction of the azimuthal perturbation, and $\lambda _m$ is the corresponding eigenvalue. Some chosen experimental conditions, including the ones shown in figure 1, are depicted here for their respective initial conditions and after $90\,\%$ of the droplet volume has evaporated. We also performed azimuthal stability analysis for higher wavenumbers $m$ , but no ${Ra}_T^c$ values were found within the considered phase space.

With an increase in ${Ma}_T$ , there is a corresponding rise in the critical ${Ra}_T^c$ , implying that thermal Marangoni effects contribute to the stabilisation of the flow in axial symmetry. This phenomenon can be linked to enhanced mixing within the droplet, a result of thermocapillary-induced convection. Interestingly, the experimental conditions all fall within the axisymmetric region, suggesting that the flow should maintain axial symmetry throughout the evaporation process. Contrary to this expectation, as discussed in § 2, the experimentally observed flows do not exhibit axial symmetry at the early stages. This discrepancy indicates, once again, the presence of an additional factor causing symmetry breaking under the given experimental conditions, potentially due to a reduction in thermocapillary-induced mixing.

5.1. Flow reversal at experimental conditions

Using the same experimental initial conditions as in § 4.1, we investigate the influence of surfactants on the flow direction. For a slight average reduction of 0.1 $\,\%$ of the static surface tension $\sigma _0$ due to the action of surfactants, the flow follows the expected thermal Marangoni direction, from the contact line towards the apex; see figure 6(a). However, if the reduction of $\sigma _0$ is 0.32 $\,\%$ , then a competing vortex in the opposite direction emerges; see figure 6(b). This opposing vortex is much larger if the surface tension reduction is 0.34 $\,\%$ ; see figure 6(c). Finally, when the reduction reaches 0.5 $\,\%$ , the flow completely reverses, moving from the apex towards the contact line; see figure 6(d). The findings from the study by van Gaalen et al. (Reference van Gaalen, Wijshoff, Kuerten and Diddens2022), which do not account for buoyancy-driven convection, indicate that a reduction of approximately 0.8 % in surface tension results in a substantial decrease of 99 % in interface velocity. It is noteworthy that these results are consistent with our own findings, despite the omission of buoyancy effects in their paper.

Figure 6. Influence of surfactants on the flow direction. The values (a) $0.1\,\%$ , (b) $0.32\,\%$ , (c) $0.34\,\%$ and (d) $0.5\,\%$ are considered for the reduction of static surface tension $| \partial _\Gamma \sigma | / \sigma _0$ . The flow direction is from the contact line towards the apex in (a). A competing vortex in the opposite direction is observed in (b) due to the increasing influence of the surfactants. There is rapid growth of the thermal Rayleigh vortex in (c). the flow is completely in the thermal Rayleigh direction in (d), with a velocity magnitude comparable to the experimental findings.

The velocity magnitude in figure 6(d) is comparable to both the experimental results of figure 1(b) and the hypothetical pure thermal Rayleigh flow depicted in figure 1(b). These results demonstrate that even a tiny amount of contaminants can cause a significant change in the flow direction. It is also evident that for a larger reduction of $\sigma _0$ , the distribution of the concentration of surfactants at the interface is more uniform. This is expected, as surfactants effectively reduce surface tension, resulting in weaker advective forces associated with thermal Marangoni flow. Consequently, surfactants diffuse along the interface, rather than concentrating at the apex of the droplet.

5.2. Phase diagram of the axisymmetric flow direction

When surfactants are taken into account, the phase diagram undergoes significant changes; see figure 7. Using the same colour code as in figure 4, we consider the ${{Ma}}_\Gamma$ values 0, 100, 500 and 5000. Even a slight $\sigma _0$ reduction of approximately 0.01 $\,\%$ , corresponding to ${{Ma}}_\Gamma = 100$ at the conditions specified for table 1, leads to a substantial portion of the phase diagram being covered by the blue area, indicating flow in the apex towards the contact line. Remarkably, for ${{Ma}}_\Gamma = 5000$ , within the considered ranges of ${Ra}_T$ and ${Ma}_T$ , the flow is almost always in the same direction as if driven by thermal buoyancy, i.e. moving from the apex towards the contact line. This corresponds to a mere $\sigma _0$ reduction of approximately 0.5 $\,\%$ . Supplementary movie 6 illustrates the change in flow direction when travelling through the different regions in the phase diagram at ${{Ma}}_\Gamma =100$ .

The influence of the contact angle on the direction of the flow is discussed in Appendix C.

Figure 7. Flow direction phase diagram for a droplet with a contact angle $\theta = 150^\circ$ , considering the influence of surfactants at the interface, for (a) ${{Ma}}_\Gamma =0$ , (b) ${{Ma}}_\Gamma =100$ , (c) ${{Ma}}_\Gamma =500$ and (d) ${{Ma}}_\Gamma =5000$ . The blue area corresponds to ${{Ra}}_T$ -dominated flows, the orange to ${{Ma}}_T$ -dominated flows, and the yellow to a combination of both. The black dashed line corresponds to 50 $\,\%$ in volume of flow in each direction, i.e. $\psi ^+=50\,\%$ .

5.2.1. Change of slope in the $\psi ^+=50\,\%$ curve

Figure 7 displays a change in the slope of the $\psi ^+=50\,\%$ curve as ${Ma}_T$ increases in the presence of surfactants. Notably, the quasi-stationary solutions along this curve feature two competing vortices – one driven by buoyancy, and the other by Marangoni forces – that occupy equal volumes within the droplet. We identify two distinct regimes of the curve corresponding to different slopes, for low and high ${Ma}_T$ values.

To explore the impact of surfactants on this phenomenon, figure 8 examines a specific case from figure 7 with ${{Ma}}_\Gamma =100$ . The surfactants’ concentration along the droplet–gas interface is depicted for both the lower (bottom) and higher (top) regimes of the curve. Additionally, figure 8 illustrates the vortices within the droplet for each corresponding ${Ma}_T$ value and the tangential velocity $u_t$ at the droplet–gas interface, plotted against the normalised arc length coordinate $s$ , where $s=0$ is at the contact line, and $s=1$ is at the apex.

Figure 8. Surfactants’ concentration along the droplet–gas interface according to the normalised droplet–gas arc length $s$ , for ${{Ma}}_\Gamma =100$ , $\theta =150^\circ$ and two ${Ma}_T$ , one below the threshold for which the flow reversal curve changes slope (bottom), and one above (top).

For low ${Ma}_T$ , the slope of the curve is approximately 1. In this regime, the vortex from the contact line towards the apex is primarily located near the droplet–gas interface, as shown at the bottom of figure 8. Thermocapillary forces counteract the effect of surfactants, resulting in a weak thermal Marangoni flow close to the interface, as depicted in the respective $u_t$ – $s$ plot.

As ${Ma}_T$ increases, surfactants are advected and compressed towards the droplet apex. In the high ${Ma}_T$ regime (top part of figure 8), the surfactants’ concentration at the contact line reaches zero, as shown in the respective $\tilde {\Gamma }$ – $s$ plot for small $s$ . Consequently, with no counteracting surfactant effect near the contact line, a strong thermal Marangoni flow emerges in that region, as indicated by the large $u_t$ in the top $u_t$ – $s$ plot for small $s$ . This flow is strong enough to reverse the flow direction not only at the interface but also near the axis, as illustrated by the vortex representations.

The change in the slope of the $\psi ^+=50\,\%$ curve is caused directly by the surfactant concentration at the contact line dropping to zero at high ${Ma}_T$ values. In the absence of surfactants near the contact line, a strong thermal Marangoni flow develops, easily overpowering the thermal Rayleigh forces. Supplementary movie 7 provides additional visual support for this explanation.

5.3. Axisymmetrybreaking

When surfactants are present ( ${{Ma}}_\Gamma \gt 0$ ), the azimuthal symmetry of the flow is broken over a wider range of ${Ra}_T$ and ${Ma}_T$ than for the ${{Ma}}_\Gamma =0$ case. Figure 9 presents the results of the azimuthal stability analysis for various values of ${{Ma}}_\Gamma$ .

Figure 9. Stability analysis at azimuthal wavenumber $m=1$ for (a) ${{Ma}}_\Gamma=0$ , (b) ${{Ma}}_\Gamma=100$ , (c) ${{Ma}}_\Gamma=500$ , (d) ${{Ma}}_\Gamma=1000$ , (e) ${{Ma}}_\Gamma=3000$ and (f) ${{Ma}}_\Gamma=5000$ of a droplet of 150 $^\circ$ contact angle in a ${{Ra}}_T$ – ${{Ma}}_T$ phase diagram. The black solid line represents the critical ${{Ra}}_T^c$ – ${{Ma}}_T$ curve for which the flow goes from axial symmetry (blue bottom part) to single-roll convection (purple top part). The black dashed lines represent chosen experimental conditions (see supplementary movies 1–5), including those shown in figure 1.

Interestingly, our calculations show that for certain values of ${{Ma}}_\Gamma$ , the flow can initially exhibit asymmetric single-roll convection, and as evaporation progresses (or as the droplet volume decreases), the flow can eventually stabilise and regain axial symmetry. In figures 9(c–f), we show one single experiment that displays this behaviour while presenting quantitative agreement between experiments and numerical simulations. Nevertheless, the numerical simulations do not accurately capture the critical parameters for the transition from single-roll convection to axial symmetry across all experimental conditions. This discrepancy may be due to simplifications in the numerical model, the absence of detailed information about the surfactants’ properties, or sporadic experimental events such as contact line pinning. Although our simplified model does not allow for precise quantitative comparison with the experiments, it does provide qualitative insight into the influence of contaminants on the stability of the axisymmetric quasi-stationary solutions. For instance, when considering ${{Ma}}_\Gamma =100$ and sufficiently large ${Ma}_T$ , as ${Ma}_T$ increases, the critical ${{Ra}}_T^c$ also rises, indicating that thermal Marangoni induced circulation stabilises the flow and restores axial symmetry. For ${Ma}_\Gamma \gtrapprox 500$ , the bifurcation curve becomes independent of ${Ma}_T$ in the limit of low ${Ma}_T$ . However, for larger ${Ma}_T$ , the critical ${Ra}_T^c$ decreases with increasing ${Ma}_T$ , suggesting that the combined effects of solutal and thermal Marangoni forces destabilise the flow and break axial symmetry.

In the next subsection, we further explore the influence of the contact angle on the azimuthal symmetry-breaking phenomenon. To simplify the problem, we isolate the thermal Rayleigh effects ( ${Ma}_T={Ma}_\Gamma =0$ ) and focus on examining how the large contact angle affects azimuthal symmetry breaking.

5.4. Contact angle influence in axisymmetrybreaking: pure Rayleigh flow

We present ${Ra}_T^c$ for a range of $\theta$ , from $85^\circ$ to $179^\circ$ , in figure 10, considering ${Ma}_T={Ma}_\Gamma =0$ . At the left of the figure, the convective roll flow field at $\theta =87^\circ$ is shown, and at the right, the flow field at $\theta =179^\circ$ . Obviously – and as the flow field representations reflects – for droplets of fixed volume, the contact radius will be much larger if the contact angle is smaller. The contact area with the substrate, which has fixed temperature and is perfectly conducting, is much smaller for larger $\theta$ , which will result in a larger temperature gradient within the droplet. As $\theta \rightarrow 180^\circ$ , the contact radius reduces to zero, approaching a heat point source. At the same time, an increasing contact angle reduces the size of the no-slip boundary condition at the substrate, thereby allowing for single-roll convection, which is almost free from viscous forces in the limit $\theta \rightarrow 180^\circ$ .

Figure 10. Critical ${Ra}_T^c$ for different contact angles $\theta$ for the $m=1$ instability to happen, for ${Ma}_\Gamma ={Ma}_T=0$ . The black solid line represents the ${Ra}_T^c$ – ${Ma}_T$ curve for which the flow goes from axial symmetry (blue bottom part) to single-roll convection (purple top part).

As the contact angle increases, ${Ra}_T^c$ decreases, indicating that the large contact angle contributes to the azimuthal instability of the flow, similar to the onset of Rayleigh–Bénard convection in a cylindrical container with aspect ratio 1. To further support this hypothesis, we considered a theoretical scenario with no thermal buoyancy effects ( ${Ra}_T={Ma}_\Gamma =0$ ), only thermal Marangoni effects. The resulting solution of the azimuthal eigenvalue problem shows that flow is always stable in the azimuthal direction, independently of ${Ma}_T$ . Surfactants thus play a key role in the complex large ${Ma}_T$ azimuthal symmetry-breaking phenomenon depicted in figures 9(c–f), as concluded from the results of this section.

For shallow droplets (small contact angle), thermal Marangoni instabilities are expected to induce highly non-axisymmetric flow patterns (Babor & Kuhlmann Reference Babor and Kuhlmann2023), also known as hydrothermal waves (Sefiane et al. Reference Sefiane, Moffat, Matar and Craster2008), which usually appear on heated substrates. Solutal Marangoni flow can exhibit even more asymmetric and even chaotic flow, as e.g. observed in evaporating ethanol–water droplets (Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b ). Since these phenomena are not observed in water droplets at ambient conditions, they are beyond the scope of this paper.