Interdiffusion of liquids of different viscosities in a microchannel

We perform a detailed study of the interdiffusion of two miscible liquids of different viscosities, water and a water/glycerol mixture, flowing side by side within a pressure-driven microfluidic flow. Using Raman imaging, we measure cartographies of the concentration in glycerol for different imposed flow rates. The analysis of these experimental data confirms a classical diffusive mixing since the extent of the diffusion layer follows a power law regime with the streamwise coordinate x. However, the interdiffusion layer is not located in the center of the microchannel, and its mean position ym evolves with x during the mixing. We also derive a theoretical model and present analytical arguments to explain the above results. We demonstrate that the displacement of the interdiffusion zone ym(x) is due to the coupling between hydrodynamics and the mixing through the dependence of the viscosity with the glycerol concentration. Eventually, we perform numerical simulations of the theoretical model and compare the solutions with the experimental data in order to estimate quantitatively an interdiffusion coefficient.


Introduction
Microfluidics offers a wide range of tools to investigate chemical reactions with an unprecedented control of the transport phenomena (for a review, see [1,2]). More specifically, coflow devices as depicted in figure 1(a) recently received a considerable amount of interest. In such experiments, two streams of miscible liquids are brought at a junction to flow side by side within a microchannel. For small channel dimensions, flows are laminar and the liquids mix by transverse diffusion and possibly react. It is thus possible to relate the distance x to the time elapsed since the two streams were put into contact, and the evolution along the flow contains the information related to the diffusion and reaction dynamics. Schematically, the distance x from the junction corresponds to time t = x/v, where v is the average velocity in the channel. Such experiments have been performed to estimate diffusion coefficients [3,4], rate constants [5], or investigate enzymatic reactions [6]. Such studies are often limited to dilute solutions for which there is no coupling between the flow and the diffusion process.
To illustrate this last point in more detail, let us first consider the interdiffusion of a dilute solution and its solvent injected at a same constant flow rate in the two arms of the device depicted in figure 1. In the dilute regime, the viscosity, the density and the diffusion coefficient D of the solution do not depend on the concentration φ in solute. In a geometry with a high aspect ratio w/ h as in figure 1(a) and in the dilute regime, the velocity profile is Poiseuille-like across the channel height and uniform over the width w [2]. At high Péclet number Pe = vh/D, where v is the mean velocity in the channel, the convective-diffusion equation for the transport of the solute simply reads v∂ x φ(x, y) = D∂ 2 y φ(x, y), where φ is assumed to be uniform over the height of the microchannel [7,8]. From equation (1), the concentration in solute φ evolves as where σ (x) 2 = Dx/v corresponds to the width of the interdiffusion layer and φ 0 is the initial concentration in solute. Measurements of the concentration field φ(x, y) in the main channel therefore lead to estimations of the diffusion coefficient D.
When the two liquids considered in the experiments have different viscosities and densities, the analysis of the interdiffusion is more complex. Actually, it is possible to describe the flow simply only when the two fluids have not significantly interdiffused, i.e. for small distances x, or at large imposed flow rates [9]. In that case, and for a geometry with a high aspect ratio as in figure 1, the two incoming streams fill different widths w A and w B (see figure 1(b)) related by: The velocity profiles in the two streams are Poiseuille-like in the height of the channel, uniform over their width, and the mean velocities are linked by: η A v A ≈ η B v B , providing w A and w B are significantly larger than h. This relation only holds in the case of an infinite aspect ratio (w/ h 1), and exact formulae can be found for rectangular microchannels [9]. For downstream distances, the two liquids mix by transverse diffusion, but it is difficult to predict the interdiffusion simply since the viscosity, the density and the interdiffusion coefficient evolve during the mixing.

4
In a previous work [10], we used Raman confocal imaging to probe the interdiffusion of various miscible liquids with chloroform in a silicon-glass Y-shaped device. Raman spectroscopy is often used to monitor chemical reactions within microfluidic devices, since it gives access to spatially resolved concentration maps of chemically reactive fluids [11]- [16]. This previous work showed first that the interdiffusion zone between two liquids of different viscosities is not centered in the main channel, and its position is roughly given by equation (3). These experiments also showed that the width σ of the interdiffusion layer widens according to a diffusive behavior σ 2 ∼ x. The main conclusion of that preliminary work was that a comprehensive description is necessary to analyze the interdiffusion process quantitatively.
Our aim is to investigate experimentally and theoretically the interdiffusion of two fluids flowing side by side in a microchannel, and for which variations of the viscosity and of the density occur during the mixing. We chose the interdiffusion of glycerol/water mixtures with water, since the dependence of the viscosity of water/glycerol mixtures with the glycerol concentration is nonlinear, and since the variation of density is well tabulated (see appendix A). The motivations of our work are also to derive a model for the mixing, and eventually to describe the interdiffusion using both numerical and analytical arguments. In the following, we first demonstrate that Raman imaging is a powerful technique to measure the glycerol concentration in aqueous solutions flowing in poly(dimethylsiloxane) (PDMS)-based microfluidic channels, with a precision level of ±1%. We then present experimental concentration maps of the interdiffusion of water/glycerol mixtures in a Y-shaped device. For all the investigated flow rates, the interdiffusion layer widens downstream according to a diffusive behavior σ 2 ∼ x, and surprisingly its mean position significantly moves along the direction y during the mixing. Then, we derive a model to describe the interdiffusion of these two liquids in such a geometry, and discuss some features of the solutions of these equations. More precisely, we provide analytical arguments to explain the displacement of the interdiffusion zone in the transverse direction y, due to the variations of the viscosity during the mixing. Eventually, we solve this model numerically and compare the numerical solutions with our experiments to quantitatively estimate an interdiffusion coefficient. Figure 2 summarizes the main features of our experimental work. Pure water and a mixture of water and glycerol at a given volumic fraction in glycerol φ 0 are injected at constant flow rates in a Y-shaped microchannel. The small length scales involved ensure a laminar flow in the main channel, and the fluids mix downstream only by interdiffusion. Microchannels are fabricated using the classical soft lithography technique [17], and dimensions of the main channel are typically h = 60-80 µm, w = 1000 µm. Such high aspect ratios ensure that velocity profiles are Poiseuille-like across the channel height [2]. The network of microchannels in the PDMS matrix is stuck onto a thin glass slide (≈ 200 µm) for microscopic observation. Spectroscopic measurements are performed in the main channel with a Raman confocal microscope (LabRam HR, Jobin-Yvon) equipped with an argon-ion laser (wavelength 514 nm, grating 600 cm −1 ).

Microfluidic experiments
Typical spectra acquired downstream of the main channel are shown in figure 2(b). Such spectra have been obtained using a high numerical aperture objective (100×, water immersion, NA = 1, Olympus) using a small confocal pinhole (200 µm) at acquisition times ranging between 1 and 2 s. These spectra correspond, respectively, to pure water (left), the initial solution of water/glycerol at φ 0 = 44.5% (right), and an intermediate spectrum at an unknown concentration φ u indicating the mixing of the two liquids (middle). We detail below a method to estimate the absolute concentration φ at a precision of 1% in such PDMS-based microchannels from such data. Figure 3(a) displays Raman spectra of different mixtures of water/glycerol at known volumic fractions φ measured in vials. The vibrational frequencies associated with hydrogen motion in liquid water lead to a broad peak at 3500 cm −1 (O-H stretching modes), whereas the glycerol contributions correspond to well-defined peaks at 2890 and 2950 cm −1 . In this figure, all the spectra are divided by their maximal signal in the 3200-3600 cm −1 spectral region, so that the water contribution is normalized to one. This figure shows that one can extract the volumic fraction φ u of a mixture with an unknown composition, by comparison of its Raman spectrum with these reference data. However, the Raman spectra recorded in the microchannel (see figure 2(b)) display a small contribution of the PDMS matrix that may interfere with precisely extracting the volumic fraction φ. Indeed, as shown in figure 3(b), reticulated PDMS displays intense Raman peaks at 2900 and 2970 cm −1 , in the spectral region of the glycerol vibrations. Despite the confocal nature of our setup, the size of the sampling volume focused in the main channel is not small enough to remove completely the intense PDMS signature, and a small contribution of 5-10% of the maximal intensity is always observed (see figure B1 in appendix B for the confocal analysis).  We show now a simple method to estimate the volumic fraction from spectra recorded in our PDMS device. Each measured spectrum I r , as shown on a typical example in figure 4, contains two contributions I r = I p + β I φ u , where I p is the contribution of the PDMS matrix, and I φ u the normalized contribution of the water/glycerol mixture with volumic fraction φ u . β is a numerical constant, and the different possible I φ are displayed in figure 3(a). Firstly, the raw spectrum I r can be normalized in the same way as the data shown in figure 3(a) since there are no PDMS vibrations in the 3300-3600 cm −1 region. The normalized spectrum now reads I n r = α I p + I φ u (α is a numerical constant). Secondly, we perform subtractions of I n r with the different reference spectra I φ . For the spectrum corresponding to the unknown composition φ u , this subtraction should only correspond to the contribution of PDMS α I p . To estimate this fraction φ u , all the subtractions are thus compared with a reference spectrum of the PDMS matrix. Such a comparison is performed using a simple least-square method on the data normalized by their maximal intensity in the 2850-2950 cm −1 region (see figure 4(b)). This method gives estimations of the best composition φ u that fits correctly the PDMS signal. Typical χ 2 are shown in the insert of figure 4(b). This fitting procedure allows estimations of the local composition φ with a precision level of about ±1% and is independent of the optical configuration of the setup (optical indexes, focusing, . . . ) since we only compare the shape of the spectra.

Interdiffusion experiments
In our experiments, pure water and a water/glycerol mixture at φ 0 = 44.5% are injected at constant flow rates Q w and Q g , respectively, in the two arms of the microdevice. The contrast of viscosity between the two incoming streams is about 6, as η w ≈ 1 and η g ≈ 6 mPas for these two liquids (see figure A1(a)). Raman spectra are acquired on different lines at typical downstream distances x ranging between 0 and 20 mm and with a transverse step y of 5-10 µm (the origin of the x-y positions is shown in figure 2). Acquisition times are about 1-2 s, and a complete map x-y is acquired in less than 15 min. Post-processing of the data leads to the concentration fields in glycerol φ(x, y) for different imposed flow rates Q w and Q g .  In order to account for these data, and following the classical results of the diffusion of a dilute solute in its solvent (see equation (2)), we tried to fit the concentration maps according to: where σ (x) and y m (x) are both free parameters corresponding, respectively, to the width of the interdiffusion layer and its position at each position x. In this approach, the diffusion region is not centered at y = 0 as in equation (2), but located at y = y m (x) to account for the observed asymmetry of the concentration profiles. This formula nicely fits all our data for different flow rates Q w and Q g as can be seen in the examples shown in figure 5(a).
For distances x for which the two fluids have not significantly interdiffused and downstream of the entrance length L e of the channel (L e w [9]), the flow is supposed to be well described by figure 1(b), and the widths of the two incoming streams related by equation (3), i.e. η w Q w /w w ≈ η g Q g /w g . The insert of figure 6(a) shows the measured width w g = w/2 + y m (x = 0) of the incoming stream of the water/glycerol mixture at x = 0 for the different experiments and normalized by w the width of the channel. The linear relation )] −1 confirms the validity of equation (3), and therefore the picture of the flow shown in figure 1 For the flow depicted in figure 1(b), the natural timescale to compare different experiments cannot be defined as usual as /w the mean velocity in the channel, but rather with a velocity that does not depend on the initial widths of the incoming streams. We thus choose the definition τ = 2x/(v A + v B ) with the velocities v A and v B defined by equation (3) to rescale our data. Figure 6 shows different sets of σ (x) and y m (x) obtained for several experiments performed at different flow rates Q w and Q g . The width σ (x) and the position y m (x) of the interdiffusion layer are plotted against the characteristic time τ defined above. Figure 6(a) shows remarkably the collapse of the widths of the interdiffusion layer σ (x) when plotted against τ . More importantly this plot demonstrates diffusive behavior σ 2 ≈ D eff τ for all the investigated flow rates with an effective interdiffusion coefficient D eff ≈ 7 × 10 −10 m 2 s −1 . Figure 6(b) also evidences a non-intuitive result: the position of the interdiffusion zone y m (x) decreases as the streamwise coordinate x increases. It means that the width of the stream which is rich in glycerol decreases during the interdiffusion. In order to understand such behavior, and also to relate quantitatively the experimental coefficient D eff to the real interdiffusion coefficient for such fluids, we need to investigate theoretically the mixing of these two liquids in the microchannel.

Theoretical model and analysis of the interdiffusion in the microchannel
We now present a model describing the mixing of two liquids with different viscosities and densities in our microfluidic device, details are given in appendix C. The basic equations describing the hydrodynamics of mixtures are based on the work of Joseph et al [18] (see also Brenner [19]). In this approach, the conservation of mass and the Stokes equation are applied to a particle of a two-fluid mixture. However, two different velocities can be used for this description: the mass-averaged velocity v and the volume-averaged velocity u. One can demonstrate for the water/glycerol case that v is not solenoidal (i.e. div v = 0) since the density of the mixture changes with mixing [18]- [20], but that u verifies div u = 0 because the water/glycerol system can be considered as a simple mixture, i.e. its volume is conserved during mixing (see figure A1(b)). When the Stokes equation is written with the mass-averaged velocity v, it involves additional terms detailed in appendix C since div v = 0 [18].

Lubrication approximation
In our experimental case, the geometry of the flow is characterized by a high aspect ratio h w (see figure 2). We thus adopt the classical lubrication approximation to derive a simpler model describing our experiments in the x-y-plane, and we demonstrate in appendix C that the steadystate equations describing the interdiffusion read: with U = (U x , U y ) the averaged velocity of u over the height of the microchannel and ∇ = (∂ x , ∂ y ) the spatial gradient in the x-y-plane. Equation (5) describes the transport of φ with the convective φ U and diffusive −D(φ)∇φ fluxes. Equation (6) relates the velocity U to the local pressure gradient ∇ . The coupling between the transport equation and the hydrodynamics comes here from the dependence of the viscosity η on the volumic fraction φ. Note that the viscosity η strongly varies with φ in a nonlinear way as shown in figure A1(a). Equation (7) is the conservation of the flow rate, and permits to compute the pressure field. In the following, we analytically derive arguments from this model to explain some of the observed experimental behavior.

Analytical arguments to describe the interdiffusion process
Solving equations (5)-(7) to obtain the concentration field φ(x, y) requires numerical methods and also the behavior of the interdiffusion coefficient D with φ. However, simple assumptions can be done in order to understand the mixing qualitatively, particularly the results displayed in figure 6. First, the velocities v in the channel have been varied in the 1-3 mm s −1 range in the experiments described above; for h = 80 µm and D = 7 × 10 −10 m 2 s −1 , the investigated Péclet numbers Pe = vh/D roughly range from 100 up to 300. We can thus in a first approximation neglect the longitudinal diffusive flux D(φ)∂ x φ compared with the convective flux φU x [7]. The conservations of φ and of the flow rate thus read: since the no-flux boundary conditions ∂ y φ = 0 and U y = 0 apply at y = ±w/2. Then, we assume that the pressure field can be approximated by (x, y) = α(x) + (x, y), where (x, y) is a first-order correction coming from the coupling between hydrodynamics and the transport of φ. Indeed, the transverse velocities U y ∼ ∂ y /η(φ) are expected to be small compared with the longitudinal ones U x ∼ ∂ x /η(φ). This assumption probably holds only for small contrast of viscosity between the two incoming streams. Using such an approximation and the incoming conditions given by equation (3), the conservation of φ and of the flow rate leads to the relation: We also found experimentally that the concentration fields φ(x, y) display a sigmoidal shape, i.e. they are well described by φ = H (s), with s = [y − y m (x)]/σ (x), and H (s) a function that varies significantly only for −s 0 < s < s 0 , and is everywhere else almost constant. Figure 7(a) displays the case H (s) = φ 0 [1 + erf(s/2)]/2 found in the experimental section to fit our data correctly (s 0 ≈ 4-5 for that particular case). We can write the above relation (9) with the natural variable s and include the limits −s 0 and s 0 in the integrals: We now assume that the interdiffusion zone lies at positions y far away from the two lateral walls. One has therefore ∂ y φ = 0 at y = ±w/2, and thus y m + s 0 σ < w/2 and y m − s 0 σ > −w/2. The functions φ/η(φ) and 1/η(φ), as shown in figure 7, only vary significantly for −s 0 < s < s 0 , and we can thus safely assume that they are almost constant for s > s 0 and s < −s 0 , the relation (10) thus becomes with η w = η(φ = 0) the viscosity of water and η g = η(φ 0 ) the viscosity of the initial water/glycerol mixture at φ 0 = 44.5%. This last equation leads to a relation between y m (x) and σ (x) that takes the following form: where y 0 is the initial position of the interface between the two incoming streams given by equation (3) and is found after calculation: depends not only on the values of η g and η w , but also on the variations of the viscosity with φ, through the two integrals ∞ −s 0 ds φ/η(φ) and s 0 −s 0 ds/η(φ). In the case φ = H (s) = φ 0 [1 − erf(s/2)]/2, we can compute the above relation, and estimate the numerical values of . We also tried different forms for H (s) such as the erf function, or the tanh case, and we found that does not depend strongly on the exact form of H (s), and also on the values of s 0 . is computed from equation (13) with H (s) = φ 0 [1 − erf(s/2)]/2 and s 0 = 5. We also take into account the small correction due to the finite aspect ratio of our microchannel to estimate y 0 [9]. The nice agreement, particularly the decrease of y m (x) along the direction of the flow, is thus a strong argument supporting the above approximations. The displacement of the interdiffusion layer during the mixing (i.e. the variations of y m with x) is due to the coupling between the mixing and the velocity field through the dependence of the viscosity η with φ.

Comparison with the experiments
It is rather difficult to find limiting cases to discuss in depth the relation (13) since the above approximations probably hold only for small contrasts of viscosity between the two coflowing liquids. However, it is clear from relation (13) that small values of are found as soon as η(φ) does not vary strongly with φ. For higher contrast of viscosity, transverse velocities U y may not be negligible, and instabilities may also occur as recently observed by Cubaud et al in diverging flows of miscible fluids [21]. Moreover, capillary effects at the walls of the microfluidic channel, or even stresses induced by high gradients of concentration (the Korteweg stresses) [18,20] may be relevant to describe the mixing for high contrast of viscosity or for strongly different miscible liquids.

Numerical simulations
We now turn to the numerical simulations of the mixing of water/glycerol mixtures in the microchannel. Our aim is to obtain estimations of interdiffusion coefficient D by comparison of the numerical solutions with the experimental data. However, full 3D calculations including the complex shape of the Y-junction, for y = ±500 µm and from x = 0 up to x ∼ 2 cm, may be too fastidious to be performed in order to extract quickly estimations of D. In the following, we will thus again assume the lubrication approximation and solve equations (5)- (7) only. Note that this simpler model cannot take into account the entire complexity due to the exact shape of the geometry close to the Y-junction.
However, complete numerical simulations of equations (5)- (7), for values of x ranging between 0 and several centimeters still remain difficult to perform rapidly. Indeed, for small distances x when the two fluids that enter the main channel have not significantly interdiffused, one expects that the flow reaches the picture given by figure 1(b) on the entrance length L e (L e w [9]), i.e. two parallel streams of widths w g and w w with velocities U g and U w related by equation (3). In this first regime, gradients of φ are large and the velocity profile evolves rapidly with possible large values of transverse velocities. A high spatial resolution is therefore needed to capture the evolution of the flow and to avoid too large a numerical diffusion. However, this high spatial resolution may lead to extremely long calculation times to solve the equations up to large values of x. However, in the second regime, i.e. for x > L e , one can probably perform simulations with initial conditions given by the theoretical picture of figure 1(b), and with a poorer spatial resolution up to large values of x in order to describe the mixing correctly. This strategy will allow us to rapidly perform accurate numerical simulations of the mixing process for different interdiffusion coefficients and to compare the solutions with our experimental data. All the details concerning the numerical methods and the discretization of the equations are given in appendix D. Figures 8(a)-(c) display the numerical concentration profile φ and the longitudinal and transverse velocities U x and U y , for the initial conditions φ(x = 0, y) = φ 0 H e (y) (H e is the Heaviside function), and for flow rates Q w = 250 µl h −1 and Q g = 150 µl h −1 . For these simulations in the x = 0−1 mm range, the diffusion coefficient is set at D = 5 × 10 −10 m 2 s −1 and is assumed to be independent of φ. In this regime of small x and for the velocities investigated, the two fluids have not significantly interdiffused as expected, and the velocity profiles present large transverse components U y . The slight diffusion that can be noticed on the concentration profile (≈20 µm) is only due to the numerical diffusion and corresponds to about 10 pixels. For x > L e ≈ w, the concentration profiles reach approximately the following form φ = φ 0 H e (y − y 0 ), i.e. two parallel streams of widths w g = w/2 − y 0 and w w = w/2 + y 0 flowing side by side in the microchannel. Downstream of L e , the velocity profile is almost along the channel direction (U x U y ), and is approximately given by U x = (U g − U w )H e (y − y 0 ) + U w . It is straightforward to show that the numerical values of U g and U e are given by Q g /(w g h) and Q e /(w g h), and that y 0 satisfies the theoretical relation (3) derived for the case of two streams of fluids flowing side by side in the same microchannel. The picture of the flow given in figure 1(b) is therefore valid downstream of L e ≈ w. initial conditions Q g and Q w , as can be seen in figure 8

Entrance regime x < L e
In the following, we will thus perform simulations of equations (5)-(7) up to large values of x, directly with the initial conditions given by equation (3), i.e. two parallel streams of widths w g = Rw and w w = w − w g with longitudinal velocities given by U g = Q g /(w g h) and U w = Q w /(w w h). These simulations can be performed safely with a poorer spatial resolution without any strong numerical diffusion to capture the interdiffusion.

5.2.
Interdiffusion regime x > L e Figure 9 displays a typical numerical solution of equations (5)-(7) in the interdiffusion regime (x > L e ) for the imposed flow rates Q w = 350 µl h −1 and Q g = 100 µl h −1 , and with an interdiffusion coefficient D = 5 × 10 −10 m 2 s −1 supposed to be independent of φ. In such simulations, we use the modified initial conditions derived from the expected theoretical picture of figure 1(b). This strategy allows us to describe correctly the mixing up to large values of x (here 20 mm) without any strong numerical diffusion and with a lower spatial resolution. The validity of this approach is due to the large difference between two different length scales: the entrance regime L e ≈ w and the mixing zone x ≈ Pe w. In this last regime, the interdiffusion layer widens downstream as expected, and its mean position evolves significantly during the mixing. It is straightforward to demonstrate that the width of the interdiffusion region σ (x) fitted from such simulations using equation (4) follows diffusive behavior σ 2 ∼ x. Figure 9 also demonstrates that the transverse velocities U y due to the coupling between hydrodynamics and the interdiffusion are small compared with the longitudinal ones (U y ≈ 20 µm s −1 at x = 2 mm), and vanish along the mixing process. These numerical results thus confirm the validity of the approximations derived in section 4.2. Again the mean position of the interdiffusion zone y m (x) estimated using equation (4) is well fitted by the theoretical relation (12) derived previously. All these results are illustrated below when comparing the simulations with the experimental data.

Comparison with the experimental data
We performed numerical simulations with different initial conditions as those displayed in figure 6, and for interdiffusion coefficients D ranging from 2 to 8 × 10 −10 m 2 s −1 . We chose for simplification interdiffusion coefficients that do not depend on the concentration φ. In all these simulations, concentration maps were fitted according to the erf formulae as in equation (4). In all cases we found diffusive behavior σ 2 = D eff τ , with the characteristic time of the flow τ defined previously in the analysis of the experimental data (see section 3). Figure 10 shows the estimated effective diffusion coefficient D eff against the interdiffusion coefficient D. We found a linear relationship as suggested by a simple dimensional analysis of the equations (5)- (7). Such a result allows us to estimate from our experimental data, an interdiffusion coefficient D ≈ 5 × 10 −10 m 2 s −1 in good agreement with other reported experimental values [22]. Figure 10 also reports the concentration maps of φ with D = 5 × 10 −10 m 2 s −1 and for the initial conditions presented previously in figure 6. The nice agreement with the numerical solutions, especially the displacement of the interdiffusion zone y m with x confirms the validity of the model, and our strategy to extract values of interdiffusion coefficients. Importantly, a φ-independent interdiffusion coefficient D is enough to correctly describe the experimental results. Nonlinearities may be relevant in the case of the interdiffusion of more concentrated water/glycerol mixtures with water.

Conclusions
We performed a detailed investigation of the interdiffusion of miscible liquids with different viscosities in a microfluidic Y-shaped channel. Our aim was to study the coupling between mixing and hydrodynamics through the dependence of the viscosity with the concentration. Using Raman imaging, we measured cartographies of the concentration in glycerol within the mixing of water and water/glycerol mixtures. The interdiffusion process is described by classical diffusive behavior (σ 2 ∼ x), but its mean position evolves in the transverse direction of the channel (y m depends on x). We explained these results thanks to an analytical model and the lubrication approximation. Particularly, we showed that the displacement y m (x) is due to the conservation of the flow rate and of the concentration. We also compared the numerical solutions of the model with the experimental data, in order to extract an interdiffusion coefficient D. It may now be interesting to vary the initial composition of the water/glycerol mixture, and test the variation of the estimated D with the glycerol concentration.
This method and the strategy described above may be useful to extract interdiffusion coefficients for different mixtures. We also believe that the observed experimental phenomena, namely the displacement of the interdiffusion layer along the mixing, may be relevant in the case of complex chemical reactions investigated using microfluidic tools.

Acknowledgments
We gratefully thank A Ajdari, F Boyer, M Colin, T Colin, J Leng and M Schindler for fruitful discussions. We also acknowledge Région Aquitaine for funding and support. verifies div u = 0, and the equation of transport reads with this velocity: with a new definition for the diffusion coefficient D.
At low Reynolds number, hydrodynamics is described using the Stokes equation: where is the pressure, η(φ) is the viscosity of the particle of fluid and D[v] is the symmetric strain tensor. The viscosity η displays nonlinear behavior with φ as shown in figure A1(a).
Since div v = 0, the Stokes equation is more complicated than for a unique fluid and involves additional terms [18,20]. We neglect in our description all possible stresses induced by high gradients of concentration (the Korteweg stresses) [18,20]. We believe that for the interdiffusion between pure water and water/glycerol mixtures up to φ = 44.5%, such stresses are negligible, but they may be relevant for higher concentration gradients [26,27]. Moreover, in our Y-junction that does not present converging or diverging flows and for the rather small contrast of viscosity investigated (≈6), we do not expect viscous instabilities such as those recently investigated by Cubaud et al in diverging flows of miscible fluids presenting large viscosity contrasts [21]. The Stokes equation (C.3) reads with the solenoidal velocity u: with ∂ φ η the derivative of η with φ. Note that the Stokes equation does not present these complex terms when the volume-averaged velocity u is used in the symmetric strain tensor in equation (C.3) [19]. We show below that the differences between these two representations are negligible in the lubrication approximation.

C.2. The case of microfluidics: lubrication approximation
In our experimental case, the geometry displays a high aspect ratio, i.e. h/w ∼ with 1, and we can adopt the classical lubrication approximation. We thus make the change of variable z → z, and the continuity equation div u = 0, imposes that the velocities u z are of the order of , i.e. u z → u z . The transport equation (C.2) reads with the above scaling: This last equation demonstrates that φ is homogeneous over the height of the channel at first order in , i.e. ∂ z φ = 0. With the same scaling z → z, u z → u z , and with ∂ z φ = 0, the first component of the Stokes equation (C.4) reads: (C.7) The same equation holds for the y-component, and along the z-direction one has 1 ∂ z = 1 ∂ x (η∂ z u x ) + ∂ x (η∂ x u z ) + 1 ∂ y (η∂ z u y ) + ∂ y (η∂ y u z ) + 1 ∂ z (2η∂ z u z ) , (C.8) since ∂ z Q(φ) = 0. In order to balance the viscous and pressure terms, one has to impose the following scaling to the pressure → / 2 , and the Stokes equations thus become at first order in which is the classical result of the lubrication theory (for a mathematical formal proof, see [28,29]). The terms due to the compressibility of the mixture, div v = 0 in equation (C.3), are negligible in geometries with a high aspect ratio, and we could have thus chosen the volumeaveraged velocity u to estimate the strain tensor in the Stokes equation [19] and obtained the same results. We now define the velocity U in the plane x-y, as the average of u over h. By integration of the previous equations along the z-direction with the classical no-slip boundary condition, one eventually demonstrates that the interdiffusion is described at steady state by the equations (5)- (7).

Appendix D. Numerical simulations
The numerical methods used in order to solve the equations (5)-(7) are based upon the following evolution model: associated with the following boundary conditions, at the inlet: φ(x = 0, y) = φ 0 H e (y − y 0 ), (D.4) with y 0 , U g and U w given by equation (3). For the two side walls, at y = ±w/2: All the numerical results shown in this paper are in fact stationary solutions of this system. Note that we have to impose at the outlet that the mixing is total through the condition ∂ x φ(x = L , y) = 0. This condition implies that L has to be long enough in order to be satisfied. This again confirms the validity of our approach to solve the system (5)-(7) on long length scales with modified initial conditions. The viscosity η(φ) is estimated from the data shown in figure A1 using a polynomial of degree 5.

D.1. Time discretization
Let us consider a sequence of time steps t n = n t. By denoting φ n (x, y) the approximation of φ(t n , x, y), a time discretization of the system (D.1)-(D.3) can be written according to: where θ is a relaxation parameter taken as 0.5. The discretization described in equation (D.10) is explicit for the convective term and semi-implicit for the diffusion term. We approximate D(φ) in an explicit way to keep the linearity of the discretized system. Equations (D.2) and (D.3) being stationary, the time discretization is straightforward.

D.2. Space discretization
Let us divide the channel into a set of subdomains i j = [x i−1/2 , x i+1/2 ] × [y j−1/2 , y j+1/2 ] as shown in figure D1. Since the velocity U is described as a gradient of the pressure , one needs to compute the pressure field and the velocity field at different nodes in order to have a precise scheme (figure D1): this is the so-called marker and cell (MAC) mesh. The space discretization of equation (D.2) is hence: (D.14) In both equations (D.1) and (D.3), an operator of the type ∇(K ∇ ) appears ( can be or φ). This diffusion operator is discretized using the 5-points scheme: x y Figure D1. Schematic view of the two-dimensional MAC mesh. •, φ and nodes; •, φ and boundary nodes; , U x nodes; and ♦, U y nodes.
In both cases, the coefficient K depends on variables which are located at the (x i , y j ) nodes. It is hence necessary to interpolate K on the edges of the mesh ((x i+1.2 , y j ) and (x i , y j+1/2 ) nodes). An harmonic mean is used in order to ensure the flux continuity across the edges: The transport term U ∇φ in equation (D.1) is discretized by a WENO-5 scheme [30]. The main advantage of such a scheme is the significantly reduced numerical diffusion due to the discontinuous inlet conditions considered in the coflow case, equations (D.4) and (D.5).
In all the simulations shown in the paper, the boundary nodes (φ, ) are taken into account to evaluate the diffusion operator ∇(K ∇ ). By writing all the discretized space operators with the suitable boundary conditions (D.4)-(D.9), one can build up a (symmetric) linear system, which is solved thanks to a conjugate gradient method.