2.2. Evolution of the instability from the interface into the cell

In one set of our experiments, we dyed the ethanol–water solutions with a trace amount of rhodamine 6G (R6G). As evaporation takes place, non-volatile R6G accumulates close to the edge, resulting in brighter regions, as can be seen in figure 2(a) (these snapshots were taken at the region marked with a red square in figure 1b). These regions are shaped as arches of approximately the same size, suggesting the presence of an instability (Linde et al. Reference Linde, Pfaff and Zirkel1964; Köllner et al. Reference Köllner, Schwarzenberger, Eckert and Boeck2015). If this were not the case, we would expect the accumulation of R6G to be distributed uniformly all along the edge (white vertical dashed line), but not in specific regions.

Figure 2. (a) Snapshots taken at the region marked with a red square in figure 1(b) at three different times. The snapshots show an evaporating 50 wt% ethanol–water solution dyed with rhodamine 6G (see movie 1 of the supplementary material for the whole time evolution). The dye accumulates in regions shaped as arches that grow and merge with each other over time. We have added the time $t_0^{d}$ because the solution can start to evaporate even before it reaches the edge. The vertical white dashed lines indicate the position of the edge of the chip and the interface between air and liquid. (b) Height $z_{F}^{d}$ of the arches (see last subpanel of panel (a) for definition) as a function of time. Different independent realizations of the experiment are shown in grey, while a single case is highlighted in red. The black dashed line corresponds to the power law $z^{d}_{F} \propto (t^{d})^{0.50\pm 0.04}$. All the realizations were fitted omitting the data for which $t^{d}<1$ s. Inset: same plot, but with linear axes.

Over time, the arches grow and merge with each other as seen in figure 2(a), where their number has gone from three to one over a few seconds. This kind of coarsening behaviour is expected for Marangoni instabilities (Schwarzenberger et al. Reference Schwarzenberger, Köllner, Linde, Boeck, Odenbach and Eckert2014). Such a phenomenon was also observed in evaporation experiments inside a Hele-Shaw cell by Linde et al. (Reference Linde, Pfaff and Zirkel1964).

We measured the distance $z_{F}^{d}(t)$ that the accumulated dye penetrated towards the inside of the cell as a function of time (see the third subpanel of figure 2a for an example of $z_{F}^{d}$). Here the superscript ${d}$ indicates a dimensional quantity. To measure $z_{F}^{d}$, we took the vertical average of the intensity field, resulting in an intensity profile as a function of $z^{d}$. Then we took the derivative of this profile with respect to $z^{d}$. Finally, we looked for the distance from the edge where the derivative was zero within noise level and we assigned this value to $z_{F}^{d}$. This process was then repeated for every time frame.

In figure 2(b), we show some examples of $z_{F}^{d}(t)$ on a double logarithmic scale. A close inspection of the plots shows that during a merging event $z_{F}^{d}$ accelerates and then slows down. However, as an overall trend, $z_{F}^{d}$ grows effectively as a power law. Fourteen independent runs were fitted to $\propto (t^{d})^{\alpha }$. An average exponent of $\alpha = 0.50 \pm 0.04$ suggests that $z_{F}^{d}$ effectively grows in a diffusive-like manner. A similar behaviour has been observed in liquid–liquid systems with Marangoni instabilities (Köllner et al. Reference Köllner, Schwarzenberger, Eckert and Boeck2015). Figure 2(b) includes cases with two different mass fractions of R6G: $10^{-4}$ and $10^{-5}$. There are no considerable changes in the exponent (see the supplementary material for a plot with all cases).

We defined the time $t^{d}=0\ \textrm {s}$ to be the moment at which the solution reached the edge within our field of view. Although evaporation can take place while the chip gets filled, we expect its effect to be negligible as long as the liquid is far away from the edge. Nevertheless, it can trigger the instability already some time before the liquid reaches the edge. For the fits, data for which $t^{d}< 1$ s was omitted to minimize the effects of our uncertainty on the initial time. The rest of the time series was included in the fit.

As the arches get bigger, interesting phenomena take place inside them. We show in figure 3(a) how the dye flowed back into the chip in a meandering way, and in figure 3(b) short-lived secondary arches. The latter travelled towards the tips of the arches in a few seconds while changing in size, as also observed by Linde et al. (Reference Linde, Pfaff and Zirkel1964) and Schwarzenberger et al. (Reference Schwarzenberger, Eckert and Odenbach2012). In movie 2 of the supplementary material the secondary arches are clearly visible for a case where the liquid was seeded with tracer particles.

Figure 3. Long-time flow dynamics for the same experiment as in figure 2. In panel (a) we highlight how the dye flows back at the centre of an arch in a meandering way (marked with a white arrow). In panel (b) secondary arches are clearly visible. They last for just a few seconds and travel towards the tip of the principal arch. For these two images we adjusted the contrast for visibility. In panel (c) the arches have grown to a size comparable to the width of the central part of the chip marked by the red rectangle in panel (d). The semiarch visible in the first subpanel of (c) is in contact with the upper wall. The white dashed line indicates the edge of the chip. The images in panel (c) were made by stitching smaller pictures (taken over 20 s).

Eventually, the lateral size of the arches becomes comparable to the width of the central part of the chip, as shown in the first subpanel of figure 3(c). The area shown corresponds to the red square in the sketch of figure 3(d). We did not observe arches spanning over this whole area. Instead the arches eventually start to reduce in size (second subpanel of figure 3c), sometimes causing the remaining ones to start travelling along the edge. Overall, after 10–25 min the instability fades away and the dye is swept to the interface (third subpanel of figure 3c). Then a narrow bright stripe parallel to the edge forms, suggesting that the gradient of surface tension is not strong enough to drive the flow anymore.

A possible reason for the instability stopping is a reduced evaporation of ethanol. Indeed, after taking the chip out from the chamber, we observed that the instability can start again. Therefore, as the chip was moved to an environment without any ethanol vapour, evaporation of ethanol can be restarted and the instability triggered again.

2.3. Visualization of the bulk flow by $\mu$PIV

To visualize the flow caused by the instability, the solution was seeded with $1\ \mathrm {\mu }\textrm {m}$ fluorescent particles at a very small volume fraction. In figure 4(a) we show the streak lines made by the particles over a period of 1 s around 1.63 s after the solution reached the edge. Clearly shown in this figure, the arches are the result of advected dye (a typical movie can be seen in movie 2 of the supplementary material).

Figure 4. Flow characterization while the instability is active. (a) Streak lines formed by tracer particles obtained by averaging the intensity of 31 frames, allowing to see the same arches as with dye. As before the white dashed line indicates the edge. (b) The $\mu$PIV measurement averaged over 1 s around $t^{d} = 1.63\ \textrm {s}$ after the solution reached the edge. The measurement was taken inside the white square marked in panel (a). The colour code shows the horizontal component $u_z$ of the velocity to highlight the periodicity of the flow. The arrows were normalized to clearly show the direction (the vertical component and the total magnitude are shown in the supplementary material). The black dotted line indicates the position of the edge. (c) Sketch of the flow around and inside the arches. The gradient in concentration of ethanol produces a gradient in surface tension that drives a flow at the interface, resulting in convective rolls which cause the arch-like shape. (d) Profiles of the velocity magnitude at different times in a semilogarithmic scale. The profiles were obtained as a spatial average over the $x$-direction of the velocity fields. Different symbols represent different realizations to show reproducibility. One of the blue lines corresponds to the case at $t=1.63$ s shown in panel (b). All cases where obtained from velocity fields averaged over 1 s around the time indicated by the legend.

The direction of the flow can be read from figure 4(b) where we have plotted the corresponding normalized velocity vector field (averaged over 30 frames of a single experiment). The colour code in this image corresponds to the horizontal component of the velocity, while the vertical component and the magnitude can be seen in the supplementary material. We did not measure the velocity close to the edge for two reasons. One is the ‘shadow’ caused by the edge of the glass plate. The other is the limited frame rate of the confocal system used for imaging, while closer to the edge the flow becomes too fast.

We would expect the flow to be as depicted in figure 4(c), where the velocity at the interface is directed towards the centre of the arches. Then the fluid recirculates back to the tips of the arches. This kind of flows is typical for Marangoni-driven systems (Linde et al. Reference Linde, Pfaff and Zirkel1964; Shi & Eckert Reference Shi and Eckert2006; Schwarzenberger et al. Reference Schwarzenberger, Köllner, Linde, Boeck, Odenbach and Eckert2014; Köllner et al. Reference Köllner, Schwarzenberger, Eckert and Boeck2015). Therefore, we can conclude that the flow changes direction at a distance smaller than the closest distance to the edge at which we can still see particles (${\sim }50\ \mathrm {\mu }\textrm {m}$ from figure 4a).

The effect of the interfacial flow is quite dramatic as can be seen from figure 4(d), where we show the mean velocity profiles (averaged over 30 frames and along the $x$-direction) for different times. From these plots we can see the huge difference in velocities close to the edge as compared with far away, going from some hundreds to approximately $20\ \mathrm {\mu }\textrm {m}\ \textrm {s}^{-1}$. In particular, this decay of the velocity seems to occur exponentially as can be seen from the approximately linear character of the velocity profiles when plotted in a semilogarithmic scale. With ongoing time the slopes get smaller and their extension larger, as result of the continuous growth of the arches. In fact the slope did not change as much at longer times, reflecting the square root dependence on time that we measured for $z_F^{d}$. At longer times the arches are bigger than our field of view, so the profiles were obtained as an average over a fraction of an arch. Additionally, for the yellow curves we have masked the region inside the arch as the density of particles was too high for a correct measurement of the velocity (an example of a frame with the mask can be found in the supplementary material).

With the measured velocity we can now calculate a Péclet number, ${Pe} = U^{d}L^{d}/D^{d}_{R6G}$, with $U^{d}$ and $L^{d}$ the typical velocity and length of the flow, and $D^{d}_{R6G}$ the diffusion coefficient of $R6G$ in ethanol ($D^{d}_{R6G} = 2.9\times 10^{-10}\ \textrm {m}^{2}\ \textrm {s}^{-1}$ (Hansen, Zhu & Harris Reference Hansen, Zhu and Harris1998)). For the typical velocity we can consider the maximum velocity inside the rolls (that we can measure) which is of the order of $10^{-4} - 10^{-3}\ \textrm {m}\ \textrm {s}^{-1}$. The typical size of the arches is of the order of $10^{-4} - 10^{-3}\ \textrm {m}$. This results in a range of ${Pe} \sim 10 - 10^{3}$, clearly meaning that advection overcomes diffusion. Therefore, the dyed arches should result from the Marangoni rolls. As a last note, if we take the velocity far away from the arches (${\sim }20\ \mathrm {\mu }\textrm {m}\ \textrm {s}^{ -1}$) as the typical velocity caused by evaporation and the size of the arches just after the stability is over, namely approximately 2 mm, then we obtain $Pe\sim 100$. The dye should indeed accumulate very close to the edge as advection overcomes the diffusion of the dye.

3.1. Model equations

In the following, the governing equations modelling the evaporation of a confined ethanol–water solution into humid air as used for numerical simulations are described. A schematic of the simulated system confined by two walls is shown in figure 5.

Figure 5. Schematic of the numerical simulation. (a) The small distance $\delta ^{d}$ between the plates of the Hele-Shaw cell allows for the assumption of a parabolic flow profile in $y$-direction, which effectively reduces the problem to two spatial dimensions $(x,z)$. (b) Sketch of considered domains, fields and boundary conditions.

The equations are made non-dimensional by using the distance between plates $\delta ^{d}$ as a length scale, the initial mass density of the liquid $\rho _0^{d}$ as a mass per volume scale, and the initial kinematic viscosity $\nu _0^{d} = \mu _0^{d}/\rho _0^{d}$ to construct the viscous time scale $(\delta ^{d})^{2}/\nu _0^{d}$. Therefore, the non-dimensional space coordinates, time and velocities are given by $\boldsymbol {r} = \boldsymbol {r}^{d}/\delta ^{d}$, $t = t^{d}\nu _0^{d}/(\delta ^{d})^{2}$ and $\boldsymbol {v} = \boldsymbol {v}^{d}\delta ^{d}/\nu _0^{d}$, respectively. Due to the presence of a binary solution in the liquid phase, the mass density $\rho = \rho ^{d}/\rho _0^{d}$, the dynamic viscosity $\mu = \mu ^{d}/\mu _0^{d}$, the diffusion coefficient $D = D^{d}/\nu _0^{d}$ and the surface tension $\gamma =\gamma ^{d}({c}_{e})/({\textrm {d}\gamma _0^{d}/\textrm {d}{c}_{e}})$ of the liquid are dependent on the local ethanol mass fraction ${c}_{e}(\boldsymbol {r},t)$. The surface tension is made non-dimensional using its derivative with respect to the mass fraction at $t=0$. It is the derivative and not the absolute value of the surface tension what determines the strength of the tangential stress. The non-dimensional mass transfer rate is given by $j_\alpha = j_\alpha ^{d}\delta ^{d}/(\nu _0^{d}\rho _0^{d})$ for $\alpha = {e},{w}$. As mentioned before, we use the superscript $d$ in case we want to make reference to a dimensional variable. The subscript 0 indicates the value of a quantity at $t=0$.

This gives rise to the 3-D Navier–Stokes equations with varying density and viscosity

(3.1)$$\begin{gather} \partial_t \rho + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\rho \boldsymbol{v}\right)=0, \end{gather}$$

(3.2)$$\begin{gather}\rho\left(\partial_t \boldsymbol{v} + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{v}\right)={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot}\left[\mu\left(\boldsymbol{\nabla}\boldsymbol{v}+\boldsymbol{\nabla}\boldsymbol{v}^{t}\right)\right], \end{gather}$$

along with the convection–diffusion equation of the ethanol mass fraction ${c}_{e}$ with the composition-dependent diffusivity $D$,

(3.3)\begin{equation} \rho \left(\partial_t {c}_{e} + \boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla}{c}_{e}\right)=\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\rho D\boldsymbol{\nabla}{c}_{e}\right) . \end{equation}

In the above equations we have neglected the effects of gravity, given the small distance between the plates $\delta ^{d}$ ($20\ \mathrm {\mu }\textrm {m}$), which is the relevant length for gravity in a horizontally oriented cell. In particular if we calculate the Galilei number ${Ga} = g(\delta ^{d})^{3}(\rho _0^{{d}}/\mu _0^{{d}})^{2}$ and the Archimedes number ${Ar}=g(\delta ^{d})^{3}\rho _0^{d}\Delta \rho _0^{d}/(\mu _0^{{d}})^{2}$ based on the acceleration of gravity g, the plate distance, the density difference between pure ethanol and water ($208.5\ \textrm {kg}\ \textrm {m}^{-3}$ at 293.15 K and $211.6\ \textrm {kg}\ \textrm {m}^{-3}$ at 298.15 K), the dynamic viscosity (2.7 mPas at 293.15 K and 2.3 mPas at 298.15 K) and the density ($914\ \textrm {kg}\ \textrm {m}^{-3}$ at 293.15 K and $910\ \textrm {kg}\ \textrm {m}^{-3}$ at 298.15 K) at 50 wt%, we obtain ${Ga} \sim 0.01$ and ${Ar} = 0.002 - 0.003$. In both cases, their magnitude is much smaller than one, meaning that we can neglect gravity effects (Yu et al. Reference Yu, Lu, Lohse and Zhang2015; Li et al. Reference Li, Diddens, Lv, Wijshoff, Versluis and Lohse2019a).

As schematically depicted in figure 5(a), the 3-D flow inside the cell is strongly confined in the $y$-direction and enclosed by no-slip boundary conditions at the top and bottom plate. To the lowest non-trivial order compatible with these boundary conditions, a parabolic velocity distribution in $y$-direction can hence be assumed, i.e. $\boldsymbol {v}(x,y,z,t)=f(y)\boldsymbol {u}(x,z,t)$ with $f(y)=(by)(1-y)$ (see figure 5a). The 2-D velocity $\boldsymbol {u}$ is defined in the $(x,z)$-plane, i.e. $\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {e}_y=0$, and the constant $b$ is set to $b=6$ in the numerics so that the lateral velocity $\boldsymbol {u}$ coincides with the height average of the 3-D velocity $\boldsymbol {v}$. With these assumptions, the 3-D Navier–Stokes equations (3.2) are simplified to two dimensions, resulting in a factor $6/5$ for the nonlinearity and a drag term $-12\mu \boldsymbol {u}$ stemming from the shear of the parabolic flow profile (Bizon et al. Reference Bizon, Werne, Predtechensky, Julien, McCormick, Swift and Swinney1997; Bratsun & De Wit Reference Bratsun and De Wit2004)

(3.4)$$\begin{gather} \partial_t \rho + \boldsymbol{\nabla}_\parallel \boldsymbol{\cdot} \left(\rho \boldsymbol{u}\right)=0, \end{gather}$$

(3.5)$$\begin{gather}\rho\left(\partial_t \boldsymbol{u} + \tfrac{6}{5}\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}_\parallel\boldsymbol{u}\right)={-}\boldsymbol{\nabla}_\parallel p +\boldsymbol{\nabla}_\parallel\boldsymbol{\cdot} \left[\mu\left(\boldsymbol{\nabla}_\parallel\boldsymbol{u}+\boldsymbol{\nabla}_\parallel\boldsymbol{u}^{t}\right)\right]-12\mu\boldsymbol{u}. \end{gather}$$

Here, $\boldsymbol {\nabla }_\parallel =(\partial _x,\partial _z)$ is the 2-D nabla operator. The liquid properties $\rho$, $\mu$ and $D$ are still allowed to depend on the local ethanol mass fraction, which is now also considered in its 2-D height-averaged projection, i.e. ${c}_{e}(x,z,t)$, which is subject to the 2-D convection–diffusion equation

(3.6)\begin{equation} \rho \left(\partial_t {c}_{e} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}_\parallel{c}_{e}\right)=\boldsymbol{\nabla}_\parallel\boldsymbol{\cdot}\left(\rho D\boldsymbol{\nabla}_\parallel{c}_{e}\right) . \end{equation}

Considering only the 2-D projection of the concentration field ${c}_{e}$ does not take into account diffusion processes in $y$-direction, which can, in cooperation with the parabolic flow profile, contribute to an effectively enhanced projected 2-D diffusivity (cf. Taylor–Aris dispersion). However, due to the non-trivial and time-dependent flow, this effect cannot be taken into account within this model.

The gas phase is modelled as a mixture of air, ethanol and water, for which we assume quasi-stationary vapour diffusion, i.e. we only consider the Laplace equations

(3.7a,b)\begin{equation} \nabla^{2} c^{g}_{e}=0 \quad\text{and}\quad\nabla^{2} c^{g}_{w}=0 \end{equation}

for the vapour mass fractions $c^{g}_{e}$ and $c^{g}_{w}$ of ethanol and water, respectively. The superscript $g$ indicates the gas phase. Then the mass fraction of air is given by $c_{a}^{g} = 1 - c^{g}_{e} - c^{g}_{w}$. The assumption of a quasi-stationary model has been validated previously for evaporating ethanol–water sessile droplets. In particular, the consideration of convection (Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b) or the full diffusion equation (Diddens et al. Reference Diddens, Kuerten, Van der Geld and Wijshoff2017a) in the gas phase resulted in negligible differences in the evaporation rates or tangential velocities at the surface of the droplets. Given that the velocities in our experiments are smaller or of the same order as in the droplet case, we have decided to follow the same approach.

Opposed to the experimental set-up, the gas phase in the numerics is considered in a 2-D projection as well. To mimic the evaporation in a 3-D space, the length of the gas domain in $z$-direction is adjusted until the typical evaporation velocity matches the experimentally observed one.

The boundary conditions for the velocity and for ${c}_{e}$ are shown in figure 5(b), namely an open stress-free inflow at the liquid boundary far away from the liquid–gas interface and no-slip boundary conditions at the sidewalls. Alternatively, when only a section of the cell with some distance to the sidewalls is of interest, periodic boundary conditions in the $x$-direction can be used (if so, this will be mentioned in the text). At the liquid–gas interface, ethanol and water evaporate with local non-dimensional mass transfer rates $j_{e}$ and $j_{w}$, respectively. These, together with the assumption of a static interface, yield the kinematic boundary condition

(3.8)\begin{equation} \rho {u_z}={-}\left(j_{e}+j_{w}\right), \end{equation}

where we have neglected the solubility of air in the binary solution, such that the air flux at the interface is zero. The negative sign arises because we have defined the fluxes as positive when directed away from the liquid.

Due to the small viscosity ratio of the gas phase with respect to the liquid phase, the varying interfacial tension induces a Marangoni shear stress only on the liquid phase, i.e.

(3.9)\begin{equation} \mu\left(\partial_z {u_x}+\partial_x {u_z}\right)= \frac{\textit{Ma}_{s}}{\textit{Sc}} \frac{\textrm{d}\gamma}{\textrm{d}{c}_{e}}\partial_x {c}_{e}, \end{equation}

where we have introduced the solutal Marangoni number (Machrafi et al. Reference Machrafi, Rednikov, Colinet and Dauby2010)

(3.10)\begin{equation} \textit{Ma}_{s} ={-}\frac{\delta^{d}}{\mu_0^{d}D_0^{d}}\frac{\textrm{d}\gamma_0^{d}}{\textrm{d}{c}_{e}} \end{equation}

and the Schmidt number

(3.11)\begin{equation} {Sc} =\frac{\nu_0^{d}}{D_0^{d}}, \end{equation}

with $D_0^{d}$ the diffusion coefficient of the liquid at time zero.

In the following section, for linear stability analysis we will introduce linear mass fraction profiles close to the edge, with slope $E$. Therefore, a characteristic change in mass fraction $\Delta c_\delta = E^{d}\delta ^{d} = E$ over a distance equal to $\delta ^{d}$ is obtained. Then a ‘modified’ Marangoni number (Machrafi et al. Reference Machrafi, Rednikov, Colinet and Dauby2010)

(3.12)\begin{equation} \textit{Ma}_{s} ^{*}= \textit{Ma}_{s}\Delta c_\delta ={-}\frac{(\delta^{d}) ^{2} E^{d}}{\mu_0^{d}D_0^{d}}\frac{\textrm{d}\gamma_0^{d}}{\textrm{d}{c}_{e}} \end{equation}

can be introduced, which determines the stability of system. Notice that dispersion can result in an enhanced effective diffusion coefficient, and reduce the value of the Marangoni number. But as mentioned before, the determination of this effect is not trivial for our flow.

Due to the preferential evaporation of ethanol, the convection–diffusion equation (3.6) is subject to the boundary condition

(3.13)\begin{equation} \rho D\partial_z {c}_{e}=((1-{c}_{e})j_{e}-{c}_{e} j_{w}). \end{equation}

The far-field boundary condition of the gas composition is set to that of ambient air, plus an optional RH. At the liquid–gas interface, vapour–liquid equilibrium is imposed, i.e.

(3.14a,b)\begin{equation} c^{g}_{e}=c^{g,VLE}_{e}({c}_{e}) \quad\text{and}\quad c^{g}_{w}=c^{g,VLE}_{w}({c}_{e}), \end{equation}

where the equilibrium is calculated by Raoult's law generalized by activity coefficients, i.e.

(3.15)\begin{equation} c^{g,VLE}_\alpha({c}_{e})=\mathcal{A}_\alpha({c}_{e}){n}_\alpha({c}_{e}) \frac{p_{{\alpha},{sat}}^{d}(T){M}_\alpha^{d}}{R^{d}T^{d}\rho^{g,d}}. \end{equation}

Here, $\mathcal {A}_\alpha$ are the activity coefficients, ${n}_\alpha$ the liquid mole fractions, ${M}_\alpha ^{d}$ the molar masses, and $p_{{\alpha },{sat}}^{d}$ the saturation pressures of the pure components, whereas $R^{d}$ is the universal gas constant, $T^{d}$ is the temperature of the system and $\rho ^{g,d}$ is the density of the gas. Notice that the equilibrium mass fraction changes in time through ${n}_\alpha ({c}_{e}(z=0,x,t))$ and $\mathcal {A}_\alpha ({c}_{e}(z=0,x,t))$. The evaporation rates are finally given by the diffusive vapour flux at the interface, i.e.

(3.16a,b)\begin{equation} j_{e}=\rho^{g}D^{g}_{e}\partial_z c^{g}_{e}/\textit{Sc} \quad\text{and}\quad j_{w}=\rho^{g}D_{w}^{g}\partial_z c_{w}^{g}/\textit{Sc}, \end{equation}

where $\rho ^{g} = \rho ^{g,d}/\rho _0^{d}$ and $D^{g}_\alpha = D^{g,d}_\alpha /D_0^{d}$ are the gas to liquid density and diffusivity ratios.

As we mentioned in the introduction, we do not consider thermal effects, which in the experiments can take place by means of evaporative cooling. In order to justify this assumption we refer to Linde et al. (Reference Linde, Pfaff and Zirkel1964) and Machrafi et al. (Reference Machrafi, Rednikov, Colinet and Dauby2010), who found that the solutal effect is dominant. In particular, Machrafi et al. (Reference Machrafi, Rednikov, Colinet and Dauby2010) showed that the ratio between solutal and thermal Marangoni numbers is large for their ethanol–water solutions. We can obtain the same ratio in our case considering the thermal Marangoni number as

(3.17)\begin{equation} {Ma}_{th} = \frac{\delta^{d} \Delta T^{d}} {\mu_0^{d}D_{0,T}^{d}} \frac{\textrm{d}\gamma_0}{\textrm{d}T^{d}}, \end{equation}

where $T^{d}$ is the temperature, $\Delta T^{d}$ is a characteristic temperature change along the interface and $D_{0,T}^{d}$ is the thermal diffusion coefficient. The specific values used are shown in table 1.

Table 1. Dimensional values used in simulations for the two phases at $t=0$ and for the thermal Marangoni number estimation.

We do not have the values of the characteristic temperature change, but when assuming the range is between 1 K and even as unrealistically high as 10 K, the ratio $\textit {Ma}_{s}^{*}/{Ma}_{th} \sim 10^{3} - 10^{4}$. Therefore, the solutal Marangoni effect is much stronger than the thermal one. This is consistent with the case of evaporating sessile ethanol–water droplets, where thermal Marangoni effects could be neglected as long as the ethanol has not yet fully evaporated (Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b).

3.2. Numerical simulation implementation

For the numerical simulation we used a finite element method implemented on the basis of the finite element package oomph-lib (Heil & Hazel Reference Heil, Hazel, Bungartz and Schäfer2006). The liquid in the Hele-Shaw cell and the adjacent gas domain are discretized by a rectangular mesh, using triangular Taylor–Hood elements for the degrees of freedom of the velocity and pressure space and linear shape function for the composition in both domains. The composition dependence of the liquid and gas properties are obtained by models or fits of experimental data (see Diddens et al. (Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b) for details). The general numerical technique is described in more detail and compared with various experiments in, for example, Diddens et al. (Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017b), Diddens (Reference Diddens2017) and Li et al. (Reference Li, Diddens, Lv, Wijshoff, Versluis and Lohse2019a,Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohseb). In the present case, however, the method is generalized by considering the drag term in the momentum equation (3.5).

In all simulations the number of triangular elements was the highest close to the interface and then reduced continuously the farther away from the interface. More details on the meshes used for different simulations can be find in the supplementary material.

For the simulation used for comparison with experiments (§ 3.3) the domain size was selected to be the same as the observation area in experiments in the vertical direction. This was already quite demanding computationally, so we used a reduced resolution such that the simulation time would not be exceedingly long. However, this came at the price that the boundary condition (3.13) was not fully satisfied due to the steep gradient at the edge which could not be captured by our linear elements. Nonetheless, we were able to observe a good agreement with experiments. More details can be seen in the supplementary material. In the case of simulations used for comparison with linear stability analysis, we made use of three different configurations of the mesh which we describe in the supplementary material.

In all simulations the initial mass fraction was ${c}_{e} = 0.5$, the RH far away was taken equal to 50 % and the temperature was set to 295.15 K. The instability was triggered by numerical noise in all cases. To process the data, we linearly interpolated the mass fraction and velocity fields from the unstructured mesh used for simulations into a regular rectangular mesh.

3.3. Numerical results and comparison with experiments

In figure 6 we show results obtained from a simulation using periodic boundary conditions and a simulation with a domain size equal to the field of view that we had for experiments. Other conditions can be found in table 1 in Appendix B and in the supplementary material. The three snapshots in figure 6(a) show the evolution of the instability over three different times. As opposed to experiments, this time we can directly observe the ethanol concentration field, and on top the normalized velocity vector field. It is clear that Marangoni rolls advect ethanol-depleted fluid into the arch-shaped regions. As in experiments, the arches merge and grow over time.

Figure 6. (a) Ethanol mass fraction field on the liquid phase ${c}_{e}$ obtained from a simulation with periodic boundary conditions. The times are the same as in figure 2(a). The normalized velocity vector field is plotted on top of the concentration field. As in experiments, with progressing time, the arches merge with each other. Close to the edge, smaller secondary arches become visible (a movie showing the whole time evolution can be seen in movie 3 of the supplementary material). (b) Averaged velocity profiles at two different times. The curves with markers show the same profiles as in figure 4(d) where different markers indicate different runs. The solid lines are results from the simulation. For the blue curves, we took the closest time obtained from simulation. (c) Height $z_F$ of the arches as a function of time. The grey dotted lines show all the experimental realizations up to $t = 6.85\times 10^{4}$ ($t^{d} = 10\ \textrm {s}$). The dark red line with squares shows the results from the simulation. (c) Surface tension $\gamma$ along the edge for the three times shown in panel (a). The small scale jumps in the curves originate from the smaller secondary rolls.

After the arches are large enough we can see secondary arches too. From movie 3 in the supplementary material it can be seen that initially only the principal arches are present. Inside the arches, ‘pockets’ of rich ethanol solution are present. Secondary rolls can appear from these pockets. We can think of this as if locally the concentration field is back to the initial state when no instability was present and the gradient of concentration just builds up. Therefore, besides the fact that now there is a flow along the interface, locally a second instability can be triggered as described by Köllner et al. (Reference Köllner, Schwarzenberger, Eckert and Boeck2013) for 3-D simulations of a liquid–liquid interface.

The length of the gas domain in the $z$-direction was selected such that the velocity far away from the interface would be close to the velocity measured in experiments. We show this matching in figure 6(b). We plot non-dimensional averaged velocity profiles obtained at two different times, both from experiments and from the simulation. The curves with markers are experimental results, corresponding to the two earliest cases that were already shown in figure 4(d). The solid lines represent the results obtained from a single simulation, a single time and averaging along the $x$-direction (for the blue curve, we took a profile at the closest time obtained from the simulation). Notice that the complete profiles agreed quite well with the experimental measurements, despite the fact that we only matched the velocity far away from the edge.

It is important to note that the merging between arches did not take place at the same times between simulations and experiments, causing a discrepancy between the experimental and numerical velocity profiles during certain periods of time. In figures 10 and 11 of the supplementary material, we make more detailed comparisons of the velocity between experiments and simulations. We have noticed that the velocity increases sharply during a merging event, but changes slowly otherwise. This is in agreement with the faster growth of $z_{F}$ during merging. Moreover, as long as single arches have similar sizes, the local velocity is approximately the same. In figure 11 of the supplementary material, we compare the experimental velocity profiles of arches of similar size at different times or from different realizations. The velocity profiles match with each other. The velocity profiles obtained from a simulated arch with a similar size compare well, too. On the other hand, the local velocities did not match between experiments and simulations if the arches had different sizes.

For $t = 6.85 \times 10^{4}$, the velocity far away from the interface is smaller than in experiments. This could be the result of the comparable size of the arch with the simulation domain. Nevertheless, the main features stay similar as in experiments with a region of exponential decay until the velocity reaches the far end value. Furthermore, we can now see that the velocity close to the edge (i.e. small values $z\lesssim 2$) can reach much larger values than those measured experimentally, approximately two orders of magnitude larger.

We tracked the height $z_{F}$ as in experiments by measuring the distance from the edge at which the vertically averaged profile of the mass fraction had reached the far-field mass fraction. The result of this process is shown in figure 6(c). Both experimental and simulation results are plotted together, again showing a good match despite our uncertainty in the initial time in the experiments. Note that there is no adjustable parameter, we only matched the initial far-field velocities in the liquid.

From the simulations we can also have access to the local surface tension $\gamma (x)$ as can be seen in figure 6(d). We show the surface tension $\gamma (x)$ corresponding to the three snapshots in figure 6(a). The curves confirm the presence of gradients going from the centre of the arches to their meeting points. Overall, the gradient in surface tension seems to stay approximately constant over the time span shown here. The small scale jumps in the curves originate from the secondary rolls.

Before a merging event, both in simulations and in experiments, the stagnation point at the summit of the arches displaces towards one side. Therefore, the flow towards the edge in-between two arches reduces. In turn, the supply of ethanol-rich solution also decreases. A similar situation was observed by Köllner et al. (Reference Köllner, Schwarzenberger, Eckert and Boeck2013), where the bulk flow caused by the Marangoni rolls hindered the flow towards the interface. Afterwards, the two adjacent rolls merged into one.

4.1. Simplified model

In this section we will perform a linear stability analysis of a simplified version of the model equations described in § 3.1. In this way we will be able to obtain analytical solutions that are very close to those obtained by Sternling & Scriven (Reference Sternling and Scriven1959), with the addition of the drag term which originates from the wall friction.

In this simplified version of our model equations we will take the density, dynamical viscosity and diffusion coefficients as constants with respect to mass fraction as has been done for different systems before, with good results (Bizon et al. Reference Bizon, Werne, Predtechensky, Julien, McCormick, Swift and Swinney1997; Bratsun & De Wit Reference Bratsun and De Wit2004; Machrafi et al. Reference Machrafi, Rednikov, Colinet and Dauby2010). Furthermore, we will assume that surface tension is a linear function of the mass fraction such that the derivative of the surface tension with respect to the mass fraction is constant. In particular we will take the values at time zero such that $\rho ^{d}=\rho _0^{d}$, $\mu ^{d} = \mu _0^{d}$, $D^{d} = D_0^{d}$ and $\textrm {d}\gamma ^{d}/\textrm {d}{c}_{e}={\textrm {d}\gamma _0^{d}/\textrm {d}{c}_{e}}$. In other words $\rho =\mu ={\textrm {d}\gamma /\textrm {d}{c}_{e}}=1$ and $D = 1/\textit {Sc}$, which results in the simplified model equations for the liquid phase

(4.1)$$\begin{gather} \boldsymbol{\nabla}_\parallel \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$

(4.2)$$\begin{gather}\partial_t \boldsymbol{u} + \tfrac{6}{5}\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}_\parallel\boldsymbol{u}={-}\boldsymbol{\nabla}_\parallel p +\boldsymbol{\nabla}_\parallel^{2}\boldsymbol{u}-12\boldsymbol{u}, \end{gather}$$

(4.3)$$\begin{gather}\partial_t {c}_{e} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}_\parallel{c}_{e}= \frac{1}{\textit{Sc}}\boldsymbol{\nabla}_\parallel^{2}{c}_{e}, \end{gather}$$

and for the gas phase

(4.4a,b)\begin{equation} \nabla^{2} c^{g}_{e}=0 \quad\text{and}\quad\nabla^{2} c^{g}_{w}=0. \end{equation}

Similarly, the boundary condition for the tangential stresses at the interface (3.9) is simplified to

(4.5)\begin{equation} \partial_z {u_x}+\partial_x {u_z}= \frac{\textit{Ma}_{s}}{\textit{Sc}} \partial_x {c}_{e}. \end{equation}

Additionally, we neglect the velocity caused by evaporation such that $u_z(x,0,t)=0$. This assumption has been validated by comparing simulations with and without making the normal velocity equal to zero at the edge. The results for both cases were virtually the same. Therefore, for the linear stability analysis we kept $u_z(x,0,t)=0$, while for the simulations $u_z(x,0,t)\neq 0$. In the experimental section we obtained Péclet numbers larger than unity in odds with this assumption, however, for the onset of instability a more proper length scale to be consider is $\delta ^{d}$ (the gap of the Hele-Shaw cell), which gives ${Pe}\sim 1$, meaning that we are in the limit where we could actually neglect this velocity.

After neglecting the normal velocity, from (3.8) we have $j_{e}+j_{w} = 0$ at the interface, then (3.13) simplifies to $j_{e} = \partial _z {c}_{e}/\textit {Sc}$. Therefore, introducing this expression for the mass transfer rate of ethanol into (3.16a,b), we get

(4.6)\begin{equation} \partial_z{c}_{e} = \rho^{g}D^{g}_{e}\partial_z c^{g}_{e}, \end{equation}

recalling that $D^{g}_{e} = D^{g,d}_{e}/D_0^{d}$.

Raoult's law is simplified by taking a Taylor expansion of the molar fraction around ${c}_{e}(t=0)\equiv c_{e,0} \ne 0$, such that ${n}_{e} = (1+(1/{c}_{e}-1)M_{e}^{d}/M_{w}^{d})^{-1} \approx a_0(c_{e,0},M_{e}^{d}/M_{w}^{d}) + a_1(c_{e,0},M_{e}^{d}/M_{w}^{d}) {c}_{e} + O({c}_{e}^{2})$. Finally, if we take the activity coefficient $\mathcal {A}_\alpha = 1$ we are left with

(4.7)\begin{equation} c^{g,VLE}_e({c}_{e}) \approx c_{sat}(a_0 + a_1 {c}_{e}), \end{equation}

where $c_{sat} =p_{e,sat}^{d}(T)M_{e}^{d}/(R^{d}T^{d}\rho ^{g,d}_0)$. An analogous expression can be found for water.

4.2. Linear stability of 2-D equations

For the base case we will consider no flow in the liquid phase, i.e. $\pmb {u}_0 = 0$, which results in a constant pressure $p_0$. While for the mass fraction linear profiles close to the edge will be taken, in the same spirit as Sternling & Scriven (Reference Sternling and Scriven1959),

(4.8)$$\begin{gather} c_{e,0} = Ez+\mathcal{C}_0, \end{gather}$$

(4.9)$$\begin{gather}c_{e,0}^{g} = E^{g}z+\mathcal{C}^{g}_0, \end{gather}$$

with $E$ and $E^{g}$ made non-dimensional with $\delta ^{d}$. A more complete model would consider a time evolving base state with an error-function-like shape, however, then the equation has to be solved numerically. Here we consider that close to the edge the profiles can be approximated as linear. We also assume that the changes in concentrations at the edge are small around the time the instability takes place. By focusing our analysis on the region close to the interface, our analysis corresponds to the small wavelength regime.

We considered small perturbations in the standard way for the velocity, pressure and mass fraction fields,

(4.10)$$\begin{gather} \pmb{u}= \epsilon \pmb {\hat {u}}_1(z) \exp({\textrm{i}\kappa x+\sigma t}), \end{gather}$$

(4.11)$$\begin{gather}p = p_0 + \epsilon \hat{p}_1(z) \exp({\textrm{i}\kappa x+\sigma t}), \end{gather}$$

(4.12)$$\begin{gather}{c}_{e} = c_{e,0}(z,t) + \epsilon \hat{c}_{e,1}(z) \exp({\textrm{i}\kappa x+\sigma t}), \end{gather}$$

(4.13)$$\begin{gather}c^{g}_{e} = c^{g}_{e,0}(z,t) + \epsilon \hat{c}^{g}_{e,1}(z) \exp({\textrm{i}\kappa x+\sigma t}), \end{gather}$$

where $\epsilon$ is a small number, $\kappa$ is the non-dimensional wavenumber of the perturbation, $\sigma$ its non-dimensional growth rate and $\pmb {\hat {u}}_1(z)$, $\hat {p}_1(z)$, $\hat {c}_{e,1}(z)$ and $\hat {c}_{e,1}^{g}(z)$, the amplitudes of the corresponding perturbations in velocity, pressure and mass fraction in the liquid and gas, respectively. These perturbed quantities can be used to linearize equations (4.1) to (4.4a,b), resulting in the following dispersion relation (see Appendix C for further details):

(4.14)\begin{align} \textit{Ma}_{s}^{*} & = \kappa^{2} \left(c_{sat}a_1\rho^{g}D^{g}_{e}+ \sqrt{1+\frac{\sigma}{\kappa^{2}}\textit{Sc}}\right)\left(1+\sqrt{1+\frac{12}{\kappa^{2}} + \frac{\sigma}{\kappa^{2}}}\right)\nonumber\\ &\quad \left(1+\sqrt{1+\frac{\sigma}{\kappa^{2}}\textit{Sc}}\right)\left(\sqrt{1+\frac{12}{\kappa^{2}} + \frac{\sigma}{\kappa^{2}}} + \sqrt{1+\frac{\sigma}{\kappa^{2}}\textit{Sc}}\right). \end{align}

Notice that this is a simplified version of equation (30) from Sternling & Scriven (Reference Sternling and Scriven1959), but with the addition of the drag term and a Marangoni number $\textit {Ma}_{s}^{*}$ based on the plate distance $\delta ^{d}$.

From this equation we can obtain the condition of marginal stability by setting $\sigma =0$, resulting in

(4.15)\begin{equation} \kappa_{\sigma=0} = \sqrt{\frac{\textit{Ma}_{s}^{*}}{8\mathcal{M}}} -3\sqrt{\frac{8\mathcal{M}}{\textit{Ma}_{s}^{*}}}, \end{equation}

where $\mathcal {M} = 1 + c_{sat}a_1\rho ^{g}D^{g}_{e}$. This relation is independent of ${Sc}$, which was expected from the results of Sternling & Scriven (Reference Sternling and Scriven1959). Their results can be recovered in the limit of large $\textit {Ma}_{s}^{*}$, i.e. $\kappa _{\sigma =0} = (\textit {Ma}_{s}^{*}/8\mathcal {M})^{0.5}$ where the dependence on $\delta ^{d}$ is not present anymore, as $\kappa _{\sigma =0}\propto \delta ^{d}$ as well as $\sqrt {\textit {Ma}_{s}^{*}}\propto \delta ^{d}$. This is equivalent to making the gap between plates large enough so that the limit of an unbounded system is reached (Martin, Rakotomalala & Salin Reference Martin, Rakotomalala and Salin2002).

We can now find a critical Marangoni number ${Ma}^{*}_c\equiv 24\mathcal {M}$ for which $\kappa _{\sigma =0} = 0$. Below this critical value, the right-hand side of (4.15) becomes negative, while $\kappa$ has to be positive. Therefore, for $0 < \textit {Ma}_{s}^{*} \leqslant Ma^{*}_c$ the system should become stable for all wavelengths. Furthermore, the Marangoni number can become negative if either the gradient $E$ or the derivative of surface tension with concentration ${\textrm {d}\gamma _0^{d}/\textrm {d}{c}_{e}}$ change in sign. Of course this is not possible as the wavenumber has to be a real positive quantity. However, this does not exclude the possibility of an oscillatory instability, which could be expected for negative Marangoni numbers considering the results from Sternling & Scriven (Reference Sternling and Scriven1959). Given that the evaporating ethanol–water system results in positive Marangoni numbers, we will restrict ourselves to this regime.

4.3. Numerical solution of the dispersion relation

Equation (4.14) was solved numerically to obtain $\sigma$ as a function of $\kappa$ and ${Ma}^{*}_{s}$ for $\textit {Sc} = 7.51\times 10^{3}$, and the properties shown in table 1 in Appendix B corresponding to an ethanol–water solution at 50 wt% and 295 K. The results can be seen in figure 7(a), where we have clamped negative $\sigma$ to zero to facilitate visualization. Equation (4.15) is plotted as a solid white line, while the purely 2-D limit, $\kappa _{\sigma =0} = (\textit {Ma}_{s}^{*}/ 8\mathcal {M})^{0.5}$ is highlighted by a dotted line. The red solid line corresponds to the locations where $\sigma$ is maximum for a given ${Ma}_{s}^{*}$, i.e. the fastest growing modes.

Figure 7. Phase diagram of the growth rate $\sigma$ for different values of $\kappa$ and (a) $\textit {Ma}_{s}^{*}$, with constant ${Sc}= 7.51\times 10^{3}$, and different values of (b) ${Sc}$, with constant $\textit {Ma}^{*}_{s}= 1.54\times 10^{5}$. The vertical white dashed line in panel (a) marks ${Ma}^{*}_{c}=121.2$, below which the system is always stable, while the dotted line corresponds to the limit $\kappa _{\sigma =0} = (\textit {Ma}_{s}^{*}/ 8\mathcal {M})^{0.5}$. In both panels (a,b) we have clamped negative values of $\sigma$ to zero to facilitate visualization (grey region). The solid white curve indicates $\sigma =0$, while the solid red line indicates $\max (\sigma )$. The black region in panel (a) indicates where we did not find solutions for the dispersion relation. (c) The growth rate $\sigma$ as a function of $\kappa$ corresponding to the red vertical dashed line shown in panels (a) and (b).

Indeed, for $\textit {Ma}_{s}^{*}<{Ma}^{*}_{c}=121.2$ the system is expected to be stable, as $\sigma <0$ for every value of $\kappa$. However, we must notice that for regions where $\sigma <-\kappa ^{2}/\textit {Sc}$, we did not find solutions (marked as a black region). This region corresponds to where three of the radicals become purely imaginary. We looked for imaginary solutions of $\sigma$, without success. While our search was not exhaustive, we are in the region where Sternling & Scriven (Reference Sternling and Scriven1959) only found real solutions for the purely 2-D case so we think it is reasonable not to find complex solutions.

Once a particular liquid–gas system has been fixed (with a given mass fraction gradient $E$), the Marangoni number $\textit {Ma}_{s}^{*}$ indicates the effect of the gap between the walls of the Hele-Shaw cell. By comparing the solid and dotted lines in figure 7 we can observe that the effect of the confinement as compared with the purely 2-D case is only visible for $\textit {Ma}_{s}^{*}\lesssim 1000$. Meaning that for larger values the purely 2-D model of Sternling & Scriven (Reference Sternling and Scriven1959) should be enough to calculate the critical wavelength. In other words, for a given set of properties, as long as $\delta ^{d}$ is large enough, but not so large that gravity starts to play a role, the use of the purely 2-D model is sufficient.

According to Sternling & Scriven (Reference Sternling and Scriven1959), if the conditions $D^{d}/D^{g,d}<1$, $\nu ^{d}/\nu ^{g,d}<1$ and $\textrm {d}\gamma ^{d}/\textrm {d}{c}_{e}<0$ are satisfied, while solute transfer occurs from the liquid to the gas, the system should be unstable. However, consideration of the drag of the walls for an evaporating binary liquid introduces the possibility of stability as long as $\textit {Ma}_{s}^{*} < \textit {Ma}^{*}_{c}$.

We also looked at the effect of the Schmidt number, ${Sc}$, as shown in figure 7(b) for a fixed $\textit {Ma}_{s}^{*}=1.54\times 10^{5}$. As mentioned by Sternling & Scriven (Reference Sternling and Scriven1959), varying ${Sc}$ shows only a minor effect which was already suggested by the fact that $\kappa _{\sigma =0}$ is not a function of this parameter. Only for ${Sc}\lesssim 50$ there is a slight increase on the value of the fastest growing wavenumber.

A plot of $\sigma$ as a function of $\kappa$ is displayed in figure 7(c). This case corresponds to the dashed red line in both panels (a) and (b) of figure 7 ($\textit {Ma}_{s}^{*} =1.54\times 10^{5}$ and $\textit {Sc} =7.51\times 10^{3}$). We have selected this particular case, because it corresponds to the simulated case in § 3.3. In order to know the correct value of $E$, we conducted simulations in a smaller domain such that the mismatch of the boundary condition was not present. If we then assume that this corresponds to the same conditions as in experiments (for which we do not know the actual gradient), it would be in accordance with the presence of an instability, as the Marangoni number is above the critical value.

With respect to the term $c_{sat}a_1\rho ^{g}D^{g}_{e}$, it affects the range over which the system becomes stable, as clearly seen from (4.15). In particular, it increases the value of ${Ma}^{*}_{c}$, while in the extreme case, in which it is equal to zero, one gets ${Ma}^{*}_{c}=24$, which corresponds to completely ignoring the gas phase and just setting a constant solute flux at the interface. Additionally, we noticed that the values of $\sigma$ decrease with increasing values of the factor $c_{sat}a_1\rho ^{g}D^{g}_{e}$.

One of the simplifications that we took for the linear stability analysis was to consider the solution as ideal by setting the activity coefficient $\mathcal {A}_{e}=1$. If this is not the case, an increase or decrease of the activity coefficient would affect the value of $c_{sat}$ correspondingly. Therefore, changes in the activity coefficient would come into effect through $c_{sat}$. In particular for ethanol–water mixtures with activity coefficients above one, consideration of the activity coefficient would cause the growth rate $\sigma$ to decrease.

4.4. Comparison of linear stability analysis with numerical simulations

Experimentally, we were not able to observe the onset of the instability since at that moment the accumulation of dye or particles was too low to create a large enough contrast. Furthermore, given that the arches merged with each other, it is possible that by the time we are able to observe them, they are already the result of different merging events, which is what we observe from simulations. Therefore, we relied on the latter for comparison with theory. To do so, we performed simulations over short times, such that merging had not taken place, and in domains at least five times larger than the fastest growing wavelength. To test different Marangoni numbers we varied the evaporation rate by changing the size of the gas domain.

In figure 8 we show the mass fraction of ethanol $c_{e}$ as a function of space and time. In particular, these are results for a case of $\textit {Ma}_{s}^{*} = 1.55\times 10^{5}$. In figures 8(a) and 8(b) we show the averaged mass fraction profiles on the liquid and gas side, respectively. In both cases the average was taken along the $x$-direction and we display the profiles over five different times. In the case of the gas profiles, the changes are so small that they lie on top of each other almost indistinguishably. Notice that we have plotted the liquid profiles on a double logarithmic scale and we have shifted them by subtracting the averaged mass fraction at the interface, i.e. $\langle {c}_{e}(z=0)\rangle _x$ (shown in figure 8c). In this way we can see that close to the interface the profiles are approximately linear, while of course far away they reach the far-field value. As in the work by Köllner et al. (Reference Köllner, Schwarzenberger, Eckert and Boeck2013) we observed a local minimum that results from the mixing zone caused by the Marangoni rolls, indicating that at those times the instability has already developed and that ethanol depleted liquid has already been advected back into the bulk.

Figure 8. (a) Ethanol mass fraction ${c}_{e}$ averaged over the $x$ coordinate as a function of the position $z$ in a log–log axis for a simulation at very early times. The curves were shifted by subtracting the mean mass fraction at $z=0$ to show that close to the edge the profile is actually linear. Different curves represent different times (see legend). (b) Mass fraction profile as a function of position $z$ on the gas phase for the same times as shown in panel (a). The changes are so small that the curves practically overlap with each other. (c) Evolution in time of the averaged mass fraction along the edge $z=0$. (d) Mass fraction profile along the edge $z=0$. All curves were shifted by the corresponding averaged mass fraction at the edge. The colours correspond to the same times as those shown in panel (a). Movie 4 of the supplementary material shows the whole concentration and velocity vector fields as a function of time.

In figure 8(d) we show the mass fraction profiles along the edge at different times. These profiles were shifted by their corresponding mean values for visibility. Clearly, the instability has not grown very much after approximately seven viscous times, while at $t=10.3$ it has dramatically increased. The presence of the instability is also suggested by the increase of the averaged mass fraction at the edge (figure 8c) and the standard deviation of the surface tension in figure 9. Furthermore, at $t = 17.1$ merging has already started, showing the very early onset of nonlinear effects.

Figure 9. Average diffusive flux $\langle \partial _z c_e\rangle_x$ at the edge and standard deviation of the surface tension $\mathrm {std}(\gamma )$ as a function of time, corresponding to the same simulations as shown in figure 8.

From the profile at $t = 27.4$ it can be seen that more than one wavelength is still present at the same time, which can be expected since the perturbation comes from numerical noise combined with the non-uniform shape of the mesh, meaning that more than one mode can be excited at the same time. Furthermore, certain modes could be excited again, given that at these early times merging takes place by the displacement (along the edge) of the formed rolls, leaving free space for new ones to be generated, as can be seen close to the end of movie 4 of the supplementary material. In fact, we tried cases where we intentionally introduced a perturbation with one particular wavelength, but other wavelengths developed too. Since the results were quite similar, here we present only the cases with random noise.

The presence of many wavelengths at the same time makes it difficult to determine which one grows the fastest just from looking at the mass fraction profiles, especially since merging starts very soon. In order to answer this question, we calculated the Fourier transform of the edge profiles (figure 8d) at each time. Then we calculated the growth rate of the power spectra for every single wavenumber present in the mass fraction profiles. This growth rate corresponds to $\sigma = \sigma(\kappa)$ under the assumption that for some time many wavelengths will be present, and not only the fastest one. A more detailed description of the procedure is presented in Appendix D.

The corresponding plot of growth rate $\sigma$ versus wavenumber $\kappa$ is shown in figure 10 for different $\textit {Ma}_{s}^{*}$. In each case, the curves are the average over different realizations and the error bars represent one standard deviation on each side. We have included cases using the different kinds of meshes described in § 3.2. We varied the size of the domain in the $x$- and $z$-directions, but always keeping the domain at least five times the size of the fastest growing wavelength in the $x$-direction. We only show cases where we did not impose any initial perturbation and the instability started from numerical error.

Figure 10. Comparison of the normalized growth rate $\sigma$ versus normalized wavenumber $\kappa$ obtained from linear instability theory (solid line) and from simulations (*). The upper horizontal axis shows the corresponding wavelength $\lambda ^{d} = 2{\rm \pi} \delta ^{d}/\kappa$. The mean $\langle \textit {Ma}_{s}^{*}\rangle$ indicated in each plot corresponds to the average value over different realizations (see the supplementary material for a table with the values of $\textit {Ma}_{s}^{*}$ of individual simulations).

Given that in simulations we control the evaporation rate through the size of the gas domain, in order to determine the modified Marangoni number $\textit {Ma}_{s}^{*}$ of each case, we obtained $E$ from the average diffusive flux at the boundary $\langle \partial _z {c}_{e} \rangle _x$. Figure 9 shows the evolution in time of $\langle \partial _z {c}_{e} \rangle _x$. For the first time steps the mass fraction gradient increases very quickly from zero up to the value determined by the boundary condition (4.6), then it increases slowly until merging takes place and the gradient slightly decreases. In all cases, we took the average value over the same period for which we fitted the growth rate of the power (see Appendix D). This average mass fraction gradient $\bar {E} \equiv \overline {\langle \partial _z {c}_{e} \rangle }_x$ was then used to calculate $\textit {Ma}_{s}^{*}$ together with material quantities corresponding to ${c}_{e} = 0.5$. The corresponding values of $\bar {E}$ and $\textit {Ma}_{s}^{*}$ obtained for each different run are shown in the supplementary material. The values of $\langle Ma^*_s \rangle$ shown in figure 10 are the average of the modified Marangoni number over all the different runs.

With the Marangoni number known, we calculated the theoretical growth rates as a function of wavenumber for each case. The results are plotted as solid red lines in figure 10. The prediction tends to underestimate $\kappa$ by a factor of approximately 0.5, and for smaller modified Marangoni numbers, the growth rate is underestimated too. Furthermore, we would expect that the inclusion of the activity coefficient would shrink the theoretical curves. However, the agreement is still fairly good considering all the simplifications done for the linear stability analysis, and the fact that in the simulations nonlinear terms and water evaporation are taken into account.

Notice that for large values of $\kappa$, the growth rate calculated from the simulations does not keep on decreasing. This could be a result of the higher resolution needed to better resolve the corresponding wavelengths. Moreover, it can be seen from figure 11(a) in Appendix D that the power of large values of $\kappa$ ($\kappa > 100$ in that example) is quite small. Therefore, above a given value of $\kappa$ we cannot trust anymore the obtained values of the growth rate.

Figure 11. (a) Power density function obtained from the fast Fourier transform of the mass fraction profile at $t = 10.3$ from figure 8(d) as a function of the normalized wavenumber $\kappa$. (b) Same power as in panel (a) but now as a function of time in a semilogarithmic scale. The four curves correspond to the vertical dashed lines in panel (a). The shadowed area represents the section over which a linear fit (in semilogarithmic scale) was used to calculate the growth rate $\sigma$.

On top of each of the plots from figure 10 we have indicated the corresponding values of the dimensional wavelength $\lambda ^{d} = 2{\rm \pi} \delta ^{d}/\kappa$. Note that for the two largest $\textit {Ma}_{s}^{*}$ the fastest growing wavelength is smaller than the gap between the plates of the cell ($\delta ^{d} = 20\ \mathrm {\mu }\textrm {m}$). Bratsun & De Wit (Reference Bratsun and De Wit2004) also found that for large Marangoni numbers their predicted fastest wavelength was smaller than the gap of their cell, which led them to disregard those solutions. A reason for this is that the Hele-Shaw model relies on the assumption that variations of the velocity in the direction parallel to the plates of the cell are much smaller than variations in the directions normal to the plates (Ruyer-Quil Reference Ruyer-Quil2001). Therefore, a wavelength smaller than the gap of the cell would result in changes of velocity that would be stronger in the plane of the cell.

In fact, from (4.15) we can obtain a value for $\textit {Ma}_{s}^{*}$ for which $\kappa = 2{\rm \pi}$, meaning that $\lambda ^{d} = \delta ^{d}$. This implies that $\textit {Ma}_{s}^{*} = 8\mathcal {M}({\rm \pi} + \sqrt {({\rm \pi} ^{2}+3)})^{2}$, meaning that for Marangoni numbers larger than this quantity, there will be at least one wavelength smaller than the gap. This is already the same regime over which the model becomes close to the purely 2-D case (from the point of view of linear stability analysis) suggesting that for large Marangoni numbers a 3-D model would be more suitable in order to look at the onset of the instability, even if the gap is as small as $20\ \mathrm {\mu }\textrm {m}$. Nevertheless, as time goes on and coarsening takes place by means of merging, this restriction is not in play anymore, explaining the good agreement between simulations and experiments at much longer times. Furthermore, if the properties of the fluids or the evaporation rate result in a low enough Marangoni number, then (4.14) could be used as a first estimation of the growth rate of long wavelengths at onset.