Taylor Dispersion in Thin Liquid Films of Volatile Mixtures: A Quantitative Model for Marangoni Contraction

The Marangoni contraction of sessile droplets occurs when a binary mixture of volatile liquids is placed on a high-energy surface. Although the surface is wetted completely by the mixture and its components, a quasi-stationary non-vanishing contact angle is observed. This seeming contradiction is caused by Marangoni flows that are driven by evaporative depletion of the volatile component near the edge of the droplet. Here we show that the composition of such droplets is governed by Taylor dispersion, a consequence of diffusion and strong internal shear flow. We demonstrate that Taylor dispersion naturally arises in a self-consistent long wave expansion for volatile liquid mixtures. Coupled to diffusion limited evaporation, this model quantitatively reproduces not only the apparent shape of Marangoni-contracted droplets, but also their internal flows.

The Marangoni contraction of sessile droplets occurs when a binary mixture of volatile liquids is placed on a high-energy surface. Although the surface is wetted completely by the mixture and its components, a quasi-stationary non-vanishing contact angle is observed. This seeming contradiction is caused by Marangoni flows that are driven by evaporative depletion of the volatile component near the edge of the droplet. Here we show that the composition of such droplets is governed by Taylor dispersion, a consequence of diffusion and strong internal shear flow. We demonstrate that Taylor dispersion naturally arises in a self-consistent long wave expansion for volatile liquid mixtures. Coupled to diffusion limited evaporation, this model quantitatively reproduces not only the apparent shape of Marangoni-contracted droplets, but also their internal flows.
Dimensional reduction in the limit of long waves is a powerful tool to analyze moving contact lines [23,33,[35][36][37][38]: The evolution of the local liquid height is derived from the net flux, treating the short (vertical) axis fully implicitly. For mixtures, however, vertical compositional gradients are unavoidably generated by shearing any horizontal gradients. This impedes the use of dimensional reduction, unless taking the limit of infinitely fast diffusion along the short axis. All existing lubrication models are formulated in this limit [39][40][41][42]. Commonly, however, the time scale of diffusion is finite and, in combination with shear flow, leads to strong dispersion. This so-called Taylor-Aris dispersion [43,44] has important consequences in many natural [45] and technological scenarios [46]. To date it remains unclear whether shear dispersion is consistent with a long-wave expansion, so no expression for the effective dispersion in general thin- film flows is available in the literature.
Here we show that small vertical compositional gradients are in agreement with the usual assumptions in a long-wave expansion. We provide a general expression for the effect of shear dispersion in thin liquid films, thus enabling lubrication theory to be consistently applied to bulk liquid mixtures. Our model is in quantitative agreement with experimental observations of Marangonicontracted droplets. We expect our analysis to be relevant well beyond droplet studies, as it offers a general route for implementing the effect of advected bulk fields in dimensional reduction problems.
Experiments.-We measured the apparent shape and the internal flows of Marangoni-contracted droplets inside an atmospheric control chamber (size ∼ 10 cm × 10 cm × 10 cm), at room temperature.
The droplets were composed of mixtures of water ("Milli-Q", resistivity 18 MΩ cm) and a carbon diol (Sigma Aldrich, ≥ 98%). Piranha cleaned microscopy glass coverslips (170 µm thick) were used as substrates (see Supplemental Material for details [47]). The humidity was set by continuously injecting a well-defined mixture of dry and moist nitrogen behind gas-permeable membranes at the side-walls of the chamber. The droplets had initial volumes of 0.5 to 1 µL. Micro particle image velocimetry (µPIV, Fig. 1 (b)) was performed with an inverted fluorescence microscope and a high-aperture water-immersion objective (20x, numerical aperture (NA) 0.95) to allow for diffraction limited imaging in the bulk droplet. Polystyrene microspheres (Thermo Fisher Scientific F8809, diameter 200 nm) were used as flow tracers, with a mass fraction of 7.8 × 10 −5 of the particle stock solution in the binary mixture. Images of the particles were captured with a high-speed camera at 600 to 1000 FPS, quickly switching between z-planes by automating the focus system of the microscope. Simultaneous side-view imaging of the drop was performed with a telecentric lens ( Fig. 1 (c)). Fig. 1 (d) shows the evolution of the apparent contact angle θ app of spreading drops of pure water (diol mass fraction φ = 0), pure 1,2-propanediol (PD, φ = 1) and of their mixtures. Pure liquids spread into complete wetting. In contrast, the binary mixtures reach a stationary non-equilibrium apparent contact angle θ app > 0 shortly after deposition. The drops stay in this contracted state for several minutes. This wetting behavior has been described previously [12][13][14], showing that θ eq depends on φ and the ambient relative humidity RH. Fig. 2 (a) shows the velocity field inside the droplet, as determined by µPIV. The arrows indicate the velocities that have been measured in different z-planes. The insets show dense velocity fields for two exemplary zplanes. Close to the free surface, the flow is directed inward, precisely balanced by an outward flow close to the substrate, leading to a quasi-stationary shape.
The surface tension gradient can be derived from the tangential stress boundary condition, using the measured shear rate and the viscosity η(φ) from the literature [50][51][52]. Fig. 3 shows the experimentally derived surface tension gradient as a function of the distance d to the contact line, for various diols, compositions, and relative humidities. All curves follow a power law ∂ r γ ∼ d −3 /2 . Despite significant differences in composition, surface activity, and ambient humidities, the curves nearly collapse in physical units.
Lubrication theory.-We consider the general case of a thin liquid film of a mixture on a flat solid surface. The film covers the entire substrate, with a continuous transition between the macroscopic droplet and a microscopically thin precursor surrounding it. The latter reflects the adsorption equilibrium of vapor molecules in the atmosphere around the droplet [53]. For water above the dew point, adsorption layers on hydrophilic surfaces are typically on the order of a few nm [54,55]. The free surface is described by h( r), where r is the location in the substrate plane (see Fig. 1 (a)). Incompressible Stokes flow without body forces is governed by where u is the fluid velocity, η is the dynamic viscosity of the fluid, and p is the fluid pressure. The evolution of the solute field φ is given by with t as time and D as the diffusion coefficient of the solute. For simplicity we limit the following analysis to isothermal, isochoric, and isoviscous cases. We consider a no-slip and no-flux boundary condition at z = 0, kinematic and stress boundary conditions at z = h, and a Stefan-type boundary condition that links composition and evaporation (see Supplemental Material for details [47]).
To derive evolution equations in terms of vertically averaged quantities, we take the limit of long waves, where the characteristic horizontal scale, r 0 , shall be much larger than the characteristic vertical scale, h 0 : h = h 0 /r 0 1. For sessile droplets, r 0 and h 0 are conveniently chosen as the footprint radius and the apical height of the droplet, respectively. We define the vertically averaged velocity u, the total hydrodynamic flux Φ, the vertically averaged composition φ, and the effective solute height Ψ through and the deviations from the average by We scale all horizontal coordinates as r = r 0 r , and all vertical coordinates as z = h r 0 z . Velocities and time are scaled with u 0 and t 0 , the characteristic velocity and time scales of the problem, respectively. Below, u 0 and t 0 will be identified with the natural scales that arise in the evolution equations. Surface tension is scaled as γ = γ 0 γ , where γ 0 = γ(φ = 0), the surface tension of the pure solvent. Diffusivity is treated similarly, D = D 0 D with D 0 = D(φ = 0). We scale j = j 0 j for the evaporation rate, where j and j 0 are determined according to Ref. [53] for the diffusion limited regime. Whether φ, the vertically averaged composition, is sufficient to describe the compositional evolution of the film, depends on the magnitude of the residual field δφ. Thus we scale δφ = φ δφ , deriving φ from the governing equation for δφ. In the following we will omit the primes for readability and work exclusively with scaled quantities.
The derivation of the evolution equation for the film height follows the standard procedure described in the reviews [23,35]. One obtains from integrating 2 along z, where E = j 0 η/( 4 h γ 0 ). Note that 6 does not involve any approximation. The long wave expansion is used only in the expressions for the evaporation rate j (see Ref. [53]) and the horizontal hydrodynamic flux [56] Here we identified the natural velocity scale u 0 = 3 h γ 0 /η, the capillary velocity for thin films, to cancel the material properties from Eq. 7. For a typical 1,2-propanediol/water droplet with φ ∼ 0.2 (i.e., γ ∼ 56 mN/m and η ∼ 2 mPa s [57]), r 0 ∼ 1.5 mm, and h 0 ∼ 0.15 mm one obtains u 0 ∼ 28 mm/s, much larger than the maximum experimentally observed velocities ∼ 2 mm/s. This is a common observation in wetting problems, where the capillary number typically remains small. The natural pressure-and time scales are p 0 = h γ 0 /r 0 and t 0 = r 0 η/( 3 h γ 0 ), respectively. The pressure p = −γ ∇ 2 h + Π(h) contains capillary and surface (disjoining) forces that stabilize the precursor film.
The capillary-( Φ C ) and Marangoni ( Φ M ) fluxes are associated with Poiseuille-and Couette-type velocity profiles, respectively [47]. These will shear any horizontal compositional gradient, such that a vertical gradient arises naturally. Combined with molecular diffusion from Eq. 3, this leads to Taylor-Aris dispersion [43,44]. In addition, the Stefan boundary condition for evaporation at the free surface requires a vertical compositional gradient [58,59]. It is commonly accepted that vertical compositional gradients are beyond the limit of the lubrication expansion [37,38,56,[60][61][62][63][64]. We challenge this paradigm, identifying three regimes, depending on aspect ratio and Péclet number: i. A regime of faint vertical compositional gradients where previous long-wave models hold [56]; ii. An intermediate regime of small but not negligible vertical gradients for which we derive a previously unknown evolution equation for φ including Taylor-Aris dispersion; and iii. A regime of large vertical gradients where the full advection-diffusion problem has to be solved [56].

Inserting 5 into 3 and integrating over the film height gives
where Pe = u 0 r 0 /D 0 is the Péclet number and D is the vertically averaged diffusion coefficient (see Supplemental Material [47] for a detailed derivation). In contrast to Eq. 6, terms with non-averaged quantities remain. These terms scale as ∼ φ , while the next-order terms in 6 with 7 scale as ∼ 2 h . Thus, the magnitude of φ relative to h determines which terms in 8 should be kept.
Limit i: φ δφ 2 h . We may ignore all terms ∼ φ and recover the previously known evolution equation [56,65,66]: Limit iii: φ δφ ∼ 1. No simplified evolution equation for vertically averaged quantities can be derived, and the full problem has to be solved [56,66].
Limit ii: φ δφ ∼ h . In this case, terms up to ∼ φ have to be retained. We consider a convection dominated problem i.e., Pe 1 and E 1. Far below the boiling point, and with typical D 0 ∼ 10 −9 m 2 /s and u 0 ∼ 10 −3 m/s, this holds for most sessile nano-to microliter droplets with small contact angles. Eq. 8 simplifies further [47]: The remaining term with δφ scales as φ . Thus a governing equation for δφ can be truncated to O(1). Inserting 5 into 3, multiplying with h, and subtracting 10 gives the leading order governing equation for δφ [47]: where we recover the natural scale of the residual field, [65,66]. φ is equivalent to a Péclet number for the characteristic vertical length-and velocity scales, h r 0 and h u 0 , respectively. Eq. 11 defines the advection-diffusion problem of the residual field in the co-moving frame of the mean flow: Diffusion along z balances the shearing due to the horizontal flow, and the residual field is stationary at leading order.
The requirements for limit ii became apparent now: The problem must be convection dominated in the horizontal direction (Pe 1), but diffusion dominated in the vertical direction, i.e., the residual field must remain small: 2 h Pe δφ 1. This is similar to the classical treatment of pipe flow by Taylor and Aris [43,44], but here the velocity field and the film height, and thus δφ vary in space. By scaling z in Eq. 11 with the local film height h, it becomes apparent that δφ ∼ h 2 δ u ∇φ if shear (δ u) dominates. For Marangoni-contracted droplets, we find δφ 1 everywhere: ∇φ and δ u, which are caused by evaporative enrichment and Marangoni convection, are strong only near the edge of the droplet where the height is small (see below for a quantitative estimate).
Eq. 11 can be integrated, and the resulting expressions for δφ and the integral in 10 can be found in the Supplemental Material [47]. In cases of axial or translational symmetry and slow evaporation (E j/h δ u), which holds for our experiments, Eq. 10 reduces to Eq. 9 with D replaced by D eff to account for Taylor-Aris dispersion:

. (12)
Numerical simulations.-We implemented the evolution equations 6 and 9 with the flux 7 and the effective diffusivity 12 in an axisymmetric finite volume scheme with convergent numerical mobilities after Refs. [67,68], and diffusion limited evaporation according to Refs. [13,53]. We used accurate material properties γ(φ), η(φ), found in the literature [50][51][52]57], and assumed D ∼ η(φ) −1 in accordance with the Stokes-Einstein relation. Simulations were initiated with a droplet of ∼ 0.7 µL volume and ∼ 30 • apparent contact angle, on top of a precursor in equilibrium with the vapor field of the droplet. See [47] for details. Fig. 1 (d) compares the apparent contact angle observed in simulations with (dashed) and without (dashdotted) Taylor-Aris dispersion. A near-quantitative agreement is observed for stationary contraction only if dispersion is taken into account. The remaining quanti-tative deviation is much smaller than the mismatch for simulations without dispersion and can be attributed to uncertainties in the material parameters, most importantly, the molecular diffusivity. We deliberately refrain from any parameter fittings. Fig. 2 (b) shows stream lines and composition for a contracted droplet. The observed difference in composition between the center and edge is merely 0.4%. The most striking feature of the simulations is a quantitative reproduction of the experimentally observed velocities (Fig. 2 a, gray, simulation, vs. black, experiments) and surface tension gradients (Fig. 3) for the case with Taylor-Aris dispersion.
The origin and importance of Taylor-Aris dispersion are highlighted in Fig. 4. Panel (a) shows the strong but compensating capillary and Marangoni fluxes in the droplet. Although the net hydrodynamic flux vanishes, the different velocity profiles of capillary and Marangoni fluxes lead to a convection roll inside the droplet, and thus strong shear dispersion. This leads to an effective diffusivity (including shear dispersion and molecular diffusion, Fig. 4 (b)) which scales quadratically with the individual flux components (Eq. 12), meaning strong dispersion in regions where the fluxes are large. Here, the effective diffusivity reaches D eff ∼ 27D around r ∼ R/2. Thus, Taylor-Aris dispersion becomes the governing phenomenon for the solute distribution in Marangonicontracted droplets.
Inserting again our characteristic quantities (r 0 ∼ 1.5 mm, h 0 ∼ 0.15 mm, φ ∼ 0.2, γ ∼ 56 mN/m, and η ∼ 2 mPa s), we obtain Pe ∼ 4 · 10 4 and 2 h Pe ∼ 4 · 10 2 . Thus the applicability of Taylor dispersion for our droplets depends on the magnitude of the scaled δφ ∼ h 2 δu ∂ r φ. The largest measured velocities were ∼ 2 mm/s ∼ 0.1u 0 , close to the contact line where the scaled height h ∼ 0.1 (see Fig. 2). The gradient of the mean composition in that region can be estimated from ∂ r γ ∼ 0.05. Thus, the deviation from the mean composition is only 2 h Pe δφ ∼ 2 · 10 −2 , which justifies the approximations.
Conclusion.-We measured the internal flow fields of Marangoni-contracted drops and derived the surface tension gradient, which follows a power law ∼ d −3 /2 . Through a systematic long wave expansion for freesurface films of mixtures, we could extend the widely used thin-film evolution equations to the case of bulk mixtures subject to Taylor-Aris dispersion. For Marangonicontracted drops, Taylor-Aris dispersion governs the composition, and our model is in quantitative agreement with the experimental findings. The theoretical analysis enables lubrication theory to be used for the very general case of advection-dominated thin free-surface films with advected bulk fields like temperature or composition.
Note added. Recently, we became aware of another study that includes Taylor dispersion in the description of Marangoni-contracted droplets [69].