2.1. Underlying dynamical equations

We focus on the evaporation of a single or two neighbouring binary droplets (or rivulets) which form a contact angle $\theta$ with the substrate. Let us consider that: only one component $B$ in the liquid phase can evaporate into the vapour phase; the surface tension monotonically decreases with the concentration of the non-volatile component; and evaporation is slow. Under these assumptions, we avoid complexities arising from Marangoni instabilities (Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b , Reference Diddens, Dekker and Lohse2024; Rocha, Lohse & Diddens Reference Rocha, Lohse and Diddens2024) and we meet the exact conditions to apply the models introduced by Diddens, Li & Lohse (Reference Diddens, Li and Lohse2021), which we describe in the following.

At ambient conditions, the evaporation is dominated by the diffusion of vapour in the surrounding air. The solution of the Laplace equation

(2.1) \begin{align} {\nabla} ^2 c_B = 0, \end{align}

with the boundary conditions

(2.2) \begin{align} c_B = c^{\textit{eq}}_B = c^{\textit{pure}}_B \gamma _B x_B \quad &\text{at liquid}{-}\text{gas interface}, \end{align}

(2.3) \begin{align} c_B = c^\infty _B \quad &\text{at far field}, \end{align}

gives the vapour concentration $c_B$ in the gas phase. Here, $c^\infty _B$ is the vapour concentration at far field, dependent of the relative humidity, $c^{\textit{eq}}_B$ is the vapour–liquid equilibrium concentration for component $B$ , dependent on the liquid’s local composition at the interface, $c^{\textit{pure}}_B$ is the saturation vapour concentration of the pure liquid $B$ , $\gamma _B$ is the activity coefficient and $x_B$ is the molar mass of component $B$ .

The diffusive mass flux of component $B$ through the liquid–gas interface due to evaporation is given by

(2.4) \begin{align} j_B = -D^{\textit{vap}}_B \boldsymbol{\nabla } c_B \boldsymbol{\cdot } \boldsymbol{n} , \end{align}

where $D^{\textit{vap}}_B$ is the diffusion coefficient of vapour and $\boldsymbol{n}$ is the normal vector at the liquid–gas interface. If only one component evaporates, then the total evaporation rate is $j = j_B$ .

Within the liquid, the mass fraction $y_A$ is governed by the advection–diffusion equation

(2.5) \begin{align} \rho \left (\partial _t y_A + \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla } y_A\right ) = \boldsymbol{\nabla } \boldsymbol{\cdot }(\rho D \boldsymbol{\nabla } y_A). \end{align}

Navier–Stokes flow without body forces in the liquid phase is described by

(2.6) \begin{align} \partial _t \rho + \boldsymbol{\nabla } \boldsymbol{\cdot } (\rho \boldsymbol{u}) &= 0, \\[-32pt] \nonumber \end{align}

(2.7) \begin{align} \rho (\partial _t \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla } \boldsymbol{u})&= - \boldsymbol{\nabla } p + {\boldsymbol{\nabla }\boldsymbol{\cdot }} \left ( \mu {\boldsymbol{\nabla } \boldsymbol{u}}\right )\!, \\[0pt] \nonumber \end{align}

where $\mu$ is the viscosity of the liquid phase and $p$ is the pressure. We neglect gravity since the droplet is small and the contact angle is low. At the substrate, the no-slip boundary condition is applied, i.e. we consider a constant contact radius (Picknett & Bexon Reference Picknett and Bexon1977). At the liquid–gas interface, the kinematic boundary condition

(2.8) \begin{align} \rho (\boldsymbol{u} - \boldsymbol{u}_I) \boldsymbol{\cdot } \boldsymbol{n} = j \end{align}

and the dynamic boundary condition

(2.9) \begin{align} \boldsymbol{n} \boldsymbol{\cdot } \left ( \mu \left(\boldsymbol{\nabla } \boldsymbol{u} + \boldsymbol{\nabla } \boldsymbol{u}^T\right) - p \boldsymbol{I} \right ) = \boldsymbol{n} \sigma \kappa + \boldsymbol{\nabla }_t \sigma \end{align}

are applied, where $\boldsymbol{u}_I$ is the velocity of the liquid–gas interface, $\sigma$ is the surface tension, $\kappa$ is the curvature of the liquid–gas interface and $\boldsymbol{\nabla }_t$ is the surface gradient operator.

Naturally, the mass fraction $y_A$ of the non-volatile component $A$ in the liquid phase increases during evaporation, which is accounted for by the Robin boundary condition in the co-moving frame of the liquid–gas interface:

(2.10) \begin{align} \rho D \boldsymbol{\nabla } y_A \boldsymbol{\cdot } \boldsymbol{n} = y_A j , \end{align}

where $\rho$ is the mass density of the liquid phase and $D$ is the diffusivity in the liquid phase.

In this paper, we neglect the influence of temperature on the flow field, for simplicity. All properties are therefore only dependent on the mass fraction $y_A$ .

2.2. Numerical simulation results for rivulets

Due to the non-axisymmetric evaporation rate of neighbouring droplets, the resulting flow field is asymmetric. Accurately computing the flow in this scenario requires solving the full three-dimensional problem, which is computationally expensive. As a starting point, we solve the full two-dimensional problem for rivulets. In later sections, we introduce simplifications to develop a model for droplets. All differential operators introduced in the previous section are expressed in Cartesian coordinates for the rivulet geometry. Note that the solution for $c$ in the gas phase within a two-dimensional domain depends on the distance $L$ from the midpoint between the neighbouring rivulets to the semicircular far-field boundary (Sneddon Reference Sneddon1966; Masoud & Felske Reference Masoud and Felske2009; Schofield et al. Reference Schofield, Wray, Pritchard and Wilson2020). We fix the far-field boundary at a distance $L / R = 100$ , where $R$ is the contact radius of each rivulet.

Figure 1. Transient direct numerical simulations, evaluated after ${2500}\,\textrm {s}$ (a), ${5000}\,\textrm {s}$ (b), ${7500}\,\textrm {s}$ (c) and ${10\,000}\,\textrm {s}$ (d), of two neighbouring rivulets, with an initial contact angle of $\theta = {50^\circ }$ , a 1,2-hexanediol mass fraction of $y_B = {5}\,\%$ and a distance $b={2.1}\,\textrm {mm}$ between their centres, evaporating in an environment with relative humidity of ${70}{\,\%}$ . The plots are magnified vertically by $\times 2.5$ , in order to better visualise the flow inside the rivulet. On the left-hand side of each subfigure, the velocity magnitude is displayed with superimposed streamlines. On the right-hand side, the 1,2-hexanediol mass fraction is shown. Additionally, the gas phase illustrates the water vapour mass fraction, while the liquid–gas interface highlights the evaporation flux of water. The $x$ position where the two vortices in the rivulet meet (which we refer to as the stagnation point $x_0$ ) increases with time.

The system of equations described in § 2 is solved using a sharp-interface Lagrangian–Eulerian finite-element method implemented in the software package pyoomph (Diddens & Rocha Reference Diddens and Rocha2024), which is based on oomph-lib (Heil & Hazel Reference Heil and Hazel2006) and GiNaC (Bauer, Frink & Kreckel Reference Bauer, Frink and Kreckel2002). All computational domains are discretised using meshes composed of triangular elements. Linear basis functions are employed for pressure $p$ fields, whereas quadratic basis functions are used for the velocity $\boldsymbol{u}$ , vapour concentration $c$ and mass fraction $y_A$ fields. Time integration is performed using a second-order implicit backward differentiation formula.

At early stages of evaporation, the streamlines in the rivulet nearly form a single vortex (see figure 1). To quantify the asymmetry in the flow within the liquid, we define the stagnation point, $x_0$ , as the $x$ coordinate where the two vortices merge – identified by a zero in the tangential interfacial velocity ( $\boldsymbol{u}_t=\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{t}$ ). For flat rivulets, this zero-velocity condition is nearly uniform along the vertical direction, which justifies using $\boldsymbol{u}_t$ to determine $x_0$ . Alternatively, we could also find $x_0$ by identifying the point where the velocity close to the substrate crosses zero, but, at least for rivulets, this would yield virtually the same result, particularly for low contact angles. We choose $x_0$ as our measure of asymmetry because it directly reports how far the internal circulation is displaced from the droplet centreline. The flow controls the redistribution of advected particles and solute near the contact line, and therefore strongly influences final deposition patterns.

As evaporation progresses, the stagnation point shifts toward the centre of the rivulet, prompting an investigation into how $x_0$ depends on the contact angle, surface tension forces and the distance between the droplets. We begin by introducing a simplified quasi-stationary system of equations, allowing us to identify the relevant dimensionless control parameters.

3.1. Quasi-stationary model

We focus on the evaporation of droplets (or rivulets) with quasi-stationary evaporation, such as mixtures of water and 1,2-hexanediol. (We ignore the slight non-monotonicity of the surface tension as a function of the 1,2-hexanediol (Diddens, Dekker & Lohse Reference Diddens, Dekker and Lohse2024). Though it leads to extremely rich interesting behaviour (Diddens et al. Reference Diddens, Dekker and Lohse2024), here it distracts from the main point of this paper. So our binary liquid can be seen as a model liquid with material parameters inspired by 1,2-hexanediol–water mixtures.) In such cases, the internal flow within the liquid remains stable and approximately quasi-steady at each instant during the drying process (Diddens et al. Reference Diddens, Li and Lohse2021), provided that the solutal Marangoni flow induces sufficient mixing. Consequently, we assume small compositional deviations around the spatially averaged mass fraction $y_{A,0}$ , and express the local mass fraction of component $A$ as $y_A = y_{A,0} + y$ , where $y$ denotes the deviation from the mean. We also consider, for simplicity, that changes in temperature do not significantly affect the flow field when compared with the effects of composition.

All liquid properties are considered constant, evaluated at the spatially averaged composition and ambient temperature, except for the liquid–gas surface tension, which varies according to the first-order Taylor expansion: $\sigma (y_A) = \sigma (y_{A,0}) + y \, \partial _{y_A} \sigma$ . Since the average composition evolves slowly – at least for droplets or rivulets with contact angles $\theta \gtrsim {5^\circ }$ – this expansion can be performed at any given time during the drying process. We introduce the following non-dimensional scales:

(3.1) \begin{align} x = R \tilde {x}, \quad t = \frac {R^2}{D_0} \tilde {t}, \quad \boldsymbol{u} = \frac {D_0}{R} \tilde {\boldsymbol{u}}, \quad p = \frac {\mu _0 D_0}{R^2} \tilde {p}, \end{align}

where $D_0$ and $\mu _0$ are the diffusivity and viscosity of the liquid, respectively. These non-dimensional scales are, in principle, time-dependent, as the liquid volume $V$ and the average composition $y_{A,0}$ evolve. However, as shown by Diddens et al. (Reference Diddens, Li and Lohse2021), the dynamics can be considered quasi-stationary at each instant during evaporation, at least after the initial short transient and as long as Marangoni flow is strong enough.

The vapour concentration $c_B$ in the gas phase is then expressed as

(3.2) \begin{align} c_B = \left(c^{\textit{eq}}_{B,0} - c^\infty _B\right) \tilde {c} + c^\infty _B, \end{align}

where $c^{\textit{eq}}_{B,0}$ is the spatially averaged vapour–liquid equilibrium concentration evaluated at specific time of interest during evaporation. Note that we neglect the effect of the variations of the local composition on the vapour concentration. As argued by Diddens et al. (Reference Diddens, Li and Lohse2021), if the local composition at the interface is close to the average composition, which happens in droplets or rivulets with strong Marangoni-induced circulation, then the equilibrium vapour concentration at the interface is approximately constant. Furthermore, for the particular case of 1,2-hexanediol and water mixtures, the equilibrium vapour concentration remains roughly constant for 1,2-hexanediol mass fractions below 0.5 (Dekker et al. Reference Dekker, van der, Marjolein and Lohse2024), based on UNIFAC activity coefficients (Constantinescu & Gmehling Reference Constantinescu and Gmehling2016). The Laplace equation (2.1) is then given by

(3.3) \begin{align} \tilde {{\nabla} ^2} \tilde {c} = 0, \end{align}

with boundary conditions $\tilde {c} = 1$ at the liquid–gas interface and $\tilde {c} = 0$ at the far field. Accordingly, the evaporation flux $j$ can be expressed as

(3.4) \begin{align} j_B = \frac {D^{\textit{vap}}_B}{R} \left(c^{\textit{eq}}_{B,0} - c^\infty _B\right) \tilde {j}, \end{align}

where $\tilde {j} = -\tilde {\partial }_n \tilde {c}$ depends solely on the liquid geometry. Neglecting higher-order terms in $y$ , the advection–diffusion equation for the deviation $y$ becomes

(3.5) \begin{align} \partial _{\tilde {t}} y_{A,0} + \partial _{\tilde {t}} y + \tilde {\boldsymbol{u}} \boldsymbol{\cdot } \tilde {\boldsymbol{\nabla }} y = \tilde {{\nabla} ^2} y, \end{align}

with the boundary condition at the liquid–gas interface:

(3.6) \begin{align} - \tilde {\boldsymbol{\nabla }} y \boldsymbol{\cdot } \boldsymbol{n} = {\textit{Ev}} j - {\textit{Ev}}_{\textit{tot}} y j. \end{align}

The dimensionless evaporation numbers are defined as

(3.7) \begin{gather} {\textit{Ev}} = \frac {- y_{A,0} D^{\textit{vap}}_B \left(c^{\textit{eq}}_{B,0} - c^\infty _B\right)}{\rho _0 D_0}, \quad {\textit{Ev}}_{\textit{tot}} = \frac {D^{\textit{vap}}_B \left(c^{\textit{eq}}_{B,0} - c^\infty _B\right)}{\rho _0 D_0}. \end{gather}

The parameter ${\textit{Ev}}$ quantifies the strength of the concentration gradients induced by the selective evaporation of component $B$ and thereby the Marangoni stresses driving the flow. It plays a central role in determining the flow field.

The parameter ${\textit{Ev}}_{\textit{tot}}$ measures the total evaporation rate and is associated with flow towards the contact line in pinned liquids, i.e. the coffee-stain effect. This contribution is generally small under strong Marangoni convection and minimal compositional deviation. However, towards the end of the drying process, when the non-volatile component dominates, local compositional effects become more significant. For simplicity, following the well-established model of Diddens et al. (Reference Diddens, Li and Lohse2021), we assume ${\textit{Ev}}_{\textit{tot}} = 0$ . The kinematic boundary condition at the liquid–gas interface is then reduced to

(3.8) \begin{align} \boldsymbol{n} \boldsymbol{\cdot } \tilde {\boldsymbol{u}} = 0 . \end{align}

Introducing the scaled variable $\xi = {y}/{{\textit{Ev}}}$ , the quasi-stationary form of (3.5) becomes

(3.9) \begin{align} \tilde {\boldsymbol{u}} \boldsymbol{\cdot } \tilde {\boldsymbol{\nabla }} \xi = \tilde {{\nabla} ^2} \xi + \dfrac {\tilde {J}}{\tilde {V}}, \end{align}

where $\tilde {J}/\tilde {V}$ is the integrated evaporation flux $\tilde {J}$ over the liquid–gas interface, averaged by the non-dimensional volume $\tilde {V}$ (or cross-sectional area for a rivulet). The Robin boundary condition (3.6) at the liquid–gas interface becomes the Neumann condition

(3.10) \begin{align} - \tilde {\boldsymbol{\nabla }} \xi \boldsymbol{\cdot } \boldsymbol{n} = \tilde {j}. \end{align}

The dynamic boundary condition (2.9) at the interface is modified to

(3.11) \begin{align} \boldsymbol{n} \boldsymbol{\cdot } \big ( \tilde {\boldsymbol{\nabla }} \tilde {\boldsymbol{u}} + \tilde {\boldsymbol{\nabla }} \tilde {\boldsymbol{u}}^T \big ) = \boldsymbol{n} \frac {\tilde {\kappa }}{{Ca}} + {\textit{Ma}} \tilde {\boldsymbol{\nabla }}_t \xi , \end{align}

where the dimensionless capillary number ${Ca}$ and Marangoni number ${\textit{Ma}}$ are defined as

(3.12) \begin{align} {Ca} = \frac {\mu _0 D_0}{\sigma _0 R}, \quad {\textit{Ma}} = \frac {R \partial _{y_A} \sigma }{\mu _0 D_0} {\textit{Ev}}. \end{align}

For small droplets or rivulets, the capillary number is small. In that case, the liquid will remain with a spherical cap.

Since the flow’s characteristic Schmidt number ${Sc} = \mu _0/(D_0 \rho _0)$ is typically large and the Reynolds number small, we can neglect the inertial terms in the Navier–Stokes equation. We therefore consider incompressible Stokes flow in the liquid phase to be given by

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

(3.14) \begin{align} \tilde {\boldsymbol{\nabla }} \tilde {p} &= \tilde {{\nabla} ^2} \tilde {\boldsymbol{u}}. \end{align}

Naturally, for droplets, the system of (3.3), (3.8)–(3.14) (hereafter termed the quasi-stationary model) is still three-dimensional and therefore computationally expensive. However, in the limit of small contact angles, the lubrication approximation, introduced in the following subsection, can be employed to effectively reduce the problem to two dimensions.

3.2. Lubrication quasi-stationary model

In the lubrication quasi-stationary model (hereafter simply lubrication model), we redefine some non-dimensional scales to make use of the short vertical length scale (as compared with the horizontal one) in small-contact-angle liquid bodies:

(3.15) \begin{align} h = \theta R \tilde {h}, \quad \boldsymbol{u}_{\Vert } = \frac {D_0}{R} \tilde {\boldsymbol{u}}_{\Vert }, \quad u_z = \theta \frac {D_0}{R} \tilde {u}_z, \quad p = \frac {\mu _0 D_0}{\theta ^2 R^2} \tilde {p}, \end{align}

where the subscript $\Vert$ indicates the horizontal component of the vector. For brevity, in the following, we drop the horizontal component subscript in the velocity, in the gradient $\boldsymbol{\nabla }$ and divergence $\boldsymbol{\nabla }\boldsymbol{\cdot }$ operators.

Following the standard derivation of the thin-film approximation (Eggers & Pismen Reference Eggers and Pismen2010; Diddens et al. Reference Diddens, Kuerten, Van der Geld and Wijshoff2017a ) and assuming negligible coffee-stain flow ( ${\textit{Ev}}_{\textit{tot}} = 0$ ), we arrive at the following evolution equation for the liquid height $h$ :

(3.16) \begin{align} \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot } \tilde {\boldsymbol{Q}} = 0, \end{align}

where

(3.17) \begin{align} \tilde {\boldsymbol{Q}} = \int _0^{\tilde {h}} \tilde {\boldsymbol{u}}\, {\rm d}\tilde {z} = - \dfrac {\tilde {h}^3}{3} \tilde {\boldsymbol{\nabla }} \tilde {p} + {\textit{Ma}} \theta \dfrac {\tilde {h}^2}{2} \tilde {\boldsymbol{\nabla }} \bar {\xi }. \end{align}

We consider the limit of zero ${Ca}$ , i.e. the droplet (or rivulet) remains a spherical cap. This means that the Laplace pressure is sufficiently strong to prevent any notable deformations. In droplets (or rivulets) with small contact angle, strong Marangoni flow can lead to, depending on its direction, a shape contraction (Karpitschka et al. Reference Karpitschka, Liebig and Riegler2017) or a pancake-like shape (Pahlavan et al. Reference Pahlavan, Yang, Bain and Stone2021).

For the compositional evolution, we consider a regime where vertical gradients of composition are small, but not negligible. We consider the evolution equation for the vertically averaged mass fraction $\bar {\xi }$ which includes Taylor–Aris dispersion (Ramírez-Soto & Karpitschka Reference Ramírez-Soto and Karpitschka2022):

(3.18) \begin{align} \tilde {\boldsymbol{Q}} \boldsymbol{\cdot } \tilde {\boldsymbol{\nabla }} \bar {\xi } = \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot } \left(D_{\textit{eff}}\, \tilde {h} \tilde {\boldsymbol{\nabla }} \bar {\xi }\right) - \dfrac {\tilde {j}}{\theta } + \dfrac {\tilde {J}\tilde {h}}{\tilde {V}}, \end{align}

where the dimensionless effective diffusivity $D_{\textit{eff}}$ is given by

(3.19) \begin{align} D_{\textit{eff}} = 1 + \theta ^2 \left ( \dfrac {2 \tilde {\boldsymbol{Q}}_C^2}{105} + \dfrac {\tilde {\boldsymbol{Q}}_C \boldsymbol{\cdot } \tilde {\boldsymbol{Q}}_M}{20} + \dfrac {\tilde {\boldsymbol{Q}}_M^2}{30} \right )\!, \end{align}

with

(3.20) \begin{align} \tilde {\boldsymbol{Q}}_C = -\dfrac {\tilde {h}^3}{3} \tilde {\boldsymbol{\nabla }} \tilde {p}, \quad \tilde {\boldsymbol{Q}}_M = \theta {\textit{Ma}} \dfrac {\tilde {h}^2}{2} \tilde {\boldsymbol{\nabla }} \bar {\xi }. \end{align}

Our numerical tests suggest that the inclusion of the Taylor–Aris dispersion correction leads to more physically meaningful results and increases numerical stability. Within the lubrication model, the local evaporative flux cannot be computed directly because doing so would require explicitly modelling the vapour concentration $c_B$ in the gas phase. Instead, this flux must be imposed as a boundary condition at the liquid–gas interface.

4.1. Droplets

For flat pure droplets, the evaporation flux is given by (Popov Reference Popov2005)

(4.1) \begin{align} j_{B,0} = \dfrac {2 D^{\textit{vap}}_B \left(c^{\textit{eq}}_{B,0} - c^\infty _B\right)}{\pi R \sqrt {1 - r^2}}. \end{align}

For the sake of simplicity, we consider that the presence of a second component in the droplet does not significantly affect the evaporation flux of the first component. In the particular case of 1,2-hexanediol and water mixtures, the equilibrium vapour concentration remains roughly constant for 1,2-hexanediol mass fractions below 0.5, which justifies the above assumption. For increased accuracy, we consider the analytical integral solution (Popov Reference Popov2005) for the evaporation flux of the droplet, considering the effect of the contact angle.

When another droplet is present, the evaporation flux is modified due to a shielding effect. The evaporation rate $j_{B,1}$ of each neighbouring droplet can be written as (Wray et al. Reference Wray, Duffy and Wilson2020)

(4.2) \begin{align} j_{B,1} = j_{B,0} \left ( 1 - \dfrac {4 \sqrt {b^2 - 1}}{\left (1 + \dfrac {2}{\pi } \arcsin {\dfrac {1}{b}}\right )\left (2\pi (r^2 + b^2 - 2r b \cos {\phi })\right )} \right )\! , \end{align}

where $b$ is the distance between the centres of the two droplets and $\phi$ the azimuthal angle in the drop. The shielding factor is independent of $\theta$ since it was derived for flat drops.

4.2. Rivulets

Figure 2. Sketch of (a) isolated and (b) neighbouring evaporating rivulets in an infinite half-space of size $L$ . In this sketch we have chosen a small value for $L$ for illustration purposes. The top row shows the numerically calculated water vapour concentration in the entire domain. The green arrows indicate the local evaporative flux. The bottom row shows a sketch of the rivulet geometry. The origin of the $x$ coordinate is chosen at the centre of the rivulet and the right rivulet in the neighbouring case. Here $R$ is the radius, or half-width, of the rivulet. The distance between rivulets can be defined as the centre-to-centre distance $b$ or the edge-to-edge distance $d$ .

Figure 2 shows the rivulet geometry, the local evaporative flux and the water vapour concentration in the semicircular infinite half-space. To the best of the authors’ knowledge, the evaporation flux of a rivulet for an arbitrary contact angle has not been derived yet. Schofield et al. (Reference Schofield, Wray, Pritchard and Wilson2020) showed the non-existence of an analytical solution for the vapour concentration in the gas phase with an infinite half-space, due to the logarithmic nature of the solution of the two-dimensional Laplace equation. When considering a finite semi-elliptic domain with radii $L$ and points at $R$ (approximately a semicircular domain for large enough $L$ ), the evaporation flux of a flat isolated rivulet is given by

(4.3) \begin{align} j_{B,0} = \dfrac {D^{\textit{vap}}_B \left(c^{\textit{eq}}_{B,0} - c^\infty _B\right)}{ \mathrm{arcsinh}{\left (L/R\right )}}\frac {1}{\sqrt {R^2 - x^2}}, \end{align}

where $x$ is the distance from the centre of the rivulet to the point of interest. Note that $\mathrm{arcsinh}{ (L/R )} \approx \ln { (2L/R )}$ for $L \gg R$ . For two neighbouring rivulets the evaporative flux is given by

(4.4) \begin{align} j_{B,1} =\frac {D^{\textit{vap}}_B \left(c^{\textit{eq}}_{B,0} - c^\infty _B\right)}{2\,\mathrm{arcsinh}\left (L/\sqrt {2 b R}\right )}\frac {1}{\sqrt {R^2 - x^2}}\frac {b+2x}{\sqrt {(x+b+R)(x+b-R)}}. \end{align}

Here, $x$ is the distance from the centre of the right rivulet, as indicated in figure 2. Note that Schofield et al. (Reference Schofield, Wray, Pritchard and Wilson2020) use the coordinate $\eta$ , which is centred between both rivulets. The first term in (4.4) is the integral flux, the second term is the local flux for an isolated rivulet and the final term is the local shielding effect. We consider (4.4) when applying the lubrication model to neighbouring rivulets.

6.1. Rivulets using the quasi-stationary model

Hereinafter, all fields are non-dimensional and the tildes are dropped for brevity. Acknowledging the accuracy of the quasi-stationary model (§ 3.1) and recognising the limitations of the lubrication model (§ 3.2) for rivulets, we now investigate, with the quasi-stationary model, how the stagnation point $x_0$ shifts with respect to its governing parameters, namely the contact angle $\theta$ , the Marangoni number ${\textit{Ma}}$ and the distance between neighbouring rivulets $b$ .

Firstly, we evaluate the stagnation point $x_0$ for a range of ${\textit{Ma}}$ values at four different contact angles – $10^\circ$ (red), $20^\circ$ (green), $30^\circ$ (orange) and $40^\circ$ (blue) – while keeping the distance between the rivulets fixed at $b = d+2 = 2.1$ , as shown in figure 5(a). It is immediately evident that, for all considered contact angles, the stagnation point $x_0$ is independent of ${\textit{Ma}}$ . An analytical explanation for this observation is developed in § 6.1.1.

Figure 5. Stagnation point in rivulets using the quasi-stationary model. (a) Stagnation point $x_0$ versus Marangoni number for rivulets at four different contact angles: $10^\circ$ (red), $20^\circ$ (green), $30^\circ$ (orange) and $40^\circ$ (blue) and $b = 2.1$ . (b) Phase space showing the variation of the displacement of the stagnation point $x_0$ (quantifying the symmetry breaking) as a function of the contact angle $\theta$ and the inter-rivulet distance $b = d+2$ at fixed ${\textit{Ma}} = 10^5$ .

Secondly, we assess the stagnation point $x_0$ within a phase space defined by $\theta$ and the distance $b$ , for a fixed ${\textit{Ma}} = 10^5$ , as illustrated in figure 5(b). Naturally, as $b$ increases, the behaviour of neighbouring rivulets converges towards that of an isolated rivulet, which has a stagnation point at $x_0 = 0$ . Interestingly, even for a large distance $b = 10$ , $x_0$ is still different from zero, indicating that the neighbouring rivulet still has an influence on the flow field. As the contact angle $\theta$ increases, the displacement of the stagnation point $x_0$ also increases. This is expected since the shielding effect near the neighbouring rivulet is amplified with a higher contact angle since the adjacent rivulets’ interfaces are closer and there is less ‘space’ for the water vapour to diffuse away, i.e. the water vapour is more confined by the neighbouring surface when the contact angle increases, thereby accentuating the asymmetry in the flow for high contact angles.

6.1.1. Analytical solutions for rivulets using the lubrication model

For rivulets we can find the shift of the stagnation point analytically when using the lubrication model, without having to solve anything numerically. To obtain an analytical solution for the $x$ shift of the stagnation point of the velocity for rivulets, we must first solve for the velocity in the drop. Similar to before, when using the lubrication theory, we assume that the evaporation is quasi-static and that coffee-stain flow is negligible. Equation (3.14) in the lubrication approximation becomes

(6.1) \begin{align} \frac {{\rm d} p}{{\rm d} x} = \frac {\partial ^2 u_x}{\partial z^2}. \end{align}

Using the boundary conditions of no slip at the substrate $u|_{z=0} = 0$ and Marangoni stress at the interface $ {\partial u_x}/{\partial z}|_{z=h}=\theta \,{\textit{Ma}}({{\rm d}\bar {\xi } }/{{\rm d} x}) $ , we find the following solution:

(6.2) \begin{align} u_x = \frac {{\rm d} p}{{\rm d} x}\left (z^2 - 2 zh\right ) + \theta \,{\textit{Ma}}\frac {{\rm d}\bar {\xi } }{{\rm d} x}z. \end{align}

Using (3.17) we can write the pressure gradient in terms of the total flow rate $\boldsymbol{Q}$ and the concentration gradient:

(6.3) \begin{align} u_x = \theta \,{\textit{Ma}}\left (\frac {3z^2}{4h} - \frac {z}{2}\right )\frac {{\rm d}\bar {\xi } }{{\rm d} x} - \frac {3}{2}\left (\frac {z^2}{h^3} - 2\frac {z}{h^2}\right )\boldsymbol{Q}\boldsymbol{\cdot } \boldsymbol{e}_x. \end{align}

Due to the no-slip condition, the velocity at the substrate must be zero. Therefore, we determine the velocity close to the surface using $ {\partial u_x}/{\partial z}$ at $z=0$ and set it to zero to find the shift of the stagnation point. When $\boldsymbol{Q}=0$ , this condition translates to

(6.4) \begin{align} \frac {\partial u_x}{\partial z}\bigg |_{z=0} = \theta \,{\textit{Ma}}\frac {{\rm d}\bar {\xi } }{{\rm d} x}=0. \end{align}

Consequently, the stagnation point coincides with the maxima in the concentration, which we can compute starting from (3.18) with $\boldsymbol{Q} = 0$ and integrating over $x$ :

(6.5) \begin{align} \frac {{\rm d}}{{\rm d} x}\left (D_{\textit{eff}} h \frac {{\rm d} \bar {\xi }}{{\rm d} x}\right ) &= \dfrac {j}{\theta } - \dfrac {Jh}{V}, \end{align}

(6.6) \begin{align} D_{\textit{eff}} \frac {{\rm d} \bar {\xi }}{{\rm d} x} &= \dfrac {1}{\theta h} \int {j \,{\rm d}x} - \dfrac {J}{Vh}\int h\, {\rm d}x. \end{align}

Knowing that $\boldsymbol{Q} = 0$ at the contact line, from (3.16), we find that $\boldsymbol{Q} = 0$ everywhere along the rivulet and therefore $\boldsymbol{Q}_C = -\boldsymbol{Q}_M$ . Equation (3.19) then simplifies to

(6.7) \begin{align} D_{\textit{eff}} = 1 + \theta ^2 \dfrac { \boldsymbol{Q}_M^2}{420} = 1 + \dfrac {1}{1680} h^4 \theta ^4 {\textit{Ma}}^2 \left ( \dfrac {{\rm d} \bar {\xi }}{{\rm d} x} \right )^2. \end{align}

We can then rewrite the vertically averaged mass fraction evolution (6.6) as

(6.8) \begin{align} \frac {{\rm d} \bar {\xi }}{{\rm d} x} + \dfrac {1}{1680} h^4 \theta ^4 {\textit{Ma}}^2 \left ( \dfrac {{\rm d} \bar {\xi }}{{\rm d} x} \right )^3 = \dfrac {1}{\theta h} \int {j \,{\rm d}x} - \dfrac {J}{Vh}\int h\, {\rm d}x + a, \end{align}

with an integration constant $a$ . In the limit approaching isolated rivulets ( $b\gg 2$ ), the concentration $\bar {\xi }$ in the rivulet must be symmetric in $x$ . Therefore, all terms in (6.8) must be odd in $x$ . Consequently, $a$ must be zero. Using the notation $\bar {\xi }^\prime (x) = {{\rm d} \bar {\xi }}/{{\rm d} x}$ , $v(x) = ({1}/{1680}) \theta ^4 {\textit{Ma}}^2 h^4$ and $w(x)=({1}/{\theta h}) \int {j \,{\rm d}x} - (J/Vh)\int h\, {\rm d}x$ , we can write (6.8) more concisely as

(6.9) \begin{align} \bar {\xi }^\prime + v\bar {\xi }^{\prime 3} = w. \end{align}

Solving for $\bar {\xi }^\prime$ gives

(6.10) \begin{align} \bar {\xi }^\prime = \frac {\psi }{v} - \frac {1}{3\psi }, \quad \mathrm{with} \ \ \psi = \left [\frac {v^2 w}{2} + \left (\frac {v^3}{3^3} + \frac {v^4w^2}{2^2}\right )^{1/2}\right ]^{1/3}. \end{align}

The $x$ shift of the maxima in the concentration can be calculated by evaluating $\bar {\xi }^\prime = 0$ . Using (6.10) we find $v^{5/2}w=0$ . Re-substituting $v$ and $w$ gives

(6.11) \begin{align} \left ( \frac {1}{1680} \theta ^4 {\textit{Ma}}^2 h^4\right )^{5/2} \left ( \dfrac {1}{\theta h} \int {j \,{\rm d}x} - \dfrac {J}{Vh}\int h\, {\rm d}x\right ) = 0. \end{align}

The first factor only depends on $x$ via the height $h$ and is positive definite everywhere except at the rim of the rivulet. Since we are interested in the non-trivial maxima/minima, we focus on the second factor, which, because of the positive-definiteness of the first factor, must be equal to zero, i.e.

(6.12) \begin{align} \int {j \,{\rm d}x} - \dfrac {\theta J}{V}\int h\, {\rm d}x= 0. \end{align}

Note that the above expression is independent of the Marangoni number. Even if we were to consider only Marangoni flow, by neglecting the $\bar {\xi }^\prime$ term in (6.9), or only diffusion, by neglecting the $\bar {\xi }^{\prime 3}$ term, we would have found the same expression for the $x$ shift. However, the $x$ shift cannot be expressed in terms of elementary functions without approximation of $j$ and $h$ . Approximating the local evaporative flux of a neighbouring rivulet, (4.4), for large rivulet separations $b\gg 2$ yields

(6.13) \begin{align} j_{1} = J \frac {1}{\pi \sqrt {1 - x^2}}\left (1-\frac {x}{b}\right ) + \mathcal{O}\left (\frac {1}{b^2}\right ). \end{align}

Assuming that the rivulet shape is a circular arc, which is analogous to a spherical cap approximation for a drop, the rivulet shape is given by

(6.14) \begin{align} h = \frac {1}{\theta }\left [\tan {\left (\frac {\theta }{2}\right )} - \frac {1-\sqrt {1-x^2\sin ^2\theta }}{\sin {\theta }}\right ]. \end{align}

For low contact angles the rivulet height can be approximated by

(6.15) \begin{align} h = \frac {1}{2}\big (1-x^2\big ) + \frac {\theta ^2}{24}\big (1 + 2x^2 - 3x^4\big ) + \mathcal{O} (\theta ^4). \end{align}

The first-order term in (6.15) is a parabola, whose shape is independent of the contact angle. Therefore, we need to include the next higher-order term to describe the contact angle dependence of the $x$ shift. Note that this is actually a third-order approximation in $\theta$ , which is obscured by the non-dimensionalisation. The volume of the rivulet is given by

(6.16) \begin{align} V = \theta \int _{-1}^{1}h\,{\rm d}x = \frac {2}{3}\theta + \frac {4}{45}\theta ^3. \end{align}

Combining the first-order approximation of the flux in $b^{-1}$ and the third-order approximation of the rivulet height in $\theta$ in (6.12) and omitting the higher-order terms, we find

(6.17) \begin{align} \int {\frac {1}{\pi \sqrt {1 - x^2}}\left (1-\frac {x}{b}\right ) \,{\rm d}x} - \dfrac {\theta }{V} \int \left [\frac {1}{2}\big(1-x^2\big) + \frac {\theta ^2}{24}\big (1 + 2x^2 - 3x^4\big )\right ]\, {\rm d}x &= 0, \end{align}

(6.18) \begin{align} \frac {\arcsin {(x)}}{\pi } - \frac {\sqrt {1-x^2}}{\pi b}-\left [\frac {x}{2} - \frac {x^3}{6} + \frac {\theta ^2}{24}\left (x - \frac {2}{3}x^3 + \frac {3}{5}x^5\right )\right ]\left (\frac {2}{3} + \frac {4}{45}\theta ^2\right )^{-1}&=0. \end{align}

Again, we cannot progress in terms of elementary functions without approximating further. Since we have already assumed that the rivulet separation is large, we expect that the $x$ shift will be small, allowing us to approximate (6.18) for small $x$ and keep only the zeroth- and first-order terms in $x$ :

(6.19) \begin{align} \frac {x}{\pi } - \frac {1}{\pi b}- x \left (\frac {1}{2} + \frac {1}{24}\theta ^2\right )\left (\frac {2}{3} + \frac {4}{45}\theta ^2\right )^{-1}=0. \end{align}

Solving for $x$ finally gives the $x$ shift of the stagnation point and the maximum concentration in the neighbouring rivulet:

(6.20) \begin{align} x_0 = \frac {1}{b}\frac {15\left (16+2\theta ^2\right )}{\left (3\pi -4\right )\left (60+5\theta ^2\right )-12\theta ^2}. \end{align}

Figure 6. Stagnation point $x_0$ as a function of inter-rivulet distance $b = d + 2$ and contact angle $\theta$ for rivulets using the lubrication model: $x_0$ for (a) constant $\theta$ and (b) constant $b$ . The analytical approximation, (6.20), valid for small $\theta$ and large $b$ , uses a third-order expansion for the rivulet shape and a first-order expansion in $1/b$ for the local flux. (c,d) Comparison of this approximation with the numerically evaluated shift obtained by solving (6.12) numerically, without simplifying the local flux or rivulet shape.

Figure 6 shows both the analytical approximation we just derived as well as the $x$ shift without any approximations for the local flux or rivulet shape. The latter is calculated by numerically solving (6.12) for $x$ . Comparing the analytical solutions with and without approximations for small values of $d$ we find that with approximations the $x$ shift is underestimated, but the overall shape is captured accurately. The dependence on $\theta$ is nearly identical for both solutions, which is to be expected since we are using the lubrication approximation where we have already assumed that $\theta$ is small.

However, when comparing the analytical solutions of $x_0$ in figure 6 with the numerical solutions of $x_0$ in figure 5, not only do we find that the analytical solutions underestimate the numerical solutions, additionally we also find that $x_0$ depends much more strongly on $\theta$ . As mentioned above, the main reason for this dependency can be attributed to the evaporative flux, which should depend on the contact angle. When an analytical expression of the local flux for neighbouring rivulets with non-zero contact angles is found, it will likely be straightforward to update (6.20) following the same steps as used here.

6.2. Droplets using the lubrication model

We employ the lubrication model to investigate the shift of the stagnation point $x_0$ for droplets. Here, we use the contact-angle-dependent evaporation rate as calculated by Popov (Reference Popov2005). Despite the limitations of applying the lubrication model to rivulets, it has been successfully used for droplets in previous studies under the same conditions (Ramírez-Soto & Karpitschka Reference Ramírez-Soto and Karpitschka2022).

Notably, although the governing equations for droplets and rivulets differ only in the coordinate system of the operators, the stagnation point $x_0$ is dependent on the Marangoni number ${\textit{Ma}}$ for droplets. Figure 7 presents $x_0$ in phase spaces defined by: ${\textit{Ma}}$ and $\theta$ for $b=2.01$ (figure 7 a); and $b$ and $\theta$ for ${\textit{Ma}} = 10^5$ (figure 7 b). For sufficiently large $\theta$ , the distance of the stagnation point $x_0$ increases with ${\textit{Ma}}$ , whereas for small contact angles it displays a slightly non-monotonic trend, first decreasing and then increasing with ${\textit{Ma}}$ .

Figure 7. Phase spaces showing the variation of the stagnation point $x_0$ measured at the substrate as a function of (a) the Marangoni number ( ${\textit{Ma}}$ ) and contact angle ( $\theta$ ) for droplets with an inter-droplet separation ( $b=d+2$ ) of 2.01 and (b) $\theta$ and $b$ for droplets with fixed ${\textit{Ma}} = 10^5$ .

In contrast to rivulets where the Marangoni flow is largely confined to the radial or vertical directions, droplets permit Marangoni flow in the azimuthal direction as well. This additional degree of freedom leads to a more intricate evolution of the stagnation point $x_0$ in droplets compared with rivulets. Figure 8 illustrates, for two values of ${\textit{Ma}}$ , namely $10^3$ (figure 8 a) and $10^6$ (figure 8 b), with a contact angle of $\theta = 40^\circ$ and $b = 2.1$ , the volumetric flow rate $\boldsymbol{Q}$ (top left of each subfigure), the vertically averaged mass fraction $\bar {y_A}$ (bottom left) and the velocity magnitude $\|\boldsymbol{u}\|$ as determined by the parabolic profile of the lubrication model, both at the interface (top right) and at $0.001\,R$ above the substrate (bottom right). In the left-hand panels, the arrows indicate the direction of the flow rate, with the colour coding distinguishing the dominance of Marangoni-driven flow (yellow) from that of pressure-driven flow (green). In the right-hand panels, the velocity directions are indicated, revealing distinct orientations near the substrate and at the interface, as expected from the vertical recirculation in droplets.

Figure 8. Flow field in droplets for two different Marangoni numbers (a) ${\textit{Ma}} = 10^3$ and (b) ${\textit{Ma}} = 10^6$ , with a contact angle of $\theta = 40^\circ$ and an inter-droplet separation of $b = 2.1$ . The top-left panels show the flow rate $\boldsymbol{Q}$ (with arrows indicating its direction, coloured in yellow if dominated by Marangoni stresses, otherwise green), while the bottom-left panels display the vertically averaged mass fraction $\bar {y}_A$ . The top- and bottom-right panels illustrate the velocity magnitude $\|\boldsymbol{u}\|$ at the liquid–gas interface and just above the substrate, respectively, highlighting the vertical recirculation induced by the Marangoni effect. The stagnation point $x_0$ , evaluated either at the interface (top) or at the substrate (bottom), is indicated by a red dot in the right-hand panels.

Close to the neighbouring droplet, the compositional gradients are more pronounced, resulting in a stronger Marangoni effect in that region. At the droplet axis, the flow rate $\boldsymbol{Q}$ is directed away from the adjacent droplet, thereby generating a backflow in the azimuthal direction.

As ${\textit{Ma}}$ increases, the enhanced mixing within the liquid leads to a more uniform composition throughout the droplet. The flow rate $\boldsymbol{Q}$ in the backflow region becomes substantially weaker than that along the droplet’s axis. Moreover, the position of the stagnation point $x_0$ exhibits a more marked change with increasing ${\textit{Ma}}$ , particularly when evaluated at the substrate in comparison with its position at the interface, revealing large vertical gradients in the velocity. We attribute this shift in $x_0$ to the nonlinear evolution of the compositional gradients. However, due to the complexity of the flow field, we do not attempt to obtain a quantitative explanation for the observed behaviour of $x_0$ in droplets.