2.1 Droplet composition

In total, we have experimentally investigated the evaporation of three different types of droplets, namely a pure water droplet, a binary water–ethanol droplet and a ternary Ouzo droplet. The liquids used for these droplets were pure Milli-Q water [produced by a Reference A+ system (Merck Millipore) at $18.2~\text{M}\unicode[STIX]{x03A9}~\text{cm}$ (at $25~^{\circ }\text{C}$ )], a mixture of 37.88 % (wt/wt) Milli-Q water and 62.12 % (wt/wt) ethanol [Sigma-Aldrich; ${\geqslant}$ 99.8 %], and a mixture of 37.24 % (wt/wt) Milli-Q water, 61.06 % (wt/wt) ethanol and 1.70 % anise oil [Sigma-Aldrich], respectively. Furthermore, for the experimental micro-PIV measurements described later on in § 2.3, tracer particles have been added to the binary water–ethanol droplet.

2.2 Set-up and image analysis

The droplets were deposited through a teflonized needle [Hamilton; 8646-02] by a motorized syringe pump [Harvard; PHD 2000] on a flat hydrophobic octadecyltrichlorosilane glass surface with a thickness of $d_{s}=0.17~\text{mm}$ . The advancing and receding contact angles of water on the surface are $112^{\circ }$ and $98^{\circ }$ , respectively. The entire evaporation process of the droplets was recorded by a CCD camera [Ximea; MD061MU-SY, 3 f.p.s. at $1372\times 1100$ pixel resolution] equipped with a long-distance microscope system [Infinity; Model K2 DistaMax] for side view recordings and a digital SLR camera [Nikon D750] equipped with a CMOS sensor [24 f.p.s. at $1920\times 1080$ pixels resolution] attached to a high-magnification zoom lens system [Thorlabs; MVL12X3Z] for top view recordings. The side view illumination was provided by a homemade collimated LED source. A universal hand-held probe [Omega; HH-USD-RP1, accuracy relative humidity is $\pm 2\,\%$ over 10 to 90 % @ $25~^{\circ }\text{C}$ and a temperature accuracy of $\pm 0.3^{\circ }~\text{K}$ @ $25~^{\circ }\text{C}$ ] was used to measure the relative humidity and temperature in the laboratory at a sampling rate of one per second. A similar sketch of the set-up is described in detail in Tan et al. (Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016).

The image analysis was performed by a custom-made MATLAB program, through which all of the geometric parameters at every frame were successfully determined. For the millimetric droplets in this study, the characteristic lengths of the droplets are smaller than the capillary length scale, which is equal to 2.7 mm for water and 1.7 mm for ethanol. Thus, gravity effects influencing the shape of the droplets can be disregarded. For pure water droplets and binary water–ethanol droplets, a spherical-cap shape assumption was used over the whole evaporation process. The program fits a circle to the contour of the droplet silhouette in side view. The contact angle $\unicode[STIX]{x1D703}$ was calculated based on the position of the intersection of the base line and the fitted circles. Since Ouzo droplets temporarily show a very characteristic deviation from the spherical-cap shape due to the appearance of an oil ring at the rim, only the top part of the contour above the oil ring was fitted by a circle (cf. Tan et al. Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016). A polynomial fit was applied to the profile of the oil ring and gave the contact angle $\unicode[STIX]{x1D703}$ . The contact angle $\unicode[STIX]{x1D703}^{\ast }$ of the spherical-cap-shaped top part was calculated based on the position of the intersection point where the fitted line of the oil ring profile and the fitted circle of the top contour cross. The droplet volume $V$ was calculated by integrating the volumes of horizontal disk layers and assuming rotational symmetry of each layer with respect to the vertical axis.

2.3 Micro-PIV

To determine the velocity in the binary water–ethanol droplet, micro-PIV measurements were performed. To that end, the binary water–ethanol mixture was seeded with fluorescent particles. Per 1 mL of the water–ethanol solution, $30~\unicode[STIX]{x03BC}\text{L}$ of aqueous tracer particle solution [microParticles GmbH; PS-FluoRed-5.0: Ex/Em 530 nm/607 nm] was used. The diameter of these tracer particles is $d_{p}=5~\unicode[STIX]{x03BC}\text{m}$ .

Figure 1. Experimental set-up of the micro-PIV system. The dual-cavity laser and the high-speed camera are synchronized to record consecutive image pairs. To apply fluorescence techniques, the inverted microscope is equipped with a dichroic cube and the droplet is seeded with fluorescent particles. The focal plane of the microscope is placed at the bottom of the droplet, i.e. close and parallel to the substrate.

An overview of the experimental micro-PIV set-up is given in figure 1. It consists of a dual-cavity Nd:YAG laser [Litron; NANO S 65-15PIV], a high-speed camera [Photron; Fastcam SA-X2 64 GB] and an inverted microscope [Olympus; GX-51] which is equipped with a dichroic cube. The laser and the high-speed camera are synchronized by a pulse/delay generator [BNC; Model-575]. The laser beam illuminates the droplet from below through a 10 $\times$ magnification objective lens with a numerical aperture of 0.30 (zoom in of figure 1). The objective focal plane was placed just above the surface at approximately $13~\unicode[STIX]{x03BC}\text{m}$ to match the thickness of the focal plane, or the depth of field, which was calculated to be $13.4~\unicode[STIX]{x03BC}\text{m}$ . Thus, we had a sharp image of the tracer particles closest to the substrate.

In this way, 10 consecutive image pairs with an interframing time of 4 ms were taken per second. To calculate the velocity field, the obtained images were first post-processed with a custom MATLAB code to enhance the contrast. Then, the image pairs were analysed with PIVlab (Thielicke Reference Thielicke2014; Thielicke & Stamhuis Reference Thielicke and Stamhuis2014), using an interrogation window of 64 $\times$ 64 pixels, corresponding to $128~\unicode[STIX]{x03BC}\text{m}\times 128~\unicode[STIX]{x03BC}\text{m}$ . An interrogation window overlap of 75 % leads to a $32~\unicode[STIX]{x03BC}\text{m}$ vector spacing.

To qualify the degree to which the tracer particles exactly follow the flow, we analysed the Stokes number and the ratio of Stokes number to a buoyancy-corrected Froude number (Mathai et al. Reference Mathai, Calzavarini, Brons, Sun and Lohse2016). The Stokes number is defined as $St\equiv t_{0}u_{max}/r_{c}$ , where $t_{0}\equiv (\unicode[STIX]{x1D70C}_{p}d_{p}^{2}/18\unicode[STIX]{x1D707}_{f})(1+\unicode[STIX]{x1D70C}_{f}/(2\unicode[STIX]{x1D70C}_{p}))$ is the characteristic time, $u_{max}$ the maximum fluid velocity ( ${\sim}10~\text{mm}~\text{s}^{-1}$ ), $\unicode[STIX]{x1D70C}_{p}$ the density of the tracer particles ( $1.04~\text{g}~\text{cm}^{-3}$ ) and $\unicode[STIX]{x1D707}_{f}$ the dynamic viscosity of the droplet liquid ( ${\sim}10^{-3}~\text{Pa}~\text{s}$ ). The factor $(1+\unicode[STIX]{x1D70C}_{f}/(2\unicode[STIX]{x1D70C}_{p}))$ accounts for the added mass force (Oliveira, van der Geld & Kuerten Reference Oliveira, van der Geld and Kuerten2015). The buoyancy-corrected Froude number is defined as $Fr\equiv u_{max}^{2}/r_{c}/(1-\unicode[STIX]{x1D70C}_{f}/\unicode[STIX]{x1D70C}_{p})g$ , where $\unicode[STIX]{x1D70C}_{f}$ is the density of the droplet liquid ( $0.977~\text{g}~\text{cm}^{-3}$ for water and $0.789~\text{g}~\text{cm}^{-3}$ for ethanol) and $g$ is the gravitational acceleration. The calculated values for $St$ and $St/Fr$ are ${\sim}10^{-5}$ and ${\sim}10^{-2}$ respectively, which indicate that the tracer particles are truly tracers for the liquid flow inside the evaporating droplets.

The side view recording set-up (§ 2.2) was coupled to the micro-PIV set-up to monitor the geometric shape evolution of the evaporating droplet synchronized with the flow inside.

2.4 Confocal microscopy

Unfortunately, the flow in ternary Ouzo droplets cannot be quantified using the micro-PIV technique, since the fluorescence signal of the tracer particles is concealed by the nucleation of oil microdroplets in its second phase (Tan et al. Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016). The nucleated microdroplets themselves cannot be used as tracer particles because of their growing size and the non-uniform dispersion in the Ouzo droplet. However, a qualitative investigation of the flow in an Ouzo droplet can be performed by visualizing the phase separation and the movement of the oil microdroplets with confocal microscopy. To that end, real-time observations of Ouzo droplets were carried out by using an inverted Nikon A1 confocal laser scanning microscope system (Nikon Corporation, Tokyo, Japan) with a 20 $\times$ magnification dry objective (CFI Plan Apochromat VC 20 $\times /0.75$ DIC, numerical aperture $=0.75$ , working distance $=1.0~\text{mm}$ ). An Ouzo droplet consisting of trans-anethole, ethanol and water was deposited on a hydrophobic glass coverslip with a thickness of $170~\unicode[STIX]{x03BC}\text{m}$ . The initial volume of the droplets varied from $0.1~\unicode[STIX]{x03BC}\text{L}$ to $1.0~\unicode[STIX]{x03BC}\text{L}$ . Trans-anethole oil was labelled by Nile Red (Sigma-Aldrich) which was excited by laser light with a wavelength of 561 nm, while the water–ethanol mixture was labelled by Fluorescein 5(6)-isothiocyanate (FITC, Sigma-Aldrich) which was excited by laser light with a wavelength of 488 nm. The scan started immediately after the droplet was deposited on the coverslip with a frame rate of 30 Hz for scanning two-dimensional images under resonant mode.

3.1 Axisymmetric finite element method

Figure 2. (a) The problem is expressed by axisymmetric cylindrical coordinates. The individual domains $\unicode[STIX]{x1D6FA}_{g}$ , $\unicode[STIX]{x1D6FA}_{l}$ , $\unicode[STIX]{x1D6FA}_{s}$ and $\unicode[STIX]{x1D6FA}_{b}$ are the gas phase, the liquid droplet, the substrate with finite thickness $d_{s}$ and the air below the substrate, respectively. (b) The model solves the coupled processes of vapour-diffusion-limited mixture evaporation, multicomponent flow with composition-dependent quantities and driven by solutal and thermal Marangoni flow, as well as the temperature field.

In the following an outline of the axisymmetric finite element method is given. A detailed description of the model can be found in Diddens (Reference Diddens2017). As depicted in figure 2, the space, represented by axisymmetric cylindrical coordinates $(r,z)$ , is separated into individual subdomains, namely the gas phase $\unicode[STIX]{x1D6FA}_{g}$ , the droplet liquid $\unicode[STIX]{x1D6FA}_{l}$ , the substrate $\unicode[STIX]{x1D6FA}_{s}$ with a finite thickness $d_{s}$ and the air below the substrate $\unicode[STIX]{x1D6FA}_{b}$ . The droplet shape, defined by the liquid–gas interface $\unicode[STIX]{x1D6E4}_{lg}$ , is assumed to be always in a spherical-cap shape. According to the ellipticity measurements of Tan et al. (Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016), this is a good approximation, at least as long as the oil ring in case of the Ouzo droplet is not taken into account.

The following derivation of the model considers the Ouzo droplet, but it can easily be simplified to the pure water droplet and the binary water–ethanol droplet. The composition dependence of the physical fluid properties, i.e. viscosity $\unicode[STIX]{x1D707}$ , diffusivity $D_{}$ , mass density $\unicode[STIX]{x1D70C}$ , surface tension $\unicode[STIX]{x1D70E}$ and thermodynamic activities, are locally considered by fitting experimental data for water–ethanol mixtures (Vazquez, Alvarez & Navaza Reference Vazquez, Alvarez and Navaza1995; González et al. Reference González, Calvar, Gómez and Domínguez2007; Pařez et al. Reference Pařez, Guevara-Carrion, Hasse and Vrabec2013) and using thermodynamic models (Zuend et al. Reference Zuend, Marcolli, Luo and Peter2008, Reference Zuend, Marcolli, Booth, Lienhard, Soonsin, Krieger, Topping, McFiggans, Peter and Seinfeld2011). Due to the small initial concentration of anise oil in the Ouzo mixture, its influence on the liquid properties is not considered. The temperature dependence of the surface tension is also considered based on the experimental data of Vazquez et al. (Reference Vazquez, Alvarez and Navaza1995) to account for thermal Marangoni flow. Plots of all these relations can be found in Tan et al. (Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016) and Diddens et al. (Reference Diddens, Kuerten, van der Geld and Wijshoff2017).

In the gas phase, the diffusion equation for the vapour concentration $c_{\unicode[STIX]{x1D6FC}}$ with $\unicode[STIX]{x1D6FC}=e$ and $\unicode[STIX]{x1D6FC}=w$ for ethanol and water, respectively, has to be solved. With the vapour diffusivity $D_{\unicode[STIX]{x1D6FC}}^{g}$ of $\unicode[STIX]{x1D6FC}$ in air, this reads

By time scale arguments, one can show that (3.1) can be approximated by its steady-state limit, i.e. by the Laplace equation $\unicode[STIX]{x1D6FB}^{2}c_{\unicode[STIX]{x1D6FC}}=0$ . This has been validated by Deegan et al. (Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000) and Hu & Larson (Reference Hu and Larson2002) for pure droplets and by Diddens et al. (Reference Diddens, Kuerten, van der Geld and Wijshoff2017) for multicomponent droplets. However, since the consideration of the time derivative has no noticeable influence on the calculation time, the full equation (3.1) is used in the model. The anise oil is assumed to be non-volatile due to the low vapour pressure of trans-anethole. As boundary condition of (3.1) for $(r,z)\rightarrow \infty$ , the ambient vapour concentrations $c_{\unicode[STIX]{x1D6FC},\infty }$ have to be imposed. While there is no ethanol vapour present far away from the droplet, i.e. $c_{e,\infty }=0$ , the water vapour concentration can be expressed by the ideal gas law:

Here, $\unicode[STIX]{x1D719}$ is the relative humidity, $M_{w}$ is the molar mass of water, $R$ is the ideal gas constant, $T_{\infty }$ is the room temperature and $p_{w,sat}$ is the saturation pressure, which temperature dependence is given by the Antoine equation. At the liquid–gas interface $\unicode[STIX]{x1D6E4}_{lg}$ , the vapour–liquid equilibrium concentration $c_{\unicode[STIX]{x1D6FC},VLE}$ is imposed, which can be calculated via Raoult’s law, i.e. by

The vapour–liquid equilibrium couples the vapour concentration at the droplet interface to the liquid mole fraction $x_{\unicode[STIX]{x1D6FC}}$ , where non-idealities are considered by the activity coefficient $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D6FC}}$ . The activity coefficients were calculated by the thermodynamic model AIOMFAC (Zuend et al. Reference Zuend, Marcolli, Luo and Peter2008, Reference Zuend, Marcolli, Booth, Lienhard, Soonsin, Krieger, Topping, McFiggans, Peter and Seinfeld2011). Plots of these coefficients can be found in Tan et al. (Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016) and Diddens et al. (Reference Diddens, Kuerten, van der Geld and Wijshoff2017). Since (3.3) has to be evaluated at the local temperature, the evaporation rates are also coupled to the temperature field $T$ . From the solution of (3.1), one can directly calculate the diffusive vapour flux at the interface:

Here, $\boldsymbol{n}_{lg}$ is the interface normal vector pointing to the gas phase. The mass transfer rates $j_{w}^{lg}$ and $j_{e}^{lg}$ can be determined from the diffusive vapour fluxes via the coupled mass transfer jump conditions

where $y_{\unicode[STIX]{x1D6FC}}^{g}$ is the mass fraction of $\unicode[STIX]{x1D6FC}$ -vapour in the gas phase. By assuming constant gas density $\unicode[STIX]{x1D70C}^{g}$ , one can approximate $y_{\unicode[STIX]{x1D6FC}}^{g}=c_{\unicode[STIX]{x1D6FC},VLE}/\unicode[STIX]{x1D70C}^{g}$ .

In the droplet, the spatio-temporal liquid composition is governed by the convection–diffusion equations for the liquid mass fractions $y_{\unicode[STIX]{x1D6FC}}$ , i.e.

For simplicity, the composition dependence of the mass density $\unicode[STIX]{x1D70C}$ is approximated by the spatially averaged composition. This assumption has been validated in Diddens (Reference Diddens2017) and simplifies (3.6) as well as the momentum equation. At the liquid–gas interface, a source/sink term is imposed via the interface delta function $\unicode[STIX]{x1D6FF}_{lg}$ . The corresponding flux $\boldsymbol{J}_{\unicode[STIX]{x1D6FC}}$ is obtained by the counterpart of (3.5) for the liquid phase, i.e.

The mass fraction of the non-volatile anise oil is not solved explicitly, but calculated via $y_{a}=1-y_{e}-y_{w}$ .

Due to the latent heat of evaporation $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D6FC}}$ , the interface is cooled down by the evaporative mass fluxes $j_{\unicode[STIX]{x1D6FC}}^{lg}$ . This effect is taken into account by considering the convection–diffusion equation for the temperature $T$ , i.e.

with $T=T_{\infty }$ for $(r,z)\rightarrow \infty$ . The mass density $\unicode[STIX]{x1D70C}$ , the specific heat capacity $c_{p}$ and the thermal conductivity $\unicode[STIX]{x1D706}$ are different on the individual subdomains. In the droplet, the values furthermore depend on the local mixture composition according to the experimental data of Grolier & Wilhelm (Reference Grolier and Wilhelm1981) and Yano, Fukuda & Hashi (Reference Yano, Fukuda and Hashi1988). While (3.8) can usually be treated in its quasi-steady limit for single-component droplets (e.g. Dunn et al. Reference Dunn, Wilson, Duffy, David and Sefiane2009), this argument does not necessarily hold in the case of multicomponent droplets, where a rapid and chaotic solutally driven Marangoni flow may be dominant, so that the consideration of the full equation (3.8) is required.

Finally, the velocity $\boldsymbol{u}$ has to be solved. For the moment, the flow in the gas phase is disregarded. The influence of convective mass and thermal transport in the gas phase is discussed later in § 3.4. Due to the micrometre-sized droplet, the Bond number $Bo$ and the capillary number $Ca$ are small, i.e. $\ll 1$ ( $Bo<0.08$ and $Ca<3\times 10^{-4}$ holds for all droplets discussed in this article). While the former allows to disregard gravitational effects, the latter ensures that the flow in the droplet is dictated by the surface tension. Thus, it is assumed that the droplet is always in an axisymmetric spherical-cap shape. Thereby, it is possible to calculate the normal interface velocity $\boldsymbol{u}_{lg}\boldsymbol{\cdot }\boldsymbol{n}_{lg}$ along the entire interface directly from the evolution of the volume $V(t)$ and the contact angle $\unicode[STIX]{x1D703}(t)$ . The volume loss $\dot{V}(t)$ can easily be calculated by integrating the evaporation rates over the liquid–gas interface and by taking the spatially averaged composition dependence of the mass density into account.

For the contact angle $\unicode[STIX]{x1D703}$ , a model could be used which predicts the contact angle based on the liquid composition at the contact line (cf. e.g. Diddens Reference Diddens2017). However, since the experimentally observed contact angle evolutions are usually more complex, we directly impose the experimentally determined contact angle evolution into the corresponding numerical simulation. To that end, we have fitted the experimental data for $\unicode[STIX]{x1D703}(V)$ . Expressing $\unicode[STIX]{x1D703}$ as a function of the volume $V$ instead of as a function of time $t$ for this purpose has the advantage that the simulated droplet follows the same geometrical evolution as in the experiment, only possibly faster or slower if the volume evolution $V(t)$ is not exactly reproduced by the model.

The normal interface velocity defines via the kinematic boundary condition the normal flow velocity $\boldsymbol{n}_{lg}\boldsymbol{\cdot }\boldsymbol{u}$ in the droplet, i.e.

With the stress tensor

the tangential velocity component is subject to the shear stress boundary condition

Due to the small viscosity ratio between gas and liquid phase, i.e. $\unicode[STIX]{x1D707}^{g}\ll \unicode[STIX]{x1D707}$ , and the continuous tangential velocity at the interface, the contribution of the shear stress in the gas phase, i.e. the term $\boldsymbol{n}_{lg}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}^{g}\boldsymbol{\cdot }\boldsymbol{t}_{lg}$ , can be neglected. This assumption is validated later in § 3.4. Since the surface tension $\unicode[STIX]{x1D70E}$ is a function of the composition and the temperature, solutally and thermally driven Marangoni flow has to be expected.

Inside the droplet, a slightly modified Stokes flow is solved:

along with the no-slip boundary condition $\boldsymbol{u}=0$ at the substrate $z=0$ .

For the implementation of the finite element method, the governing equations are converted to corresponding weak forms (Diddens Reference Diddens2017) and solved with the help of the finite element package FEniCS (Logg Reference Logg2012).

The time stepping of the implicit integration scheme is dynamically adjusted delimited by a maximum Courant–Friedrichs–Lewy (CFL) number of $1/2$ . The dynamic mesh resolves the initial droplet domain with approximately 30 000 triangles, whereas the initial gas domain and the region below $z=0$ , i.e. the substrate and the air below the substrate, consist of approximately 10 000 triangles each. The circumradius of the total mesh is approximately eight times the initial contact line radius. The choice of these discretization parameters is validated in Diddens (Reference Diddens2017).

Figure 3. (a,b) Snapshots of a pure water droplet with velocity field in the droplet and vapour concentration in the gas phase (left) and temperature field (right) at two different times $t=100~\text{s}$ (a) and $t=1000~\text{s}$ (b). The white lines on the right half represent the isotherms, which are furthermore divided by equidistant sublevels inside the droplet (dashed lines). A corresponding movie is available as supplementary material (movie 1 available at https://doi.org/10.1017/jfm.2017.312). (c) Volume evolution of the experiment and according to simulations with and without considering thermal effects. (d) Contact angle of the experiment and corresponding fit, which is used as input parameter for the simulation. (e) Averaged tangential fluid velocity at the interface.

3.2 Pure water droplet

To validate the model on the basis of the simplest case, the experimental set-up was used to measure the evolution of a $0.63~\unicode[STIX]{x03BC}\text{L}$ pure water droplet at ambient temperature $T_{\infty }=20.4~^{\circ }\text{C}$ and relative humidity $\unicode[STIX]{x1D719}=54\,\%$ . The case of pure water can easily be simulated by setting $y_{w}=1$ in the model equations. The experimental data are depicted along with the corresponding simulation data in figure 3. The right side of panels (a) and (b) show the temperature field of the droplet and its surrounding. Due to the thin substrate, the latent heat of evaporation leads to an intense cooling, which can be more than $3~\text{K}$ . From the volume evolution $V(t)$ depicted in (c), it is apparent how the cooling influences the evaporation rate: If the latent heat and the corresponding reduction of the vapour pressure are considered, the numerical simulation (black solid line) perfectly agrees with the experimental data (red crosses). For comparison, also the corresponding isothermal simulation, i.e. without considering thermal effects, is indicated (blue dashed line). This simulation, which resembles the result of the model equation by Popov (Reference Popov2005), predicts a noticeably faster drying than the experiment. The role of evaporative cooling and its influence on the evaporation rate has already been investigated in several other studies (Girard et al. Reference Girard, Antoni, Faure and Steinchen2006; Girard, Antoni & Sefiane Reference Girard, Antoni and Sefiane2008; Dunn et al. Reference Dunn, Wilson, Duffy, David and Sefiane2009; Sefiane & Bennacer Reference Sefiane and Bennacer2011). Here, we extend this knowledge by pointing out the importance of the substrate thickness: although the used quartz glass is a rather good conductor ( $\unicode[STIX]{x1D706}\approx 1.4~\text{W}~(\text{mK})^{-1}$ ), thermal effects cannot be disregarded since the substrate is thin, i.e. smaller than the initial droplet size, and air is below.

The contact angle evolution in figure 3(d) shows the typical stick-slip behaviour, i.e. an initial decrease of $\unicode[STIX]{x1D703}(t)$ at constant contact line radius $r_{c}$ , followed by a constant contact angle with a receding contact line $r_{c}(t)$ . While this behaviour could be incorporated into the simulation by introducing a receding contact angle, the contact line dynamics of mixture droplets is more rich and complicated to account for with an accurate model. To that end, we have simply fitted the contact angle measured experimentally and used it as an input parameter throughout the simulation.

As can be seen from the left side of figure 3(a,b), the simulation predicts a strong thermal Marangoni flow from the contact line along the liquid–gas interface towards the apex. The corresponding temporal evolution can be inferred from figure 3(e), where the average tangential fluid velocity at the interface

is plotted. The tangent $\boldsymbol{t}_{lg}$ is defined to point along the interface towards the apex. One can infer that the simulation predicts a persistent thermal Marangoni flow throughout the entire lifetime of the droplet. In experiments, the observed thermal Marangoni flow in water droplets is usually much slower than the theoretical predictions (Hu & Larson Reference Hu and Larson2005). We discuss this issue later in § 4.2.

3.3 Binary water–ethanol droplet

As a next step, the model is validated by a binary water–ethanol droplet at ambient temperature $T_{\infty }=21.4~^{\circ }\text{C}$ and relative humidity $\unicode[STIX]{x1D719}=55\,\%$ . Due to practical reasons in the experiment, it took some seconds between the deposition of the droplet and the first snapshot at $t=0~\text{s}$ . Since the average composition changes within this offset time, the initial composition in the simulation has a slightly lower ethanol content (57.7 wt %) than specified in § 2.1. This particular correction was determined by extrapolating the numerically obtained initial compositional rate of change over the offset time.

Figure 4. (a–c) Snapshots of a binary water–ethanol droplet with ethanol mass fraction $y_{e}$ in the liquid phase and water vapour concentration in the gas phase (left) and temperature field (right) at three different times $t=20~\text{s}$ (a), $t=150~\text{s}$ (b) and $t=250~\text{s}$ (c). The arrows at the interface indicate the evaporation rates of water (left) and ethanol (right). A corresponding movie is available as supplementary material (movie 2). (d) Volume evolution according to experiment, full simulation and isothermal simulation. (e) Partial masses of both components in the drop. (f) Experimentally obtained contact angle and corresponding fit for the model. (g) Mean tangential fluid velocity at the interface and its root mean square average. (h) Average interface temperature.

Representative snapshots of the simulation are shown in figure 4(a–c). Initially, the droplet is rather flat and both components evaporate with the typical singularity near the contact line. In combination with the fact that water (evaporation rate by arrows on the left side) has a lower volatility than ethanol (arrows on the right side), this leads to an enhanced water concentration near the contact line. The resulting surface tension gradient drives an intense solutal Marangoni flow, which breaks up into chaos due the Marangoni instability of water–ethanol mixtures (Machrafi et al. Reference Machrafi, Rednikov, Colinet and Dauby2010). Since the chaotic bulk flow leads to a rapidly changing mixture composition and temperature at the interface, also the local evaporation rate is subject to chaotic fluctuations. This effect can be seen in movie 2 in the supplementary material.

In the initial regime, thermal Marangoni flow is almost irrelevant, but the incorporation of the latent heat of evaporation is still important, as it affects again the volume evolution, which is depicted in figure 4(d). As for the case of pure water, the volume $V(t)$ of the simulation agrees perfectly with the experimental data, but only if the evaporative cooling is considered. In the case of a binary mixture, the evaporative cooling has another effect, namely the suppression of the water evaporation at initial stages: Due to the highly volatile ethanol, the interface is cooled down to a temperature, at which the vapour–liquid equilibrium for the water vapour concentration at the interface almost reaches the level of the vapour concentration in the ambient. At a higher relative humidity, this coupling between the evaporation rates and the temperature would even induce an initial condensation of water (Diddens Reference Diddens2017).

The two distinct slopes in the volume curve $V(t)$ , which we already observed and explained in Tan et al. (Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016), result from the preferential evaporation of ethanol in the initial regime followed by the slower evaporation of the remaining water residual at later times. This interpretation is validated by resolving the evolution into the partial masses of water and ethanol as shown in figure 4(e). It is apparent that the evaporation rate of water, i.e. the slope of the partial mass ${\dot{m}}_{w}(t)$ , is slightly reduced as long as ethanol is present. This is due to Raoult’s law (3.3), but is also caused by the reduction of the vapour pressure due to the additional cooling stemming from the evaporation of ethanol.

The contact angle was again fitted from experimental data and used as input parameter for the simulation. The data for $\unicode[STIX]{x1D703}(t)$ depicted in figure 4(f) reveal an initial decrease due to a pinned contact line followed by an intense increase, which can be understood by considering the composition dependence of the surface tension in the Young–Laplace equation. The non-monotonic contact angle evolution qualitatively agrees with other reported data on evaporating water–ethanol droplets (Sefiane et al. Reference Sefiane, Tadrist and Douglas2003; Cheng et al. Reference Cheng, Soolaman and Yu2006; Wang et al. Reference Wang, Peng, Mujumdar, Su and Lee2008; Shi et al. Reference Shi, Shen, Zhang, Lin and Jiang2009).

At later times, when virtually all ethanol has evaporated, a flow transition occurs (cf. figure 4 c): the chaotic solutal Marangoni flow suddenly switches over to regular thermal Marangoni flow as for the pure water droplet of figure 3. The abrupt transition can also be inferred from the plot of the average tangential velocity $\bar{u}_{t}(t)$ in figure 4(g). While the chaotic solutal flow induces a net flow towards the contact line, $\bar{u}_{t}<0$ , it suddenly converges to the regular thermal Marangoni flow towards the apex, i.e. $\bar{u}_{t}>0$ . Additionally, to estimate the average flow speed at the interface without respecting the local direction, the root mean square averaged tangential fluid velocity

is plotted in figure 4(g). It is apparent that the flow is most intense in the limit of rather dilute ethanol. In this region, the compositional surface tension gradient, i.e. $\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}/\unicode[STIX]{x2202}y_{e}$ , has its maximum and hence drives the most intense solutal Marangoni flow.

Finally, in figure 4(h), the average temperature of the liquid–gas interface $\bar{T}_{lg}$ is depicted. Again, cooling can be up to 3.5 K and is most pronounced during the initial regime, where both components evaporate. Although water has a $2.7$ times higher specific latent heat than ethanol, most of the initial cooling is induced by the up to $9$ times faster evaporation rate of ethanol.

3.4 Influence of convective transport in the gas phase

Before addressing the ternary Ouzo droplet, the flow in the gas phase is investigated. Until now, this flow has been disregarded, but since it generates convective vapour and thermal transport, it can influence the evaporation rates and the temperature field. The flow in the gas phase can most generally be driven by four mechanisms, namely forced convection, natural convection, Marangoni convection and Stefan flow. While forced convection can be ruled out due to the geometry of the experimental set-up, the mass density gradient required for natural convection can be caused by a thermal gradient as well as by vapour concentration gradients. Since the temperature at the droplet is lower than the ambient temperature, thermally driven natural convection can be ruled out. To estimate the influence of solutally induced natural convection, we define the solutal Rayleigh number

where $g$ is the acceleration due to gravity, $\unicode[STIX]{x1D708}$ the kinematic viscosity of the gas phase (assumed to be independent of the vapour concentration) and $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D6FC}}=(\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{}^{g}/\unicode[STIX]{x2202}c_{\unicode[STIX]{x1D6FC}})/\unicode[STIX]{x1D70C}_{}^{g}$ is the solutal expansion coefficient, which can be calculated by the ideal gas law. The solutal Rayleigh numbers with the largest modulus of the simulations discussed so far read $Ra_{w}\sim 1$ and $Ra_{e}\sim -10$ for water vapour and ethanol vapour, respectively. Note that the different signs originate from the fact that water vapour is less dense than dry air, whereas ethanol vapour is more dense. However, according to Dietrich et al. (Reference Dietrich, Wildeman, Visser, Hofhuis, Kooij, Zandvliet and Lohse2016b ), an influence of solutally driven natural convection on the evaporation rate can be ruled out at these Rayleigh numbers.

Since the tangential velocity is continuous at the interface, the fast Marangoni flow will be also present in the gas phase. Furthermore, the density difference of liquid and vapour leads to a discontinuous jump in the normal velocity component, which constitutes the so-called Stefan flow. To investigate the influence of these effects, the Stokes flow is solved in the gas phase as well, where the normal component of the gas velocity is imposed analogously to (3.9) and the tangential component is imposed to be continuous with its counterpart in the liquid. Furthermore, the convection–diffusion equations (3.1) and (3.8) are augmented with the corresponding convective term in the gas phase and the additional convective vapour transport at the interface, i.e. $c_{\unicode[STIX]{x1D6FC},VLE}(\boldsymbol{u}_{g}-\boldsymbol{u}_{lg})\boldsymbol{\cdot }\boldsymbol{n}_{lg}$ , is added to the evaporative mass fluxes.

Figure 5. (a,b) Same as in figure 4 at two different times $t=25~\text{s}$ (a) and $t=250~\text{s}$ (b), but now considering in the calculation and showing in the picture the flow in the gas phase (depicted on the right side). (c–e) Comparison of relevant quantities obtained by simulations with and without flow in the gas phase, revealing hardly any difference.

In figure 5(a,b), the simulation of the binary water–ethanol droplet of figure 4 with consideration of Marangoni and Stefan flow in the gas phase is depicted. From the velocity field on the right side of the snapshots, it is apparent how the Marangoni flow creates vortices in the gas phase, which are initially chaotic and become regular once the ethanol has evaporated. The contribution of the Stefan flow, i.e. the normal velocity jump, is barely visible, since the typical Stefan flow velocity is only approximately $1~\text{mm}~\text{s}^{-1}$ . As can be inferred from figure 5(c–e), the additional convective transport of vapour and energy has virtually no influence on the volume evolution, the fluid velocity and the interfacial temperature. Hence, disregarding the flow in the gas phase is justified for the discussed cases.

Furthermore, the simulations with flow in the gas phase allow us to validate the assumption of omitting the shear stress term $\boldsymbol{n}_{lg}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}^{g}\boldsymbol{\cdot }\boldsymbol{t}_{lg}$ in (3.11): the ratio of the shear stresses in the gas phase and in the liquid phase is indeed below 2 %.

3.5 Ternary Ouzo droplet

Now that the model has been validated based on a pure water and a binary water–ethanol droplet, we focus on a ternary Ouzo droplet at $T_{\infty }=22.5~^{\circ }\text{C}$ and $\unicode[STIX]{x1D719}=40\,\%$ in the following. To that end, we took the experimental data of our previous article (Tan et al. Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016) and performed corresponding simulations with the finite element method described in § 3.1. Different from the lubrication theory model used in our previous article (Tan et al. Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016), the present model takes thermal effects into account and is not subject to the limitation of the lubrication approximation where the contact angle has to be small. In the previous article, we furthermore had to adjust the relative humidity to reproduce the experimental data with the isothermal lubrication theory model in Tan et al. (Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016).

Figure 6. (a–c) Snapshots of a ternary Ouzo drop with ethanol mass fraction $y_{e}$ in the liquid phase and water vapour concentration in the gas phase (left) and temperature field (right) at three different times $t=20~\text{s}$ (a), $t=22.5~\text{s}$ (b) and $t=25~\text{s}$ (c). Regions where the Ouzo effect occurs are shaded white on the left side. The arrows at the interface indicate the evaporation rates of water (left) and ethanol (right). A corresponding movie is available as supplementary material (movie 3). (d) Volume evolution according to experiment and the full simulation. (e) Partial masses of all components in the drop. (f) During the evolution, the droplet shape temporarily deviates from a spherical-cap shape due to the presence of the oil ring. During this period, the contact angle in the model was set to the contact angle $\unicode[STIX]{x1D703}^{\ast }$ of the spherical-cap-shaped upper part. Otherwise, the usual contact angle $\unicode[STIX]{x1D703}$ was imposed. (g) Mean tangential fluid velocity at the interface and its root mean square average. (h) Average interface temperature.

In figure 6(a–c), the numerically determined onset of the Ouzo effect is depicted (white regions on the left side). As in the experiment, the nucleation of oil microdroplets starts at about $t=20~\text{s}$ at the contact line. Different from the predictions of the isothermal lubrication theory model (Tan et al. Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016), the flow in the droplet is chaotic, which leads to a chaotic propagation of the phase-separation front through the droplet. At $t\approx 30~\text{s}$ , nucleation can occur everywhere in the droplet.

As it is apparent from the volume evolution curve in figure 6(d), the initial slope matches the experimental results perfectly, but shows deviations at later times. Since the volume evolution of the pure water droplet and the binary water–ethanol droplet perfectly match the experiments, the disagreement in the ternary case can definitely be attributed to the presence of the third component, namely the anise oil. While it is initially miscible, which is an assumption of the model, the simulation shows perfect agreement in the first regime. At later times, however, the oil ring starts to form, which changes the geometry from the typical spherical-cap shape to a more complicated shape, i.e. a water droplet sitting on an oil ring. This leads to a reduction of the water–air interface area and thereby reduces the evaporation rate. The model cannot account for this change in geometry since it always assumes a perfect spherical-cap shape and hence predicts a faster evaporation rate. Another possible factor that could contribute to the disagreement is the fact that the anise oil used is not a pure grade of trans-anethole and consists of approximately 10 % of unknown components (Rodrigues et al. Reference Rodrigues, Rosa, Marques, Petenate, Meireles and Angela2003). These could also influence the evaporation rate, e.g. by forming a shielding monolayer suppressing evaporation (Langmuir & Schaefer Reference Langmuir and Schaefer1943).

The other quantities depicted in figure 6(e–h), e.g. the fluid velocity and temperature, are qualitatively similar to the case of the binary water–ethanol droplet of figure 4. While there are temporarily two contact angles in the experiment, i.e. the contact angle $\unicode[STIX]{x1D703}$ between the oil ring and the substrate and the contact angle $\unicode[STIX]{x1D703}^{\ast }$ of the spherical-cap-shaped upper part of the droplet, the present simulation can only take a single contact angle. Here, $\unicode[STIX]{x1D703}^{\ast }$ was imposed as the contact angle in the simulation during the time the oil ring was present in the experiment.

4.1 Three-dimensional model

To investigate the axial symmetry breaking within the framework of the model, it has to be generalized to three dimensions. Mathematically, this generalization from the two-dimensional axisymmetric model to a full three-dimensional variant is an easy step. One just has to transfer the weak forms (cf. Diddens Reference Diddens2017) to their corresponding three-dimensional equivalents. However, from the perspective of implementation, one faces the challenge of having to treat a free boundary problem. In the axisymmetric model, the movement of the sharp liquid–gas interface was made possible by shifting interfacial mesh nodes and performing a local reconstruction of the mesh whenever necessary to prevent the elements from collapsing. For the interpolation from the previous mesh to the reconstructed mesh, the supermesh method introduced by Farrell et al. (Reference Farrell, Piggott, Pain, Gorman and Wilson2009) was utilized. Both the local mesh reconstruction and the supermesh method constitute a fundamental challenge to generalize to three dimensions. To circumvent these problems, we focus on a static three-dimensional mesh instead. This approach does not allow us to perform a single simulation over the entire lifetime of the droplet, but can still capture the three-dimensional flow statistics for a specific droplet composition. The use of a static mesh is possible, since the typical flow velocity in the droplet is several orders of magnitude larger than the movement of the interface.

Figure 7. (a) View into a mesh for the three-dimensional model. (b) Axial symmetry breaking for snapshot time $t=150~\text{s}$ at different three-dimensional simulation times $t_{3d}$ . The droplet is depicted in top view showing the ethanol mass fraction at the interface. The composition and the tangential fluid velocity at the interface (white arrows) quickly break up into a chaotic, highly non-axisymmetric behaviour with a cellular structure. (c) Same as (b), but for the snapshot time $t=200~\text{s}$ . Due to the low ethanol concentration, the breakup is slower. At later times, a cellular pattern can be observed along the rim, which is pulled in filaments towards the apex by thermal Marangoni flow. (d) Evolution of the root mean square averaged tangential fluid velocity components $\bar{u}_{t,rms}^{apex}$ and $\bar{u}_{t,rms}^{circum}$ for the simulation of (c). After a transient regime, the dynamics end up in a chaotic long-time regime with converged mean values and standard deviations.

The procedure is as follows: for selected individual times $t$ of the axisymmetric simulation depicted in figure 4, we project the axisymmetric data onto a corresponding three-dimensional mesh (cf. figure 7 a). This is used as initial condition for the generalized three-dimensional model, which is integrated over another time $t_{3d}$ to obtain the characteristic three-dimensional flow at time instant $t$ . Again, the vapour diffusion equations for the volatile components are solved in the gas phase, the convection–diffusion equations for the composition and the Stokes flow are solved in the droplet, and the temperature field is solved in all domains, including the substrate and the air below it. Since the mesh is static, the kinetic boundary condition has to be modified to $\boldsymbol{n}_{lg}\boldsymbol{\cdot }\boldsymbol{u}=0$ . Furthermore, to capture the characteristic three-dimensional flow at time instant $t$ , one has to ensure that the average composition is not changing within the simulation time $t_{3d}$ of the three-dimensional model. This can be achieved by augmenting the convection–diffusion equation (3.6) by a spatially uniform correction term which compensates for the evaporation, i.e.

The axisymmetry of the initial condition from the axisymmetric model can break up by the Marangoni instability. Azimuthal perturbations of the composition can arise, yielding to non-axisymmetric solutal Marangoni flow which in turn feeds back to the coupled dynamics of multicomponent flow, evaporation and thermal effects. A representative example is shown in figure 7(b) for the snapshot time $t=150~\text{s}$ , which corresponds to the axisymmetric initial condition of figure 4(b). Within a very short three-dimensional simulation time $t_{3d}$ , the axial symmetry of mixture composition and Marangoni flow is broken and the system exhibits highly non-axisymmetric dynamics with a cellular structure.

When the ethanol concentration is lower and thermal Marangoni flow becomes relevant, the symmetry breakup is slower. This can be seen from the simulation at snapshot time $t=200~\text{s}$ depicted in figure 7(c). In an initial transient regime, perturbations arise at the apex and near the contact line, which coalesce to four large cells. With ongoing simulation time, however, new perturbations arise at the contact line, which form small spots with enhanced ethanol concentration. Due to thermal Marangoni flow, filaments of enhanced ethanol concentration are pulled out of these spots towards the apex. The spots near the contact line have a remarkable similarity to the hydrothermal waves in ethanol droplets on heated substrates as reported by Sefiane et al. (Reference Sefiane, Moffat, Matar and Craster2008b ) and Sobac & Brutin (Reference Sobac and Brutin2012). The fundamental difference is, however, that the hydrothermal waves are a result of an intense thermal Marangoni effect alone, whereas here a combination of thermal and solutal Marangoni flow induces these structures.

For even lower ethanol concentrations, i.e. for snapshot times $t>225~\text{s}$ , the initial axisymmetry persists, i.e. a breaking of the axial symmetry cannot be observed.

To determine the typical three-dimensional flow characteristics, the simulation time $t_{3d}$ has to be sufficiently long. In particular, we are interested in the long-time regime only, i.e. not in the transient regime as e.g. depicted in the central picture of figure 7(c). To determine whether the simulation has already reached its long-time behaviour, the root mean square averaged tangential fluid velocity at the interface is investigated. In order to distinguish between axisymmetric behaviour and broken symmetry later on, the root mean square is split into its two components, namely

where $\boldsymbol{t}_{lg}^{apex}$ is the interface tangent pointing towards the apex and $\boldsymbol{t}_{lg}^{circum}$ is the tangent pointing along the circumference of the droplet. Hence, a non-vanishing value of $\bar{u}_{t,rms}^{circum}$ with respect to $\bar{u}_{t,rms}^{apex}$ indicates the absence of axial symmetry. A typical evolution of these quantities is depicted in figure 7(d). Initially, the Marangoni instability breaks the axial symmetry, which results in an increasing value of $\bar{u}_{t,rms}^{circum}$ . After that, it can take several seconds of simulation time $t_{3d}$ until the system enters its long-time behaviour. The long-time behaviour regime can be identified when both $\bar{u}_{t,rms}^{apex}$ and $\bar{u}_{t,rms}^{circum}$ show a converged mean value and variance.

Figure 8. Evolution of the root mean square averaged tangential velocity components based on the long-time behaviour of individual three-dimensional simulations corresponding to the axisymmetric simulation of figure 4. The error bars indicate the standard deviation and the insets show typical snapshots of the three-dimensional simulations at the indicated times. The numbers below the snapshot are the time and the averaged ethanol mass fraction. The scale bar is the same for all snapshots.

The characteristic flow statistics were finally extracted based on at least $t_{3d}=30~\text{s}$ of long-time behaviour for different snapshot times $t$ . The resulting data for the tangential fluid velocity components are depicted in figure 8. Initially, the droplet is highly non-axisymmetric with $\bar{u}_{t,rms}^{apex}\approx \bar{u}_{t,rms}^{circum}$ . With ongoing time, ethanol evaporates and since the composition dependence of the surface tension has its steepest gradient in the limit of dilute ethanol, the flow becomes faster. After a maximum at approximately $t\approx 130~\text{s}$ , the average velocity decreases again. Until $t=175~\text{s}$ , the characteristic flow pattern is predominantly driven by solutal Marangoni flow. The ethanol concentration at the interface exhibits spatio-temporally chaotic cellular patterns. Moreover, it is possible that a net flow from one side of the rim to the opposing side can build up (e.g. from right to left in the inset at $t=175~\text{s}$ ). After $t=175~\text{s}$ , i.e. after the local minimum of $\bar{u}_{t,rms}^{apex}$ , a transition can be observed. The flow is mainly directed towards the apex due to the thermal Marangoni effect, but the dilute ethanol still causes symmetry-breaking perturbations (cf. inset at $t=185~\text{s}$ ). This characteristic flow pattern has already been discussed in figure 7(c). With ongoing time and vanishing ethanol, the flow becomes more and more axisymmetric and regular, which can be inferred by the vanishing value of $\bar{u}_{t,rms}^{circum}$ .

4.2 Comparison with experiment

As described in § 2.3, micro-PIV measurements were performed during the evaporation of the binary water–ethanol droplet. Thereby, experimental data for the flow in the focal plane, i.e. close to the substrate became available. Since the simulation solves the entire flow field in the bulk, it is possible to extract the corresponding velocity by slicing the three-dimensional mesh at the focal plane and projecting the local velocity vectors onto it. This allows for a direct comparison of the experimentally and numerically determined velocity.

Figure 9. (a–f) Experimental snapshots of the velocity field in the focal plane at different instants. The arrows display the local velocity and the vorticity is colour coded. The yellow circle indicates the intersection of the liquid–air interface and the focal plane. (a) Small vortices are present close to the liquid–air interface, breaking the axial symmetry in the entire droplet. (b) Small vortices can coalesce into a pairs of bigger vortices. (c,d) The overall flow velocity and vorticity increase. (e,f) After a while, the flow has an inward flow pattern (but not perfectly axisymmetric) and finally calms down. (g–l) Velocity in the focal plane as extracted from the numerics. While the simulations shows initially good agreement with the experiment, deviations are present at later times. Note that the colour code is the same for both experiments (a–f) and numerics (g–l). Corresponding movies are available as supplementary material (movies 4 and 5).

Representative snapshots of the experimental micro-PIV measurement are depicted in figure 9(a–f). It is apparent that the axisymmetry is initially broken by the presence of multiple vortices near the liquid–air interface (a,b). With ongoing time, the flow gets more intense and more chaotic up to a maximum at approximately $t=100~\text{s}$ (c,d). Directly after that, a net flow towards the centre of the focal plane is building up, which is not perfectly axisymmetric but globally it is directed radially inward (e). The flow is slowing down until it completely stops at approximately $t=140~\text{s}$ (f).

While the initial chaotic stage in figure 9(a–d) qualitatively agrees with the data of Christy et al. (Reference Christy, Hamamoto and Sefiane2011) and Bennacer & Sefiane (Reference Bennacer and Sefiane2014), we could not reproduce the remarkable peak in the radial outward flow during the transition regime, which is mentioned in these articles. Instead, we find a nearly axisymmetric inward flow (e) as also reported by Zhong & Duan (Reference Zhong and Duan2016). This difference can presumably be attributed to different set-ups, e.g. different contact angles, but it points out the importance of further studies. Since we have not measured the flow until the end of the droplet lifetime, figure 9(f) does not show the usual radial outward capillary flow reported in the previous studies as last regime.

Figure 10. Averaged velocity (a) and vorticity (b) in the focal plane based on the experiment and the numerical simulation. The error bars indicate the standard deviation of the numerically obtained quantities.

The data extracted from the numerical simulation are depicted in figure 9(g–l) for comparison. It is apparent how the typical flow is very similar to the micro-PIV data at early times (g,h), but in the numerics the maximum averaged flow intensity can only be found at $t\approx 150~\text{s}$ (i,j), i.e. when the flow has already ceased in the experiments. After this maximum, the flow does not turn inward, but outward instead, which is caused by the thermal Marangoni effect. As long ethanol resides in the droplet, deviations from perfect axisymmetric thermal Marangoni flow can be seen (k). These eventually decrease over time until the flow finally converges to a perfectly axisymmetric intense thermal Marangoni flow. The initial agreement of experiment and simulation can also be inferred from figure 10(a), where the averaged velocities in the focal plane are plotted following both methods. At later times, however, the disagreement is obvious.

Since the disagreement is most pronounced at later times, it cannot be attributed to the presence of a binary mixture. In fact, for $t\gtrsim 250~\text{s}$ the droplet is virtually a pure water droplet, which should show an intense thermal Marangoni flow according to the simulation, which is not observed in experiment. This disagreement between prediction and observation of thermal Marangoni flow in water droplets is well known (Hu & Larson Reference Hu and Larson2005). At $t=250~\text{s}$ , the numerically obtained temperature difference between contact line and apex reads $\unicode[STIX]{x0394}T\approx 0.57~\text{K}$ , which corresponds to a surface tension difference of $\unicode[STIX]{x0394}\unicode[STIX]{x1D70E}\approx -0.084~\text{mN}~\text{m}^{-1}$ , i.e. 0.1 % of the averaged surface tension. As already discussed by Hu & Larson (Reference Hu and Larson2005), this small thermally induced difference can be counteracted by a very small amount of surface-active contaminants, which are unavoidable during the experiment.

Therefore, the present scenario, i.e. good agreement at initial stages and large disagreement in the limit of pure water, can presumably be attributed to the presence of contaminants. Since the initial solutal Marangoni instability is intense and chaotic, the contaminants are unable to arrange to a counteracting distribution along the interface. With ongoing time, three mechanisms can contribute to an increase of the overall contaminant concentration at the interface. These are the fact that the droplet and thus the interface are shrinking, a possible ongoing adsorption of contaminants from the gas phase during the evaporation process, and a possibly increasing interface affinity of the contaminants with decreasing ethanol concentration. Once the contaminant concentration is sufficiently high, they start to stabilize the solutal Marangoni flow leading to the almost symmetric radially inward flow as seen in figure 9(e). Eventually, the thermal Marangoni flow is completely suppressed by the contaminants (cf. figure 9 f). However, to substantiate this scenario, more elaborate experiments, in particular measurements of the contaminant concentration at the interface, or a generalization of the model considering surface-active contaminants would be necessary. It can also not be ruled out that the hydrophobic polystyrene tracer particles influence the interface dynamics.