2. Problem configuration

We consider the configuration depicted in figure 1, in which an axisymmetric droplet of initial volume $V^{\ast}_{0}$ evaporates from a solid substrate. Here and hereafter, an asterisk denotes a dimensional variable. We let $(r^{\ast},\theta,z^{\ast})$ be cylindrical polar coordinates centred along the line of symmetry of the droplet with the substrate lying in the plane $z^{\ast} = 0$: by axisymmetry, we shall assume that all the variables are independent of $\theta$. The droplet contact line is thus circular and we assume that it is pinned throughout the drying process, which is observed in practice for a wide range of liquids for the majority of the drying time (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten1997; Hu & Larson Reference Hu and Larson2002; Kajiya, Kaneko & Doi Reference Kajiya, Kaneko and Doi2008; Howard et al. Reference Howard, Archer, Sibley, Southee and Wijayantha2023). We let $r^{\ast} = R^{\ast}$ be the radius of the contact line. Throughout this analysis, we shall assume that the droplet is thin, which reduces to the assumption that

(2.1)\begin{equation} 0<\delta = \frac{V^{\ast}_{0}}{R^{*3}} \ll1. \end{equation}

As we discuss presently, the thin-droplet assumption allows us to greatly simplify the flow and solute transport models; the assumption has been extensively validated and has shown to be reasonable even for droplets that should realistically fall outside of this regime (Larsson & Kumar Reference Larsson and Kumar2022).

Figure 1.

A side-on view of a solute-laden droplet evaporating under an evaporative flux $E^*(r^*)$ from a solid substrate that lies in the plane $z^* = 0$. The droplet is axisymmetric and the contact line is assumed to be pinned on the substrate at $r^* = R^*$. The droplet free surface is denoted by $h^*(r^*,t^*)$. The solute is assumed to be inert and sufficiently dilute that the flow of liquid in the droplet is decoupled from the solute transport.

The droplet consists of a liquid of constant density and viscosity denoted by $\rho ^{\ast}$ and $\mu ^{\ast}$, respectively. The droplet free surface is denoted by $z^{\ast} = h(r^{\ast},t^{\ast})$ and the air–water surface tension coefficient, $\sigma ^*$, is assumed to be constant.

The liquid evaporates into the surrounding gas and we assume that the evaporative process is quasi-steady, which is a reasonable assumption for a wide range of liquid–substrate configurations (Hu & Larson Reference Hu and Larson2002). We denote the evaporative flux by $E^*(r^*)$. The exact form of the flux will depend upon the dominant evaporative processes, which can change significantly with the properties of the droplet, ambient gas and substrate. Although it is not the goal of the present study to determine the correct form for $E^*(r^*)$, the differences do merit further discussion, which we pursue shortly in § 2.2.

The droplet contains an inert solute of initially uniform concentration $\phi ^{\ast}_{0}$. The solute is assumed to be sufficiently dilute that the flow and transport problems completely decouple. We shall discuss the validity of the dilute assumption further in § 6.

2.1. Flow model

The droplet is assumed to be sufficiently thin and the evaporation-induced flow sufficiently slow that the flow is governed by the lubrication equations

(2.2)$$\begin{gather} \frac{\partial h^{\ast}}{\partial t^{\ast}} + \frac{1}{r^{\ast}}\frac{\partial}{\partial r^{\ast}}(r^{\ast} h^{\ast}u^{\ast}) =- \frac{E^{\ast}}{\rho^{\ast}}, \end{gather}$$

(2.3)$$\begin{gather}u^{\ast} =- \frac{h^{*2}}{3\mu^{\ast}}\frac{\partial p^{\ast}}{\partial r^{\ast}}, \end{gather}$$

(2.4)$$\begin{gather}p^{\ast} = p_{atm}^{\ast}-\rho^{\ast}g^{\ast}(z^{\ast}-h^{\ast}) - \sigma^{\ast}\frac{1}{r^{\ast}}\frac{\partial}{\partial r^{\ast}}\left(r^{\ast}\frac{\partial h^{\ast}}{\partial r^{\ast}}\right), \end{gather}$$

for $0< r^{\ast}< R^{\ast}$, $t^{\ast}>0$, where $u^{\ast}(r^{\ast},t^{\ast})$ is the depth-averaged radial fluid velocity, $p^{\ast}(r^{\ast},z^{\ast},t^{\ast})$ is the liquid pressure and $p^{\ast}_{atm}$ denotes atmospheric pressure (Hocking Reference Hocking1983; Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000; Oliver et al. Reference Oliver, Whiteley, Saxton, Vella, Zubkov and King2015).

We note here that we assume throughout the analysis that the lubrication equations (2.2)–(2.4) remain applicable in the regions of interest, most notably in the region close to the contact line where solutal diffusion becomes important. This will introduce restrictions on the size of $\delta$ in cases where the effect of gravity dominates over surface tension, as we discuss in detail in § 3.

Equations (2.2)–(2.4) must be solved subject to the symmetry conditions

(2.5a,b)\begin{equation} r^{\ast} h^{\ast}u^{\ast} = \frac{\partial h^{\ast}}{\partial r^{\ast}} = 0 \,\quad \text{at} \ r^{\ast} = 0, \end{equation}

and the fact that the free surface touches down at, and we require no-flux of liquid through, the pinned contact line, that is,

(2.6a,b)\begin{equation} h^{\ast} = r^{\ast} h^{\ast}u^{\ast} = 0 \,\quad \text{at} \ r^{\ast} = R^{\ast}. \end{equation}

We close the problem by specifying the initial droplet profile, that is,

(2.7)\begin{equation} h^{\ast}(r^{\ast},0) = h^{\ast}_{0}(r^{\ast})\quad \text{for} \ 0< r^{\ast}< R^{\ast}. \end{equation}

It is worth noting at this stage that, while this initial condition is needed to fully specify the mathematical problem, in our analysis we do not explicitly use (2.7). In what follows, it is assumed that the rate of evaporation is sufficiently slow that the droplet quickly relaxes under capillary action to the quasi-steady profile found in § 3 (see, for example, Lacey Reference Lacey1982; De Gennes Reference De Gennes1985; Oliver et al. Reference Oliver, Whiteley, Saxton, Vella, Zubkov and King2015). Thus, we shall, for simplicity, assume that $h_{0}^*(r^*)$ is of the same functional form of the free surface we find in § 3. While this assumption is reasonable for a wide range of applications, for extremely rapid evaporation (for example, laser-induced evaporation, Volkov & Strizhak Reference Volkov and Strizhak2019), a more careful consideration of the evolution after deposition would be needed.

2.2. Evaporation model

As alluded to above, there are a number of different viable evaporation models depending on the physical and chemical characteristics of the problem (see, for example, Hu & Larson Reference Hu and Larson2002; Shahidzadeh-Bonn et al. Reference Shahidzadeh-Bonn, Rafai, Azouni and Bonn2006; Kelly-Zion et al. Reference Kelly-Zion, Pursell, Vaidya and Batra2011; Murisic & Kondic Reference Murisic and Kondic2011, and references therein). We shall outline a few of the more common cases.

A common evaporation model for small droplets (typically up to around a millimetre) with a pinned contact line is a diffusion-limited model, in which evaporation is limited by how quickly vapour is transported away from the droplet surface by diffusion and exhibits a singular flux at the contact line. Such a model has been shown to accurately predict both the total evaporation rate (Hu & Larson Reference Hu and Larson2002) and the pointwise evaporative flux (see, for example, Sáenz et al. Reference Sáenz, Wray, Che, Matar, Valluri, Kim and Sefiane2017; Wray & Moore Reference Wray and Moore2023) for a range of different droplet geometries.

When the ambient gas consists solely of vapour or when the diffusion of vapour is rapid, evaporation may instead be limited by kinetic effects at the liquid–vapour interface (see, for example, Murisic & Kondic Reference Murisic and Kondic2011; Jambon-Puillet et al. Reference Jambon-Puillet, Carrier, Shahidzadeh, Brutin, Eggers and Bonn2018). These effects are governed by the Hertz–Knudsen relation and may be shown for thin droplets to give an approximately constant evaporative flux. A constant flux may also be shown to be relevant in situations where a droplet evaporates above a hydrogel bath (Boulogne et al. Reference Boulogne, Ingremeau and Stone2016).

For larger droplets evaporating into a mixture of the liquid vapour and another gas, natural convection may be shown to play an important role in the evaporation rate for both pinned (Dollet & Boulogne Reference Dollet and Boulogne2017) and non-pinned (see, for example, Shahidzadeh-Bonn et al. Reference Shahidzadeh-Bonn, Rafai, Azouni and Bonn2006; Kelly-Zion et al. Reference Kelly-Zion, Pursell, Vaidya and Batra2011) droplets, although its role is greatly diminished when the density difference between the vapour and ambient gas is small (Radhakrishnan, Anand & Bakshi Reference Radhakrishnan, Anand and Bakshi2019). The relative importance of natural convection to vapour diffusion is measured by the Grashof number,

(2.8)\begin{equation} {Gr} = \left|\frac{(\rho_{s}^*-\rho_\infty^*)}{\rho_{\infty}^*}\right|\frac{g^*R^{*3}}{\nu_{air}^*}, \end{equation}

where $\rho ^*_s$ is the vapour saturation density at the droplet free surface, $\rho ^*_\infty$ is the ambient vapour density and $\nu _{\text {air}}^*$ is the kinematic viscosity of air (Kelly-Zion et al. Reference Kelly-Zion, Pursell, Vaidya and Batra2011; Dollet & Boulogne Reference Dollet and Boulogne2017). When ${Gr}\ll 1$, diffusion dominates the evaporation, while for moderate and large ${Gr}$, convection effects are non-negligible and often dominant.

For the particular example of a large ($R^* \geq 4$ mm) circular disk of water evaporating into ambient air, Dollet & Boulogne (Reference Dollet and Boulogne2017) show that the total evaporative flux $F^* = 2{\rm \pi} \int _{0}^{R^*} r^*E^*(r^*)\,dr^*$ may be approximated by the form

(2.9)\begin{equation} F^* \approx 2{\rm \pi} D^* R^* (c_s^*-c_\infty^*)\left[a_1{Gr}^{\beta} + a_2\right], \end{equation}

where $D^{\ast}$ is the vapour diffusion coefficient, $c^*_s$ is the saturation concentration at the droplet free surface and $c^*_\infty$ is the ambient vapour concentration. The first term in the brackets corresponds to the importance of convection, while the second term is akin to the importance of diffusion. The exponent $\beta$ is estimated analytically to be $1/5$, which is shown to be very close to experimental data, where $\beta \approx 0.18$ with $a_1 \approx 0.31$ and $a_2\approx 0.48$. They show that for ${Gr}\approx 20$ – corresponding to a droplet radius of $R^*\approx 5$ mm – the contribution of convection is around $50\,\%$ larger than diffusion. The difference grows wider as $R^*$ is increased further. Kelly-Zion et al. (Reference Kelly-Zion, Pursell, Vaidya and Batra2011) find a similar scaling law, although with slightly different coefficients, for a range of different evaporating liquids.

Once convection is significant, it becomes prohibitively difficult to find explicit expressions of the evaporative flux valid for all $r^*$. One can make scaling arguments based on boundary layer theory for large Grashof numbers, but these are unable to capture the flux near the droplet edge or at the centre of the droplet where a plume of vapour rises (Dollet & Boulogne Reference Dollet and Boulogne2017). Since these regions respectively play a crucial role in both the early time evolution of the coffee-ring (Moore et al. Reference Moore, Vella and Oliver2021, Reference Moore, Vella and Oliver2022) and the late time ‘fadeout’ profile of the ring (Witten Reference Witten2009), this severely hinders any prospect of analytic progress in studying solute transport in such regimes, and such approaches would necessarily be numerical.

It is not the goal of the present study to determine the correct model for evaporation for a given configuration and nor do we seek to cover the manifold possible models herein. Instead, our aim is to concentrate on the interplay between gravity, surface tension, solutal advection and solutal diffusion in solute transport, so here we shall choose an illustrative model for evaporation to use throughout the analysis, namely diffusive evaporation. Our choice is partly due to its predominance for smaller droplets (Hu & Larson Reference Hu and Larson2002) and for droplets evaporating in an ambient gas whose density is close to the vapour density (Radhakrishnan et al. Reference Radhakrishnan, Anand and Bakshi2019), but also due to access to an explicit form of the evaporative flux, as discussed shortly. While we recognize that there are regimes where this evaporation model will be lacking, the principles of the asymptotic analysis we pursue will be similar for a given flux – provided that we are in a regime where a coffee ring forms – but the sizes of the different regions may change, as is seen in the surface-tension-dominated regime discussed in Moore et al. (Reference Moore, Vella and Oliver2021, Reference Moore, Vella and Oliver2022) wherein both a diffusive and a kinetic evaporative flux are considered in detail. We shall return to this discussion further in § 6.

2.2.1. Diffusive evaporation

In a diffusive evaporation model, the vapour concentration $c^*(r^*,z^*)$ is determined by solving

(2.10)\begin{equation} \nabla^2 c^* = 0 \quad \text{in the gas phase}, \end{equation}

subject to

(2.11)\begin{equation} c^* = c^*_s \ \text{on the droplet free surface}, \quad \frac{\partial c^*}{\partial z^*} = 0\ \text{on the solid substrate}, \end{equation}

and such that

(2.12)\begin{equation} c^*\rightarrow c_\infty^* \quad \text{as} \ r^{*2}+z^{*2}\rightarrow \infty. \end{equation}

In the limit in which the droplet is thin, this boundary value problem may be linearized onto $z^* = 0$ and is thus equivalent to the classical problem for finding the electrostatic potential outside of a charged disk of radius $R^*$. The evaporative flux $E^*(r^*)$ may be shown to be given by

(2.13)\begin{equation} E^{\ast}(r^{\ast}) = \left.-D^* M_\ell^*\frac{\partial c^*}{\partial z^*}\right|_{z^* = 0} = \frac{2D^{\ast}M_\ell^*(c_{s}^{\ast}-c_{\infty}^{\ast})}{{\rm \pi}\sqrt{R^{*2}-r^{*2}}}, \end{equation}

where $M_\ell ^*$ is the molar mass of the liquid vapour (see, for example, Sneddon Reference Sneddon1966).

Assuming the contact line is pinned, the volume of the droplet $V^{\ast}(t^{\ast})$ is given by

(2.14)\begin{equation} V^{\ast}(t^{\ast}) = 2{\rm \pi}\int_{0}^{R^{\ast}}r^{\ast} h^{\ast}(r^{\ast},t^{\ast})\,\text{d}r^{\ast}, \quad V^{\ast}(0) = V^{\ast}_{0}. \end{equation}

The total mass loss due to evaporation $F^{\ast}$ is given by

(2.15)\begin{equation} F^{\ast} = 2{\rm \pi}\int_{0}^{R^{\ast}}r^{\ast} E^{\ast}(r^{\ast})\,\text{d}r^{\ast} = 4D^{\ast} M_\ell^* (c_{s}^{\ast}-c_{\infty}^{\ast})R^{\ast}. \end{equation}

Thus, conservation of mass in the liquid phase is

(2.16)\begin{equation} \frac{\text{d}V^{\ast}}{\text{d} t^{\ast}} =-\frac{F^{\ast}}{\rho^*} =-\frac{4D^{\ast} M_\ell^* (c_{s}^{\ast}-c_{\infty}^{\ast})R^{\ast}}{\rho^*} \end{equation}

so that

(2.17)\begin{equation} V^{\ast}(t^{\ast}) = V^{\ast}_{0} - \frac{4D^{\ast} M_\ell^* (c_{s}^{\ast}-c_{\infty}^{\ast})R^{\ast}t^{\ast}}{\rho^*}. \end{equation}

In particular, the dryout time, that is the time when the drop has fully evaporated, is

(2.18)\begin{equation} t_{f}^{\ast} = \frac{\rho^*V^{\ast}_{0}}{4D^{\ast} M_\ell^* (c_{s}^{\ast}-c_{\infty}^{\ast})R^{\ast}}. \end{equation}

2.3. Solute model

The droplet is assumed to be sufficiently thin that the transport of the solute is governed by the depth-averaged advection–diffusion equation

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

for $0< r^{\ast}< R^{\ast}$, $0< t^{\ast}< t^{\ast}_f$, where $\phi ^{\ast}(r^{\ast},t^{\ast})$ is the depth-averaged solute concentration and $D^{\ast}_{\phi }$ is the solutal diffusion coefficient (Wray et al. Reference Wray, Papageorgiou, Craster, Sefiane and Matar2014; Pham & Kumar Reference Pham and Kumar2017; Moore et al. Reference Moore, Vella and Oliver2021).

While there is an acknowledged effect of the solute particles eventually being trapped at and transported along the free surface (Maki & Kumar Reference Maki and Kumar2011; Kang et al. Reference Kang, Vandadi, Felske and Masoud2016; D'Ambrosio Reference D'Ambrosio2022), this effect is less pronounced for thin droplets, where the capture tends to occur closer to the contact line due to the stronger outward radial flow, and for droplets where Marangoni effects may be neglected. Thus, we shall neglect its effects here, as our study concerns the interplay between gravity, surface tension and solutal advection/diffusion. A more focused analysis on the final deposit profile would certainly need to account for such effects.

Equation (2.19) must be solved subject to the symmetry condition

(2.20)\begin{equation} \frac{\partial\phi^{\ast}}{\partial r^{\ast}} = 0 \,\quad \text{at} \ r^{\ast} = 0, \end{equation}

and the condition that there can be no flux of solute particles through the pinned contact line,

(2.21)\begin{equation} r^{\ast}\left(h^{\ast}u^{\ast}\phi^* - D^{\ast}_{\phi}h^{\ast}\frac{\partial\phi^{\ast}}{\partial r^{\ast}}\right) = 0 \,\quad \text{at} \ r^{\ast} = R^{\ast}. \end{equation}

Finally, we impose an initially uniform distribution of solute throughout the droplet, so that

(2.22)\begin{equation} \phi^{\ast}(r^{\ast},0) = \phi^{\ast}_{0} \quad \text{for} \ 0< r^{\ast}< R^{\ast}. \end{equation}

2.4. Non-dimensionalization

We assume that the fluid velocity is driven by evaporation and, for now, we retain both gravity and surface tension, so that the pertinent scalings are

(2.23)\begin{equation} \left. \begin{aligned} & (r^{\ast}, z^{\ast}) = R^{\ast}( r,\delta z), \quad u^{\ast} = \frac{D^{\ast}M_\ell^* (c_{s}^{\ast}-c_{\infty}^{\ast})}{\delta\rho^{\ast}R^{\ast}}u, \quad t^{\ast} = t_{f}^{\ast}t, \quad \phi^{\ast} = \phi^{\ast}_{0}\phi,\\ & (h^{\ast},h_{0}^{\ast}) = \delta R^{\ast} (h,h_{0}), \quad p^{\ast} = p_{atm}^{\ast} + \frac{\mu^{\ast}D^{\ast}M_\ell^* (c_{s}^{\ast}-c_{\infty}^{\ast})}{\delta^{3}\rho^{\ast} R^{*2}} p, \quad V^{\ast} = V_{0}^{\ast}V. \end{aligned} \right\} \end{equation}

Note, in particular, that the choice of time scale fixes the dimensionless dryout time to be $t = 1$.

Upon substituting the scalings (2.23) into (2.2)–(2.4), we see that

(2.24)$$\begin{gather} \frac{\partial h}{\partial t} + \frac{1}{4r}\frac{\partial}{\partial r}\left(rhu\right) =- \frac{1}{2{\rm \pi}\sqrt{1-r^{2}}}, \end{gather}$$

(2.25)$$\begin{gather}u = \frac{h^{2}}{3{Ca}}\frac{\partial}{\partial r}\left[-{Bo} h + \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial h}{\partial r}\right)\right], \end{gather}$$

for $0< r,t<1$, where the Capillary and Bond numbers are defined by

(2.26)\begin{equation} {Ca} = \frac{\mu^{\ast}D^{\ast}M_\ell^* (c_{s}^{\ast}-c_{\infty}^{\ast})}{\delta^{4}\rho^{\ast}R^{\ast}\sigma^{\ast}} \quad \text{and} \quad {Bo} = \frac{\rho^{\ast}g^{\ast}R^{*2}}{\sigma^{\ast}}, \end{equation}

respectively.

Under scalings (2.23), the symmetry conditions (2.5) become

(2.27a,b)\begin{equation} rhu = \frac{\partial h}{\partial r} = 0 \,\quad \text{at} \ r = 0, \end{equation}

while the contact line conditions (2.6) are

(2.28a,b)\begin{equation} h = rhu = 0 \,\quad \text{at} \ r = 1. \end{equation}

The initial condition (2.7) becomes

(2.29)\begin{equation} h(r,0) = h_{0}(r) \quad \text{for} \ 0< r<1. \end{equation}

Finally, the dimensionless form of the conservation of liquid volume conditions (2.14) and (2.17) is

(2.30)\begin{equation} 1-t = 2{\rm \pi}\int_{0}^{1}rh(r,t)\,\text{d}r. \end{equation}

After scaling, the solute transport equation (2.19) becomes

(2.31)\begin{equation} \frac{\partial}{\partial t}\left(h\phi\right) + \frac{1}{4r}\frac{\partial}{\partial r}\left[r\left(hu\phi - \frac{h}{{Pe}}\frac{\partial\phi}{\partial r}\right)\right] = 0 \end{equation}

for $0< r,t<1$, where the solutal Péclet number is

(2.32)\begin{equation} {Pe} = \frac{D^{\ast}M_\ell^*(c_{s}^{\ast}-c_{\infty}^{\ast})}{\delta\rho^{\ast}D_{\phi}^{\ast}}. \end{equation}

The symmetry (2.20) and boundary (2.21) conditions become

(2.33)\begin{equation} \frac{\partial \phi}{\partial r} = 0 \,\quad \text{at} \ r = 0 \end{equation}

and

(2.34)\begin{equation} r\left(hu\phi - \frac{h}{{Pe}}\frac{\partial \phi}{\partial r}\right) = 0 \,\quad \text{at} \ r = 1, \end{equation}

respectively. Finally, the initial condition (2.22) becomes

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

2.5. Integrated mass variable formulation

The assumption that the solute is dilute decouples the flow and solute transport problems, so that we may solve for $h$ and $u$ from (2.24)–(2.30) independently of the solute concentration, $\phi$. We shall discuss the resulting flow solution shortly in § 3.

First, however, we present a reformulation of the solute transport problem (2.31)–(2.35), which will greatly aid us in our asymptotic and numerical investigations. In this, we follow Moore et al. (Reference Moore, Vella and Oliver2021, Reference Moore, Vella and Oliver2022) by introducing the integrated mass variable

(2.36)\begin{equation} \mathcal{M}(r,t) = \int_{0}^{r}sh(s,t)\phi(s,t) \,\text{d}s. \end{equation}

By integrating the advection–diffusion equation (2.31) from $0$ to $r$ and applying the symmetry conditions (2.27a), (2.33), we find that

(2.37)\begin{equation} \frac{\partial\mathcal{M}}{\partial t} + \left[\frac{u}{4} + \frac{1}{4{Pe}}\left(\frac{1}{r} + \frac{1}{h}\frac{\partial h}{\partial r}\right)\right]\frac{\partial \mathcal{M}}{\partial r} - \frac{1}{4{Pe}}\frac{\partial^{2}\mathcal{M}}{\partial r^{2}} = 0 \quad \text{for} \ 0< r,t<1. \end{equation}

This must be solved subject to the boundary conditions

(2.38a,b)\begin{equation} \mathcal{M}(0,t) = 0, \quad \mathcal{M}(1,t) = \frac{1}{2{\rm \pi}}\quad \text{for} \ t>0, \end{equation}

where the latter condition dictates that mass is conserved along a radial ray, which replaces the no-flux condition (2.34). The initial condition (2.35) becomes

(2.39)\begin{equation} \mathcal{M}(r,0) = \int_{0}^{r}sh(s,0)\,\text{d}s \quad \text{for} \ 0< r<1. \end{equation}

Finally, we note that, once we have determined the integrated mass variable from (2.37)–(2.39), the solute mass $m = \phi h$ can then be retrieved from

(2.40)\begin{equation} m = \frac{1}{r}\frac{\partial\mathcal{M}}{\partial r}. \end{equation}

3. Flow solution in the small-${Ca}$ limit

We now suppose that surface tension dominates viscosity in the flow problem, that is, ${Ca} \ll 1$. Importantly, this means that the problems for the free surface profile and the flow velocity decouple, an assumption that is valid for a wide range of different liquids and evaporation models in practice (Moore et al. Reference Moore, Vella and Oliver2021, Reference Moore, Vella and Oliver2022). Unlike these previous studies, however, we shall retain gravity in (2.25) to investigate what role it plays in the formation of the nascent coffee ring.

To this end, we neglect the left-hand side of (2.25), so that upon integrating and applying the symmetry condition (2.27b), the contact line condition (2.28a) and the conservation of liquid volume condition (2.30), we deduce that

(3.1)\begin{equation} h(r,t) = \frac{(1-t)}{\rm \pi}\frac{{\rm I}_{0}(\sqrt{{Bo}})}{{\rm I}_{2}(\sqrt{{Bo}})}\left(1-\frac{{\rm I}_{0}(\sqrt{{Bo}}\,r)}{{\rm I}_{0}(\sqrt{{Bo}})}\right), \end{equation}

where $\textrm {I}_{\nu }(z)$ is the modified Bessel function of the first kind of order $\nu$. We note that this result has been reported previously for non-evaporating sessile drops by, for example, Allen (Reference Allen2003).

With the free surface found, the velocity is determined immediately from (2.24) and the no-flux condition (2.28b) to be

(3.2) \begin{align} u(r,t) &= \frac{1}{rh}\Biggl[\frac{2}{\rm \pi}\sqrt{1-r^{2}} + \frac{4{\rm I}_{0}(\sqrt{{Bo}})}{{\rm \pi} {\rm I}_{2}(\sqrt{{Bo}})}\Biggl(\frac{r^{2}-1}{2}\nonumber\\ &\quad + \frac{1}{\sqrt{{Bo}}{\rm I}_{0}(\sqrt{{Bo}})}({\rm I}_{1}(\sqrt{{Bo}})-rI_{1}(\sqrt{{Bo}}\,r))\Biggr)\Biggr]. \end{align}

Notably, as in the surface-tension-dominated regime where ${Bo} \rightarrow 0$, time is separable in both the free surface and fluid velocity profiles, and so merely acts to scale the functional form. In particular, this means that the streamlines and pathlines coincide, which we shall exploit when considering the regime in which solutal diffusion is negligible in § 5.3.

We display the scaled forms of the free surface and fluid velocity for various values of the Bond number in figure 2(a,b). For the droplet free surface profile, we see the expected transition from the spherical cap as ${Bo}\rightarrow 0$ (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000) to the flat ‘pancake’ (or ‘puddle’) droplet as ${Bo}\rightarrow \infty$ (Rienstra Reference Rienstra1990). For each Bond number, the velocity is singular at the contact line – as expected for a diffusive evaporative flux (as discussed in § 2.2). We see that as the effect of gravity increases, the sharp increase in $u$ occurs closer to the contact line, corresponding to the progressively smaller region in which surface tension effects are important.

Figure 2.

(a) The quasi-steady droplet free surface, (b) the fluid velocity and (c) the divergence of the velocity displayed for ${Bo} = 0.1$ (black), $1$ (dark purple), $10$ (blue), $20$ (cyan), $50$ (green) and $100$ (yellow). Notably, we see the transition from the spherical cap to the ‘pancake’ droplet profile as the effect of gravity increases. The divergence of the fluid velocity also shows a transition from a monotonic to a non-monotonic profile as the Bond number increases.

Finally, since this will be important in our discussions of the secondary peaks seen in the solute mass profile in § 5.3, we show the divergence of the fluid velocity in figure 2(c). For small Bond numbers, the divergence is monotonically increasing with $r$ and, as with the velocity, singular at the contact line. However, for moderate and large Bond numbers $\gtrsim 15$, we see a clear change of behaviour, with a region of non-monotonic behaviour in the droplet interior. This behaviour is accentuated as ${Bo}\rightarrow \infty$.

For future reference, the asymptotic behaviours of the free surface and fluid velocity as $r\rightarrow 1^-$ for ${Bo} = O(1)$ are given by

(3.3)$$\begin{gather} h = \theta_{c}(t;{Bo})(1-r) + O((1-r)^2), \end{gather}$$

(3.4)$$\begin{gather}u = \frac{2\chi}{\theta_{c}(t;{Bo})}(1-r)^{-1/2} + O((1-r)^{1/2}), \end{gather}$$

where

(3.5)\begin{equation} \theta_{c}(t;{Bo}) =-\lim_{r\rightarrow 1^-}\frac{\partial h}{\partial r} = (1-t)\psi(Bo), \quad \psi({Bo}) = \frac{\sqrt{{Bo}}{\rm I}_{1}(\sqrt{{Bo}})}{{\rm \pi} {\rm I}_{2}(\sqrt{{Bo}})} \end{equation}

is the leading-order contact angle in the thin-droplet limit – again, as reported by Allen (Reference Allen2003) – and

(3.6)\begin{equation} \chi = \frac{\sqrt{2}}{\rm \pi} \end{equation}

is the dimensionless coefficient of the inverse square root singularity at the contact line in the evaporative flux (2.13). Note that we have chosen this notation to highlight the similarities with the previous analysis of Moore et al. (Reference Moore, Vella and Oliver2022), who consider a surface-tension-dominated droplet of arbitrary contact set.

On the other hand, if we take $1-r = O(1)$ and consider the large-${Bo}$ limit of (3.1), (3.2), we find that

(3.7)$$\begin{gather} h = h_{0}(t) + {Bo}^{-1/2}h_{1}(t) + O({Bo}^{-1}), \end{gather}$$

(3.8)$$\begin{gather}u = u_{0}(r,t) + {Bo}^{-1/2}u_{1}(r,t) + O({Bo}^{-1}), \end{gather}$$

as ${Bo}\rightarrow \infty$, where

(3.9a,b)\begin{equation} h_{0}(t) = \frac{(1-t)}{\rm \pi}, \quad h_{1}(t) = \frac{2(1-t)}{\rm \pi}, \end{equation}

and

(3.10a,b)\begin{equation} u_{0}(r,t) = \frac{2\sqrt{1-r^{2}}}{r(1-t)}(1-\sqrt{1-r^{2}}), \quad u_{1}(r,t) = \frac{4}{r(1-t)}(1-\sqrt{1-r^{2}}). \end{equation}

Notably, in the droplet bulk, the droplet free surface $h$ is flat to all orders: the aforementioned characteristic of ‘pancake’ droplets associated with large Bond numbers (Rienstra Reference Rienstra1990). These expansions break down close to the contact line where surface tension effects become important. We find that for $1-r = {Bo}^{-1/2}\bar {r}$, we have

(3.11)$$\begin{gather} h = \bar{h}_{0}(\bar{r},t) + {Bo}^{-1/2}\bar{h}_{1}(\bar{r},t) + O({Bo}^{-1}), \end{gather}$$

(3.12)$$\begin{gather}u = {Bo}^{-1/4}\left[\bar{u}_{0}(\bar{r},t) + {Bo}^{-1/4}\bar{u}_{1}(\bar{r},t) + {Bo}^{-1/2}\bar{u}_{2}(\bar{r},t) + O({Bo}^{-3/4})\right] \end{gather}$$

as ${Bo}\rightarrow \infty$, where

(3.13)$$\begin{gather} \bar{h}_{0}(\bar{r},t) = \frac{(1-t)}{\rm \pi}(1-\text{e}^{-\bar{r}}), \end{gather}$$

(3.14)$$\begin{gather}\bar{h}_{1}(\bar{r},t) = \frac{(1-t)}{2{\rm \pi}}(4(1-\text{e}^{-\bar{r}}) -\bar{r}\text{e}^{-\bar{r}}), \end{gather}$$

and

(3.15)$$\begin{gather} \bar{u}_{0}(\bar{r},t) = \frac{2\sqrt{2\bar{r}}}{(1-t)(1-\text{e}^{-\bar{r}})}, \end{gather}$$

(3.16)$$\begin{gather}\bar{u}_{1}(\bar{r},t) = \frac{4}{(1-t)}\left(1 - \frac{\bar{r}}{(1-\text{e}^{-\bar{r}})}\right), \end{gather}$$

(3.17)$$\begin{gather}\bar{u}_{2}(\bar{r},t) = \frac{\sqrt{\bar{r}}}{\sqrt{2}(1-t)(1-\text{e}^{-\bar{r}})}\left(3\bar{r} - 8 + \frac{2\bar{r}\text{e}^{-\bar{r}}}{(1-\text{e}^{-\bar{r}})}\right). \end{gather}$$

We note here that as $\bar {r}\rightarrow 0$, we retrieve the expected inverse square root singularity in the fluid velocity.

It is worth stressing here that the above expansions when the Bond number is large assume that we remain in a regime where the lubrication approximation (2.24)–(2.25) is still valid. Notably, this requires the restriction that variations in the free surface are small when $1-r = {Bo}^{-1/2}\bar {r}$, which holds provided

(3.18)\begin{equation} {Bo} \ll \frac{1}{\delta^{2}}. \end{equation}

Throughout our analysis, we shall assume that $\delta$ is sufficiently small that (3.18) holds. In problems where ${Bo} = O(\delta ^{-2})$, we would need to reformulate the analysis to include vertical variation in the Stokes equations close to the contact line, which is likely to significantly change the resulting behaviours in the large-${Bo}$ asymptotic regime. This is beyond the scope of the present study, so we do not pursue this further here.

4. Solute transport in the large-${Pe}$ limit

Having fully determined the leading-order flow, we now seek to understand the transport of solute within the drop and to make predictions about the early stages of coffee-ring formation. We follow the analyses of Moore et al. (Reference Moore, Vella and Oliver2021, Reference Moore, Vella and Oliver2022) by considering the physically relevant regime in which ${Pe} \gg 1$. In this regime, in the bulk of the droplet, advection dominates solutal diffusion, with the latter only being relevant close to the contact line.

Previous studies of this problem have concentrated on surface-tension-dominated drops (i.e. ${Bo}\rightarrow 0$) and have shown how the competition between solutal advection and diffusion near the contact line leads to the early stages of coffee-ring formation in drying droplets. In this analysis, we wish to investigate how this behaviour changes as we allow ${Bo}$ to vary, which we pursue using a hybrid asymptotic-numerical approach.

There are two main asymptotic regimes depending on the relative sizes of ${Bo}$ and ${Pe}$:

1. (i) moderate Bond number, ${Bo} = O(1)$ as ${Pe}\rightarrow \infty$, where the asymptotic structure of the solute transport depends solely on the large Péclet number;

2. (ii) large Bond number, ${Bo} = O({Pe}^{4/3})$ as ${Pe}\rightarrow \infty$, where the asymptotic structure of the solute transport now depends on the relative sizes of ${Bo}$ and ${Pe}$.

Note that there are further possibilities when the Bond number is large, namely, $1\ll {Bo}\ll {Pe}^{4/3}$ or ${Bo}\gg {Pe}^{4/3}$ as ${Pe}\rightarrow \infty$. These may be approached by taking the appropriate limit in the second regime, as we shall see presently.

In the first regime where ${Bo} = O(1)$ as ${Pe}\rightarrow \infty$, the asymptotic structure of the flow is a natural extension of the surface-tension-dominated case considered in Moore et al. (Reference Moore, Vella and Oliver2021). In the droplet bulk where $1-r = O(1)$, solute advection dominates diffusion. However, close to the contact line, a balance between solutal advection and diffusion occurs when

(4.1)\begin{equation} rhu\phi \sim \frac{rh}{{Pe}}\frac{\partial\phi}{\partial r} \implies 1-r = O({Pe}^{-2}). \end{equation}

We discuss the asymptotic solution for this regime in § 4.1.

In the second regime, the relative sizes of the boundary layer where surface tension enters the flow profile and the solutal diffusion boundary layer are comparable. As detailed in § 3, for a large Bond number the free surface is flat in the bulk of the droplet, with the effect of surface tension restricted to a boundary layer at the contact line of size $1 - r = O({Bo}^{-1/2})$, where $h = O(1)$ and $u = O({Bo}^{-1/4})$. Turning to the solute transport equation (2.31), since $h$ is order unity and $u$ is square root bounded in this region, advection and diffusion are comparable when

(4.2)\begin{equation} 1 - r = O({Pe}^{-2/3}). \end{equation}

Hence, the size of the two boundary layers are comparable when ${Bo} = O({Pe}^{4/3})$ as ${Pe}\rightarrow \infty$. We introduce the parameter

(4.3)\begin{equation} \alpha = {Bo}^{-1/2}{Pe}^{2/3}, \end{equation}

and note that $\alpha = O(1)$ as ${Pe}\rightarrow \infty$ in this regime. The asymptotic analysis in this regime is somewhat more involved and is presented in § 4.2.

4.1. Asymptotic solution when ${Bo} = O(1)$ as ${Pe}\rightarrow \infty$

In this section, we present the asymptotic solution of the solute transport problem when ${Bo} = O(1)$ as ${Pe}\rightarrow \infty$. The analysis herein is a natural extension of Moore et al. (Reference Moore, Vella and Oliver2021). For the purposes of this section, we shall use the concentration form of the advection–diffusion equation (2.31)–(2.35) and, in particular, find the solution in terms of the solute mass $m = \phi h$, where $h$ is given by (3.1).

4.1.1. Outer region

In the droplet bulk away from the contact line, we seek a solution of the form $m = m_{0} + O({Pe}^{-1})$ as ${Pe}\rightarrow \infty$. Substituting into (2.31), (2.35), we find that

(4.4)\begin{equation} \frac{\partial m_{0}}{\partial t} + \frac{1}{4r}\frac{\partial}{\partial r}(rm_{0}u) = 0 \quad \text{for} \ 0< r,t<1, \end{equation}

where $u$ is given by (3.2), subject to $m_0(r,0) = h(r,0)$. This is the usual advection equation, with solution given by

(4.5)\begin{equation} m_{0}(r,t) = \frac{h(R,0)}{J(R,t)}, \end{equation}

where $R$ is the initial location of the point that is at $r$ at time $t$ and $J(R,t)$ is the Jacobian of the Eulerian–Lagrangian transformation, which satisfies Euler's identity,

(4.6)\begin{equation} \frac{\text{D}}{\text{D}t}(\log{J}) = \frac{1}{4r}\frac{\partial}{\partial r}(ru), \quad J(R,0) = 1, \end{equation}

where $\textrm {D}/\textrm {D}t$ is the convective derivative (see, for example, Moore et al. Reference Moore, Vella and Oliver2022).

A straightforward asymptotic analysis of (4.4) reveals that

(4.7)\begin{equation} u\frac{\partial m_0}{\partial r} \sim \frac{m_0}{r}\frac{\partial }{\partial r}(ru) \end{equation}

as $r\rightarrow 1^-$, so that $m_{0} = O(\sqrt {1-r})$ as $r\rightarrow 1^-$, and hence, the concentration $\phi _0$ is square root singular. This sharp local concentration increase necessitates the inclusion of a diffusive boundary layer.

4.1.2. Inner region

Close to the contact line, we set

(4.8)\begin{equation} r = 1 - {Pe}^{-2}\hat{r}, \quad h = {Pe}^{-2}\hat{h}, \quad u = {Pe}\hat{u}, \quad m = {Pe}^{2}\hat{m}, \end{equation}

where the last scaling on the mass comes from global conservation of solute considerations (Moore et al. Reference Moore, Vella and Oliver2021). We seek an asymptotic solution of the form $\hat {m} = \hat {m}_{0} + O({Pe}^{-1})$ as ${Pe}\rightarrow \infty$ and find to leading order

(4.9)\begin{equation} \frac{\partial}{\partial \hat{r}}\left[\left( \frac{2\chi}{\theta_{c}(t;{Bo})\sqrt{\hat{r}}} - \frac{1}{\hat{r}}\right)\hat{m}_{0} + \frac{\partial\hat{m}_{0}}{\partial\hat{r}}\right] = 0 \quad\ \text{in}\ \hat{r}>0, \, 0< t<1, \end{equation}

such that

(4.10)\begin{equation} \left( \frac{2\chi}{\theta_{c}(t;{Bo})\sqrt{\hat{r}}} - \frac{1}{\hat{r}}\right)\hat{m}_{0} + \frac{\partial\hat{m}_{0}}{\partial\hat{r}} = 0 \quad \text{on} \ \hat{r} = 0. \end{equation}

It is straightforward to show that the solution to (4.9)–(4.10) is given by

(4.11)\begin{equation} \hat{m}_{0}(\hat{r},t) = C(t;{Bo})\hat{r}\exp\left(-\frac{4\chi}{\theta_{c}(t;{Bo})}\sqrt{\hat{r}}\right), \end{equation}

where, by pursuing a similar matching process to Moore et al. (Reference Moore, Vella and Oliver2022), we find that the coefficient $C(t;{Bo})$ is given by

(4.12)\begin{equation} C(t;{Bo}) = \frac{64\chi^4}{3\theta_{c}(t;{Bo})^4}\mathcal{N}(t;{Bo}), \end{equation}

where $\mathcal {N}(t;{Bo})$ is the leading-order accumulated mass advected into the contact line region up to time $t$, viz.

(4.13)\begin{equation} \mathcal{N}(t;{Bo}) = \tfrac{1}{4}\int_{0}^{t} m_{0}(1^-,\tau)u( 1^-,\tau)\,\text{d}\tau. \end{equation}

It is worth noting that this solution follows directly from the ${Bo} = 0$ regime discussed in Moore et al. (Reference Moore, Vella and Oliver2021, Reference Moore, Vella and Oliver2022), with the alterations due to gravity entering into the accumulated mass flux into the contact line and the leading-order contact angle. In particular, we note that in the limit ${Bo}\rightarrow 0$, since $\psi = 4/{\rm \pi} + O({Bo})$, this yields the expected form found in the surface-tension-dominated problem in Moore et al. (Reference Moore, Vella and Oliver2022) (see § 3.7.2 therein). We display the accumulated mass flux and the local contact angle for a wide range of Bond numbers in figure 3. We see that as the influence of gravity increases, the accumulated mass flux into the contact line at a fixed percentage of the evaporation time is reduced from the surface-tension-dominated regime. On the other hand, the local contact angle increases, commensurate with the droplet profile transitioning from a spherical cap to a ‘pancake’ droplet. We note that this combined behaviour leads to $C(t;{Bo})$ decreasing as ${Bo}$ increases. We discuss how these findings impact coffee-ring formation in more detail in § 5.2.1.

Figure 3.

(a) The accumulated mass flux, $\mathcal {N}(t;{Bo})$, as defined by (4.13) and (b) the leading-order local contact angle $\theta _{c}(t;{Bo})$ as defined by (3.5), for ${Bo} = 10^{-2}$ (purple), ${Bo} = 10^{-1}$ (purple) (dark blue), ${Bo} = 1$ (light blue), ${Bo} = 10$ (green) and ${Bo} = 10^{2}$ (yellow).

4.1.3. Composite solution

We may use van Dyke's rule (Van Dyke Reference Van Dyke1964) to formulate a leading-order composite solution for the solute mass that is valid throughout the drop by combining the leading-order-outer solution (4.5) and the leading-order-inner solution (4.11), finding

(4.14)\begin{equation} m_{comp}(r,t) = m_{0}(r,t) + {Pe}^{2}\hat{m}_0\left({Pe}^{2}(1-r),t\right). \end{equation}

4.2. Asymptotic solution when ${Bo} = O({Pe}^{4/3})$ as ${Pe}\rightarrow \infty$

In this section, we present the asymptotic solution of the solute transport problem in the limit in which ${Bo} = O({Pe}^{4/3})$ as ${Pe}\rightarrow \infty$ so that $\alpha$ defined by (4.3) is order unity. For convenience, we choose to use ${Pe}^{-2/3}$ as our small parameter in the asymptotic expansions. Moreover, it transpires that it is easier to analyse the integrated mass variable formulation of the problem (2.37)–(2.40).

4.2.1. Outer region

In the droplet bulk, we recall from (3.9) that the droplet free surface $h$ is flat to all orders and that the velocity is given by (3.10). Upon substituting these expressions into (2.37) and (2.39), and then expanding $\mathcal {M} = \mathcal {M}_{0} + {Pe}^{-2/3}\mathcal {M}_{1} + O({Pe}^{-4/3})$ as ${Pe}\rightarrow \infty$, we find to leading order

(4.15a,b)\begin{equation} \frac{\partial\mathcal{M}_{0}}{\partial t} + \frac{u_{0}}{4}\frac{\partial\mathcal{M}_{0}}{\partial r}= 0 \quad \text{for} \ 0< r,t<1, \quad \mathcal{M}_{0}(r,0) = \frac{r^{2}}{2{\rm \pi}}\quad \text{for} \ 0< r<1. \end{equation}

This may be solved using the method of characteristics, yielding

(4.16)\begin{equation} \mathcal{M}_{0}(r,t) = \frac{(1-t)r^{2}}{2{\rm \pi}} + \frac{\sqrt{1-t}(1-\sqrt{1-t})}{\rm \pi}(1-\sqrt{1-r^{2}}). \end{equation}

We see that this solution automatically satisfies the boundary condition (2.38a).

At $O({Pe}^{-2/3})$, the problem for $\mathcal {M}_{1}(r,t)$ is given by

(4.17a,b)\begin{align} \frac{\partial\mathcal{M}_{1}}{\partial t} + \frac{u_{0}}{4}\frac{\partial\mathcal{M}_{1}}{\partial r}=-\frac{\alpha u_{1}}{4}\frac{\partial\mathcal{M}_{0}}{\partial r} \quad \text{for}\ 0< r,t<1, \quad \mathcal{M}_{1}(r,0) = \frac{\alpha r^2}{\rm \pi} \quad \text{for} \ 0< r<1. \end{align}

This may be solved in a similar manner, yielding

(4.18)\begin{align} \mathcal{M}_{1}(r,t) = \frac{2\alpha\kappa(r,t)}{\rm \pi}(1-\kappa(r,t))\log\left(\frac{\sqrt{1-t}-\kappa(r,t)}{1-\kappa(r,t)}\right) + \frac{\alpha}{\rm \pi}(1-(1-\kappa(r,t))^{2}), \end{align}

where $\kappa (r,t) = \sqrt {1-t}(1-\sqrt {1-r^{2}})$.

Expanding the leading-order solution (4.16) as we approach the contact line, we have

(4.19) \begin{align} \mathcal{M}_{0}(r,t) &\sim \frac{\sqrt{1-t}}{\rm \pi} - \frac{(1-t)}{2{\rm \pi}} - \frac{\sqrt{2(1-t)}}{\rm \pi}(1-\sqrt{1-t})\sqrt{1-r}\notag\\ &\quad - \frac{(1-t)}{\rm \pi}(1-r) + O((1-r)^{3/2}) \end{align}

as $r\rightarrow 1^-$. A similar expansion of (4.18), yields

(4.20) \begin{align} \mathcal{M}_{1}(r,t) &\sim \frac{\alpha\sqrt{1-t}(1-\sqrt{1-t})}{\rm \pi}\log(1-r)\notag\\ &\quad + \frac{\alpha\sqrt{1-t}(1-\sqrt{1-t})}{\rm \pi}\log\left(\frac{2(1-t)}{(1-\sqrt{1-t})^{2}}\right)\nonumber\\ &\quad + \frac{\alpha}{\rm \pi}(1-(1-\sqrt{1-t})^{2}) + O(\sqrt{1-r}\log(1-r)) \end{align}

as $r\rightarrow 1^-$. We can clearly see this will necessitate an inner expansion that contains logarithmic terms; a similar behaviour is displayed for surface-tension-dominated drops under different evaporative fluxes (Moore et al. Reference Moore, Vella and Oliver2021).

Finally, if we expand the solute mass $m = m_{0} + O({Pe}^{-2/3})$ as ${Pe}\rightarrow \infty$ in (2.40), we find that

(4.21)\begin{equation} m_{0}(r,t) = \frac{\sqrt{1-t}}{{\rm \pi}\sqrt{1-r^{2}}}\left[1-\sqrt{1-t}(1-\sqrt{1-r^{2}})\right]. \end{equation}

Notably, this means that the leading-order-outer solute mass $m_{0}$ is singular at the contact line, which gives a strong indication of the importance of diffusive effects local to the edge of the droplet. This is in stark contrast to the ${Bo} = O(1)$ solution, where the outer solute mass was square root bounded as $r\rightarrow 1^-$. Lastly, we note that whilst we could proceed to $O({Pe}^{-2/3})$ in the solute mass expansion in the outer region, we shall not require this when constructing a composite profile that is valid to $O(1)$ throughout the droplet, so we do not present this here.

4.2.2. Inner region

Recalling (3.11)–(3.12), (4.2) and (4.3), in order to retain a balance between the advective and diffusive effects in (2.37) close to the contact line, we set

(4.22)\begin{equation} r = 1 - {Pe}^{-2/3} \tilde{r}, \quad u = {Pe}^{-1/3} \tilde{u}, \quad h = \tilde{h}, \quad \mathcal{M} = \tilde{\mathcal{M}}, \quad m = {Pe}^{2/3}\tilde{m} \end{equation}

in (2.37)–(2.40). Note that we therefore have

(4.23)\begin{equation} \tilde{h} = \tilde{h}_0 + {Pe}^{-2/3}\tilde{h}_1 + O({Pe}^{-4/3}), \quad \tilde{u} = \tilde{u}_0 + {Pe}^{-1/3}\tilde{u}_1 + {Pe}^{-2/3}\tilde{u}_2 + O({Pe}^{-1}) \end{equation}

as ${Pe}\rightarrow \infty$ where

(4.24)\begin{equation} \tilde{h}_0(\tilde{r},t) = \bar{h}_0(\tilde{r}/\alpha,t), \quad \tilde{h}_1(\tilde{r},t) = \alpha\bar{h}_1(\tilde{r}/\alpha,t) \end{equation}

and $\bar {h}_{0}$, $\bar {h}_{1}$ are given by (3.13)–(3.14), and

(4.25)\begin{equation} \tilde{u}_0(\tilde{r},t) = \sqrt{\alpha}\bar{u}_{0}(\tilde{r}/\alpha,t), \quad \tilde{u}_1(\tilde{r},t) = \alpha\bar{u}_{1}(\tilde{r}/\alpha,t), \quad \tilde{u}_2(\tilde{r},t) = \alpha^{3/2}\overline{u_{2}}(\tilde{r}/\alpha,t) \end{equation}

and $\bar {u}_{0}$, $\bar {u}_{1}$, $\bar {u}_{2}$ are given by (3.15)–(3.17).

Seeking an asymptotic expansion of the integrated mass of the form $\tilde {\mathcal {M}} = \tilde {\mathcal {M}}_{0} + {Pe}^{-1/3}\tilde {\mathcal {M}}_{1} + {Pe}^{-2/3}\log {Pe}^{-2/3}\tilde {\mathcal {M}}_{2} + {Pe}^{-2/3}\tilde {\mathcal {M}}_{3} + o({Pe}^{-2/3})$ as ${Pe}\rightarrow \infty$, we find that the leading-order-inner problem is given by

(4.26)\begin{equation} \frac{\partial^{2}\tilde{\mathcal{M}}_{0}}{\partial \tilde{r}^{2}} + \left(\tilde{u}_{0} - \frac{1}{\tilde{h}_0}\frac{\partial \tilde{h}_0}{\partial \tilde{r}}\right)\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial r} = 0\quad \text{for} \ \tilde{r}>0, 0< t<1, \end{equation}

subject to the boundary condition $\tilde {\mathcal {M}}_{0}(0,t) = 1/2{\rm \pi}$ for $0< t<1$ and, in order to match with the local expansion of the leading-order-outer solution at the contact line (4.19), we must have

(4.27)\begin{equation} \tilde{\mathcal{M}}_{0}\rightarrow \frac{\sqrt{1-t}}{\rm \pi} - \frac{(1-t)}{2{\rm \pi}} \quad \text{as}\ \tilde{r}\rightarrow\infty. \end{equation}

Defining the integrating factor

(4.28)\begin{equation} I(\tilde{r},t) = \left(\frac{1}{1-\text{e}^{-\tilde{r}/\alpha}}\right)\exp\left(\frac{2\sqrt{2}}{(1-t)}\int_{0}^{\tilde{r}}\frac{\sqrt{\xi}}{1-\text{e}^{-\xi/\alpha}}\,\text{d}\xi\right), \end{equation}

we find that the solution is given by

(4.29)\begin{equation} \tilde{\mathcal{M}}_{0}(\tilde{r},t) = \frac{1}{2{\rm \pi}} + B_{0}(t)\int_{0}^{\tilde{r}}\frac{1}{I(s,t)}\,\text{d}s, \end{equation}

where

(4.30)\begin{equation} B_{0}(t) =-\frac{1}{\rm \pi}\left(1-\sqrt{1-t}-\frac{t}{2}\right)\left(\int_{0}^{\infty} \frac{1}{I(s,t)}\,\text{d}s\right)^{-1}. \end{equation}

We note here that the first term on the right-hand side of $B_{0}(t)$ is simply the leading-order accumulated mass at the contact line as a function of time, $\mathcal {N}(t)$, that is,

(4.31)\begin{equation} \mathcal{N}(t) = \frac{1}{4}\int_{0}^{t} (m_{0} u_{0})(1^-,\tau)\,\text{d}\tau = \frac{1}{\rm \pi}\left(1 - \sqrt{1-t} - \frac{t}{2}\right). \end{equation}

It is worth noting the similarities between (4.31) and the equivalent expression for a surface-tension-dominated drop evaporating under a constant evaporative flux (Freed-Brown Reference Freed-Brown2015; Moore et al. Reference Moore, Vella and Oliver2021).

At $O({Pe}^{-1/3})$, we have

(4.32)\begin{equation} \frac{\partial^{2}\tilde{\mathcal{M}}_{1}}{\partial \tilde{r}^{2}} + \left(\tilde{u}_{0} - \frac{1}{\tilde{h}_0}\frac{\partial \tilde{h}_0}{\partial \tilde{r}}\right)\frac{\partial\tilde{\mathcal{M}}_{1}}{\partial r} = 4\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial t} - \tilde{u}_{1}\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial\tilde{r}} \quad \text{for} \ \tilde{r}>0, 0< t<1, \end{equation}

subject to $\tilde {\mathcal {M}}_{1}(0,t) = 0$ for $0< t<1$ and the far-field matching condition

(4.33)\begin{equation} \tilde{\mathcal{M}}_{1}\rightarrow -\frac{\sqrt{2(1-t)}}{\rm \pi}(1-\sqrt{1-t})\sqrt{\tilde{r}} \quad \text{as}\ \tilde{r}\rightarrow\infty. \end{equation}

While in practice it may be easier to find $\tilde {\mathcal {M}}_{1}(\tilde {r},t)$ from (4.32)–(4.33) numerically, for posterity, we state that this boundary value problem has the solution

(4.34)\begin{equation} \tilde{\mathcal{M}}_{1}(\tilde{r},t) = \int_{0}^{\tilde{r}}\frac{1}{I(s,t)} \left(\int_{0}^{s}\left(4\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial t} - \tilde{u}_{1}\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial \tilde{r}}\right)I(\sigma,t)\,\text{d}\sigma\right)\,\text{d}s + B_{1}(t)\int_{0}^{\tilde{r}}\frac{1}{I(s,t)}\,\text{d}s, \end{equation}

where

(4.35)\begin{align} B_{1}(t) &=-\left(\int_{0}^{\infty}\left\{\frac{1}{I(s,t)} \left(\int_{0}^{s}\left(\frac{4\partial\tilde{\mathcal{M}}_{0}}{\partial t} - \tilde{u}_{1}\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial \tilde{r}}\right)I(\sigma,t)\,\text{d}\sigma\right)- \sqrt{2}(1-t)\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial t}\frac{1}{\sqrt{s}}\right\}\,\text{d}s \right)\nonumber\\ &\quad\times \left(\int_{0}^{\infty}\frac{1}{I(s,t)}\,\text{d}s\right)^{-1} \end{align}

is chosen to remove the $O(1)$ term in the far-field expansion of $\tilde {\mathcal {M}}_{1}(\tilde {r},t)$.

The $O({Pe}^{-2/3}\log {Pe}^{-2/3})$ problem is given by

(4.36)\begin{equation} \frac{\partial^{2}\tilde{\mathcal{M}}_{2}}{\partial \tilde{r}^{2}} + \left(\tilde{u}_{0} - \frac{1}{\tilde{h}_0}\frac{\partial \tilde{h}_0}{\partial \tilde{r}}\right)\frac{\partial\tilde{\mathcal{M}}_{2}}{\partial r} = 0 \quad \text{for} \ \tilde{r}>0, 0< t<1, \end{equation}

subject to $\tilde {\mathcal {M}}_{2}(0,t) = 0$ for $0< t<1$ and the far-field matching condition

(4.37)\begin{equation} \tilde{\mathcal{M}}_{2}\rightarrow \frac{\alpha\sqrt{1-t}(1-\sqrt{1-t})}{\rm \pi} \quad \text{as}\ \tilde{r}\rightarrow\infty. \end{equation}

The solution may be found in a similar manner to the leading-order problem, yielding

(4.38)\begin{equation} \tilde{\mathcal{M}}_{2}(\tilde{r},t) = B_{2}(t)\int_{0}^{\tilde{r}}\frac{1}{I(s,t)}\,\text{d}s, \end{equation}

where

(4.39)\begin{equation} B_{2}(t) = \frac{\alpha\sqrt{1-t}(1-\sqrt{1-t})}{\rm \pi}\left(\int_{0}^{\infty}\frac{1}{I(s,t)}\,\text{d}s\right)^{-1}. \end{equation}

Lastly, at $O({Pe}^{-2/3})$, we have

(4.40)$$\begin{gather} \frac{\partial^{2}\tilde{\mathcal{M}}_{3}}{\partial \tilde{r}^{2}} + \left(\tilde{u}_{0} - \frac{1}{\tilde{h}_{0}}\frac{\partial \tilde{h}_{0}}{\partial \tilde{r}}\right)\frac{\partial\tilde{\mathcal{M}}_{3}}{\partial r} = 4\frac{\partial\tilde{\mathcal{M}}_{1}}{\partial t} - \tilde{u}_{1}\frac{\partial\tilde{\mathcal{M}}_{1}}{\partial\tilde{r}} - \tilde{u}_{2}\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial\tilde{r}} \nonumber\\ -\, \frac{1}{\tilde{h}_{0}}\left(\frac{\tilde{h}_{1}}{\tilde{h}_{0}}\frac{\partial\tilde{h}_{0}}{\partial\tilde{r}} - \frac{\partial\tilde{h}_{1}}{\partial\tilde{r}}\right)\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial\tilde{r}} - \frac{\partial\tilde{\mathcal{M}}_{0}}{\partial\tilde{r}}=: \mathcal{V}(\tilde{r},t) \end{gather}$$

for $\tilde {r}>0$, $0< t<1$, subject to $\tilde {\mathcal {M}}_{3}(0,t) = 0$ for $0< t<1$ and the far-field condition

(4.41)$$\begin{align} &\tilde{\mathcal{M}}_{3}\rightarrow - \frac{(1-t)}{\rm \pi}\tilde{r} +\left(\frac{\alpha\sqrt{1-t}(1-\sqrt{1-t})}{\rm \pi}\right)\left(\log{\tilde{r}} + \log\left(\frac{2(1-t)}{(1-\sqrt{1-t})^{2}}\right)\right) \nonumber\\ &\quad\qquad+\, \frac{\alpha}{\rm \pi}(1-(1-\sqrt{1-t})^{2}) \quad \text{as}\ \tilde{r}\rightarrow\infty. \end{align}$$

The solution is given by

(4.42)\begin{equation} \tilde{\mathcal{M}}_{3}(\tilde{r},t) = \int_{0}^{\tilde{r}}\frac{1}{I(s,t)} \left(\int_{0}^{s} \mathcal{V}(\sigma,t)I(\sigma,t)\,\text{d}\sigma\right)\,\text{d}s + B_{3}(t)\int_{0}^{\tilde{r}}\frac{1}{I(s,t)}\,\text{d}s, \end{equation}

where

(4.43)\begin{align} B_{3}(t) &= \left[- \int_{1}^{\infty}\left\{\frac{1}{I(s,t)} \left(\int_{0}^{s}\mathcal{V}(\sigma,t)I(\sigma,t)\,\text{d}\sigma\right) + \frac{(1-t)}{\rm \pi} - \frac{\alpha\sqrt{1-t}(1-\sqrt{1-t})}{{\rm \pi} s}\right\}\,\text{d}s \right.\nonumber\\ &\quad \left. -\int_{0}^{1}\frac{1}{I(s,t)} \left(\int_{0}^{s}\mathcal{V}(\sigma,t)I(\sigma,t)\,\text{d}\sigma\right)\,\text{d}s - \frac{(1-t)}{\rm \pi} + \frac{\alpha}{\rm \pi}(1-(1-\sqrt{1-t})^{2}) \right.\nonumber\\ & \quad \left. +\, \frac{\alpha\sqrt{1-t}(1-\sqrt{1-t})}{\rm \pi}\log\left(\frac{2(1-t)}{(1-\sqrt{1-t})^{2}}\right)\right]\left(\int_{0}^{\infty}\frac{1}{I(s,t)}\,\text{d}s\right)^{-1} \end{align}

has been chosen to give the correct far-field behaviour.

We are now in a position to find the inner solution for the solute mass. By substituting the scalings (4.22) into (2.40), we see that

(4.44)\begin{equation} \tilde{m} =-\frac{1}{1 - {Pe}^{-2/3}\tilde{r}}\frac{\partial\tilde{\mathcal{M}}}{\partial\tilde{r}}, \end{equation}

so that expanding $\tilde {m} = \tilde {m}_{0} + {Pe}^{-1/3}\tilde {m}_{1} + {Pe}^{-2/3}\log {Pe}^{-2/3}\tilde {m}_{2} +{Pe}^{-2/3}\tilde {m}_{3} + O({Pe}^{-1}\log {Pe}^{-2/3})$ as ${Pe}\rightarrow \infty$, we have

(4.45a–d)\begin{align} \tilde{m}_{0} =-\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial\tilde{r}}, \quad \tilde{m}_{1} =-\frac{\partial\tilde{\mathcal{M}}_{1}}{\partial\tilde{r}}, \quad \tilde{m}_{2} =-\frac{\partial\tilde{\mathcal{M}}_{2}}{\partial\tilde{r}}, \quad \tilde{m}_{3} =-\frac{\partial\tilde{\mathcal{M}}_{3}}{\partial\tilde{r}} - \tilde{r}\frac{\partial\tilde{\mathcal{M}}_{0}}{\partial\tilde{r}}. \end{align}

4.2.3. Composite solutions

We now have all of the necessary components needed to construct (additive) composite solutions for comparison to the numerical results. To construct a composite solution for the integrated mass variable, we combine the first two outer solutions (4.16) and (4.18), the first four inner solutions (4.29), (4.34), (4.38) and (4.42), and the overlap terms given by (4.19)–(4.20) using van Dyke's matching rule (Van Dyke Reference Van Dyke1964), which yields

(4.46) \begin{align} \mathcal{M}_{comp}(r,t)& = \mathcal{M}_{0}(r,t) + {Pe}^{-2/3}\mathcal{M}_{1}(r,t) \nonumber\\ &\quad + \tilde{\mathcal{M}}_{0}({Pe}^{2/3}(1-r),t) + {Pe}^{-1/3}\tilde{\mathcal{M}}_{1}({Pe}^{2/3}(1-r),t) \nonumber\\ & \quad+ {Pe}^{-2/3}\log{Pe}^{-2/3}\tilde{\mathcal{M}}_{2}({Pe}^{2/3}(1-r),t) + {Pe}^{-2/3}\tilde{\mathcal{M}}_{3}({Pe}^{2/3}(1-r),t) \nonumber\\ & \quad- \left[\frac{\sqrt{1{\,-\,}t}}{\rm \pi} {\,-\,} \frac{(1{\,-\,}t)}{2{\rm \pi}} {\,-\,} \frac{\sqrt{2(1{\,-\,}t)}}{\rm \pi}(1-\sqrt{1{\,-\,}t})\sqrt{1{\,-\,}r} {\,-\,} \frac{(1{\,-\,}t)}{\rm \pi}(1{\,-\,}r) \right.\nonumber\\ & \quad+ {Pe}^{-2/3}\left(\frac{\alpha}{\rm \pi}(1-(1-\sqrt{1-t})^{2}) + \frac{\alpha\sqrt{1-t}(1-\sqrt{1-t})}{\rm \pi}\right.\nonumber\\ &\left.\left.\quad \times \left(\log(1-r) + \log\left(\frac{2(1-t)}{(1-\sqrt{1-t})^{2}}\right)\right)\right) \right]. \end{align}

This composite solution is valid up to and including $O({Pe}^{-2/3})$ throughout the whole of the droplet.

Similarly, for the solute mass, the equivalent composite profile is compiled by taking the first outer solution (4.21) as well as the first four inner solutions given by (4.45), so that, accounting for the overlap contributions,

(4.47) \begin{align} m_{comp}(r,t) &= m_{0}(r,t) + {Pe}^{2/3}\tilde{m}_{0}({Pe}^{2/3}(1-r),t) + {Pe}^{1/3}\tilde{m}_{1}({Pe}^{2/3}(1-r),t) \nonumber\\ &\quad + \log{Pe}^{-2/3}\tilde{m}_{2}({Pe}^{2/3}(1-r),t) + \tilde{m}_{3}({Pe}^{2/3}(1-r),t) \notag\\ &\quad - \frac{\sqrt{(1-t)}(1-\sqrt{1-t})}{\sqrt{2}{\rm \pi}\sqrt{1-r}} -\frac{(1-t)}{\rm \pi}. \end{align}

We note that this composite solution is valid up to and including $O(1)$ throughout the droplet.

6. Summary and discussion

In this paper, we have performed a detailed asymptotic and numerical analysis into the effect of gravity on the famous coffee-ring phenomenon observed in solute-laden droplets. For droplets evaporating in a diffusion-limited evaporation regime, in the physically relevant limit of small droplet capillary number, ${Ca}\ll 1$, and large solutal Péclet number, ${Pe}\gg 1$, we identified two asymptotic regimes based on the size of the Bond number, ${Bo}$:

1. (i) a moderate-Bond-number regime, where ${Bo} = O(1)$ as ${Pe}\rightarrow \infty$;

2. (ii) a large-Bond-number regime, ${Bo} = O({Pe}^{4/3})$ as ${Pe}\rightarrow \infty$.

In the first of these regimes, gravity acts to flatten the droplet profile from the spherical cap of the zero-gravity problem, while reducing the liquid velocity. Moreover, the asymptotic structure of the solute transport follows exactly that presented by Moore et al. (Reference Moore, Vella and Oliver2021) for surface-tension-dominated droplets, with advection dominating in the droplet bulk, while the competition between advection and diffusion in a boundary layer of width of $O({Pe}^{-2})$ near the pinned contact line drives the nascent coffee ring. Gravity acts to modify the surface-tension-dominated solution both through the accumulated mass flux of solute into the contact line and a parameter dependent on the local contact angle. In particular, as the Bond number increases, the height of the nascent coffee ring is reduced – which is consistent with the reduced flow velocity as ${Bo}$ is increased. Moreover, the peak is situated further from the contact line.

To categorize the role of gravity more explicitly, we derived an approximate similarity profile, $\hat {m}_{0}$, for the nascent coffee-ring profile, given by

(6.1)\begin{equation} \frac{\hat{m}_{0}(R,t)}{{Pe}_{t}^{2}\mathcal{N}(t;{Bo})} = \frac{2\chi}{3\psi({Bo})}f\left(\sqrt{R},3,\frac{4\chi}{\psi({Bo})}\right), \quad R = {Pe}_{t}^{2}(1-r), \end{equation}

where ${Pe}_{t} = {Pe}/(1-t)$ is the time-dependent Péclet number, $\mathcal {N}(t;{Bo})$ is the accumulated mass flux of solute at the contact line from the droplet bulk, $\chi$ is the coefficient of the inverse square root singularity in the evaporative flux at the contact line; $\psi ({Bo})$ is the leading-order initial local contact angle and $f(x,k,l) = l^kx^{k-1}\textrm {e}^{-lx}/\varGamma (k)$ is the probability density function of a gamma distribution. Clearly, the Bond number acts to scale the coffee-ring profile through the accumulated mass flux, while it acts to change the shape of the profile through the initial contact angle $\psi ({Bo})$. The nascent coffee ring has a peak height that scales with ${Pe}_t^2$ and peak location that is situated a distance from the contact line that scales with ${Pe}_t^{-2}$.

In the second regime, the Bond number is large, so that the droplet is approximately flat, with surface tension confined to a narrow region near the pinned contact line – a ‘pancake’ droplet. Thus, the asymptotic analysis discussed above is no longer valid, since there are two competing boundary layers near the edge of the droplet – the diffusion boundary layer in the solute transport and the surface tension boundary layer in the droplet free surface profile (and, hence, the liquid velocity). We derived the resulting solute distribution in the most general regime in which the two boundary layers are comparable, which reduces to the assumption that $\alpha = {Bo}^{-1/2}{Pe}^{2/3} = O(1)$ as ${Pe}\rightarrow \infty$. Under this assumption, diffusion and advection balance in a region of size ${Pe}^{-2/3}$ near the contact line, noticeably larger than in the moderate gravity regime. This is a further indication of gravity acting to shift the coffee ring further from the contact line and, moreover, tends to cause shallower solute profiles in the boundary layer region.

The asymptotic analysis in the large-Bond-number regime is more challenging than that in the moderate-Bond-number regime and, in particular, the nascent coffee ring no longer collapses onto an approximate similarity profile. However, we were able to derive expressions for the location (5.12) and height (5.13) of the peak, demonstrating that it still strongly depends on the accumulated mass flux of solute into the contact line alongside the parameter $\alpha$. In particular, the peak height scales with ${Pe}^{2/3}$ and peak location scales with ${Pe}^{-2/3}$ when $\alpha = O(1)$. Moreover, increasing $\alpha$ leads to higher coffee rings that are located closer to the contact line. This latter behaviour is indicative of the overlap between the two regimes when $1\ll {Bo}\ll {Pe}^{4/3}$, which we demonstrated explicitly in § 5.2. In particular, for a fixed Péclet number, the transition between the regimes occurs linearly in ${Bo}$.

In each regime, we demonstrated that our asymptotic predictions were in excellent agreement with numerical simulations of the advection–diffusion problem for the solute mass distribution.

Alongside the anticipated nascent coffee ring driven by the competition between advection and diffusion of the solute, our asymptotic and numerical analysis also revealed a novel phenomenon: that the solute profile may have a secondary peak. The secondary peak was characterized by being situated upstream of and significantly smaller than the primary coffee ring. Moreover, the presence of this peak strongly depends on the Bond number, Péclet number and evaporation time.

Further investigation revealed that, for a fixed Péclet number, there exists a band in $({Bo},t)$ space at which two peaks are present in the profile. We demonstrated that the onset of this band is independent of the Péclet number and is caused by the critical point at the centre of the droplet changing in type from a maximum (as in the spherical cap droplet in the ${Bo} = 0$ regime) to a minimum. When the critical point at the droplet centre changes type, an internal maximum forms downstream of the centre and it is this that corresponds to the secondary peak. This behaviour only occurs above a critical Bond number, ${Bo}_{c} \approx 15.21$, and then only after a given drying time, given by

(6.2)\begin{equation} t_{c}({Bo}) = 1 - \left(\frac{2\mathcal{U}_{1}({Bo})\mathcal{H}_{0}({Bo})}{2\mathcal{U}_{1}({Bo})\mathcal{H}_{0}({Bo})+\mathcal{H}_{1}({Bo})\mathcal{U}_{0}({Bo})}\right)^{2/\mathcal{U}_{0}({{Bo}})}. \end{equation}

In particular, as ${Bo}$ increases, $t_c$ decreases, so the secondary peak emerges earlier in the evaporative process. These predictions were shown to be in excellent agreement with the numerical results and, remarkably, are independent of the Péclet number.

However, the above analysis suggests that, for all ${Bo} > {Bo}_c$ and $t>t_c$, a secondary peak exists – something that we did not find in our analysis. The reason for this discrepancy was shown to be due to the presence of the primary peak. In particular, as time increases, the secondary peak is located further from the droplet centre so that it may get subsumed in the tail of the primary peak. For a fixed Bond number, this possibility was shown to depend strongly on both the Péclet number and the evaporation time; this is due to the fact that the size of the primary peak increases with both $t$ and ${Pe}$, while the size of the secondary peak only varies with $t$.

Beyond this subsuming effect, however, we were able to demonstrate that the Péclet number plays a negligible role in the size and location of the secondary peak for a range of Bond numbers, suggesting that this feature may be reliably controlled simply by altering ${Bo}$.

In previous studies of coffee-ring formation (e.g. Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000; Popov Reference Popov2005; Moore et al. Reference Moore, Vella and Oliver2021), gravity has frequently been neglected under the assumption of small Bond number, which is a reasonable assumption for sufficiently small droplets. However, given that the Bond number may be increased in an experimental or industrial setting by steadily increasing the droplet radius, the influence of gravity may be of fundamental interest in applications that utilize droplet drying to adaptively control the shape of the residual deposit, such as colloidal patterning (Harris et al. Reference Harris, Hu, Conrad and Lewis2007; Choi et al. Reference Choi, Stassi, Pisano and Zohdi2010) and fabrication techniques using inkjet printing (Layani et al. Reference Layani, Gruchko, Milo, Balberg, Azulay and Magdassi2009). Our analysis thus plays a dual role in the field. First, we have presented the first formal categorization of the role of gravity in the early stages of coffee-ring formation and given a quantitative prediction of the resulting features of the solute profile. The emergence of two distinct regimes based on the size of the Bond number leading to different scalings for the size of the nascent ring is a key finding of our analysis. Second, we have found a novel phenomenon – the secondary peak – which may also be exploited in such processes, particularly when the size of the primary peak can be carefully controlled or in analysing a post-jamming regime, when finite particle size effects become important (Popov Reference Popov2005; Kaplan & Mahadevan Reference Kaplan and Mahadevan2015). This is particularly relevant given that the secondary peak emerges at a moderate critical Bond number.

There are, naturally, limitations to our analysis. Throughout, we have assumed that the contact line remains pinned as the droplet evaporates. This has been shown to be a reasonable assumption for many configurations (see, for example, the experiments in Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000; Kajiya et al. Reference Kajiya, Kaneko and Doi2008; Howard et al. Reference Howard, Archer, Sibley, Southee and Wijayantha2023) and may further be enhanced by solute aggregation near the edge of the droplet (Orejon, Sefiane & Shanahan Reference Orejon, Sefiane and Shanahan2011; Weon & Je Reference Weon and Je2013; Larson Reference Larson2014). However, at late stages of the evaporation, the contact line may depin and become mobile, moving inwards towards the droplet centre. The contact line may then become pinned at a new location and the process may repeat. This behaviour is known as ‘stick-slip’ evaporation and represents an important class within the field that is beyond the scope of the present study, but may represent an interesting future direction in terms of the effect of gravity, particularly with the presence of the secondary peak and its associated increased solute mass, which may promote re-pinning.

As discussed previously, in this analysis, we have presented results for a diffusive evaporative flux, neglecting any effects of convection in the vapour transport. Depending on the configuration, this assumption may be invalid for larger droplets, so that the functional form of the evaporative flux changes. Similarly, there are other situations where different evaporative models may be appropriate. Examples include water evaporating on glass, which may more appropriately be modelled using a kinetic evaporative model (Murisic & Kondic Reference Murisic and Kondic2011), droplets evaporating above a bath of the same liquid, where the evaporation is effectively constant (Boulogne et al. Reference Boulogne, Ingremeau and Stone2016) and situations where a mask is used to control the evaporative flux so that it is stronger towards the centre (Vodolazskaya et al. Reference Vodolazskaya and Tarasevich2017). In each case, provided that the evaporation-driven flow is predominantly outwards for the majority of the drying time, we would still expect a coffee ring to form. However, the sizes of the asymptotic regimes will likely differ from the diffusive evaporative problem (see, for example, the differences between kinetic and diffusive evaporative fluxes for surface-tension-dominated droplets in Moore et al. Reference Moore, Vella and Oliver2021); nevertheless, the over-arching methodology will extend in a natural way, with the interaction between solutal advection and diffusion driving the formation of the ring.

Another effect that we have neglected in the present analysis is the possibility of solute becoming trapped at the free surface of the droplet. If this occurs, the solute is then transported to the contact line along the free surface, and has been suggested as an alternative mechanism for coffee-ring formation (Kang et al. Reference Kang, Vandadi, Felske and Masoud2016). This behaviour has been demonstrated to occur for a wide variety of droplets but is more pronounced for droplets with large contact angles (Kang et al. Reference Kang, Vandadi, Felske and Masoud2016) or when vertical diffusion happens over a longer time scale than evaporation (D'Ambrosio Reference D'Ambrosio2022). Since we deal with the opposite case of a thin droplet with fast vertical diffusion (i.e. so that the solute concentration may be assumed to be uniform across the droplet to leading order), we have not considered this phenomenon here. It would be interesting to see how such behaviour impacted the solute profile in the current case, although it should be noted that the aforementioned studies neglect gravity entirely.

A further aspect that would form the basis of an exciting future study surrounds the assumption that the solute is dilute in the droplet. Naturally, the build up of the solute in the coffee ring means that the concentration rapidly approaches levels where finite particle size effects can no longer be ignored. This has been analysed in detail for surface-tension-dominated droplets in Moore et al. (Reference Moore, Vella and Oliver2021, Reference Moore, Vella and Oliver2022) and a similar analysis would follow here with the inclusion of gravity. Indeed, the reduction of the peak height with increasing gravity is likely to extend the validity of the dilute assumption for larger droplets. Moreover, a further possible aspect that would differentiate droplets where gravity is included is in the vicinity of the secondary peak. It is an interesting open question as to whether the dilute assumption may also break down in the vicinity of the secondary peak as well as the primary. Once finite particle size effects become important, there are a number of different approaches that can be followed to continue the analysis, such as a sharp transition between a dilute and jammed region (Popov Reference Popov2005), using a viscosity and solute diffusivity that vary with concentration (Kaplan & Mahadevan Reference Kaplan and Mahadevan2015) or through more complicated two-phase suspension models (see, for example, Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018).

A future direction of interest would be to extend the analysis herein to non-axisymmetric droplets. Such droplets occur widely in applications, particularly in printing OLED/AMOLED screens (see, for example, Mai & Richerzhagen Reference Mai and Richerzhagen2007; Huo et al. Reference Huo, Shao, Dong, Liang, Bi, He, Li, Gao and Song2020). It is well known that droplet geometry plays a strong role in the behaviour of the evaporative flux (Sáenz et al. Reference Sáenz, Wray, Che, Matar, Valluri, Kim and Sefiane2017; Wray & Moore Reference Wray and Moore2023) and the transient and final deposit profiles (Freed-Brown Reference Freed-Brown2015; Sáenz et al. Reference Sáenz, Wray, Che, Matar, Valluri, Kim and Sefiane2017; Moore et al. Reference Moore, Vella and Oliver2022). It would be of significant theoretical and practical interest to explore the influence of gravity in such problems as well.

Finally, we note that another context in which gravity may play an important role is that of binary/multi-component droplets, particularly in situations where the different fluids have different densities. Multi-component droplets occur widely, from commercial alcohols such as whiskey and ouzo (Tan et al. Reference Tan, Wooh, Butt, Zhang and Lohse2019; Carrithers et al. Reference Carrithers, Brown, Rashed, Islam, Velev and Williams2020) to various inks (Shargaieva et al. Reference Shargaieva, Näsström, Smith, Többens, Munir and Unger2020). While it would be certainly of interest to extend the analysis presented here to such droplets, a careful treatment of the internal flow would be needed, as the multi-component nature of the droplet significantly complicates the dynamics (Li et al. Reference Li, Diddens, Lv, Wijshoff, Versluis and Lohse2019).