2.1. The dimensional problem

We consider the configuration in figure 1 in which a droplet of liquid of volume $V^{*}$ lies on a rigid, planar substrate. The substrate lies along $z^{*} = 0$, where $(r^{*},\theta ,z^{*})$ are cylindrical polar coordinates defined with respect to the centre of the droplet footprint, which is taken to be a circle of radius $R^{*}$. Here and hereafter, an asterisk denotes a dimensional variable. The droplet contact line is assumed to remain pinned throughout the motion, which is a reasonable assumption for the majority of the drying time (see Hu & Larson Reference Hu and Larson2002) and certainly while the solute remains dilute.

Following the symmetry of the problem, we make the assumption that the dynamics are independent of the polar angle $\theta$. The free surface delineating the droplet from the surrounding air is thus denoted by $z^{*} = h^{*}(r^{*},t^{*})$. We shall make the assumption that the droplet is thin so that $H^{*} = h^{*}(0,0) \ll R^{*}$.

The liquid contains a non-volatile solute of initial concentration $\phi _{{init}}^{*}$, which we shall assume to be evenly distributed throughout the droplet. We also assume throughout that the solute is sufficiently dilute that the flow within the drop is unaffected by its presence. This means that, crucially, we can decouple the flow in the liquid drop from solute transport.

2.1.1. Flow problem

The liquid has density $\rho ^{*}$ and viscosity $\mu ^{*}$, while the surface tension of the air–liquid interface is denoted by $\sigma ^{*}$; all of these material parameters are taken to be constant. For the purposes of our analysis, we shall neglect the effect of gravity; that is, we assume the Bond number $Bo = \rho ^{*} g^{*} R^{*2}/\sigma ^{*}$ is small, where $g^{*}$ is the gravitational acceleration. The liquid velocity and pressure are denoted by $\boldsymbol {u}^{*} = u_{r}^{*}\boldsymbol {e}_{r} + u_{z}^{*}\boldsymbol {e}_{z}$ and $p^{*}$, respectively, where $\boldsymbol {e}_{r}, \boldsymbol {e}_{z}$ are the unit vectors in the $r^{*}$- and $z^{*}$-directions.

Evaporation of the liquid into the surrounding air induces a flux of vapour $E^{*}$ at the droplet surface. We shall assume that the evaporation is a quasi-steady process, which is reasonable for a wide range of applications, including the evaporation of water on glass (Hu & Larson Reference Hu and Larson2002). The evaporative flux in turn drives a flow within the droplet whose size $U^{*}$ depends upon the dominant evaporative process, as discussed in more detail in § 2.2.

In each of our evaporation models, the induced velocity will be assumed to be sufficiently small and the droplet sufficiently thin that, to leading order in $\delta = H^{*}/R^{*}$, the equations of fluid motion within the drop are well-approximated by the lubrication equations (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000; Freed-Brown Reference Freed-Brown2015), i.e.

(2.1a–c)\begin{align} \frac{\partial h^{*}}{\partial t^{*}} + \frac{1}{r^{*}}\frac{\partial}{\partial r^{*}}\left(r^{*}h^{*}\bar{u}_{r}^{*}\right) ={-}\frac{E^{*}}{\rho^{*}}, \quad \bar{u}_{r}^{*} ={-}\frac{h^{*2}}{3\mu^{*}}\frac{\partial p^{*}}{\partial r^{*}}, \quad p^{*} = p_{{atm}}^{*}-\frac{\sigma^{*}}{r^{*}}\frac{\partial}{\partial r^{*}}\left(r^{*}\frac{\partial h^{*}}{\partial r^{*}}\right) \end{align}

for $0< r^{*}< R^{*}$, $t^{*}>0$, where $p_{{atm}}^{*}$ denotes atmospheric pressure and $\bar {u}_{r}^{*}$ is the depth-averaged radial velocity. (Note that, in a slight abuse of language, we shall often refer to $\bar {u}_{r}^{*}$ as simply the radial velocity for brevity.) Equations (2.1) must be solved subject to the symmetry conditions

(2.2a,b)\begin{equation} \frac{\partial h^{*}}{\partial r^{*}} = r^{*}h^{*}\bar{u}_{r}^{*} = 0 \quad \text{at} \ r^{*} = 0, \ t^{*}>0, \end{equation}

and the zero-thickness and no-flux conditions at the contact line

(2.3a,b)\begin{equation} h^{*} = r^{*}h^{*}\bar{u}_{r}^{*} = 0 \quad \text{at} \ r^{*} = R^{*} \ t^{*}>0. \end{equation}

While the initial droplet profile should capture its shape shortly after being deposited on the substrate and must be prescribed to close the problem (2.1)–(2.3), after an initial transient on the time scale of capillary action (in which $t^{*} = \mu ^{*}R^{*}/\delta ^{3}\sigma ^{*}$ by balancing the first two terms on the left-hand side of (2.1a)), the free surface rapidly approaches a spherical cap, see Lacey (Reference Lacey1982). Hence, for simplicity, we shall impose the initial profile

(2.4)\begin{equation} h^{*}(r^{*},0) = 2V^{*}(R^{*2}-r^{*2})/{\rm \pi} R^{*4} \quad \text{for} \ 0< r^{*}< R^{*}. \end{equation}

This fixes $H^{*} = 2V^{*}/{\rm \pi} R^{*2}$ and hence we require that $\delta = 2V^{*}/{\rm \pi} R^{*3}\ll 1$.

Assuming $h^{*}>0$ for $0< r^{*}< R^{*}$ for $0< t^{*}< t_{f}^{*}$, where $t_{f}^{*}$ is the evaporation or dry out time of the drop (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000), integrating (2.1a) from $r^{*} = 0$ to $r^{*} = R^{*}$ and applying the no-flux conditions (2.2b) and (2.3b), we obtain an expression representing global conservation of mass of the liquid phase, namely

(2.5)\begin{equation} \frac{\text{d}}{\text{d}t^{*}}\int_{0}^{R^{*}}r^{*}h^{*}(r^{*},t^{*})\,\text{d}r^{*} ={-} \int_{0}^{R^{*}}\frac{r^{*}E^{*}(r^{*})}{\rho^{*}}\,\text{d}r^{*}. \end{equation}

Integrating and applying the initial condition (2.4) then gives

(2.6)\begin{equation} \int_{0}^{R^{*}}r^{*}h^{*}(r^{*},t^{*})\,\text{d}r^{*} = \int_{0}^{R^{*}} r^{*}h^{*}(r^{*},0)\,\text{d}r^{*} - \left(\int_{0}^{R^{*}}\frac{r^{*}E^{*}(r^{*})}{\rho^{*}}\,\text{d}r^{*}\right)t^{*}, \end{equation}

so that

(2.7)\begin{equation} t_{f}^{*} = \frac{V^{*}}{2{\rm \pi}}\left(\int_{0}^{R^{*}}\frac{r^{*}E^{*}(r^{*})}{\rho^{*}}\,\text{d}r^{*}\right)^{{-}1}. \end{equation}

For a given evaporative flux, (2.1)–(2.4) and (2.7) fully specify the flow problem for $h^{*}(r^{*},t^{*})$, $\bar {u}_{r}^{*}(r^{*},t^{*})$, $p^{*}(r^{*},t^{*})$ and $t_{f}^{*}$.

2.1.2. Solute problem

Since we are assuming that the dilute solute has no effect on the liquid flow, the solute concentration $\phi ^{*}$ simply satisfies an advection–diffusion equation. In the limit in which $\delta \ll 1$ and $\delta ^{2}U^{*}R^{*}/D_{\phi }^{*}\ll 1$ where $D_{\phi }^{*}$ is the solutal diffusion coefficient, it is straightforward to show that $\phi ^{*}$ is independent of $z^{*}$ at leading order. Hence averaging the advection–diffusion equation across the droplet thickness yields

(2.8)\begin{equation} \frac{\partial}{\partial t^{*}}\left(h^{*}\phi^{*}\right) + \frac{1}{r^{*}}\frac{\partial}{\partial r^{*}}\left[r^{*}h^{*}\phi^{*}\bar{u}_{r}^{*} - D_{\phi}^{*}r^{*}h^{*}\frac{\partial\phi^{*}}{\partial r^{*}}\right] = 0 \end{equation}

for $0< r^{*}< R^{*}$, $t^{*}>0$ (Wray et al. Reference Wray, Papageorgiou, Craster, Sefiane and Matar2014; Pham & Kumar Reference Pham and Kumar2017). There are two clear competing physical effects in (2.8). The first is an advective flux of solute mass $\phi ^{*}h^{*}\bar {u}_{r}^{*}$, where the evaporation-induced flow carries solute to the contact line. The solute concentration thus increases local to the contact line, in turn driving a competing diffusive flux $-D_{\phi }^{*}h^{*}\partial \phi ^{*}/\partial r^{*}$ towards the droplet bulk. The relative importance of these effects is what drives the formation of the nascent coffee ring in our analysis. We have schematically illustrated these competing fluxes in figure 1.

By symmetry, we have

(2.9)\begin{equation} \frac{\partial\phi^{*}}{\partial r^{*}} = 0 \quad \text{at} \ r^{*} = 0, \ t^{*}>0, \end{equation}

while there can be no flux of particles through the contact line, so that

(2.10)\begin{equation} r^{*}h^{*}\phi^{*}\bar{u}_{r}^{*} - D_{\phi}^{*}r^{*}h^{*}\frac{\partial\phi^{*}}{\partial r^{*}} = 0 \quad \text{at} \ r^{*} = R^{*}, \ t^{*}>0. \end{equation}

Finally, the initial solute distribution is taken to be uniform and given by

(2.11)\begin{equation} \phi^{*}(r^{*},0) = \phi_{{init}}^{*} \quad \text{for} \ 0< r^{*}< R^{*}. \end{equation}

For a given flow, (2.8)–(2.11) completely specify the solute transport problem for $\phi ^{*}(r^{*},t^{*})$.

We emphasize that, in the thin-droplet, dilute-solute limit, the solute problem decouples from the flow problem, so that we can solve (2.1)–(2.4), (2.7) for $\bar {u}_{r}^{*}$, $p^{*}$, $h^{*}$ and $t_{f}^{*}$, before solving for the solute concentration $\phi ^{*}$ from (2.8)–(2.11).

To close the problem, we require the velocity field $\bar {u}_{r}^{*}$, which depends on the evaporative flux via (2.1). The evaporative flux (and hence the induced flow and structure) is dependent on the dominant evaporative process. However, we stress that, in our analysis, we are interested in the structure of the solute concentration profile near the contact line rather than determining which evaporative model is the most appropriate for a given problem. To that end, we shall consider two well-established evaporative models in this paper: a kinetic evaporation model and a diffusive evaporation model, as we shall now describe.

2.1.3. Evaporation models

In a kinetic evaporation model, it is often the case that either the surrounding gas consists entirely of the droplet vapour, or that diffusion away from the droplet surface happens sufficiently quickly that evaporation is limited by the liquid phase alone. As described in, for example, Murisic & Kondic (Reference Murisic and Kondic2011), the evaporative flux in this regime is well-approximated by the expression $E^{*} = \mathcal {A}_{k}^{*} / (h^{*}+\mathcal {B}_{k}^{*})$, where $\mathcal {A}_{k}^{*}$ and $\mathcal {B}_{k}^{*}$ depend on the properties of the droplet and thermodynamic properties of the system. For a thin droplet, Murisic & Kondic (Reference Murisic and Kondic2011), note that there are liquid/substrate systems (such as water/silicon) for which $\mathcal {B}_{k}^{*} \gg h^{*}$, so that $E^{*}\approx \mathcal {A}_{k}^{*}/\mathcal {B}_{k}^{*}$, a constant. For the purposes of this paper we therefore take the evaporative flux in this regime to be

(2.12)\begin{equation} E^{*}(r^{*}) = \mathcal{E}_{k}^{*} \end{equation}

for $0< r^{*}< R^{*}$, where $\mathcal {E}_{k}^{*}$ is a constant.

A constant evaporative flux may also be a reasonable approximation in other regimes for which kinetic evaporation is not the dominant effect. For example, Boulogne et al. (Reference Boulogne, Ingremeau and Stone2016) consider water droplets evaporating on a glass substrate resting in a large hydrogel bath. Since the hydrogel dries at a similar speed to the water, the vapour concentration in the surrounding gas is greatly affected by the bath. For sufficiently large baths, the flux around the droplet is well-approximated by a constant flux, although Boulogne et al. (Reference Boulogne, Ingremeau and Stone2016) note that convective effects may also be important in the air.

In a diffusive evaporation model, the dominant transport of the liquid vapour away from the droplet–air interface is diffusion, and for a wide range of problems, the vapour Péclet number is sufficiently large that this process is quasi-steady (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000; Hu & Larson Reference Hu and Larson2002). The vapour concentration thus satisfies a mixed boundary value problem for Laplace's equation in the air. Since the droplet is thin, this problem is equivalent to solving for the potential outside a charged disk in classical electrostatics (Sneddon Reference Sneddon1966), with the resulting evaporative flux given by

(2.13)\begin{equation} E^{*}(r^{*}) = \frac{2R^{*}\mathcal{E}_{d}^{*}}{\rm \pi}\frac{1}{\sqrt{R^{*2}-r^{*2}}}, \end{equation}

for $0< r^{*}< R^{*}$, where again $\mathcal {E}_{d}^{*}$ depends on the properties of the system (Murisic & Kondic Reference Murisic and Kondic2011).

2.2. Non-dimensionalization

The evaporative flux $E^{*}$ has a typical scale $\mathcal {E}^{*}$ that is given by $\mathcal {E}_{k}^{*}$ and $\mathcal {E}_{d}^{*}$ in the kinetic and diffusive evaporation models, respectively. The induced radial velocity is $U^{*} = \mathcal {E}^{*}/\rho ^{*}\delta$. Hence, we non-dimensionalize (2.1)–(2.4), (2.7), (2.8)–(2.11) and (2.12)–(2.13) by setting

(2.14)\begin{equation} \left. \begin{gathered} (r^{*},\ z^{*}) = R^{*}( r,\delta z), \quad \bar{u}_{r}^{*} = \frac{\mathcal{E}^{*}}{\delta\rho^{*}}\bar{u}_{r}, \quad t^{*}_{f} = \frac{\delta\rho^{*} R^{*}}{\mathcal{E}^{*}}t_{f},\quad t^{*} = \frac{\delta\rho^{*} R^{*}t_{f}}{\mathcal{E}^{*}}t ,\\ h^{*} = \delta R^{*} h, \quad p^{*} = p_{{atm}}^{*} + \frac{\mu^{*}\mathcal{E}^{*}}{\rho^{*} R^{*}\delta^{3}} p, \quad E^{*} = \mathcal{E}^{*} E, \quad \phi^{*} = \phi_{{init}}^{*}\phi. \end{gathered} \right\} \end{equation}

We note that in (2.14), the dimensionless dry out time $t_{f}$ is given by

(2.15–2.16) \begin{equation}{t_{f} = \frac{1}{4}\left(\int_{0}^{1} rE(r) \,\text{d}r\right)^{-1} =}\begin{cases} \displaystyle{\frac{1}{2}} & \text{in the kinetic regime,} \\ \displaystyle{\frac{\rm \pi}{8}} & \text{in the diffusive regime,} \end{cases}\end{equation}

and we have rescaled time to fix its domain to be $0< t<1$, which will simplify substantially numerous expressions in the sequel.

Under the scalings (2.14), the flow problem (2.1)–(2.4) becomes

(2.17a–c)\begin{equation} \frac{1}{t_{f}}\frac{\partial h}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}\left(rh\bar{u}_{r}\right) ={-}E, \quad \bar{u}_{r} ={-}\frac{h^{2}}{3}\frac{\partial p}{\partial r}, \quad p ={-}\frac{1}{{\textit{Ca}}}\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial h}{\partial r}\right) \end{equation}

for $0< r<1$, $t>0$, such that

(2.18a,b)\begin{equation} \frac{\partial h}{\partial r} = rh\bar{u}_{r} = 0 \quad \text{at} \ r = 0, \ t>0 \end{equation}

and

(2.19a,b)\begin{equation} h = rh\bar{u}_{r} = 0 \quad \text{at} \ r = 1, \ t>0, \end{equation}

along with the initial condition

(2.20)\begin{equation} h(r,0) = 1-r^{2} \quad \text{for} \ 0< r<1, \end{equation}

where

(2.21)\begin{equation} {\textit{Ca}} = \frac{\mathcal{E}^{*}\mu^{*}}{\sigma^{*}\rho^{*}\delta^{4}} \end{equation}

is an enhanced droplet capillary number.

Similarly, the solute problem (2.8)–(2.11) becomes

(2.22)\begin{equation} \frac{1}{t_{f}}\frac{\partial}{\partial t}\left(h\phi\right) + \frac{1}{r}\frac{\partial}{\partial r}\left[rh\phi\bar{u}_{r} - \frac{rh}{{\textit{Pe}}}\frac{\partial\phi}{\partial r}\right] = 0\end{equation}

for $0< r<1$, $t>0$, such that

(2.23a,b)\begin{equation} \frac{\partial\phi}{\partial r} = 0, \quad \text{at} \ r = 0, \ t>0, \quad rh\phi\bar{u}_{r} - \frac{rh}{{\textit{Pe}}}\frac{\partial\phi}{\partial r} = 0 \quad \text{at} \ r = 1, \ t>0, \end{equation}

and initially

(2.24)\begin{equation} \phi(r,0) = 1 \quad \text{for} \ 0< r<1. \end{equation}

Here, the solutal Péclet number, ${\textit {Pe}}$, is defined by

(2.25)\begin{equation} {\textit{Pe}} = \frac{\mathcal{E}^{*}R^{*}}{\rho^{*} D_{\phi}^{*}\delta}. \end{equation}

Finally, the dimensionless evaporative flux is given for $0< r<1$ by

(2.26–2.27)\begin{equation}{E(r) =} \begin{cases}1 & \text{in the kinetic regime,} \\ \displaystyle{\frac{2}{\rm \pi}\frac{1}{\sqrt{1-r^{2}}}} & \text{in the diffusive regime.} \end{cases} \end{equation}

2.3. Small capillary number limit

There are two dimensionless parameters in the problem, ${\textit {Ca}}$ and ${\textit {Pe}}$. We shall consider the regime in which ${\textit {Ca}}\ll 1$, so that the fluid motion is dominated by surface tension. We shall then proceed to consider the large-${\textit {Pe}}$ sublimit of this model, so that the importance of solute diffusion is confined to a region near the pinned contact line. Firstly, however, we will briefly illustrate the validity and applicability of these limits by considering several experimental studies from the literature.

For a wide range of physical problems, the relevant limit for the flow model is that in which the droplet capillary number is small, so that surface tension dominates viscous forces in determining the droplet profile. To illustrate the pertinence of this assumption, consider a typical example from Hu & Larson (Reference Hu and Larson2002): a water droplet of volume and radius $V^{*} = 0.5\ \mathrm {\mu }\textrm {l}$ and $R^{*} = 1 \ \textrm {mm}$, respectively, sits on a glass substrate and evaporates diffusively (cf. § 4) into the surrounding air. The typical flow velocities induced by the evaporation are reported to be of the order of $U^{*} = 1\ \mathrm {\mu }\textrm {ms}^{-1}$. Thus, $\delta \approx 0.3$ and ${\textit {Ca}} \approx 4.3\times 10^{-7}\ll 1$.

We see similar orders of magnitude for the droplet aspect ratio and the capillary number for other common liquids as well. Kajiya et al. (Reference Kajiya, Kaneko and Doi2008) study the coffee-ring structure for evaporating droplets of anisole, for which $\mu ^{*} = 1.03\times 10^{-3}\ \textrm {Pa}\ \textrm {s}$ and $\sigma^{\ast} = 3.5\times 10^{-2}\ \textrm {Nm}^{-1}$. Again, the anisole droplet has volume and radius $V^{*} = 0.5\ \mathrm {\mu }\textrm {l}$ and $R^{*} = 1 \ \textrm {mm}$, respectively, while the authors report a drying time of $\approx 8\ \textrm {min}$, which gives a fluid velocity $U^{*}\sim 0.8 \ \mathrm {\mu }\textrm {ms}^{-1}$. Hence, for this example, we find that $\delta \approx 0.3$ and ${\textit {Ca}} \approx 7.47\times 10^{-7}$.

As a final example, we consider the experimental and numerical analysis of evaporating ethanol droplets in Sáenz et al. (Reference Sáenz, Wray, Che, Matar, Valluri, Kim and Sefiane2017). The authors present experiments for a wide variety of droplet geometries, with the circular geometry of interest to the present study concerning droplets of radius $R^{*} \approx 2\ \textrm {mm}$ and volume $V^{*} = 7\ \mathrm {\mu }\textrm {l}$. Although the authors do not explicitly estimate the induced flow velocity, they state that the final evaporation time of the droplets is $O(100)\ \textrm {s}$, which leads to $U^{*}\approx 10\ \mathrm {\mu }\textrm {ms}^{-1}$. Thus, since ethanol has viscosity $\mu ^{*} = 1.15\times 10^{-3}\ \textrm {Pa}\ \textrm {s}$ and surface tension $\sigma^{\ast} = 2.1\times 10^{-2}\ \textrm {Nm}^{-1}$, we find that $\delta \approx 0.6$ and the droplet capillary number is given by ${\textit {Ca}}\approx 2.49\times 10^{-6}$.

Hence, as we see, for a wide variety of different liquids, we are very much entrenched in the small-${\textit {Ca}}$ regime, so that in what follows, we shall assume that surface tension dominates the fluid motion.

For all of these examples, we have seen that $R^{*}\approx 1\ \textrm {mm}$ and $U^{*}\approx 1 - 10\ \mathrm {\mu }\textrm {ms}^{-1}$. Hence, perhaps unsurprisingly, the size of the solutal Péclet number is most dependent on the size of the diffusion coefficient, $D^{*}_{\phi }$. If we return to the evaporating anisole droplets in Kajiya et al. (Reference Kajiya, Kaneko and Doi2008), in which the solute is a fluorescent polystyrene for which the authors note $D_{\phi }^{*} \approx 2\times 10^{-11}\ \textrm {m}^{2}\ \textrm {s}^{-1}$, we find that ${\textit {Pe}} \approx 40$.

More generally, if we assume that the solute particles are approximately spherical and hence that $D_{\phi }^{*}$ is well-modelled by the Stokes–Einstein equation, it is the size of the solute particles that is the dominant factor in determining the size of the solutal Péclet number (note that for all three fluids considered above, the viscosities were approximately $10^{-3}\ \textrm {Pa}\ \textrm {s}$). Thus, if we consider solute particle radii varying from $10\ \mathrm{\mu}$m (e.g. polystyrene microspheres, as in Deegan et al. (Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000)) to $5\times 10^{-4}\ \mathrm {\mu }\textrm {m}$ (e.g. a sugar molecule), the Péclet number varies from ${\textit {Pe}}\sim 3$ to ${\textit {Pe}}\sim 5\times 10^{4}$. Thus, it is clear our assumption that ${\textit {Pe}}\gg 1$ is applicable to a wide range of problems.

It is worth noting that, in the examples considered above, $\delta$ is not particularly small initially, so that we must be careful in the assumption of small-reduced Péclet number, $\delta ^{2}{\textit {Pe}}$, that is made in writing down (2.8). However, during evaporation the droplet necessarily becomes more slender, so that this assumption becomes more reasonable.

The above estimates motivate us to pursue a small-${\textit {Ca}}$, large-${\textit {Pe}}$ solution of (2.17)–(2.27). Firstly, let us consider the flow in the droplet. As ${\textit {Ca}}\rightarrow 0$, it is evident that

(2.28)\begin{gather} h \sim \left(1-t\right)(1-r^{2}), \end{gather}

(2.29)\begin{gather}p \sim \frac{4}{{\textit{Ca}}}(1-t), \end{gather}

(2.30)\begin{gather}\bar{u}_{r} \sim{-}\frac{1}{4t_{f}r}\frac{(1-r^{2})}{(1-t)} + \frac{1}{r(1-t)(1-r^{2})}\int_{r}^{1}sE(s)\,\text{d}s, \end{gather}

to leading order (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000). Inspecting the final term in (2.30) reveals that, depending on the exact form of the evaporative flux, the depth-averaged radial velocity may be bounded as $r\rightarrow 1$ (for (2.26)) or singular as $r\rightarrow 1$ (for (2.27)). The behaviour of $\bar {u}_{r}$ as $r\rightarrow 1$ will play an instrumental role in determining the structure of the diffusive boundary layer governed by (2.22)–(2.24) close to the contact line, further justifying our consideration of these two possible behaviours in §§ 3 and 4.

2.4. Formulation in terms of the solute mass

Having determined the leading-order flow in the droplet, we turn to the solute problem. We assume that $\bar {u}_{r}$ and $h$ are well-approximated by their leading-order forms in the small-${\textit {Ca}}$ expansion and now concentrate on determining the asymptotic solution to (2.22)–(2.24) in the limit in which $\varepsilon := 1/{\textit {Pe}} \rightarrow 0^{+}$. Note, to be asymptotically consistent with our derivation of (2.22)–(2.24), we are therefore considering the regime in which $\delta ^{2}\ll \varepsilon \ll 1$. While the concentration, $\phi$, is important in determining when the dilute regime breaks down within the boundary layer (as we discuss in detail in § 5), we find it is convenient to proceed with our asymptotic analysis by introducing the solute mass per unit area, $m = \phi h$, which satisfies

(2.31)\begin{equation} \frac{1}{t_{f}}\frac{\partial m}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}\left[r\left(\bar{u}_{r} + \frac{\varepsilon}{h}\frac{\partial h}{\partial r}\right)m - \varepsilon r\frac{\partial m}{\partial r}\right] = 0 \quad \text{for}\ \ 0< r<1, \ t>0, \end{equation}

subject to

(2.32a,b)\begin{align} \frac{\partial m}{\partial r} = 0, \quad \text{at} \ r = 0, \ t>0, \quad r\left(\bar{u}_{r} + \frac{\varepsilon}{h}\frac{\partial h}{\partial r}\right)m - \varepsilon r\frac{\partial m}{\partial r} = 0 \quad \text{at} \ r = 1, \ t>0 \end{align}

and

(2.33)\begin{equation} m(r,0) = h(r,0) = 1-r^{2} \quad \text{for} \ 0< r<1. \end{equation}

Finally, since it will be useful in what follows, we note that global conservation of solute dictates that

(2.34)\begin{equation} \int_{0}^{1} rm(r,t)\,\text{d}r = \int_{0}^{1} rm(r,0)\,\text{d}r = \frac{1}{4}. \end{equation}

We also note that, in addition to being mathematically convenient, using the solute mass has an advantage over $\phi$ because it is related to the absorbance of the deposit via the Beer–Lambert law (Swinehart Reference Swinehart1962); it may therefore be easier to compare predictions of $m$ directly with experimental data.

As mentioned previously, the asymptotic structure is sensitive to the behaviour of the radial velocity, which is given by (2.30), close to the contact line. For a kinetic evaporative flux, the velocity is bounded at the contact line, which we consider in § 3. On the other hand, for a diffusive evaporative flux, the velocity is singular at the contact line, and we consider this case in § 4. Our goal is to formulate a composite expansion for the solute mass valid everywhere in the droplet to leading order and to use the asymptotic results to establish thereby the dynamics of the nascent coffee ring for each flux law. We will compare these predictions with numerical simulations of the full system (2.31)–(2.33) (aided by a further reformulation that is motivated, described and assessed in Appendix B).

3.1. Outer region

In the outer region, we naïvely expand $m = m_{0} + \varepsilon m_{1} + O(\varepsilon ^{2})$ as $\varepsilon \rightarrow 0$, where, as we shall shortly see, we must proceed to $O(\varepsilon )$ in the outer region in order to be able to construct a composite mass solution that vanishes at the contact line. To leading order in $\varepsilon$, (2.31) and (2.33) become

(3.2)\begin{equation} \frac{\partial m_{0}}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}\left(\frac{r^{2}}{4(1-t)}m_{0}\right) = 0 \quad \text{for} \ 0< r<1, \ t>0, \end{equation}

with $m_{0}(r,0) = 1-r^{2}$ for $0< r<1$. Hence, as expected, the leading-order solute mass is simply advected to the contact line by the radial flow. We can solve (3.2) using the method of characteristics, finding

(3.3)\begin{equation} m_{0}(r,t) = a(t)\left(1-r^{2}a(t)\right), \end{equation}

where we have introduced the function $a(t) = \sqrt {1-t}$. We note that (3.3) can be used to determine, to leading order in $\varepsilon$, the total amount of solute swept into $r = 1$ by time $t$, denoted by $\mathcal {M}(t)$, namely

(3.4)\begin{equation} \mathcal{M}(t) = t_{f}\int_{0}^{t} \bar{u}_{r}(1^{-},\tau)m_{0}(1^{-},\tau)\,\text{d}\tau = \frac{1}{2}\left(1-\sqrt{1-t}-\frac{t}{2}\right), \end{equation}

which was previously reported for a kinetic evaporative flux by Freed-Brown (Reference Freed-Brown2015) and Boulogne et al. (Reference Boulogne, Ingremeau and Stone2016).

At $O(\varepsilon )$, we have

(3.5)\begin{equation} \frac{\partial m_{1}}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}\left(\frac{r^{2}}{4(1-t)}m_{1}\right) = \frac{2a(t)(1-a(t))}{(1-r^{2})^{2}} \quad \text{for} \ 0< r<1, \ t>0, \end{equation}

where $m_{1}(r,0) = 0$ for $0< r<1$. This can be solved in a similar manner, yielding

(3.6)\begin{gather} m_{1}(r,t) = \frac{1}{r}\frac{\partial}{\partial r}\left[r^{2}a(t)\left((1-r^{2}a(t))t - \frac{2}{3}\left(1-a(t)^{3}\right) \right. \right. \nonumber\\ - \left. \left. 2a(t)r^{2}\left(1-r^{2}a(t)\right)\left(a(t) - 1 + a(t)r^{2}\log\left(\frac{a(t)(1-r^{2})}{1-a(t)r^{2}}\right)\right)\right)\right]. \end{gather}

Note that, as $r\rightarrow 1$, $m_{0} \rightarrow a(t)(1-a(t))$ and $m_{1} = O(1/(1-r))$, both of which are physically unreasonable, as the mass should vanish at the contact line because the droplet thickness vanishes there. This is a clear indication of the need to consider the behaviour close to $r = 1$, where the effects of solute diffusion become relevant.

3.2. Inner region

We scale into the inner region by setting

(3.7)\begin{equation} r = 1 - \varepsilon R, \ m = \varepsilon^{{-}1}M \end{equation}

in (2.31)–(2.33), where the scaling for the mass has been determined from the global conservation of mass condition, (2.34). Then, to account for the logarithmic terms in the local expansion of (3.6) at the contact line, we seek an asymptotic series of the form

(3.8)\begin{equation} M = M_{0} + (\varepsilon\log\varepsilon) M_{1} + \varepsilon M_{2} + O(\varepsilon^{2}\log\varepsilon) \end{equation}

as $\varepsilon \rightarrow 0$. The leading-order inner problem is given by

(3.9)\begin{equation} \frac{\partial}{\partial R}\left[\left(\alpha(t)-\frac{1}{R}\right) M_{0} + \frac{\partial M_{0}}{\partial R}\right]=0 \quad \text{for} \ R>0, \ t>0, \end{equation}

such that

(3.10)\begin{equation} \left(\alpha(t)-\frac{1}{R}\right)M_{0} +\frac{\partial M_{0}}{\partial R} = 0 \quad \text{at} \ R = 0, \ t>0, \end{equation}

where we have introduced the function $\alpha (t) = 1/(2(1-t))$. Equations (3.9)–(3.10) are readily solved, yielding

(3.11)\begin{equation} M_{0} = A(t) R\exp{\left(-\alpha(t) R\right)}, \end{equation}

where $A(t)$ is a function of time that we shall determine shortly. The functional form of $M_{0}$ – which may be viewed as the probability density function for a gamma distribution – drives the characteristic nascent coffee-ring profile in the boundary layer. We shall discuss this further in § 3.5.

At $O(\varepsilon \log \varepsilon )$, the inner problem is the same as that at leading order, so that

(3.12)\begin{equation} M_{1} = B(t) R\exp{\left(-\alpha(t) R\right)}, \end{equation}

where $B(t)$ is to be determined.

At $O(\varepsilon )$, we find that

(3.13)\begin{equation} \frac{\partial}{\partial R}\left[\left(\alpha(t)-\frac{1}{R}\right)M_{2} + \frac{\partial M_{2}}{\partial R}\right] = 2\frac{\partial M_{0}}{\partial t} - \frac{\partial}{\partial R}\left[\left(\frac{3}{2}-2\alpha(t) R\right) M_{0} - R\frac{\partial M_{0}}{\partial R}\right] \end{equation}

in $R>0$, $t>0$, such that

(3.14)\begin{equation} \left(\alpha(t)-\frac{1}{R}\right)M_{2} + \frac{\partial M_{2}}{\partial R} + \left(\frac{3}{2}-2\alpha(t) R\right) M_{0} - R\frac{\partial M_{0}}{\partial R} = 0 \quad \text{at} \ R = 0, \ t>0. \end{equation}

The solution to this problem is given by

(3.15)\begin{align} M_{2} = \left[C_{1}(t)R + C_{2}(t)R^{2} + C_{3}(t)R^{3} + C_{4}(t)R\left(\text{Ei}(\alpha(t) R) - \log(\alpha(t) R)\right) \right]\exp\left(-\alpha(t) R\right), \end{align}

where $\text {Ei}(x)$ is the exponential integral, $C_{1}(t)$ is to be determined and

(3.16a–c)\begin{align} C_{2}(t) = \frac{15}{2}A(t) - \frac{2}{\alpha(t)}\dot{A}(t), \quad C_{3}(t) = \frac{5\alpha(t)}{2} A(t), \quad C_{4}(t) = \frac{2}{\alpha(t)^{2}}\dot{A}(t) - \frac{8}{\alpha(t)}A(t). \end{align}

3.3. Determining the unknown coefficients

It remains to determine the unknown coefficients $A(t)$, $B(t)$ and $C_{1}(t)$ in (3.11), (3.12) and (3.15), respectively. In principle, we can obtain these by matching between the inner and outer solutions, although this is notably challenging as (3.11) and (3.12) decay exponentially in the far-field, so we must proceed to higher order to complete the matching. We already have enough information to match for $A(t)$ in this manner, but we would need to go to even higher order in the inner region to match for $B(t)$ and $C_{1}(t)$, which becomes algebraically taxing. To avoid this work, we can instead appeal to global conservation of solute, (2.34). The details of this approach are given in Appendix A, and show that

(3.17a,b)\begin{equation} A(t) = \frac{1}{4(1-t)^{2}}\left[\frac{1}{2}-\frac{\sqrt{1-t}}{2}-\frac{t}{4}\right], \quad B(t) = \frac{1-\sqrt{1-t}}{2\sqrt{1-t}} \end{equation}

and

(3.18)\begin{gather} C_{1}(t) ={-}\alpha(t)^{2}a(t)\left[(1-a(t))t - \frac{2}{3}\left(1-a(t)^{3}\right) \right. \nonumber\\ - \left. 2a(t)\left(1-a(t)\right)\left(a(t) - 1 + a(t)\log\left(\frac{2a(t)}{1-a(t)}\right)\right)\right] \nonumber\\ - \frac{2(C_{2}(t)-A(t))}{\alpha(t)} - \frac{6C_{3}(t)}{\alpha(t)^{2}} - 2C_{4}(t)(\gamma-1) - C_{4}(t)\log\alpha(t), \end{gather}

where $\gamma $ is the Euler–Mascheroni constant.

Note that the term in square brackets in $A(t)$ is simply the total solute mass advected into the boundary layer from the outer region at leading order, (3.4).

3.4. Composite expansion

Having determined the unknown coefficients in the inner and outer expansions, we are left to form a composite mass profile. A natural first attempt would be to simply write down

(3.19)\begin{equation} m_{{comp}} = m_{0}(r,t) + \frac{1}{\varepsilon} M_{0}\left(\frac{1-r}{\varepsilon},t\right), \end{equation}

where $m_{0}$ and $M_{0}$ are given by (3.3) and (3.11), respectively. As we shall see, while (3.19) does a fine job in approximating the mass profile in the droplet bulk, close to the contact line, there is an $O(1)$ error in the mass profile caused by the fact that $m_{0}$ is finite as $r\rightarrow 1$.

As forecast by our consideration of the higher-order outer and inner problems, to remedy this we can write down a second composite expansion that is valid up to $O(1)$ for all $r$. This entails using $m_{0}, M_{0}, M_{1}$ and $M_{2}$. To form the composite, we turn to Van Dyke's matching rule (Van Dyke Reference Van Dyke1964), which requires a knowledge of the common contribution from the outer and inner solutions. We introduce the intermediate variable $r = 1 - \varepsilon ^{k}\bar {r} = 1-\varepsilon R$, where $0< k<1$. Substituting this scaling into the inner and outer solutions and expanding as $\varepsilon \rightarrow 0$, we see that

(3.20a,b)\begin{align} m_{{inner}} = \frac{C_{4}(t)}{\alpha(t)^{2}} + o(1) = a(t)(1-a(t)) + o(1), \quad m_{{outer}} = a(t)(1-a(t)) + o(1), \end{align}

where we have exploited the known value of $A(t)$ in the inner expansion to evaluate the expression in (3.20a).

Therefore, the appropriate composite expansion valid to $O(1)$ for all $r$ is

(3.21)\begin{gather} m_{{comp}} = m_{0}(r,t) + \frac{1}{\varepsilon} M_{0}\left(\frac{1-r}{\varepsilon},t\right) \nonumber\\ + \,(\log\varepsilon)M_{1}\left(\frac{1-r}{\varepsilon},t\right) + M_{2}\left(\frac{1-r}{\varepsilon},t\right) - a(t)(1-a(t)), \end{gather}

where $M_{1}$ and $M_{2}$ are given by (3.12) and (3.15), respectively, and we recall $a(t) = \sqrt {1-t}$. We note that, although we do not need the precise form of $m_{1}$ to write down this composite solution, we did need to ascertain its behaviour as $r\rightarrow 1$ to perform the correct asymptotic expansion in the inner region, which is why we considered the $O(\varepsilon )$-outer problem in § 3.1. Moreover, the logarithmic singularity in $m_{1}$ given by (3.6) makes it challenging to write down a higher-order composite expansion that is valid everywhere.

3.5. Profiles and properties of the coffee ring

Now we have determined two composite solutions for the solute mass, we are able to demonstrate the formation and evolution of the nascent coffee-ring effect as solute is advected from the outer region into the diffusive boundary layer, which grows in a characteristic gamma-distribution-type profile.

In the following, an important quantity is the modified time-dependent Péclet number (henceforth the modified Péclet number),

(3.22)\begin{equation} {\textit{Pe}}_{t} = \frac{{\textit{Pe}}}{1-t}, \end{equation}

which is the ratio of advective to diffusive transport modified to take into account the time dependence of the evaporation-induced velocity, the latter scaling with the reciprocal of the time to dry out, i.e. $(1-t)^{-1}$. Since we shall see that at leading order the thickness and height of the diffusive boundary layer scale with functions of the modified Péclet number, it captures in a single quantity the scaling behaviour of the nascent coffee ring with respect to both the Péclet number, ${\textit {Pe}}$, and the time remaining until dry out.

Of particular note is that our asymptotic analysis suggests that the form of the nascent coffee ring is dominated by the leading-order inner solution (3.11) for ${\textit {Pe}}_{t}\gg 1$. In particular, we expect the nascent coffee ring to tend to the similarity form given by

(3.23)\begin{equation} \frac{M_{0}(R,t)}{\mathcal{M}(t){\textit{Pe}}_{t}} = \frac{R}{4}\text{e}^{{-}R/2} = f\left(R;2,\frac{1}{2}\right), \; R = {\textit{Pe}}_{t}(1-r), \end{equation}

where $\mathcal {M}(t)$ is given by (3.4) and $f(x,k,l) = l^{k}x^{k-1}\text {e}^{-l x}/\varGamma (k)$ is the probability density function of a gamma distribution.

To validate our asymptotic predictions, we have also solved (2.31)–(2.33) numerically. This solution is hindered by the boundary layer thickness and height scaling with $1/{\textit {Pe}}_{t}$ and ${\textit {Pe}}_{t}$ as ${\textit {Pe}}_{t}\rightarrow \infty$; we have ensured that this boundary layer is resolved by discretizing in a suitable manner, as described in Appendix B.

We plot both the numerical and asymptotic predictions of the solute mass profile in figure 2 for ${\textit {Pe}} = 100$. In the top row, we compare the composite mass solution given by (3.19) with the numerical solution, while in the bottom row we show the comparison for the composite mass profile given by (3.21). In each plot, we see the evolution of the solute mass from the initial profile (bold, black curve). As the droplet evaporates, the geometry of the droplet induces a radial flow towards the contact line, which advects solute outwards. In the boundary layer, diffusion leads to the characteristic coffee-ring spike, which is clearly seen in the insets to each plot. While it is clear that both composite profiles do an excellent job of capturing the solute in the droplet bulk, there is a clear deviation of the lower-order composite (3.19) from the expected solution in the boundary layer, leading to $O(1)$ errors. The higher-order inner solutions accounted for in (3.21) rectify this deficiency, and we see excellent agreement between the asymptotics and the numerics, particularly as $t$ increases. We note in particular that, at 90 % of the dry out time (the final curve in each plot), the peak mass in the boundary layer is approximately $22$ times the initial peak at the centre of the drop (for which $m_{{max}} = 1$).

Figure 2.

Profiles of the solute mass as the droplet evaporates under a kinetic flux, $E = 1$, with ${\textit {Pe}} = 100$. In each panel, the bold, black curve represents the initial mass profile (which is identical to the droplet profile since $\phi (r,0) = 1$). Also shown are plots at time intervals of 0.1 up to $t = 0.9$ in which solid, blue curves represent the numerical results and the dashed, red curves show the composite mass profiles predicted by (3.19) (a,b) and (3.21) (c,d) for large ${\textit {Pe}}$. Panels (b,d) display a doubly logarithmic plot of the mass profile near the contact line, where we can clearly see the formation of the concentration spike that becomes the nascent coffee ring as $t$ increases.

We depict in figure 3 the anticipated collapse of the nascent coffee ring onto the similarity form (3.23). For a wide range of ${\textit {Pe}}$, we see clear evidence of this similarity form emerging after an initial transient, demonstrating the universality of the gamma distribution profile (3.23) in the diffusive boundary layer.

Figure 3.

The scaled solute mass $m/\mathcal {M}(t){\textit {Pe}}_{t}$ against the radial similarity coordinate ${\textit {Pe}}_{t}(1-r)$ as a function of time. The dashed, black curve represents the scaled leading-order inner solute mass profile (3.23). The coloured curves are results from the numerical simulations at $t = 0.3$ (dark purple), $t = 0.6$ (blue) and $t = 0.9$ (light green) for ${\textit {Pe}} = 10$ (solid), ${\textit {Pe}} = 100$ (dots) and ${\textit {Pe}} = 1000$ (dash-dots). As $t\rightarrow 1$ and ${\textit {Pe}}\rightarrow \infty$, the results collapse onto the similarity form given by (3.23), as expected.

Since we are now armed with sufficient evidence of the validity of our asymptotic results, we can use them to investigate various aspects of the solutal flow dynamics and the development of the nascent coffee ring. In particular, this analysis was motivated by establishing the competition between the advective mass flux $J_{{adv}} = \bar {u}_{r} h\phi$ and the diffusive mass flux $J_{{diff}} = -(h/{\textit {Pe}})\partial \phi /\partial r$ in different regions of the droplet. We display both of these leading-order fluxes and the resulting leading-order mass flux at various times for ${\textit {Pe}} = 100$ in figure 4. As is clearly seen, the advective flux carries solute particles towards the contact line, while the diffusive flux acts to take particles away from the high concentration at the contact line. The magnitude of the maximum advective and diffusive fluxes is comparable at each instant, which is expected from the asymptotic structure of the model: in the boundary layer, the dominant balance is between these fluxes at the expense of the time derivative term in (2.31). In each of the profiles shown in figure 4, it is clear that the advective flux is larger, so that the resultant movement of solute is towards the pinned contact line, which leads to the development of the coffee ring. At each time instant, there is a slow increase of the overall mass flux towards the contact line, which rapidly falls as we enter the boundary region where the diffusive flux is important. It should be noted that, as we can see in figure 4, the peak mass flux moves towards the contact line as $t\rightarrow 1$; this is indicative of the majority of mass having been advected radially out from the droplet centre as it dries. Moreover, it is clear that, while the overall mass flux towards the contact line generally increases in magnitude as t increases, the relative mass flux has fallen: at time $t = 0.9$, the maximum advective and diffusive fluxes are roughly three orders of magnitude larger than those at $t = 0.3$, while the maximum overall mass flux has increased by just a factor of $4$.

Figure 4.

Transient profiles of the advective particle flux $J_{{adv}}$ (a), the diffusive particle flux $J_{{diff}}$ (b) and the overall mass flux (c) for a droplet evaporating in the kinetic evaporative flux regime with ${\textit {Pe}} = 100$ at times $t = 0.3$ (dark purple), $t = 0.6$ (blue) and $t = 0.9$ (light green). These profiles have been calculated from the composite solution (3.21).

We conclude by using our asymptotic results to predict characteristics of the nascent coffee ring that may be of use in applications and, in particular, measurable in an experimental setting. As it forms, the maximum height, $m_{{max}}$, of the solute profile and the radial position of the maximum, $r_{{max}}$, can be well-approximated by maximizing the leading-order inner solution, (3.11), yielding

(3.24a,b)\begin{equation} m_{{max}}(t) = \frac{\mathcal{M}(t){\textit{Pe}}_{t}}{2\text{e}}, \quad r_{{max}}(t) = 1 - \frac{2}{{\textit{Pe}}_{t}}, \end{equation}

where $\mathcal {M}(t)$ is given by (3.4). In particular, we note that, as the droplet evaporates, the peak moves closer to the contact line and the maximum solute profile height diverges as $t\rightarrow 1$. This is illustrated for three different Péclet numbers in figure 5. The figure also demonstrates that the asymptotics provide a good approximation of the numerical results at an earlier time for larger Péclet numbers, as expected. We note that the sharp changes in the figure correspond to the maximum in the solute mass profile moving from the droplet centre to the nascent coffee ring in the boundary layer. This happens at earlier times for larger Péclet numbers.

Figure 5.

The maximum height of the solute mass profile $m_{{max}}$ (a), its location $1-r_{{max}}$ (b) and the profile full-width at half-maximum $w_{1/2}$ (c) as a function of time in the kinetic evaporative flux regime. In this figure, we consider the numerical results for three different Péclet numbers, ${\textit {Pe}} = 10^{2}$ (dark purple), ${\textit {Pe}} = 10^{3}$ (blue), ${\textit {Pe}} = 10^{4}$ (light green). For each case, the dashed, black curve corresponds to the asymptotic predictions (3.24a,b) and (3.25). Note that $w_{1/2}$ is only well-defined for $t \gtrsim 0.2$ ($t \gtrsim 0.5$ for ${\textit {Pe}} = 100$); see main text for details.

Finally, we consider a measure of the size of the coffee ring by considering the width of the profile at half its peak height – the full-width at half-maximum – which we denote by $w_{1/2}$. Asymptotically, this is given by

(3.25)\begin{equation} w_{1/2}(t) = \frac{2}{{\textit{Pe}}_{t}}\left[W_{0}\left(\frac{-1}{2{e}}\right) - W_{{-}1}\left(\frac{-1}{2{e}}\right)\right], \end{equation}

where $W_{0}(x)$ and $W_{-1}(x)$ are the Lambert-W functions (i.e. solutions to $w\text {e}^{w} = x$), (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010). Clearly, the width of the nascent coffee ring shrinks as the drop evaporates, with the majority of the mass confined to a sharp, narrow peak – note that (3.11) tends to a delta function, $\delta (R)$, as $t\rightarrow 1$. We plot the asymptotic prediction (3.25) in figure 5 alongside the evolution of $w_{1/2}$ extracted from the numerical results. We note that results are only shown for $t$ such that this definition makes sense: for small times, the minimum in the solute mass profile between the coffee ring and the droplet bulk lies above half the maximum height, making $w_{1/2}$ ill-defined. As the Péclet number increases, $w_{1/2}$ is well-defined for a longer time period. As with the comparisons for $m_{{max}}$ and $r_{{max}}$, we see that the leading-order asymptotic prediction holds over a longer time interval as ${\textit {Pe}}$ increases, while even for the smallest Péclet number depicted, the asymptotic prediction captures the behaviour of the half-width $w_{1/2}$ well as $t\rightarrow 1$.

4.1. Outer region

In the outer region, we expand $m = m_{0} + O(\varepsilon )$ as $\varepsilon \rightarrow 0$ in (2.31)–(2.33). Then, to leading order in $\varepsilon$, we obtain the advection equation

(4.2)\begin{equation} \frac{\partial m_{0}}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}\left[\frac{1}{4(1-t)}\left(\frac{1}{\sqrt{1-r^{2}}} - (1-r^{2})\right)m_{0}\right] = 0 \quad \text{for} \ 0< r<1, \ t>0, \end{equation}

with $m_{0}(r,0) = (1-r^{2})$ for $0< r<1$.

Equation (4.2) can again be solved by the method of characteristics; we find that

(4.3)\begin{equation} m_{0} = \sqrt{1-r^{2}}(1-t)^{3/4}\left[1-(1-t)^{3/4}(1-(1-r^{2})^{3/2})\right]^{1/3}. \end{equation}

In particular, we note that $m_{0}$ vanishes as $r\rightarrow 1$, in contrast to the kinetic regime for which $m_{0}$ was bounded, but finite at the contact line. We will show that it is for this reason that we do not need to proceed to higher order here: the leading-order outer solution is sufficient to construct a composite mass profile that is accurate enough for our purposes. We also note that we can use (4.3) to recover (to leading order in $\varepsilon$) the total mass of solute advected into $r = 1$ by the radial outward flow as a function of time, viz.

(4.4)\begin{equation} \mathcal{M}(t) = t_{f}\int_{0}^{t} \bar{u}_{r}(1^{-},\tau)m_{0}(1^{-},\tau)\,\text{d}\tau = \frac{1}{4}\left(1-(1-t)^{3/4}\right)^{4/3}, \end{equation}

a result previously reported in Deegan et al. (Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000) and Popov (Reference Popov2003).

4.2. Inner region

In the diffusive regime, we scale into the contact line region by setting

(4.5a,b)\begin{equation} r = 1 - \varepsilon^{2}R, \quad m = \varepsilon^{{-}2}M, \end{equation}

where, as previously, the scaling for $m$ can be determined from the conservation of solute condition, (2.34). We can see immediately that the diffusive boundary layer is an order of magnitude thinner, while the solute mass profile is an order of magnitude larger for this evaporation model than in the kinetic evaporation model (in which the width and height are of $O(\varepsilon )$ and $O(\varepsilon ^{-1})$, respectively). This is consistent with the experimental results of Kajiya et al. (Reference Kajiya, Kaneko and Doi2008), in which droplets allowed to evaporate naturally (i.e. diffusively) produced thinner coffee rings than those constrained to evaporate in a box (i.e. closer to the kinetic regime).

Upon substituting the scalings (4.5a,b) into (2.31)–(2.33) and expanding the mass in an asymptotic series of the form $M = M_{0} + O(\varepsilon )$ as $\varepsilon \rightarrow 0$, we find to leading order that

(4.6)\begin{equation} \frac{\partial}{\partial R}\left[\left(\frac{1}{4\sqrt{2}(1-t)\sqrt{R}} - \frac{\rm \pi}{8R}\right)M_{0} + \frac{\rm \pi}{8}\frac{\partial M_{0}}{\partial R}\right] = 0 \quad \text{for} \quad R>0, \ t>0, \end{equation}

such that

(4.7)\begin{equation} \left(\frac{1}{4\sqrt{2}(1-t)\sqrt{R}} - \frac{\rm \pi}{8R}\right) M_{0} + \frac{\rm \pi}{8}\frac{\partial M_{0}}{\partial R} = 0 \quad \text{at} \ R=0, \ t>0. \end{equation}

The leading-order inner solute mass profile is therefore given by

(4.8)\begin{equation} M_{0} = F(t)R\text{exp}\left(-\frac{2\sqrt{2}}{{\rm \pi}(1-t)}\sqrt{R}\right), \end{equation}

where $F(t)$ is an arbitrary function of $t$.

4.3. Determining the unknown coefficient $F(t)$

Since the leading-order inner solution is exponentially small as $R\rightarrow \infty$, we again appeal to conservation of solute to find an expression for $F(t)$. Although the evaporative flux has changed, conservation of solute (2.34) must still hold. We can then determine $F(t)$ in a similar manner to the kinetic regime, see Appendix A; we find that

(4.9)\begin{equation} F(t) = \frac{16}{3{\rm \pi}^{4}}\frac{1}{(1-t)^{4}}\left[\frac{1}{4}\left(1-(1-t)^{3/4}\right)^{4/3}\right]. \end{equation}

We have written this in a slightly unusual manner in order to highlight that, as with $A(t)$ in the kinetic regime, the term in square brackets is the total solute mass advected into the contact line region up to time $t$, (4.4).

4.4. Composite expansion

We can write the leading-order composite mass solution for the diffusive evaporative flux model as

(4.10)\begin{equation} m_{comp} = m_{0}(r,t) + \frac{1}{\varepsilon^{2}}M_{0}\left(\frac{1-r}{\varepsilon^{2}},t\right), \end{equation}

where $m_{0}$ and $M_{0}$ are given by (4.3) and (4.8), respectively. To obtain a higher-order composite approximation for the solute mass in this regime requires tackling the second-order outer problem, which will typically need to be done numerically. However, as we shall see, the leading-order composite expansion (4.10) is in excellent agreement with our numerical simulations.

4.5. Profiles and properties of the coffee ring

With a diffusive evaporative model, the stronger evaporative flux at the contact line leads to a sharper, thinner solute profile in the boundary layer. In particular, as the modified Péclet number, ${\textit {Pe}}_{t}$, increases, we expect the nascent coffee ring to approach the similarity form given by

(4.11)\begin{equation} \frac{M_{0}(R,t)}{\mathcal{M}(t){\textit{Pe}}_{t}^{2}} = \frac{16R}{3{\rm \pi}^{4}}\text{e}^{{-}2\sqrt{2}\sqrt{R}/{\rm \pi}} = \frac{\sqrt{2}}{3{\rm \pi}}f\left(\sqrt{R};3,\frac{2\sqrt{2}}{\rm \pi}\right), \quad R = {\textit{Pe}}_{t}^{2}(1-r), \end{equation}

where $\mathcal {M}(t)$ is given by (4.4) and we again see the emergence of a gamma distribution probability density function, $f$, although with an increased shape parameter ($3$ rather than $2$ in the kinetic case).

To check the veracity of our asymptotic predictions, we again compare them with numerical simulations of the full system (2.31)–(2.33). The nature of the extremely thin boundary layer with diffusive evaporation – recall, an order of magnitude thinner than the kinetic regime – means more care has to be taken in the range of Péclet numbers considered to balance numerical convergence with remaining in the large-${\textit {Pe}}$ asymptotic regime. In figure 6, we consider both numerical solutions (solid curves) and the asymptotic predictions (dashed curves) of the solute mass profiles as $t$ increases up to 90 % of the drying time for ${\textit {Pe}} = 100$. Again, we can clearly see the formation of the peak in the solute profile that should ultimately become the coffee ring as the drop evaporates. We note that the stronger evaporative flux in this regime induces a much sharper increase in the coffee-ring height: in the inset to figure 6, we see that by $t = 0.9$, the solute mass is $\approx 7000$, which is two orders of magnitude larger than the height of the mass profile in the kinetic regime at the same stage of the dynamics. Nevertheless, we again see excellent agreement between the asymptotic predictions and the results of our numerical simulations.

Figure 6.

Profiles of the solute mass as the droplet evaporates under a diffusive flux given by (2.27) with ${\textit {Pe}} = 100$. In each figure, the bold, black curve represents the initial droplet solute profile (which mimics the initial droplet profile). Plots at time intervals of $0.1$ up to $t = 0.9$ show the results of numerical simulations (solid, blue curves) and the composite mass profile (dashed, red curves) predicted by our large-${\textit {Pe}}$ asymptotics (4.10). Panel (b) is a doubly logarithmic plot of the mass profile near the contact line, where we can clearly see the formation of the nascent coffee ring as $t$ increases.

In figure 7, we demonstrate the collapse of the nascent coffee ring to (4.11) as ${\textit {Pe}}_{t}$ increases. The convergence to the similarity profile is very rapid in the diffusive regime and gives us confidence, alongside the excellent comparisons of the composite mass profile in figure 6, in using our model to predict features of the nascent coffee ring.

Figure 7.

The scaled solute mass $m/\mathcal {M}(t){\textit {Pe}}_{t}^{2}$ against the radial similarity coordinate ${\textit {Pe}}_{t}^{2}(1-r)$ as a function of time. The dashed, black curve represents the scaled leading-order inner solute mass profile (4.11). The coloured curves are results from the numerical simulations at $t = 0.3$ (dark purple), $t = 0.6$ (blue) and $t = 0.9$ (light green) for ${\textit {Pe}} = 10$ (solid), ${\textit {Pe}} = 100$ (dots) and ${\textit {Pe}} = 200$ (dash-dots). We see the expected collapse to the similarity form (4.11) as ${\textit {Pe}}_{t}$ increases.

The maximum height, $m_{{max}}$, and its radial position, $r_{{max}}$, can again be well-approximated by considering the leading-order inner solution $M_{0}$, yielding

(4.12a,b)\begin{equation} m_{{max}} = \frac{8\mathcal{M}(t){\textit{Pe}}_{t}^{2}}{3{\rm \pi}^{2}\text{e}^{2}}, \quad r_{{max}} = 1 - \frac{{\rm \pi}^{2}}{2{\textit{Pe}}_{t}^{2}}. \end{equation}

As with the kinetic regime (3.24a,b), $m_{{max}}$ diverges and $r_{{max}}$ approaches the contact line as $t\rightarrow 1$, although for diffusive evaporation, these effects are much more pronounced. This can clearly be seen in figure 8, where we plot the numerical and asymptotic predictions for ${\textit {Pe}} = 50, 100$ and $200$. In comparison with, in particular, the ${\textit {Pe}} = 100$ curve in figure 5, we see that the concentration peak forms much sooner in the diffusive regime than in the kinetic regime due to the enhanced effects of evaporation. We can also see that, even compared with the ${\textit {Pe}} = 1000$ case in figure 5, the concentration peak is much closer to the contact line in this regime, driven by the stronger velocity profile.

Figure 8.

The maximum height of the solute mass profile $m_{{max}}$ (a), its location $1-r_{{max}}$ (b) and the full-width at half-maximum $w_{1/2}$ (c) as a function of time in the diffusive regime. In this figure we consider three different Péclet numbers, ${\textit {Pe}} = 50$ (dark purple), ${\textit {Pe}} = 100$ (blue), ${\textit {Pe}} = 200$ (light green). For each case, the dashed, black curve corresponds to the asymptotic predictions (4.12a,b) and (4.13).

The width of the concentration peak at half of its maximum height in the diffusive regime can be found analytically to be

(4.13)\begin{equation} w_{1/2} = \frac{{\rm \pi}^{2}}{2{\textit{Pe}}_{t}^{2}}\left[W_{{-}1}\left(\frac{-1}{\sqrt{2}\text{e}}\right)^{2} - W_{0}\left(\frac{-1}{\sqrt{2}\text{e}}\right)^{2}\right], \end{equation}

which, consistent with the behaviours of $m_{{max}}$ and $r_{{max}}$, shrinks more rapidly than its equivalent in the kinetic regime (3.25) as the droplet evaporates. We plot (4.13) alongside the corresponding numerical results in figure 8. Commensurate with the rapidity with which the nascent coffee ring forms in this regime, we are able to capture $w_{1/2}$ much sooner than in the kinetic regime, and we see that even for Péclet numbers as small as $50$, there is excellent agreement between the asymptotic prediction (4.13) and the numerical results.

Finally, we note that it is more difficult to display the competing fluxes in the diffusive evaporative regime than it was for the kinetic evaporative regime (cf. figure 4). This is primarily due to the fact that, while the composite asymptotic expansion (4.10) is asymptotically consistent for the solute mass, $m$, it loses this consistency if we attempt to differentiate it, which is necessary when computing $J_{{diff}} = -(h/{\textit {Pe}})\partial \phi /\partial r$. The inconsistency arises from the fact that, although the leading-order outer solute mass (4.3) is square root bounded at the contact line, $\partial \phi /\partial r \sim (1-r)^{-3/2}$ as $r\rightarrow 1$. Hence, $J_{{diff}}$ is inverse square root singular at the contact line. In order to address this singularity, we would need to proceed to higher order in both the inner and outer regions of our asymptotic analysis, which is a challenging procedure and beyond the scope of the present paper. However, as we shall now discuss, there are even more significant challenges facing the diffusive evaporative flux regime.