Marangoni instability triggered by selective evaporation of a binary liquid inside a Hele-Shaw cell

Abstract

Interfacial stability is important for many processes involving heat and mass transfer across two immiscible phases. For the evaporation of a binary solution with one component more volatile than the other, gradients in surface tension can arise. These gradients can ultimately destabilize the liquid–gas interface. We study the evaporation of an ethanol–water solution, for which ethanol is more volatile. The solution is contained in a horizontal Hele-Shaw cell which is open from one end to allow for evaporation into air. A Marangoni instability is then triggered at the liquid–air interface. We study the temporal evolution of the instability through its effects on the bulk of the liquid. More specifically, the growth of convective cells is measured with confocal microscopy and the velocity field with microparticle image velocimetry. The results of numerical simulations based on quasi-two-dimensional equations satisfactorily compare with the experimental observations, even without consideration of evaporative cooling, although this cooling can play an extra role in experiments. Furthermore, a linear stability analysis applied to a simplified version of the quasi-two-dimensional equations showed reasonably good agreement with the results from simulations at early times, when the instability has just been triggered and before coarsening. In particular, we find a critical Marangoni number below which a regime of stability is predicted.

1. Introduction

The stability of an interface is important for processes where mass or heat transfer between two phases take place, e.g. drying of paint, coating, distillation, absorption and extraction. In the particular case of a Marangoni instability, variations of temperature or concentration along the interface give rise to surface tension gradients that trigger convection. The study of interfacial instability goes back to the pioneering work of Bénard (1901). Although in that work natural convection was considered to be the cause of the instability, Pearson (1958) showed by means of linear stability analysis that some of the observations could be explained by gradients in surface tension caused by variations in temperature (thermal Marangoni instability). Soon after, Sternling & Scriven (1959) showed, by a linear stability analysis too, that the same effect could be driven by gradients of concentration (solutal Marangoni instability).

Since then, a large body of studies has extended and generalized the works by Sternling & Scriven (1959), for example, to three dimensions with deformable interfaces (Scriven & Sternling 1964; Hennenberg, Sørensen & Sanfeld 1977) and gravity waves (Smith 1966). Also the effects of chemical reactions (Ruckenstein & Berbente 1964), finite domains (Reichenbach & Linde 1981), finite domains with binary mixtures and the Soret effect (Bergeon et al. 1998) and modified geometries (such as spherically symmetric (Sørensen 1980)) have been considered. Exponential concentration profiles and surfactant accumulation at the interface were studied by Sørensen et al. (1977) and Sørensen, Hennenberg & Hansen (1978). More recently, Bratsun & De Wit (2004) studied the instability between two liquids inside a Hele-Shaw cell where a chemical reaction takes place at the interface. Picardo, Radhakrishna & Pushpavanam (2016) considered the case of evaporation coupled with a Poiseuille flow, and Machrafi et al. (2010) looked at the stability of an evaporating ethanol–water mixture in two dimensions taking evaporative cooling and the Soret effect into account. This list is by no means extensive, the reviews by Levich & Krylov (1969), Sanfeld & Steinchen (1984) and Kovalchuk & Vollhardt (2006) present a more comprehensive overview, including cases where heat transfer is also considered. We refer to the reviews of Linde, Schwarzenberger & Eckert (2013) and Schwarzenberger et al. (2014) for further experimental and numerical work.

While a considerable part of the experimental work has focused on extended interfaces (Zhang, Oron & Behringer 2011; Linde et al. 2013), there has also been quite some work on binary systems confined within a Hele-Shaw cell. Employing this kind of cell has the advantage that it allows us to have visual access to the effects that the instability has in the bulk of the two phases. Furthermore, the system can be considered as quasi-two-dimensional (quasi-2-D). Many of the papers that make use of a Hele-Shaw cell have focused on liquid–liquid interfaces where chemical reactions take place, including experimental (Eckert, Acker & Shi 2004; Shi & Eckert 2006, 2007, 2008; Eckert et al. 2012; Schwarzenberger, Eckert & Odenbach 2012), numerical (Grahn 2006; Mokbel et al. 2017) and theoretical (Bratsun & De Wit 2004) approaches. In particular, Köllner et al. (2015) performed both experiments and numerical simulations, showing good qualitative agreement between them, though the Marangoni rolls in experiments grew faster than in simulations.

Following the criteria established by Sternling & Scriven (1959), Linde, Pfaff & Zirkel (1964) studied the stability of liquid–gas interfaces of binary systems, where either selective evaporation or adsorption triggered a Marangoni instability. Linde et al. (1964) pointed out that both thermal and solutal instabilities were at play simultaneously. By looking at the evaporation of pure liquids to isolate the thermal instability, they observed that the flow was weaker than when the solutal Marangoni instability was also present. In this way they were able to suggest that solutal Marangoni was the strongest mechanism. However, in that paper, the quantitative description of the phenomena was limited.

More generally, great efforts have been devoted to the understanding of the physicochemical hydrodynamics of multicomponent systems (Levich 1962), especially of multicomponent droplets, due to their importance in many fields such as chemical diagnosis, inkjet printing or nanotechnology, to mention a few. For a recent review on this subject we refer to Lohse & Zhang (2020). In many of these systems, mass transfer through interfaces causes gradients in chemical concentration that in turn result in a Marangoni flow. Such a flow can then lead to quite rich phenomena; for example, it can participate in the segregation of evaporating binary droplets (Li et al. 2018), make a droplet inside a bath jump against the pull of gravity (Li et al. 2019b, 2021), and more generally cause the self-propulsion of droplets (Izri et al. 2014; Maass et al. 2016; Lohse & Zhang 2020). These Marangoni flows can also assist the emulsification of evaporating ternary droplets (Tan et al. 2016, 2017), or drive the so-called Marangoni bursting (Keiser et al. 2017). In some cases Marangoni forces compete with gravity, leading to unexpected flow directions inside evaporating sessile and pendant droplets (Li et al. 2019a; Diddens, Li & Lohse 2021). Clearly, these phenomena are closely connected with the presence of instabilities at interfaces where mass transfer and evaporative cooling take place. However, the three-dimensional (3-D) shape of the droplet which changes in time makes the geometry complicated. Therefore, the Hele-Shaw geometry is a good candidate to further understand the Marangoni convection present when selective evaporation takes place. Moreover, the very common use of microfluidic devices in research nowadays makes it important to understand the influence that the confinement has on the Marangoni convection. Indeed, Marangoni flow in microfluidic devices can be used to enhance mixing (Michelin et al. 2020).

The aim of this work is to revisit one of the systems studied by Linde et al. (1964), namely a binary solution of ethanol–water confined by a Hele-Shaw cell evaporating into air. In particular, we study the evolution of the instability in time by means of confocal microscopy and microparticle image velocimetry ($\mu$PIV) (§ 2). In § 3, this system is modelled and numerical simulations of the model equations are performed. This model is quasi-2-D and takes into account the drag caused by the walls of the cell (Bizon et al. 1997; Bratsun & De Wit 2004; Köllner et al. 2015). The simulation results are able to reproduce our experimental observations quite well, even though we do not consider the effects of evaporative cooling. In § 4, inspired by a previous work (Chen et al. 2021), we present a linear stability analysis of a simplified version of the quasi-2-D model equations. We compare the results of the linear stability analysis with those from numerical simulations at early times, showing not perfect, but close agreement. Finally, in § 5 we conclude and give an outlook on future work with ternary systems.

2. Experimental methods and results

2.1. Materials and methods

2.1.1. Design of Hele-Shaw cell

A microfluidic chip (Beijing First Mems Co., Ltd) was used as a Hele-Shaw cell. In figure 1(a) a cross-section of the chip is shown. The black region was made of an edged silicon wafer, while the light blue plate was made of glass to allow for optical visualization. The gap between the two plates was of $20 \ \mathrm {\mu }\textrm {m}$ and kept constant by means of a series of pillars depicted as thin vertical lines in figure 1(a). Further details can be seen in Appendix A.1. At the back of the chip (right-hand edge in figure 1b,c), a channel $50 \ \mathrm {\mu }\textrm {m}$ in depth was used to fill the chip with the binary liquid. Three sides of the chip were closed and one was open to allow for evaporation (left-hand side in figure 1b,c). The surface of the chip was made hydrophobic as described in Appendix A.2.

Figure 1. Sketch of the experimental set-up from (a) frontal, (b) bottom and (c) side views. The main component was a microfluidic chip made of an edged silicon wafer and a glass plate for observation. The gap between both plates was $20\ \mathrm {\mu }\textrm {m}$. To maintain a homogeneous gap, a matrix of pillars was placed between the wafer and the glass, represented as thin vertical lines in panel (a) and as circles in panel (b) (not to scale). Additionally, two walls divided the cell in three different parts. At the back of the chip a channel of $50\ \mathrm {\mu }\textrm {m}$ in depth was used to continuously refill the cell and compensate for mass loses due to evaporation. As depicted in panel (a), the whole chip was surrounded by a chamber where the humidity was monitored and controlled with a flow of dry and humid nitrogen. In panel (b), the gap region is represented in grey and the red rectangle indicates the region of observation. In panel (c), the meniscus is drawn (solid green line) to indicate a slightly wetting liquid. The dimensions are not to scale for better visualization. The silicon wafer is $400\ \mathrm {\mu }\textrm {m}$ thick without edging and the glass plate is $500\ {\mathrm {\mu }\textrm {m}}$ thick. The holder cap was approximately $170 \ {\mathrm {\mu }}\textrm {m}$ above the chip.

2.1.2. Experimental set-up

We visualized the evaporation process at a region around the open edge (red square in figure 1b) with a confocal microscope (Nikon Confocal Microscopes A1 system). As opposed to the liquid, the gas phase was not contained by the Hele-Shaw cell; instead it was a 3-D volume of approximately $3000\ \textrm {cm}^{3}$. This volume was created by covering part of the microscope with plastic wrap (see supplementary material available at https://doi.org/10.1017/jfm.2021.555). The relative humidity (RH) in this chamber was maintained within a desired range using a negative feedback. Every time the RH was outside the desired range, a flow of either humid or dry nitrogen was gently blown in. As long as the RH was within the range, the nitrogen flow was turned off. Humidity and temperature sensors (HIH6130, Honeywell) were used to monitor these two quantities. The readings were taken with an Arduino board (Arduino Nano) which also controlled the valves of the humid and dry nitrogen, accordingly. The range over which there was no flow was either $50 \pm 1\,\%$ or $51 \pm 1\,\%$, depending on the RH in the room. This resulted in an overall range between 48 %–52 % for all our experiments. Examples of the time series of the RH in the chamber and extra details on the chamber can be found in the supplementary material.

In the case of the temperature of the air in the chamber, we only measured it during each experiment, but did not control it. Due to the heating caused by the electronics, the temperature was constantly rising at approximately $2\ \textrm {K}\ \textrm {h}^{-1}$. In the majority of our experiments the focus was on the first 100 s, during which the variations in temperature stayed within 0.3 K. The initial temperature of different experiments ranged from 294 to 298 K, with the majority starting at approximately 297 K. We did not observe a strong effect on the main evolution of the instability and the velocity field. Examples of the temperature time series can be found in the supplementary material.

We assumed that at the beginning of each experiment the concentration of ethanol in the gas phase was negligible. We did not introduce ethanol into the chamber intentionally besides that coming from evaporation of the liquid inside the cell. Between each experiment, the chamber was opened to change the chip, allowing the ethanol that had evaporated during the experiment to escape.

In each experiment, the solution was delivered from a syringe connected through a capillary tube to the back channel of the chip (see figure 1b). The flow rate was controlled with a motorized syringe pump (Harvard Instruments, PHD 2000). In each experiment the chip was filled within a minute at a rate of $15 \ {\mathrm {\mu }}\textrm {l}\ \min ^{-1}$. After the chip was filled, we had a delay of approximately 30 s to 60 s before a replenishing flow was started to compensate the mass loss due to evaporation. This flow was kept constant throughout the rest of the experiment (the direction of the flow is depicted with blue arrows in figure 1b). Despite the hydrophobicity of the chip, the ethanol–water mixture slightly wets the chip. The contact angle with the glass wall is within $75^{\circ }$ to $80^{\circ }$ and $75^{\circ }$ to $95^{\circ }$ with the silicon wall (see the supplementary material for more details on the estimations of the contact angle). Therefore, as long as there was liquid at the back channel, the Hele-Shaw cell would remain filled. The replenishing flow would then prevent the back channel from emptying. We found the right volumetric flow by slightly underfilling the chip such that a meniscus was visible at the back channel (figure 1b). A volumetric flux of approximately $0.43\ \mathrm {\mu }\textrm {l}\ \min ^{-1}$ kept the meniscus stationary. In experiments with dye, we used a slightly larger flow, $0.5 \ \mathrm {\mu }\textrm {l}\ \min ^{-1}$, which resulted in the chip eventually getting fully filled. However, there was no noticeable overflow at the observation region. Furthermore, the edge of the chip helped to keep the interface pinned. The pressure at the inlet of the chip was estimated to be approximately 19 Pa higher than at the evaporating edge. This pressure is smaller than the pressure drop across the meniscus at the evaporating edge ($10^{2}$ to $10^{3}$ Pa), in accordance with no overflow (see the supplementary materials for details on the calculation of the pressure drop).

We used rhodamine 6G (SigmaAldrich, dye content 99 %) and $1\ {\mathrm {\mu }}\textrm {m}$ particles (FluoSpheres, carboxilate modified, Life Technologies) for visualization and for $\mu$PIV (30 f.p.s.), respectively. Further details on the $\mu$PIV experiments can be found in Appendix A.3.

After each experiment the chip was cleaned by immersing it in a bath of acetone, isopropyl alcohol and ethanol, sequentially. In-between each bath we used a flow of nitrogen for drying.

The solutions were made at a 50 wt% ratio between Milli-Q water (produced by a Reference A$+$ system (Merck Millipore) at $18.2\ \textrm {M}\Omega \ \textrm {cm}$ (at $25\,^{\circ }\textrm {C}$)) and ethanol (Boom; 100 %(v/v), technical grade), prepared by measuring the mass of each component using a balance (Secura 224-1S, Sartorius). The solution was not degassed before experiments.

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. 1964; Köllner et al. 2015). 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. 2014). Such a phenomenon was also observed in evaporation experiments inside a Hele-Shaw cell by Linde et al. (1964).

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. 2015). 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. (1964) and Schwarzenberger et al. (2012). 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. 1964; Shi & Eckert 2006; Schwarzenberger et al. 2014; Köllner et al. 2015). 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 1998)). 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. Model equations and their numerical simulations

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. 2015; Li et al. 2019a).

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. 1997; Bratsun & De Wit 2004)

(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. 2017b) or the full diffusion equation (Diddens et al. 2017a) 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. 2010)

(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. 2010)

(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. (1964) and Machrafi et al. (2010), who found that the solutal effect is dominant. In particular, Machrafi et al. (2010) 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. 2017b).

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 2006). 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. (2017b) for details). The general numerical technique is described in more detail and compared with various experiments in, for example, Diddens et al. (2017b), Diddens (2017) and Li et al. (2019a,b). 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. (2013) 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. (2013), where the bulk flow caused by the Marangoni rolls hindered the flow towards the interface. Afterwards, the two adjacent rolls merged into one.

4. Linear stability analysis

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 (1959), 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. 1997; Bratsun & De Wit 2004; Machrafi et al. 2010). 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 (1959),

(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 (1959), 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 (1959). 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 2002).

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 (1959). 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 (1959) 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 (1959) 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 (1959), 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 (1959), 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. (2013) 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 (2004) 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 2001). 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.

5. Summary, conclusions and outlook

We have investigated by means of experiments, numerical simulations and linear stability analysis the evaporation of a solution of ethanol and water confined in a Hele-Shaw cell. The selective nature of the evaporation process, due to the larger volatility of ethanol, resulted in the appearance of a Marangoni instability at the liquid–air interface. This instability results in the appearance of convective rolls in the bulk of the liquid phase, which caused the accumulation of dye in arch-like regions. With ongoing time, we observed the rolls to grow and merge until they reached the lateral size of the cell. Eventually, the instability stopped and the flow ceased.

Qualitatively, such phenomena were observed before in experiments by Linde et al. (1964) for different evaporating binary mixtures confined inside a Hele-Shaw cell, however, here we measured and analysed the growth of the rolls in a quantitative way by means of confocal microscopy, $\mu$PIV and numerical simulations. We found that the size of the arches increases following an overall square root of time dependence, similarly to other multicomponent liquid–liquid systems (Köllner et al. 2015). Furthermore, by looking at averaged velocity profiles, we found that the velocity in the rolls decays exponentially towards the bulk of the liquid.

For numerical simulations and linear stability analysis we considered a quasi-2-D model which results from an average over the direction normal to the walls, assuming a parabolic profile for the velocity and a constant mass fraction in that direction. In this way we took into account the drag caused by the walls without having to use the full 3-D equations. In spite of these simplifications, simulations and experiments showed a good quantitative agreement supporting the applicability of the quasi-2-D model at long enough times.

Due to the merging process that characterizes the evolution of the instability, which also slows down with time, we cannot be sure if the first rolls that we observe experimentally represent the onset of the instability. In fact, from simulations we were able to corroborate that the merging process can start at shorter times than those accessible experimentally. Therefore, for comparison with linear stability theory, we relied only on the results from simulations.

By looking at the perturbation of the mass fraction at the interface at early times, we were able to calculate the growth rate of the perturbation as a function of the different wavenumbers present in it. This allowed us to compare these calculations with the predictions that we made by linear stability theory. The results are in reasonable agreement, taking into account that simulations consider the complete nonlinear equations with material properties that depend on mass fraction, while for linear stability analysis we (i) made use of a simplified model with constant material properties, (ii) disregarded the velocity caused by the evaporation of both ethanol and water, (iii) did not include the evaporation of water and (iv) did not consider a time evolving base state.

The consideration of the drag of the walls of the Hele-Shaw cell in our analysis resulted in the possibility of stability for conditions that would be considered unstable under the purely 2-D model of Sternling & Scriven (1959). More specifically, for $\textit {Ma}_{s}^{*} < \textit {Ma}^{*}_c=24\mathcal {M}$ the system is expected to be stable for all values of the wavenumber $\kappa$. For large Marangoni numbers the quasi-2-D model tends to the purely 2-D model, making it unnecessary to consider the drag of the walls, at least in terms of the linear stability prediction. However, at the same time at large Marangoni numbers, the model predicts wavelengths that can be smaller than the gap of the cell, suggesting that the full 3-D equations need to be considered in order to account for the onset of the instability. Nevertheless, as long as $\textit {Ma}_{s}^{*} \leqslant 8\mathcal {M}({\rm \pi} + \sqrt {({\rm \pi} ^{2}+3)})^{2}$, the wavelengths predicted by the Hele-Shaw model become smaller than the cell thickness. Therefore, the dispersion relation given by (4.14) could be used as a first estimate for the initial size of the instability, as long as thermal effects stay small. Moreover, as time goes on and coarsening takes place, the restriction which the gap size imposes is overcome, and, as we mentioned before, the numerical results show good agreement with experiments, even at large values of the Marangoni number. It would be interesting to further test the predictions of the stability analysis using systems with lower solutal Marangoni numbers.

In future work we plan to apply our results to evaporating ternary systems, also confined in a Hele-Shaw cell. In those cases the Marangoni instability will determine the location of spontaneous emulsification (ouzo effect) and could compete with other mechanisms of pattern formation. In particular, this competition could preclude the formation of branch-like structures, composed of emulsion droplets, that were observed for a ternary miscible liquid–liquid system (Lu et al. 2017).