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

Interfacial stability is important for many processes involving heat and mass transfer across two immiscible phases. When this transfer takes place in the form of evaporation of a binary solution with one component being more volatile than the other, gradients in surface tension can arise. These gradients can ultimately destabilise the liquid-gas interface. In the present work, we study the evaporation of an ethanol-water solution, for which ethanol has a larger volatility. 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 this instability by observing the effects that it has on the bulk of the liquid. More specifically, the growth of convective cells is visualized with confocal microscopy and the velocity field close to the interface is measured with micro-particle-image-velocimetry. The results of numerical simulations based on quasi 2D 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 2D equations showed reasonably good agreement with the results from simulations at early times, when the instability has just been triggered and no coarsening has taken place. In particular, we find a critical Marangoni number below which a regime of stability is predicted.


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 et al. 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. (1977Sørensen et al. ( , 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 et al. (2016) considered the case of evaporation coupled with a Poiseulle 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 et al. (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 (Linde et al. 2013;Zhang et al. 2011), 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 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 2D. 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 et al. 2004;Shi & Eckert 2006, 2007Eckert et al. 2012;Schwarzenberger et al. 2012), numerical (Mokbel et al. 2017;Grahn 2006), 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 et al. (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 like chemical diagnosis, injekt 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 results in a Marangoni flow. Such a flow can then lead to quite rich phenomena, for example, it can make a droplet inside a bath jump against the pull of gravity (Li et al. 2019b, 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(Tan et al. , 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 et al. 2021). Clearly, these phenomena are closely connected with the presence of instabilities at interfaces where mass transfer and evaporative cooling take place. However, the threedimensional 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. Morever, 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 µ-PIV (section 2). In section 3, this system is modelled and numerical simulations of the model equations are performed. This model is quasi-2D and takes into account the drag caused by the walls of the cell (Köllner et al. 2015;Bizon et al. 1997;Bratsun & De Wit 2004). The simulation results are able to reproduce our experimental observations quite well, even though we do not consider the effects of evaporative cooling. In section 4 we present a linear stability analysis of a simplified version of the quasi-2D model equations. In section 5 we compare its results with those from numerical simulations at early times, showing not perfect, but close agreement. Finally, in section 6 we conclude and give an outlook on future work with ternary systems.

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 µ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 edge in figure 1(b) and (c)), a channel 50 µ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 side in figure 1(b) and (c)). The surface of the chip was made hydrophobic as described in appendix A.2.  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 µm. To maintain a homogeneous gap, a matrix of pillars was placed between the wafer and the glass, represented as thin vertical lines in (a) and as circles in (b) (not to scale). Additionally, two walls divided the cell in three different parts. At the back of the chip a channel of 50 µm in depth was used to continuously refill the cell and compensate for mass loses due to evaporation. As depicted in (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 (b), the gap region is represented in grey and the red rectangle indicates the region of observation. In (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 µm thick without edging and the glass plate is 500 µm thick. The holder cap was about 170 µm above the chip.

Experimental setup
We visualized the evaporation process at a region around the open edge (red square in figure 1(b)) 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 three-dimensional volume of about 3000 cm 3 . This volume was created by covering part of the microscope with plastic wrap (see SI). 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 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 ± 1% or 51 ± 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 relative humidity in the chamber and extra details on the chamber can be found in the SI.
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 about 2 K 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 around 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 SI.
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 experi-ments, the chamber was opened to change the chip, allowing the ethanol that 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 1(b)). 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 µL min −1 . After the chip was filled, we had a delay of about 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  1(b)). Despite the hydrophobicity of the chip, the ethanol-water mixture slightly wets the chip. The contact angle with the glass wall is within 75°to 80°and 75°to 95°with the silicon wall (see the SI 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 under-filling the chip such that a meniscus was visible at the back channel (figure 1b). A volumetric flux of about 0.43 µL min −1 kept the meniscus stationary. In experiments with dye, we used a slightly larger flow 0.5 µ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 about 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 SI for details on the calculation of the pressure drop). Furthermore, the edge of the chip helped to keep the interface pinned.
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 MΩcm (at 25°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.

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 1(b)). These regions are shaped as arches of approximately same sizes, suggesting the presence of an instability (Köllner et al. 2015;Linde et al. 1964). 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.
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 phenomenon was also observed in evaporation experiments inside a Hele-Shaw cell by Linde et al. (1964).   We measured the distance z d F (t) that the accumulated dye penetrated towards the inside of the cell as function of time (see third panel of figure 2(a) for an example of z d F ). Here the superscript d indicates a dimensional quantity. To measure z d F , we took the vertical average of the intensity field, resulting in an intensity profile as 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 d F . This process was then repeated for every time frame. In figure 2(b), we show some examples of z d F (t) on a double logarithmic scale. A close inspection of the plots shows that during a merging event z d F accelerates and then slows down. However, as an overall trend, z d F grows effectively as a power law. Fourteen independent runs were fitted to ∝ (t d ) α . An average exponent of α = 0.50 ± 0.04 suggests that z d F effectively grows in a diffusive-like manner. A similar behavior 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 SI for a plot with all cases).
We defined the time t d = 0 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 SI the secondary arches are clearly visible for a case where the liquid was seeded with tracer particles.
Eventually, the lateral size of the arches becomes comparable with the width of the central part of the chip, as shown in the first panel 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 panel of figure 3(c)), sometimes causing the remaining ones to start traveling along the edge. Overall, after 10-25 min the instability fades away and the dye is swept to the interface (third panel of figure 3(c)). 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 to stop 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.

Visualization of the bulk flow by micro particle image velocimetry
To visualize the flow caused by the instability, the solution was seeded with 1 µ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 video can be seen in movie 2 of the SI).
The direction of the flow can be read from panel (b) of figure 4 where we have plotted the corresponding normalized velocity vector field (averaged over 30 frames and a single experiment). The color code in this image corresponds to the horizontal component of the velocity, while the vertical component and the magnitude can be seen in the SI. We x d hju  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 semi-log 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 (b). All cases where obtained from velocity fields averaged over 1 s around the time indicated by the legend.
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 s Figure 5. Schematic of the numerical simulation. (a) The small distance δ 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.
interface is directed towards the center 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;Köllner et al. 2015;Schwarzenberger et al. 2014;Shi & Eckert 2006). 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 (∼ 50 µm from figure 4(a)). 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 to far away, going from some hundreds to about 20 µm 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 semi-logarithmic 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 d F . 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 SI).
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 × 10 −10 m 2 /s (Hansen et al. 1998)). For the typical velocity we can consider the maximum velocity inside the rolls (that we can measure) which is in the order of 10 −4 − 10 −3 m s −1 . The typical size of the arches is in the order of 10 −4 − 10 −3 m. This results in a range of Pe ∼ 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 (∼20 µm s −1 ) as the typical velocity caused by evaporation and the size of the arches just after the stability is over, namely about 2 mm, then we obtain Pe ∼ 100. The dye should indeed accumulate very close to the edge as advection overcomes the diffusion of the dye.

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.
The equations are made non-dimensional by using the distance between plates δ d as a length scale, the initial mass density of the liquid ρ d 0 as a mass per volume scale, and the initial kinematic viscosity ν d Therefore, the non-dimensional space coordinates, time, and velocities are given by r = respectively. Due to the presence of a binary solution in the liquid phase, the mass density of the liquid are dependent on the local ethanol mass fraction c e (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 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 three-dimensional Navier-Stokes equations with varying density and viscosity: along with the convection-diffusion equation of the ethanol mass fraction c e with the composition-dependent diffusivity D: In the above equations we have neglected the effects of gravity, given the small distance between the plates δ d (20 µm), which is the relevant length for gravity in a horizontally oriented cell. In particular if we calculate the Galilei number Ga = g(δ d ) 3 (ρ d 0 /µ d 0 ) 2 and the Archimedes number Ar = g(δ d ) 3 ρ d 0 ∆ρ d 0 /(µ d 0 ) 2 based on the plate distance, the density difference between pure ethanol and water (207.7 kg/m 3 at 293.15 K and 211.3 kg/m 3 at 298.15 K), the dynamic viscosity (3.1 mPas at 293.15 K and 2.4 mPas at 298.15 K), and the density (0.918 kg/m 3 at 293.15 K and 0.911 kg/m 3 at 298.15 K) at 50 wt%, we obtain Ga ∼ 0.01 and Ar = 0.0016 − 0.0026. 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 three-dimensional 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.
v(x, y, z, t) = f (y)u(x, z, t) with f (y) = (by)(1 − y) (see 5(a)). The two-dimensional velocity u is defined in the (x, z)-plane, i.e. u · e y = 0, and the constant b is set to b = 6 in the numerics so that the lateral velocity u coincides with the height average of the three-dimensional velocity v. With these assumptions, the three-dimensional Navier-Stokes equations (3.2) are simplified to two dimensions, resulting in a factor 6/5 for the non-linearity and a drag term −12µu stemming from the shear of the parabolic flow profile (Bizon et al. 1997;Bratsun & De Wit 2004): Here, ∇ = (∂ x , ∂ z ) is the two-dimensional nabla operator. The liquid properties ρ, µ and D are still allowed to depend on the local ethanol mass fraction, which is now also considered in its two-dimensional height-averaged projection, i.e. c e (x, z, t), which is subject to the two-dimensional convection-diffusion equation Considering only the two-dimensional 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 twodimensional 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 modeled as a mixture of air, ethanol, and water, for which we assume quasi-stationary vapour diffusion, i.e. we only consider the Laplace equations 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 g a = 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 in 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 two-dimensional projection as well. To mimic the evaporation in a three-dimensional 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 side walls. Alternatively, when only a section of the cell with some distance to the side walls is of interest, periodic boundary conditions in 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. This, together with the assumption of a static interface, yields the kinematic boundary condition 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 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.
where we have introduced the solutal Marangoni number (Machrafi et al. 2010) (3.10) and the Schmidt number with D d 0 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 ∆c δ = E d δ d = E over a distance equal to δ d is obtained. Then a "modified" Marangoni number (Machrafi et al. 2010) 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) The far-field boundary condition of the gas composition is set to that of ambient air, plus an optional relative humidity. At the liquid-gas interface, vapour-liquid equilibrium is imposed, i.e. where the equilibrium is calculated by Raoult's law generalised by activity coefficients, i.e.
( 3.15) Here, A α are the activity coefficients, n α the liquid mole fractions, M d α the molar masses, and p d α,sat the saturation pressures of the pure components, whereas R d is the universal gas constant, T d is the temperature of the system, and ρ g,d is the density of the gas. Notice that the equilibrium mass fraction changes in time through n α (c e (z = 0, x, t)) and A α (c e (z = 0, x, t)). The evaporation rates are finally given by the diffusive vapour flux at the interface, i.e. j e = ρ g D g e ∂ z c g e /Sc and where ρ g = ρ g,d /ρ d 0 and D g α = D g,d α /D d 0 are the gas to liquid density and diffusivity ratios.
As we mentioned in our 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 where T d is the temperature, ∆T d is a characteristic temperature change along the interface and D d 0,T is the thermal diffusion coefficient. The specific values used are shown in table 1.
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 Ma * s /Ma th ∼ 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).

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 discretised 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 e.g. Diddens et al. (2017b); Diddens (2017); Li et al. (2019a,b). In the present case, however, the method is generalised 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 further away from the interface. More details on the meshes used for different simulations can be find in the SI.
For the simulation used for comparison with experiments (section 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 SI. 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 SI.
In all simulations the initial mass fraction was c e = 0.5, the relative humidity 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.

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 SI. The three snapshots of 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.
After the arches are large enough we can see secondary arches too. From movie 3 in the SI 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 builts 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 3D 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 SI, 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 SI, 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 × 10 4 , the velocity far away from the interface is smaller than in experiments. This could be a 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 2) can reach much larger values than those measured experimentally, about 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 γ(x) as can be seen in figure 6(d ). We show the surface tension γ(x) corresponding to the three snapshots in figure 6(a). The curves confirm the presence of gradients going from the center 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.

Simplified model
In this section we will perform a linear stability analysis of a simplified version of the model equations described in section 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 (Bratsun & De Wit 2004;Machrafi et al. 2010;Bizon et al. 1997). 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 ρ In other words ρ = µ = dγ/dc e = 1, and D = 1/Sc, which results in the simplified model equations for the liquid phase ∇ · u = 0 (4.1) and for the gas phase ∇ 2 c g e = 0 and ∇ 2 c g w = 0. (4.4) Similarly, the boundary condition for the tangential stresses at the interface (equation (3.9)) is simplified to (4.5) 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) = 0. In the experimental section we obtained Péclet numbers larger than 1 in odds with this assumption, however for the onset of instability a more proper length scale to be consider is δ d (the gap of the Hele-Shaw cell), which gives Pe ∼ 1 meaning that we are in the limit where we could actually neglect this velocity.
After neglecting the normal velocity, from equation (3.8) we have j e + j w = 0 at the interface, then equation (3.13) simplifies to j e = ∂ z c e /Sc. Therefore, introducing this expression for the mass transfer rate of ethanol into equation (3.16), we get where c sat = p d e,sat (T )M d e /(R d T d ρ g,d 0 ). An analogous expression can be find for water.

Linear stability of 2D equations
For the base case we will consider no flow in the liquid phase, i.e. u u 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), c g e,0 = E g z + C g 0 , (4.9) with E and E g made non-dimensional with δ 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 to 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, u u u = ǫû û u 1 (z)e iκx+σt (4.10) p = p 0 + ǫp 1 (z)e iκx+σt (4.11) c e = c e,0 (z, t) + ǫĉ e,1 (z)e iκx+σt (4.12) c g e = c g e,0 (z, t) + ǫĉ g e,1 (z)e iκx+σt , (4.13) where ǫ is a small number, κ is the non-dimensional wavenumber of the perturbation, σ its non-dimensional growth rate, andû û u 1 (z),p 1 (z),ĉ e,1 (z), andĉ g e,1 (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.4), resulting in the following dispersion relation (see appendix C for further details): (4.14) 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 Ma * s based on the plate distance δ d .
From this equation we can obtain the condition of marginal stability by setting σ = 0, resulting in (4.15) where M = 1 + c sat a 1 ρ 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 Ma * s , i.e. κ σ=0 = (Ma * s /8M) 0.5 where the dependence on δ d is not present anymore, as κ σ=0 ∝ δ d as well as Ma * s ∝ δ d . This is equivalent to making the gap between plates large enough so that the limit of an unbounded system is reached (Martin et al. 2002).
We can now find a critical Marangoni number Ma * c ≡ 24M for which κ σ=0 = 0. Below this critical value, the right hand side of equation (4.15) becomes negative, while κ has to be positive. Therefore for 0 < Ma * s 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 dγ d 0 /dc 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.

Numerical solution of the dispersion relation
Equation (4.14) was solved numerically to obtain σ as function of κ and Ma * s for Sc = 7.51 × 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 σ to zero to facilitate visualization. Equation 4.15 is plotted as a solid white line, while the purely 2D limit, κ σ=0 = (Ma * s /8M) 0.5 is highlighted by a dotted line. The red solid line corresponds to the locations where σ is maximum for a given Ma * , i.e. the fastest growing modes.
Indeed, for Ma * s < Ma * c = 121.2 the system is expected to be stable, as σ < 0 for every value of κ. However, we must notice that for regions where σ < −κ 2 /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 σ, without success. While our search was not exhaustive, we are in the region where Sternling & Scriven (1959) only found real solutions for the purely 2D 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 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 to the purely 2D case is only visible for Ma * s 1000. Meaning that for larger values the purely two dimensional 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 δ d is large enough, but not so large that gravity starts to play a role, the use of the purely 2D model is sufficient.
According to Sternling & Scriven (1959), if the conditions D d /D g,d < 1, ν d /ν g,d < 1, and dγ d /dc 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 Ma * s < Ma * c . We also looked at the effect of the Schmidt number, Sc, as shown in figure 7(b) for a fixed Ma * s = 1.54 × 10 5 . As mentioned by Sternling & Scriven (1959), varying Sc shows only a minor effect which was already suggested by the fact that κ σ=0 is not a function of this parameter. Only for Sc 50 there is a slight increase on the value of the fastest growing wavenumber. A plot of σ as function of κ is displayed in figure 7(c). This case corresponds to the dashed red line in both panels (a) and (b) of figure 7 (Ma * s = 1.54 × 10 5 and Sc = 7.51 × 10 3 ). We have selected this particular case, because it corresponds to the simulated case in section 3.3. In order to know the right 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 ρ g D g e , it affects the range over which the systems becomes stable, as clearly seen from equation (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 σ decrease with increasing values of the factor c sat a 1 ρ 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 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 σ to decrease.

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 5 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 function of space and time. In particular, these are results for a case of Ma * s = 1.55 × 10 5 . In figures 8(a) and (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 5 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. c e (z = 0) x (shown in figure 8(c)). 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.
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 about 7 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 8(c)) 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 non-linear effects.
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 SI. 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, specially since merging starts very soon. In order to answer this question, we calculated the Fourier transform of the edge profiles (figure 8(d )) at each time. As an example, the profile at t = 6.8 is plotted again in frequency space in figure 11(a) of appendix D. Then we looked at the power as function of time as shown in figure 11(b). There is a range where the power growth is approximately exponential, and after some time, it saturates and stays at an approximately constant value (at least for this small range of time).
By fitting a line (in semi-logarithmic scale) to the region shadowed in figure 11(b), we obtained the growth rate of every single wave number present in the mass fraction profiles. This relies on the assumption that for some time many wavelengths will be present, and not only the fastest one. The corresponding plot of growth rate σ versus wave number κ is shown in figure 10 for different 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 section 3.2. We varied the size of the domain in the x and z directions, but always keeping the domain at least 5 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.
Given that in simulations we control the evaporation rate through the size of the gas domain, in order to determine the modified Marangoni number Ma * s of each case, we obtained E from the average diffusive flux at the boundary ∂ z c e x . Figure 9 shows the evolution in time of ∂ z c e 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 E ≡ ∂ z c e x was then used to calculate Ma * s together with material quantities corresponding to c e = 0.5. The corresponding values of E and Ma * s obtained for each different run are shown in the SI. The values 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 function of wavenumber for each case. The results are plotted as solid red lines in figure  10. The prediction tends to underestimate κ by a factor of about 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 non-linear terms and water evaporation are taken into account.
Notice that for large values of κ, 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) that the power of large values of κ (κ > 100 in that example) is quite small. Therefore, above a given value of κ we cannot trust anymore the obtained values of the growth rate. On top of each of the plots from figure 10 we have indicated the corresponding values of the dimensional wavelength λ d = 2πδ d /κ. Note that for the two largest Ma * s the fastest growing wavelength is smaller than the gap between the plates of the cell (δ d = 20 µ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 lead 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 equation (4.15), we can obtain a value for Ma * s for which κ = 2π, meaning that λ d = δ d . This implies that Ma * s = 8M(π+ (π 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 2D case (from the point of view of linear stability analysis) suggesting that for large Marangoni numbers a 3D model would be more suitable in order to look at the onset of the instability, even if the gap is as small as 20 µ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 equation (4.14) could be used as a first estimation of the growth rate of long wavelengths at onset.

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, µ-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-2D 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 3D equations. In spite of these simplifications, simulations and experiments showed a good quantitative agreement supporting the applicability of the quasi-2D 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 its growth rate as 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 non-linear 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, we (ii) disregarded the velocity caused by the evaporation of both ethanol and water, we (iii) did not include the evaporation of water, and we (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 2D model of Sternling & Scriven (1959). More specifically, for Ma * s < Ma * c = 24M the system is expected to be stable for all values of the wavenumber κ. For large Marangoni numbers the quasi-2D model tends to the purely 2D 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 3D equations need to be considered in order to account for the onset of the instability. Nevertheless, as long as the Ma * s 8M(π + (π 2 + 3)) 2 , the wavelengths predicted by the Hele-Shaw model become smaller than the cell thickness. Therefore, the dispersion relation given by equation (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).
as in Zhang et al. (2006). To activate the surface, we dipped the chip in a bath of piranha solution (30% of hydrogen peroxide and 70% sulfuric acid) for 40 min at 75°C. We then washed it with plenty of Milli-Q water and ethanol. The chip was kept in an oven at 120°C for 2 h to eliminate traces of water. After cooling down, the chip was dipped in a solution of 0.5% OTS in toluene for 17 h. Finally the chip was sonicated for 10 min in chloroform and 15 min in each of the following solvents one after the other: toluene, acetone, isopropyl alcohol, ethanol and water. Finally, we dried the chip under a nitrogen flow.

A.3. µ-PIV
A suspension of fluorescent particles (FluoSpheres, Life Technologies, 1 µm, carboxilate-modified, 2% solid in aqueous solution, with 2 mM of sodium azide) was added at 0.9 wt% to the ethanol-water solution. This resulted in an initial volume fraction of particles of 0.016 % and mass fraction of sodium azide of 1.2×10 −4 wt%, meaning that we can disregard any effects from the latter even after 30 min of evaporation.
We recorded images with a confocal microscope (Nikon Confocal Microscopes A1 system) at 30 fps using a 488 nm laser and a pinhole size of 26.8 µm. This settings resulted in a vertical resolution of 8 µm. We made measurements at a single position in the y direction, approximately in the middle of the cell. By translating the position of the focal plane along the gap of the cell, we observed that the visible particles were always the same, but with different intensities. Therefore, despite the vertical resolution of 8 µm all the particles inside the gap were visible at the same time. Different particles moved at different speeds, which we attribute to the Poiseuille flow. Therefore, the velocity results shown in this work represent a weighted average of the velocity, where the dominant peak during the correlation step in the PIV process determines the velocity at a given position.
We calculated the relaxation time according to (Tan et al. 2019;Li et al. 2019b) where ρ d sol ≈ 910.73 kg/m 3 and µ d sol ≈ 2.5 mPa s are the density and dynamic viscosity of the ethanol-water solution, and ρ d p = 1050.00 kg/m 3 and d d p = 1 µm are the density and diameter of the tracer particles. Then the Stokes number is estimated to be St = t d 0 u d max /z d F , with u d max the maximum velocity at a given time, and z d F the roll height at that time. The velocity ranges from ∼ 10 −4 − 10 −3 m/s while the roll height goes from 0.1 − 1 mm, resulting in a Stokes number of at most St ∼ 10 −7 . If we consider that from simulations we observed that the velocity at the edge reaches up to 0.01 m/s and we combine this with the smallest roll that we can observe experimentally, then we get St ∼ 10 −6 , meaning that in any case the particles should follow the flow properly.

Appendix B. Values of different properties used for simulations and theoretical calculations
In by the thermodynamic model AIOMFAC (Zuend et al. 2008(Zuend et al. , 2011. On the other side, as also mentioned before, a Taylor expansion around an initial mass fraction of c e,0 = 0.5 was used for the theoretical calculations. The corresponding coefficients were calculated according to a 0 (c e,0 = 0.  Yano et al. (1988), corresponding to c e = 0.494 at 288 K. The gradient of surface tension with respect to temperature was obtained by fitting data from Vazquez et al. (1995).
For the case of E, as discussed in the main text, we did not used the value at t = 0 as there is no gradient initially. Instead we used the mean value E calculated over the same period used for calculating the growth rate σ as described in section D. The corresponding values for each simulation can be seen in the SI. equation (C 2) leads to the solutionŝ u z,1 = − κf p κ 2 − β 2 e −κz − e −βz , (C 6) in the z direction, andû in the x direction, with β 2 = κ 2 + 12 + σ. In these solutions we have already considered that the velocity should stay finite, assumed u z,1 (x, 0, t) = 0 as described in the main text, and took into account continuity ∇ ∇ ∇ · u u u = 0. Introduction of equation (C 6) into equation (C 3) leads tô where ζ 2 = κ 2 + σSc, B 0 = −EScf p κ/(κ 2 − β 2 ), and B 1 is a constant. Again we took a finite perturbation far away from the edge. The solution of equation (C 4) is given bŷ c g e,1 = B g 1 e κz , (C 9) with B g 1 a constant. Once more, we have consider that the perturbation tends to zero far away from the edge.
The expanded Raoult's law (4.7) for the base case results in C g 0 = c sat (a 0 + a 1 C 0 ), (C 10) meaning that for the perturbed mass fraction the condition readŝ c g e,1 (0) = c sat a 1ĉe,1 (0), (C 11) thus B g 1 = c sat a 1 (I(0) + B 1 ), (C 12) where I(0) = B0 p − B0 q , with p = κ 2 − ζ 2 and q = β 2 − ζ 2 . On the other side, the mass conservation of solute at the interface (4.6) for the base case results in E = D g e ρ g E g , (C 13) while for the perturbed mass fraction we have I ′ (0) + ζB 1 = −D g e ρ g κB g 1 , (C 14) with I ′ (0) = κB0 p − βB0 q , and D g e = D g,d e /D d 0 . In this way we can obtain B 1 and B g 1 in terms of B 0 , but more importantly the mass fraction at the edge, which in the liquid side readŝ c e,1 (0) = ζI(0) − I ′ (0) κc sat a 1 D g ρ g + ζ .
(C 15) As a last step we apply boundary condition (4.5), dû x,1 (0) dz = iκ M a Scĉ e,1 (0), (C 16) where again we have taken into account that u z,1 (x, 0, t) = 0. In this way, the dispersion relation is expressed as EMa s = κ 2 c sat a 1 D g ρ g + ζ κ 1 + β κ 1 + ζ κ β κ + ζ κ , (C 17) which becomes equivalent to equation (30) from Sternling & Scriven (1959), provided that we neglect surface viscosity, and consider that µ g,d ≪ µ l,d , ρ g,d ≪ ρ l,d , and D l,d ≪ D g,d , with the superscripts g and l referring to the gas and liquid phases, respectively. We solved the complete equation from Sternling & Scriven (1959) for a few cases and compared it with its simplified counterpart (C 17) resulting in negligible differences when taking into account the material properties of ethanol and water and keeping the surface viscosity term equal to zero. Finally, if we notice that Ma * s = EMa s , then equation (C 17) leads to equation (4.14) in the main text.

Appendix D. Calculation of the growth rate
Here we show typical plots used for the calculation of the growth rate σ as function of wavenumber κ. In figure 11(a) we show the power density function (PDF) as function of the normalized wavenumber κ. This PDF corresponds to the mass fraction profile at t = 10.3 from figure 8(d ). By calculating the PDF over different times, we can look at the evolution in time of the power for every single wavenumber (not only those for which the power has a peak). In figure 11(b) we show four examples corresponding to the vertical dashed lines in figure 11(a).
If plotted in a semilog scale, there is a region where the power can be approximately fitted by a straight line (shadowed region in figure 11(b)), reflecting exponential growth. The growth rate was obtained from a straight line fit in the shadowed region. We did not apply the fit for earlier times to avoid considering the period over which the concentration gradient is building up (transients, see figure 9), but also not at later times, because eventually either non-linear coupling or spectral bleeding causes all the wavelengths to grow at the same rate as the fastest wavelength. One can see that both the blue line with circles and the green line with right triangles in figure 11(b) eventually start growing at a faster rate before reaching a saturation value. However, in some cases we included part of these regions of non-linear growth (for the largest wave numbers), in order to have enough fitting points for the smallest wave numbers, which is the reason for the noisy behavior observed in figure 10 at larger values of κ. Figure 10 shows the average over different runs, where we have changed the mesh kind (see SI for a description of the different kinds of meshes), the resolution and the size of the simulation domain. The use of different resolutions and different domain sizes, caused that the specific values of κ obtained for the Fourier transform to be different between simulations. In order to obtain an average curve of σ(κ), we interpolated the values of σ over a set of values of κ common to every different run.
From figure 11(b) it is evident that the power is quite small in the region over which we made the fit (notice that the power represents variations on the mass fraction). However, by the time the value becomes large, either merging has already taken place or the fastest wavelength has already shadowed the growth rate of the others. Therefore, we made the fit at early times when different wavelengths have different growth rates. In movie 5 of the SI we show the evolution in time of the mass fraction variation showing that merging already takes place before the variations in mass fraction becomes large. This early merging makes it hard to observe the onset of the instability in experiments.
In some of our simulations the instability was triggered first in certain regions of the interface, while it started later in others. We discarded these simulations.