Competing Marangoni and Rayleigh convection in evaporating binary droplets

Abstract For a small sessile or pendant droplet it is generally assumed that gravity does not play any role once the Bond number is small. This is even assumed for evaporating binary sessile or pendant droplets, in which convective flows can be driven due to selective evaporation of one component and the resulting concentration and thus surface tension differences at the air–liquid interface. However, recent studies have shown that in such droplets gravity indeed can play a role and that natural convection can be the dominant driving mechanism for the flow inside evaporating binary droplets (Edwards et al., Phys. Rev. Lett., vol. 121, 2018, 184501; Li et al., Phys. Rev. Lett., vol. 122, 2019, 114501). In this study, we derive and validate a quasi-stationary model for the flow inside evaporating binary sessile and pendant droplets, which successfully allows one to predict the prevalence and the intriguing interaction of Rayleigh and/or Marangoni convection on the basis of a phase diagram for the flow field expressed in terms of the Rayleigh and Marangoni numbers.


Introduction
Evaporating droplets frequently occur in nature and applications, be it a rain droplet evaporating on a leaf, a droplet on a hot surface in spray cooling, a droplet of insecticides sprayed on a leaf or an inkjet-printed ink droplet on paper. Many of such droplets are multicomponent, i.e. consisting of a mixture of liquids. From the physical point of view, an evaporating multicomponent droplet in a host gas is paradigmatic for combined multiphase and multi-component flow including a phase transition. Scientifically, this process encompasses the various fields of fluid mechanics, thermodynamics and also aspects from the field of chemistry. The evaporation dynamics is also relevant for the deposit left behind the evaporation of a particle-laden droplet. Here, pioneering work was done by Deegan et al. (1997) around 20 years ago, when they identified the coffee-stain effecti.e. the phenomenon of finding a typical ring structure of deposited particles after the evaporation of a coffee droplet -and successfully explained it by the combination of a non-uniform evaporation rate along the droplet interface and a pinned contact line. In applications, one usually wants to prevent such coffee-stain effect, e.g. for obtaining a homogeneous deposition pattern in inkjet printing (Park & Moon 2006;Kuang et al. 2014;Hoath 2016;Sefiane 2014). For reviews on evaporating pure droplets we refer to Cazabat & Guéna (2010) and Erbil (2012).
Preventing the coffee-stain effect can be achieved by altering the flow inside the droplet during the drying process by inducing gradients in the acting forces. Focussing on the interfacial forces first, a tangential gradient of the surface tension along the liquid-gas interface leads to the well-known Marangoni effect, i.e. a tangential traction that drives the liquid towards positions of higher surface tension (Scriven & Sternling 1960;Pearson 1958). By that, the entire flow in the droplet can be altered from the typical outwards flow towards the contact line to a recirculating flow driven by a persistent Marangoni effect (Hu & Larson 2006). For the case of a pure droplet, the necessary gradient in surface tension can be generated by thermal effects, e.g. either self-induced by latent heat of evaporation or externally imposed by heating or cooling the substrate (Girard et al. 2006;Sodtke et al. 2008;Dunn et al. 2009;Tam et al. 2009).
The other mechanism to induce Marangoni flow is known as solutal Marangoni effect, which is usually much stronger. For solutal Marangoni flow, the droplet must consist of more than one component, e.g. a solvent and one or more surfactants (Still et al. 2012;Marin et al. 2016;Kwieciński et al. 2019) or a solvent and possibly multiple co-solvents (Sefiane et al. 2003;Christy et al. 2011;Tan et al. 2016;Li et al. 2018) or dissolved salts (Soulié et al. 2015;Marin et al. 2019). For a recent perspective review on droplets consisting of more than one component, we refer to Lohse & Zhang (2020). The difference in the volatilities of the individual constituents leads to preferential evaporation of one or the other component and thereby compositional gradients are induced. Since the surface tension is a function of the composition and due to the non-uniform evaporation profile, a surface tension gradient along the liquid-gas interface can build up and result in a similar Marangoni circulation as in the thermally-driven case. The nature of the resulting flow can be quite different, mostly depending on whether the evaporation process leads to an overall decreasing or increasing surface tension, i.e. whether the more volatile component has a higher or lower surface tension than the less volatile component.
In a binary droplet consisting e.g. of water and glycerol, with water being more volatile and having the higher surface tension, the overall surface tension decreases during the preferential evaporation of water and the resulting Marangoni flow is usually regular, axisymmetric and directed towards the position of the lowest evaporation rate of water, i.e. towards the contact line for contact angles above 90 • and towards the apex for contact angles below 90 • (Diddens et al. 2017a;Diddens 2017).
On the contrary, e.g. in case of a binary droplet consisting of water and ethanol, where the overall surface tension increases due to the predominant evaporation of ethanol, the typical Marangoni effect is way more violent and chaotic (Christy et al. 2011;Bennacer & Sefiane 2014). Here, in particular, the axial symmetry of the droplet is usually broken, leading to a complicated scenario of initially chaotic flow driven by the solutal Marangoni effect and followed by either thermal Marangoni flow or the typical coffee-stain flow, when almost only water is left (Diddens et al. 2017b). Remarkably, the presence of a strong Marangoni effect can also have a significant influence on the shape and wetting behaviour of droplets (Tsoumpas et al. 2015;Karpitschka et al. 2017). Finally, the evaporation of mixture droplets can show a variety of additional intriguing phenomena, e.g. multiple phase changes and microdroplet nucleation in ternary droplets like ouzo (Tan et al. 2016(Tan et al. , 2017, and phase segregation in binary droplets (Li et al. 2018) or rather homogeneous deposition patterns by an interplay of Marangoni flow, surfactants and polymers (Kim et al. 2016).
As highlighted above, besides the gradient in the surface tension, i.e. in the interfacial forces, also gradients in the mass density, i.e. in the bulk force due to gravity, can influence the flow by natural convection. Similar to the surface tension, the mass density is a function of the temperature and, in the case of mixtures, of the composition, so that thermally and solutally driven natural convection can be realized in evaporating droplets. Flow driven by natural convection is one of the most important fields of fluid mechanics, as e.g. in Rayleigh-Bénard systems, however, these are usually investigated at large spatial dimensions. For small droplets, on the other hand, one would naively expect that the flow in case of thermal or solutal gradients is predominantly driven by the Marangoni effect, since a small droplet is associated with a small Bond number and hence surface tension effects would dominate over gravity. As a consequence, most studies on droplet evaporation focus on the Marangoni effect, but disregard the presence of natural convection by this argument.
Recent studies, however, showed that even for small droplets with small Bond numbers, the internal flow can be decisively determined by natural convection and not by Marangoni flow (Edwards et al. 2018;Li et al. 2019a). This has been even found at the later stages of water-ethanol droplets, which initially show a very intense chaotic Marangoni flow (Edwards et al. 2018). Obviously, these findings give rise to the following question: Under what circumstances which kind of flow pattern can be found in an evaporating binary droplet, i.e. when is the flow dominated by the Marangoni effect and when by natural convection?
In this manuscript, we answer this questions by carefully investigating both kinds of driving forces and their mutual interaction. The corresponding effects can be quantified by non-dimensional numbers, namely the Marangoni number for flow due to surface tension gradients and the Rayleigh (or Archimedes/Grashof) number for the natural convection. By considering quasi-stationary instants during the drying process, these numbers successfully allow to predict the flow inside the droplets on the basis of phase diagrams in the Ra-Ma parameter space. We also validated these phase diagrams with full simulations and corresponding experiments.
The paper is organised as follows: We will first present the complete set of dynamical equations describing the evaporation of a binary mixture droplet. In section 3, these equations are solved to discuss an illustrative example case. We will then introduce the quasi-stationary approximation in section 4 and discuss the phase diagrams obtained by this model in section 5. The paper ends with a conclusion and a comparison with experimental data in the appendix.

Governing Equations
The evaporation of a mixture droplet is a multi-phase and free interface problem with multi-component dynamics in both the liquid and gas phase. For a binary droplet, the liquid is constituted by two components, α = A, B, whereas the gas phase is in general a ternary gas mixture of the ambient medium, e.g. air, and the vapours of the components A and B. When the droplet evaporates at a temperature T far below the boiling point and in absence of forced or strong natural convection, the diffusive dynamics in the gas phase can be approximated by the quasi-stationary Laplace equation (Deegan et al. 2000;Hu & Larson 2002;Popov 2005;Diddens et al. 2017a) for the vapour concentrations c α , i.e. the partial mass densities. The corresponding boundary conditions are given by the vapour-liquid equilibrium according to Raoult's law at the liquid-gas interface and the ambient vapour concentration at infinity, i.e. interface kinematic boundary condition with mass transfer (2.9), Marangoni shear stress (2.11), preferential evaporation (2.6) gravity sessile pendant figure 1. Schematic of the model. The problem is considered to be axisymmetric and isothermal. The flow and the advection-diffusion equation for the composition in the droplet is solved with consideration of gravity and the composition-dependence of the liquid mass density, dynamic viscosity and diffusivity. The transport in the gas phase is assumed to be diffusion-limited. At the interface, Raoult's law is used to enforce the vapour-liquid equilibrium, mass transfer due to evaporation is considered and the Marangoni shear stress due to a composition-dependent surface tension is taken into account.
where c pure α is the saturation vapour concentration in case of the pure liquid α and x α is the liquid mole fraction. The activity coefficients γ α account for thermodynamic non-idealities and are functions of the composition, i.e. of x α . Neglecting the small contribution of the Stefan flow at temperatures sufficiently below the boiling point, the evaporation rates j α are given by the diffusive fluxes at the liquid-gas interface, i.e.
While the dynamics in the gas phase can be considered in the diffusive and quasistationary limit, convection can be dominant in the liquid phase, which can be attributed to the typical diffusion coefficients, namely D vap α ∼ 10 −5 m 2 /s in the gas phase and D ∼ 10 −9 m 2 /s in the liquid phase. Therefore, the liquid phase has to be described by the full convection-diffusion equation for the liquid mass fraction y α , which is expressed for the component A only due to the identity y A + y B = 1, i.e. ρ (∂ t y A + u · ∇y A ) = ∇ · (ρD∇y A ) . (2.5) The liquid density ρ and the mutual diffusivity D are in general functions of the composition, i.e. of y A . The mass transfer rates j α due to evaporation induce a change in composition at the liquid-gas interface. Using the mass transfer expression j α = ρy α (u α −u I )·n with u α denoting the velocity of component α and u I and n denoting the interface velocity and normal, respectively, this compositional change can be expressed by a Robin boundary condition for Eq. (2.5) in the frame co-moving with the interface in normal direction, namely (2.6) Finally, the flow in the droplet is given by the Navier-Stokes equations Here, we have chosen the z-axis to point towards the apex of the droplet, i.e. a sessile droplet and a pendant droplet can be realized by negative and positive values for g, respectively. Note that the viscosity µ and the mass density ρ are in general functions of the composition y A . A dependency on the temperature is disregarded in the following due to the fact that thermal effects at lower temperatures are usually considerably inferior to the impact of solutal gradients. For the contact line dynamics, we are focussing here on a pinned contact line, i.e. evaporation in the constant radius mode (CR-mode, Picknett & Bexon (1977); Stauber et al. (2014)). To resolve the incompatibility of a no-slip boundary condition at the substrate and the evaporative mass loss at the contact line, we impose a Navier-slip boundary condition with a small slip-length in the nanometre scale instead. This effectively resembles a no-slip boundary condition in the main part of the dropletsubstrate interface, but still allows for a consistent mass transfer at the contact line. The free liquid-gas interface is subject to the kinematic boundary condition considering the mass transfer, i.e.
and furthermore to the Laplace pressure in normal direction − p + µn · ∇u + (∇u) t · n = σκ , (2.10) where the traction in the gas phase has been neglected due to the viscosity ratio. Here, σ is the local surface tension, κ the curvature of the interface and p denotes the pressure difference with respect to the ambient gas pressure. Finally, also the Marangoni shear stress in tangential direction has to be considered: (2.11) Here, ∇ t = (1 − nn) · ∇ is the surface gradient operator. A sketch of the model is depicted in figure 1.

Numerical solution of the dynamical equations for an instructive example
In order to solve the given set of equations numerically, we have generalised the sharpinterface arbitrary Lagrangian-Eulerian finite element method described in Diddens (2017) by considering the gravitational force, and also validated it by a more general reimplementation of the same model with the finite element package oomph-lib † (Heil & Hazel 2006), which allows for interface deformations and considers the general continuity equation (2.8). The latter method has been successfully validated against various experiments (Li et al. 2018(Li et al. , 2019aGauthier et al. 2019;Li et al. 2019b).
In figure 2, a simulation of a sessile glycerol-water droplet (initially 5 wt.% glycerol) with an initial volume of 1 µL and an initial contact angle of 120 • evaporating at a constant temperature of 22 • C and a relative humidity of 20 % is shown. The contact line remains pinned during drying and glycerol (liquid B) is assumed to be non-volatile due to its low volatility compared to water (liquid A), i.e. c B = c ∞ B = 0 and j B = 0. For more details about these kind of simulations we refer to Diddens et al. (2017b);Diddens (2017), where however we did not consider of the influence of gravity.
Initially, in figure 2(a), one can see a single vortex in the entire droplet, directed from the apex towards the contact line. This vortex is generated for two reasons, namely Marangoni convection and natural convection (Rayleigh convection). Due to the enhanced  evaporation rate of water at the apex at the high contact angle, the water content is predominantly reduced at the top of the droplet, resulting in a lower surface tension compared to the region near the contact line. This drives a Marangoni flow towards the contact line. Since glycerol is more dense than water, the glycerol-rich outer shell of the droplet also sinks down due to natural convection, which also results in a flow from the apex to the contact line due to the spherical geometry. Hence, both mechanisms support recirculating flow in the same direction.
Remarkably, in figure 2(b), the situation changes. The contact angle is still above 90 • , i.e. still having the highest water evaporation rate at the top of the droplet. According to the afore-mentioned discussion, one would still naively expect the same kind of single vortex flow. However, the simulation clearly shows two vortices, one in the bulk driven by natural convection (white) and another one close to the interface, which is driven by Marangoni flow in the opposite direction (black). The reason why the Marangoni flow is reversed, i.e. why there is more water at the top of the droplet although the evaporation rate of water is still dominant at the apex, is the fact that there is enhanced water replenishment by diffusion at the apex, which compensates for the rather small difference in the evaporation rates at the top and near the contact line. This can be seen by the rather steep concentration gradient in normal direction at the apex compared to the region near the contact line. The reason of the steep concentration gradient in normal direction close at the apex is the upward directed convective water replenishment from the bulk, which is governed by the internal vortex driven by natural convection. This means that sufficiently strong natural convection in the bulk can reverse the Marangoni flow at the interface, although one would not anticipate this by just considering the profile of the evaporation rate at this contact angle.
Upon further evaporation, in figure 2(c), the contact angle falls below 90 • , resulting in a higher water evaporation rate near the contact line. Hence, less water is present at the contact line as compared to the apex due two facts, namely the effect of preferential evaporation at the contact line and the lower replenishing diffusive flux of water from the bulk at the contact line. Thereby, the Marangoni flow gets enhanced compared to the situation in figure 2(b) and the relative size of the Marangoni-induced vortex at the interface grows at the expense of the counter-rotating bulk vortex by natural convection.
Finally, in figure 2(d), the contact angle becomes rather small so that the Marangoni flow at the interface is even stronger due to the enhanced non-uniformity of the evaporation rate. Furthermore, the influence of natural convection also diminishes rather quickly, i.e. with cubic power in terms of the length scale according to the Rayleigh number (see later for its definition), due to the reduced volume of the droplet. This results in the depicted situation, i.e. that the flow direction within the entire droplet is completely determined by the Marangoni effect.
In a nutshell, one can infer from the direct numerical simulation results in figure 2 that there can be multiple flow scenarios during the drying of a single binary droplet, driven by an interplay of natural (i.e. Rayleigh) convection and Marangoni convection. One also clearly sees that, for a particle laden droplet, the coffee-stain effect would not occur as there is no noticeable flow towards the contact line (which for pure evaporating droplets transports the suspended particles to the rim of the pinned droplet) as compared to the strongly recirculating flow due to Marangoni flow and gravity.

Quasi-stationary approximation of the dynamical equations
After discussing some possible flow scenarios by considering a representative numerical example in the previous section, we will now focus on a simplification of the full model described in section 2. We generalise again from the particular case of a water-glycerol droplet to the general case of a binary droplet, where both liquids A and B are allowed to evaporate. The goal is to find the simplest model possible that allows to predict the expected flow scenario in the droplet by a minimum number of non-dimensional quantities.

Evaporation Numbers
As shown in the example simulation, the liquid recirculates multiple times during the evaporation process due to the fast flow in the droplet. Hence, the typical liquid velocity u is much larger than the normal interface movement velocity u I . Moreover, this leads to a rather well-mixed droplet, i.e. with typical compositional deviations of about a few percent in terms of mass fractions. These observations allow for several simplifications of the model. First of all, the liquid composition is expanded into two terms, i.e. y A (x, t) = y A,0 + y, namely the spatially averaged composition y A,0 , which slowly evolves over the entire drying time, and the small local composition deviations y(x, t). Since the composition-dependent liquid properties are usually rather smooth functions of the composition, this separation can be transferred to a first order Taylor expansion of the liquid properties, i.e.
Since the averaged composition y A,0 evolves slowly, this expansion can be done at any specific time of interest during the evaporation process. In particular, this means that the coefficients of the Taylor expansions (4.1) can be treated as constants during some time close to the considered instant. This allows us to introduce the following nondimensionalized scales where the spatial scale is chosen in that way, that the nondimensionalized droplet volumẽ V becomes unity. In a next step, the vapour fields are decomposed in a similar manner as (4.1), namely in a normalized contributionc 0 which is one at the interface and zero at infinity and a contributionc ∆ accounting for the effect of local composition variations on the vapour concentration via Raoult's law to the first order, i.e.
The Laplace equation ( Thereby, the evaporation rates (2.4) separate in the same way, i.e.
wherej 0 = −∂ nc0 only depends on the shape of the droplet, i.e. resembles the normalized evaporation profile of a homogeneous droplet, andj y = −∂ nc∆ is a linear functional of y, i.e. the Dirichlet-to-Neumann map, accounting for deviations in the evaporation rate due to a varying interfacial composition via the composition-dependent vapour-liquid equilibrium, i.e. Raoult's law. When dropping terms of quadratic order in y, the convection-diffusion equation (2.5) within the droplet becomes and the corresponding interface boundary condition (2.6) reads −∇y · n = Ev yj0 + Ev vapjy − Ev tot yj 0 (4.8) with the non-dimensional evaporation numbers (4.11) The number Ev y quantifies the intensity of the concentration gradient induced in the liquid by preferential evaporation of one of the components, i.e. it compares the differences of the two diffusive vapour transports in the gas phase with the mutual diffusion in the binary liquid. Since the resulting composition gradient along the interface and in the bulk is the driving mechanism for Marangoni flow and natural convection, this number will become important to quantify these processes later on. Note that dependent on the volatilities of the components and their mass fractions in the liquid, Ev y may be positive or negative. Ev vap is an estimate for the influence of local variations in the liquid concentration on the preferential evaporation, i.e. the linear feedback due to the quasi-stationary diffusion in the gas phase. If the composition is rather uniform in the droplet, which is usually by fast recirculating convection, the term Ev vapjy provides only a minor contribution in (4.8), meaning that the profile of the evaporation rates is similar to the one of a pure droplet. Since ∂ yA c eq A > 0 and ∂ yA c eq B < 0, i.e. the vapour concentration of A increases and B decreases for an increasing fraction of A in the liquid, Ev vap is always positive. Large numbers of Ev vap can actually arise towards the end of the drying of a glycerol-water droplet, as discussed later on in section 6.
Finally, Ev tot is a measure for the total evaporation speed, i.e. for the typical interface speedũ I and the volume evolution. Note that the total evaporation speed and volume evolution is a measure for the flow towards a pinned contact line, i.e. the flow leading to the coffee-stain effect. If none of the components condensates, i.e. both either evaporate or are non-volatile, the modulus of Ev y is smaller than Ev tot . Nevertheless, since the deviation from the average composition is small, i.e. y 1, and the Marangoni convection and/or natural convection are sufficiently large, the contribution of the latter to the flow can still be dominant compared to Ev tot yj 0 in (4.8).

Nondimensionalized flow
For the Navier-Stokes equations, we employ the established Boussinesq approximation, which is valid as long as y∂ yA ρ is small compared to ρ 0 (Gray & Giorgini 1976). Due to the usually small composition gradients, this assumption is valid here. Therefore, except for the body force term ρge z , only the zeroth order terms proportional to ρ 0 are kept, whereas y∂ yA ρ-terms are disregarded. With the same argument, terms proportional to y∂ yA µ and y∂ yA σ can be disregarded, whenever there is a dominant term proportional to µ 0 and σ 0 , respectively. Following this argument, the Navier-Stokes equations can be written as Here, the shifted nondimensionalized pressure, the Schmidt number and the incomplete Rayleigh number read (4.14) The Schmidt number for liquids is usually Sc > 10 3 which suggests that the lhs of (4.12) can be disregarded. However, since the chosen velocity and time scale in (4.2) does not necessarily coincide with the actual present scales, this argument is only valid for small Reynolds numbers. In small droplets with rather low volatilities and regular Marangoni flow, however, this assumption is surely met, e.g. Re < 0.05 in the case of the simulation in figure 2. The incomplete Rayleigh number Ra * deviates from the conventional definition of the Rayleigh number just by the lack of an estimate for the composition difference, i.e. a term like ∆y. The dynamic boundary conditions at the interface, (2.10) and (2.11), read in the Boussinesq approximation (4.16) Here, the non-dimensional number Ca * , the Bond number and the incomplete Marangoni number read Note that the definition of Ca * does not coincide with the capillary number, i.e. it does not consider the actually present typical velocity scale, i.e. the intensity of the capillary shape relaxations during evaporation cannot be inferred from Ca * . However, both the real capillary number Ca = µ 0 U/σ 0 and Ca * are small in the systems considered here (Ca < 1 × 10 −6 and Ca * < 1 × 10 −7 in the simulation depicted in figure 2). Similar to Ra * , the incomplete Marangoni number Ma * lacks in an estimate for the composition difference, i.e. ∆y, as compared to the conventional definition. Finally, the kinematic boundary condition (2.9) becomes where is the analogue of Ev vap for the total evaporation rate, i.e. the effect of a change in the saturation pressure due to a locally deviating composition on the total evaporation rate.

Estimation of the outwards flow
Before focusing on natural convection and Marangoni flow in the droplet, it is beneficial to obtain an estimate for the velocity in the droplet in absence of these mechanisms, i.e. Ma * = Ra * = 0. This case, exemplified e.g. by a pure isothermal droplet, combined with a pinned contact line represents the purely capillary-driven outward flow, which causes the coffee-stain effect.
For small droplets, the capillary number Ca is small, so that the surface tension forces according to (4.15) lead to an intense relaxing traction, whenever the droplet deviates from the equilibrium shape. Since the Bond number Bo is small as well for small droplets, the hydrostatic term in (4.15) can be neglected, leading to a spherical cap with a homogeneous curvatureκ as equilibrium shape. Hence, the shape evolution and thereby the interface velocityũ I is solely given by the evaporation rate and the contact line kinetics, which is assumed to be pinned here. Since the term Ev tot,vapjy in (4.18) is proportional to y, it can be disregarded with respect to Ev totj0 in accordance with the Boussinesq approximation. As a consequence, one ends up with a linear Stokes flow problem, where the entire bulk velocity is given by the instantaneous shape relaxation, which is proportional to the rate of evaporation, i.e. to Ev tot .
By integrating the evaporation rate Ev totj0 one obtains the volume loss and thereby one can reconstruct the normal velocity of the interfaceũ I · n. The flow in the bulkũ is subsequently given by solving the Stokes flow with the normal boundary conditioñ u · n =ũ I · n + Ev totj0 . In figure 3, representative solutions for the bulk flowũ with unity evaporation number, i.e. Ev tot = 1, are depicted. It is apparent that the typical bulk velocity is of the order unity, i.e. ũ ∼ 1. Since the flow is proportional to Ev tot , the typical velocity corresponding to an arbitrary evaporation number Ev tot is hence ũ ∼ Ev tot . This holds also for the typical interface movement, i.e. ũ I ∼ Ev tot .

Quasi-stationary limit
Knowing the fact that the capillary flow due to the volume loss is on the order of Ev tot , we now focus on the contributions to the flow by Marangoni forces and natural convection. In a first step, one can consider the case where Ev tot = 0, i.e. no total mass transfer and hence a constant volume and shape of the droplet. This scenario can be realized by tuning the ambient humidities of A and B so that the evaporative mass loss of component A is balanced by the condensation of component B. In this case Ev tot = 0 and Ev y > 0 holds. Again, due to the small capillary number and the small Bond number, one can assume a spherical cap shape with volumeṼ = 1 and contact angle θ, which are both constant now. Furthermore, there is no interface movement,ũ I · n = 0, and no total mass transfer,ũ · n = 0.
By averaging (4.7) over the droplet volumeṼ = 1, defining the integrated evaporation rateJ = j 0 dÃ (4.20) and considering only the zeroth order term in the boundary condition (4.8) in accordance with the Boussinesq approximation, one can separate the average composition y A,0 and the deviation y as follows: ∂ty +ũ ·∇y =∇ 2 y + Ev yJ0 . Here, the term Ev yJ0 assures that y A,0 is indeed the average composition and that the average of y remains zero, i.e. the term compensates for the imposed composition gradient at the liquid-gas interface. As already stated in section 4.1, this splitting holds only for limited time, since a variation in y A,0 leads to a change in the liquid properties which where used for the nondimensionalization. Usually, however, the coupled dynamics of flowũ and compositional differences y due to Marangoni and natural convection is considerable faster than Ev yJ0 , This was already apparent from the simulations depicted in figure 2 and it will be validated later on in section 6. Furthermore, this observation allows to focus on stationary solutions. Finally, upon introducing ξ = y Ev y , (4.23) one ends up at the following set of coupled equations: u ·∇ξ =∇ 2 ξ +J 0 (4.24) −∇p +∇ · ∇ũ + (∇ũ) t + Ra ξe z = 0 (4.25) ∇ ·ũ = 0 (4.26) subject to the following boundary conditions −∇ξ · n =j 0 (4.27) u · n = 0 (4.28) n · ∇ũ + (∇ũ) t · t = Ma∇ t ξ · t (4.29) at the liquid-gas interface and∇ ξ · n = 0 (4.30) u = 0 (4.31) at the liquid-substrate interface. Note that the simplified kinematic boundary condition (4.28) is now compatible with the no-slip boundary condition (4.31) at the contact line, i.e. a slip length is not required. Besides the contact angle θ, only two parameters enter the system, namely the Marangoni number and the Rayleigh number, which read (4.32) Note that the characteristic numbers for both mechanisms are proportional to the induced composition gradient due to mass transfer, i.e. Ev y . Of course, in particular the tangential gradient along the interface is also strongly dependent on the contact angle θ, since this determines the profile of the evaporation ratej 0 . These equations are not only valid for the specific assumed case of Ev tot = 0, but also when the combination of Marangoni and Rayleigh flowũ predicted by this model is considerably faster than the capillary flow, i.e. a flow situation with outwards flow leading to the coffee-stain effect. According to the estimations in section 4.3, this is the case ifũ Ev tot .

Procedure
Unfortunately, the analytical treatment of the model equations (4.24)-(4.26) is hampered by the geometry, which demands rather complicated toroidal coordinates, and by the very strong nonlinear coupling ofũ and ξ due to the advection term. Therefore we investigate the system by numerical means. Our analysis is limited to axisymmetric solutions and we only consider the case Ev vap = 0, i.e. neglecting the feedback of the altered gas composition due to the liquid-vapour equilibrium on the local evaporation rate. Finally, we will focus on the case Ma 0, for which evaporation leads to an overall reduction of the surface tension, as in the case of the water-glycerol depicted in figure 2. This results in a regular flow, i.e. no chaotic behaviour can be found, at least not for moderate flow conditions. In the case of negative Marangoni numbers, chaotic flow patterns cannot be excluded due to the Marangoni instability (Christy et al. 2011;Bennacer & Sefiane 2014;Machrafi et al. 2010). Of course, this spatio-temporal evolving type of flow cannot be captured within the assumption of a quasi-stationary process.
One can, on the other hand, test the linear stability of the quasi-stationary solutions in the case of negative Marangoni numbers to find the transition to chaotic flow, but since also the axial symmetry is usually broken in case of negative Marangoni numbers, one also would have to generalize the entire solution procedure from axisymmetric cylindrical coordinates to the full three-dimensional problem, as done by Diddens et al. (2017b). In order to find solutions of the system, we employed a finite element method on an axisymmetric mesh with triangular elements. We used linear basis functions for ξ and p and quadratic basis functions forũ, i.e. typical Taylor-Hood elements. The equations have been implemented in both FEniCS † (Logg et al. 2012) and oomph-lib for mutual validation. The condition of zero velocity in normal direction, i.e. Eq. (4.28), has been implemented by Lagrange multipliers. For an enhanced stability in the Newton method during the solution process, it has been found beneficial to replaceJ 0 in (4.24) by a Lagrange multiplier which ensures that the average of ξ is zero. This removes the null space with respect to a constant shift in ξ and a corresponding adjustment of the pressurẽ p.
Due to the nonlinear advection term, it is in general possible that multiple solutions exist for a given parameter combination (θ, Ma, Ra). For the parameter ranges considered in the following, however, we are confident that we found the generic solutions due to the following strategy: For every considered contact angle θ, we performed adiabatic scans along Ra in increasing and decreasing direction for fixed Ma and vice versa. During that, no hysteresis, i.e. bistable regions, have been found. Furthermore, by tracing these parameter paths with continuation, we have not detected any unstable branches. This has been furthermore validated by investigating the eigenvalues with a shift-inverted Arnoldi method ‡. Finally, for each parameter combination, we performed temporal integrations of the unsteady model equations starting from a homogeneous state ξ = 0. Since these runs converged to the same solutions as obtained by the steady parameter scans, we are sure that all solutions discussed in the following are indeed generic and stable. Note, however, that this is in general not true outside the considered parameter ranges.

Phase diagrams
The phase diagrams for small and large contact angles, i.e. for θ = 60 • and θ = 120 • , are depicted in figures figure 4 and figure 5, respectively. Here we have assumed, as in the case of the glycerol-water droplet, that the blue liquid A (e.g. water) is more volatile, less dense but associated with a higher surface tension than the red liquid B (e.g. glycerol). This means by definition that Ma > 0 holds and that a sessile droplet is described by Ra > 0 whereas a pendant droplet is given by Ra < 0 .
Depending on the Marangoni number Ma, the Rayleigh number Ra and the contact angle θ, different qualitative flow scenarios can be found. For high Marangoni numbers and small Rayleigh numbers, the Marangoni flow dominates (Ma dominant) and vice versa (Ra dominant). In between, however, for sessile droplets with a contact angle below 90 • and for pendant droplets with a contact angle above 90 • , there is a region † https://fenicsproject.org/ ‡ using Spectra https://spectralib.org In between these regions, however, there is a regime where the bulk flow is driven by natural convection whereas the flow close to the interface is dominated by the Marangoni effect (Ma vs. Ra). Here, two counter-rotating vortices can be seen. For pendant droplets (Ra < 0), both mechanisms driving the flow in the same direction (Ma & Ra same dir.). Hence, one cannot identify the main driving mechanism from the direction of the flow, so that the streamlines are coloured grey. If both mechanisms are sufficiently strong, however, the bulk flow due to natural convection can become so intense that the composition gradient along the interface changes direction and a flow reversal due to the Marangoni effect can arise in the vicinity of the interface (Ra reverses Ma).
where the Marangoni effect determines the flow direction at the interface, whereas the bulk flow is driven by natural convection (Ma vs. Ra). In the opposite cases, i.e. for pendant droplets with θ < 90 • and sessile droplet with θ > 90 • , both mechanisms drive a flow in the same direction, so that one cannot directly distinguish between the two mechanisms driving the flow (Ma & Ra same dir.). In the limit of very strong driving of both mechanisms, however, natural convection can become so intense, that the surface tension gradient is reversed, leading to a Marangoni-induced reversal of the flow at the interface (Ra reverses Ma). This effect can be explained by the distortion of the internal composition field due to natural convection. For pendant droplets with θ < 90 • , the composition gradient in the bulk in normal direction is much more pronounced near the contact line as opposed to the apex. As a consequence, the diffusive replenishment of the blue liquid at the interface is enhanced near the contact line so that in fact more blue liquid, i.e. the one with higher surface tension, can be found near the contact line instead of at the apex -despite of its higher volatility and the pronounced evaporation rate at the contact line. The resulting Marangoni flow is therefore reversed as anticipated by considering the profile of the evaporation rate alone. The same explanation holds for the case θ > 90 • , except that one finds a more pronounced normal composition gradient near the apex as compared to the region close to the contact line, and the situation is reversed.
All transitions between the afore-mentioned regimes are continuous. The drawn phase boundaries are defined by the emergence or disappearance of a second vortex. There is no bifurcation and/or hysteresis present at the boundaries of the regimes. In supplementary movie 2, a path through the parameter space is traversed and the corresponding stationary solution is shown, which illustrates the behaviour of the flow upon crossing the phase region boundaries, i.e. how the stationary solution gradually changes between single and two-vortex solutions.
Finally, we also investigate the contact angle dependence of the phase diagrams by showing the corresponding regions for θ = 40 • , 60 • and 80 • in figure 6(a) and for θ = 100 • , 120 • and 140 • in figure 6(b). Obviously, the phase boundaries are shifted, but qualitative differences in the phase diagrams cannot be found.
Note again that we have assumed in the phase diagrams that the more volatile liquid (blue) is less dense in the insets in the phase diagrams. In the other case, the droplets depicted in the insets are required to be mirrored vertically, as the Rayleigh number is then negative for sessile droplets and positive for pendant droplets. Furthermore, it is noteworthy that these diagrams are for pinned droplets (CR-mode) and droplets with a constant contact angle (CA-mode), as only stationary solutions are considered anyhow. As long as the dominant velocity contribution is given by recirculating flow due to Marangoni and/or natural convection, any capillary flow due to shape relaxations can be disregarded. Finally, the diagrams can also be used for condensation instead of evaporation, as long as the Ma > 0 holds. For condensation of component A, Ev y < 0 holds, so that Ma > 0 is true if component A has the lower surface tension, i.e. ∂ yA σ < 0. We therefore anticipate that the diagrams can predict the flow when ethanol condenses on a pure water droplet, whereas it would fail to predict the flow when water condenses on a pure glycerol droplet (Ma < 0). In fact, the latter case has been investigated experimentally, showing indeed chaotic cellular flow structures (Shin et al. 2016).

Validation of the quasi-stationary approximation against the full numerical simulation
Since there were a number of assumptions made in the simplification of the problem, it is necessary to validate the predicted flow by comparing it with results of the full numerical simulation, i.e. with the full set of equations as described in section 2. We focus on the representative simulation depicted in figure 2. At each instant in time, we have extracted the spatially averaged water mass fraction y A,0 , the volume V to determine the spatial scale 3 √ V and the contact angle θ from the simulation. From y A,0 , we obtain ρ 0 and ∂ yA ρ, σ 0 and ∂ yA σ, as well as µ 0 , D 0 and c eq A,0 from the composition-dependent properties of binary glycerol-water mixtures. This allows to calculate the normalized evaporation-induced composition gradient Ev y and the characteristic numbers Ra and Ma. On the basis of these numbers and the contact angle, we solve the simplified quasistationary model and re-dimensionalise the resulting velocity and composition field as   well as the evaporation rate using the scales (4.2). This procedure allows to compare the full unsteady evolution of the droplet with the corresponding predictions at each instant by the simplified quasi-stationary model. The results are depicted for several instants in figure 7, where the full simulation is shown on the left and the corresponding prediction of the quasi-stationary model is depicted on the right. Initially, i.e. in figure 7(a), the full simulation has not attained the quasi-stationary limit. Hence, the quasi-stationary model slightly overpredicts the composition variations, i.e. it shows more glycerol (red) near the interface and more water (blue) in the bulk. Therefore, the flow field also slightly differs, i.e. the transient full simulation shows a single vortex, whereas the quasi-stationary model predicts the presence of a small counter-rotating vortex near the apex. Furthermore, a very gentle deviation in the spherical cap shape due to the gravitational effect in the full simulation can be seen at the apex as well (regime Ra reverses Ma). At later time steps, i.e. in figure 7(b-d), however, the flow and the composition predicted by the quasi-stationary model match the results of the full simulation almost perfectly, be it in terms of the composition field, the flow pattern, the shape or the evaporation rate. This result substantiates the fact that the capillary outwards flow, which has been disregarded in the quasi-stationary model, can indeed be neglected as long as there is a prominent recirculating flow due to Marangoni and/or Rayleigh convection. To assess the quality of the quasi-stationary model in more detail, we have extracted some characteristic quantities of both simulations, i.e. from the full simulation of figure 2 and the corresponding quasi-stationary limit at each instant. In figure 8(a), the time evolution of the three key parameters, namely the Rayleigh number Ra, the Marangoni number Ma and the contact angle θ is shown. These numbers were used as input for the quasi-stationary model. The rms of the velocity inside the droplet is depicted in figure 8(b). Again one can see an initial disagreement due to the fact that the full simulation has not yet attained its quasi-stationary limit. After that, i.e. after about 50 s to 100 s, the rms-velocity is well predicted until it shows again a disagreement towards the end of the drying time. The reason for the overpredicted velocity in the quasi-stationary model can be found in figure 8(c), where the minimum and maximum glycerol concentration in the droplet according to both simulations are plotted against time. While it shows good agreement in the main part of the drying, the quasi-stationary model shows an enhanced maximum glycerol concentration towards the end of the drying, i.e. when almost only glycerol is left in the droplet. In fact, the glycerol concentration predicted by the quasi-stationary model even exceeds the physically realistic threshold of 100 %. Obviously, this overprediction of the composition differences explains the elevated prediction of the rms velocity. The reason of the overpredicted composition difference can finally be seen in figure 8(d), where the evaporation numbers are depicted. When the droplet almost entirely consists of glycerol, the evaporation number Ev vap , quantifying the reduction of the water vapour pressure for vanishing water at the interface on the evaporation dynamics (cf. Eq. (4.10)), becomes very large (Ev vap → 18 at t = 1000 s). This effect is not considered in the quasi-stationary model, since it assumes the averaged composition y A,0 to predict the vapour-liquid equilibrium, not the local composition at the interface. Thereby, the amount of water vapour is strongly overestimated which results in a high evaporation rate and thereby in an unrealistically high composition difference. Obviously, the quasi-stationary model loses validity when Ev vap becomes too large, meaning that the dependence of the vapour-liquid equilibrium on the local interface composition cannot be neglected any more. For a more detailed model, this effect can easily be incorporated into the quasi-stationary model, but it would introduce a fourth parameter besides Ma, Ra and θ into the set of equations, which is beyond the scope of this article.

Conclusion
During the evaporation of a binary droplet, multiple flow scenarios can be found, which is a result of an interplay of differences in the volatilities, mass densities and surface tensions of the two constituents. The difference in the volatilities induces compositional gradients in the bulk and also, due to the in general non-homogeneous evaporation rate, along the interface. Due to the composition-dependent mass density and surface tension, natural convection and Marangoni flow can set in, leading to a recirculating flow in the droplet that is usually much faster than the typical capillary outwards flow towards the contact line, which can be seen in pure droplets and leads to the coffee-stain effect in particle-laden droplets.
Based on justified assumptions, we simplified the full model equations to a quasistationary model that only requires three parameters, namely the contact angle, the Rayleigh and the Marangoni number. Both, the Rayleigh and Marangoni number, linearly scale with a non-dimensional evaporation number, Ev y , which is a measure for the induced composition gradient by the preferential evaporation of one or the other component. By numerically solving for stationary solutions of the simplified model, we have explored the phase space in terms of these three quantities. The obtained phase diagrams allow for the prediction of the flow types in sessile and pendant binary droplets, with contact angles below and above 90 • .
We found in total five different flow patterns: If one of both mechanisms, i.e. either natural convection or Marangoni flow, gets sufficiently strong as compared to the other one, it can dominate and control the flow direction in the entire droplet. This scenario can usually be seen in the case when the corresponding number, namely the Rayleigh or Marangoni number, respectively, is much larger than the other. In these cases, a single vortex can be seen in the droplet. If both mechanisms drive the flow into a different direction and are comparably strong in terms of their non-dimensional numbers, one can find two vortices, one in the bulk driven by natural convection, and a counterrotating vortex at the interface due to Marangoni flow. The fourth flow type is the case, when both mechanisms act in the same direction, so that one cannot distinguish the particular cause of the driving and only a single vortex is present. Remarkably, however, in particular regimes in the phase space, the Marangoni flow can be reversed due to the natural convection in the bulk, leading to the fifth solution, where again two vortices can be found. In this situation, the bulk flow driven by natural convection deforms the internal composition field so that the diffusion dynamics in the liquid are altered, which eventually reverses the composition gradient at the interface and hence the Marangoni flow.
To use the phase diagrams presented in this article, several requirements have to be fulfilled: First of all, the influence of thermal effects must be negligible compared to the solutal ones. The two liquids must be miscible and the droplet must not be too large, so that the capillary number and the Bond number are small in order to guarantee a spherical cap shape during the evaporation. Furthermore, the Reynolds number must be small and the spatial variations in the composition must be small enough in order to allow for a first order Taylor expansion of the composition-dependent liquid properties according to (4.1). Also the requirements for the Bousinessq approximation must hold. The Marangoni number, as defined in (4.32), must be positive, i.e. evaporation leads to an overall decrease of the surface tension, so that Marangoni-unstable chaotic flow can be excluded and the recirculating flow must be sufficiently faster than the movement of the interface. Finally, the influence of a change in the local composition on the vapourliquid equilibrium may not be too strong, as it has been discussed on the basis on the evaporation number Ev vap describing the feedback of local composition changes on the evaporation rate in section 6.
If all these requirements are fulfilled and the composition-dependence of the required physical properties are known, the phase diagrams of this article allow for a prediction of the qualitative flow pattern in an evaporating binary droplet, probably with exception of a short initial transient phase.
The method described in this article can be directly transferred to thermally driven Marangoni flow and natural convection in a pure droplet. Instead of a convection-diffusion equation for one component, one would have to consider the convection-diffusion equation for the temperature field. The boundary conditions will be different, e.g. a Dirichlet boundary condition of constant temperature at a highly conducting substrate and nondimensional evaporative cooling instead of the number Ev y , but the methodological principle can remain the same. Also a generalization to negative Marangoni numbers could be interesting, but it would require the consideration of the problem in three dimensions. This would allow to predict axial symmetry breaking and also bifurcations into chaotic Marangoni flow regimes by performing a linear stability analysis of the quasistationary solutions. with h being the height of the droplet and µ the averaged viscosity of both liquids, was used in our previous publication (Li et al. 2019a) as an indicator whether the flow in the droplet is dominated by natural convection (Gr 1) or not (Gr 1). Compared to the non-dimensional numbers presented in this manuscript, i.e. the evaporation number Ev y and the Rayleigh number Ra, this number is independent of the current droplet composition. Instead, it just takes the pure densities of both fluids and the averaged density and viscosity into account. This means that the Grashof number Gr is easily accessible, whereas the non-dimensional numbers used throughout this manuscript require the knowledge of the instantaneous average composition and the full compositiondependence of all properties, which is not always possible in an experimental setup. Therefore, we are interested to substantiate the argumentation by Li et al. (2019a) that the much simpler Grashof number Gr can be used as indicator whether to expect natural convection (Gr 1) or not (Gr 1). To investigate the validity, we replace the Rayleigh number by the Grashof number in the following. If one assumes that ∂ yA ρ is independent of the composition, i.e. a linear dependence of the mass density on the mass fractions, one can obtain the Grashof number via the relation Gr = 3 π 1 − cos θ 2 + cos θ Ra Ev y Sc , (7.2) where the factor depending on the contact angle θ is a consequence of the different characteristic length scales, i.e. 3 √ V for Ra and h for Gr. While the Schmidt number Sc = µ/(ρD) for liquids is typically Sc ∼ O(10 4 − 10 5 ), for moderately volatile liquids like e.g. water at typical ambient conditions, Ev y ∼ O(10 −1 − 10 0 ) holds. In order to obtain diagrams independent of these quantities, we set the factor Ev y Sc = 1000 in (7.2) for the determination of the boundaries in the phase diagrams.
The phase diagrams rescaled to the Grashof number via this way are depicted in figure 9. One can infer from these diagrams that even in competition with a strong Marangoni effect, the onset of gravity-driven bulk flow happens approximately at a Grashof number of Gr ∼ O(1) for a contact angle of θ = 70 • in (a) and θ = 100 • in (b). Furthermore, the experimental results of Li et al. (2019a) are indicated as dots. The 1,2-propanediol-water droplets with a contact angle of θ = 70 • discussed in Li et al. (2019a) clearly show the effect of natural convection for an apex height of h = 800 µm, whereas is was not visible for h = 410 µm. This clearly coincides with the prediction of the phase diagram in figure 9(a). The experiments on glycerol-water droplets with θ = 100 • , as discussed in the supplementary information of Li et al. (2019a), reveal an absence of observable natural convection for h = 154 µm, whereas the presence of natural convection was found at heights h 320 µm, with increasing velocity for elevated heights. Also this can be inferred from the Ma-Gr-diagram depicted in figure 9(b). Thus we conclude that the Grashof number Gr is indeed an indicator for the presence or absence of decisive natural convection in a binary droplet. The Ma-Ra-diagrams presented in this manuscript, however, provide a much more detailed prediction of the possible flow scenarios.