1. Introduction
Micro- and nanodroplet generation is of tremendous interest due to its wide range of applications in drug delivery (Gursoy & Benita Reference Gursoy and Benita2004; Attama & Nkemnele Reference Attama and Nkemnele2005; Devarajan & Ravichandran Reference Devarajan and Ravichandran2011), liquid–liquid extraction (Rezaee et al. Reference Rezaee, Assadi, Hosseini, Aghaee, Ahmadi and Berijani2006; Rezaee, Yamini & Faraji Reference Rezaee, Yamini and Faraji2010; Yu et al. Reference Yu, Chin, Chen, Zhang, Lo, Ayi, Yap, Kwong and Liu2010; Jain & Verma Reference Jain and Verma2011), (nano)material synthesis (Liff, Kumar & McKinley Reference Liff, Kumar and McKinley2007; Kumar, Liff & McKinley Reference Kumar, Liff and McKinley2008; Duraiswamy & Khan Reference Duraiswamy and Khan2009), catalytic reactions (Yabushita et al. Reference Yabushita, Enami, Sakamoto, Kawasaki, Hoffmann and Colussi2009; Shen et al. Reference Shen, Zou, Wang, Dryfe, Huang, Dou, Dai and Wang2014) and cosmetics (Xu et al. Reference Xu, Nie, Seo, Lewis, Kumacheva, Stone, Garstecki, Weibel, Gitlin and Whitesides2005; Lee et al. Reference Lee, Yoo, Cheng and Jeong2008; Yeh et al. Reference Yeh, Zhao, Lee and Lin2009; Kuehne & Weitz Reference Kuehne and Weitz2011), etc. One way to generate microdroplets is to utilize microfluidic devices such as T-junctions (Yeh et al. Reference Yeh, Zhao, Lee and Lin2009), flow focusing setups (Anna, Bontoux & Stone Reference Anna, Bontoux and Stone2003; Teh et al. Reference Teh, Lin, Hung and Lee2008; Seemann et al. Reference Seemann, Brinkmann, Pfohl and Herminghaus2011) or co-flowing devices (Utada et al. Reference Utada, Lorenceau, Link, Kaplan, Stone and Weitz2005; Serra & Chang Reference Serra and Chang2008; Shah et al. Reference Shah2008), where monodispersed microdroplets with well-defined properties could be generated successively. All these devices and methods utilize a top-down approach, in which a liquid jet or drop is split into smaller parts. This limits the smallest droplet size which can be achieved.
This limitation can be overcome in a bottom-up approach such as solvent exchange (Lou et al. Reference Lou, Ouyang, Zhang, Li, Hu, Li and Yang2000; Lohse & Zhang Reference Lohse and Zhang2015, Reference Lohse and Zhang2020; Zhang et al. Reference Zhang, Lu, Tan, Bao, He, Sun and Lohse2015), where a large number of micro- and nanodroplets are generated by nucleation out of an oversaturated solution. This method, also called nanoprecipitation or solvent shifting (Fessi et al. Reference Fessi, Puisieux, Devissaguet, Ammoury and Benita1989; Galindo-Rodriguez et al. Reference Galindo-Rodriguez, Allemann, Fessi and Doelker2004; Aubry et al. Reference Aubry, Ganachaud, Cohen Addad and Cabane2009; Lepeltier, Bourgaux & Couvreur Reference Lepeltier, Bourgaux and Couvreur2014; Hajian & Hardt Reference Hajian and Hardt2015), though commonly used, is much less well understood.
In solvent exchange, a good solvent of the target droplet component (the solute) is replaced by a poor solvent, where the two solvents are miscible. A typical example is an oil saturated aqueous ethanol solution being replaced by oil saturated water. Upon contact of the ethanol and water solution, the two solvents start to mix with each other. Due to the addition of water, the solubility of oil is lowered, and the subsequent oversaturation leads to droplet nucleation and growth. The micro- and nanodroplets can nucleate in the bulk (Vitale & Katz Reference Vitale and Katz2003) or on a hydrophobic surface. For solvent exchange in a microchannel, it has been found that droplets nucleated in the bulk tends to migrate to and then stay in the centre in a co-flowing device (Hajian & Hardt Reference Hajian and Hardt2015), where the droplet movement is controlled by solutal Marangoni flow and the composition of the mixture. On the other hand, for droplets that nucleated on the surface, their average volume is found to increases with the Péclet number 
$Pe$ as 
$\propto h^{3}Pe^{3/4}$ (Zhang et al. Reference Zhang, Lu, Tan, Bao, He, Sun and Lohse2015), where the flow rate 
$Q$ is included in 
$Pe$ and 
$h$ the channel height. Later, the effect of flow geometry (Yu et al. Reference Yu, Maheshwari, Zhu, Lohse and Zhang2017) and solution composition (Lu et al. Reference Lu, Xu, Zeng and Zhang2015; Lu, Peng & Zhang Reference Lu, Peng and Zhang2016) on the average oil droplet size were also qualitatively investigated. The mutual interaction between a multitude of surface droplets and the resulting effect on the growth dynamics was also studied (Xu et al. Reference Xu, Yu, Peng, Lu, Lei, Lohse and Zhang2017; Dyett et al. Reference Dyett, Kiyama, Rump, Tagawa, Lohse and Zhang2018). Despite all of these studies, a thorough understanding of the oversaturation pulse – which is crucial to droplet nucleation and growth – is still lacking, because of its transient nature and a lack of means to directly measure it.
To have a quantitative understanding of the oversaturation pulse, we study the total amount of oversaturated oil inside a Hele-Shaw-like microfluidic channel, by measuring the total volume of the oil droplets that nucleate out from it using confocal microscopy. In this paper, we are neither interested in the nucleation process itself nor in the droplet morphology, since they do not help to quantify the oversaturation pulse. A theoretical model for the total nucleated oil volume is developed, based on the ternary phase diagram and Taylor–Aris dispersion. The model accurately predicts the scaling behaviour of the total volume 
$V$ of oil with respect to the Péclet number, 
$V\propto Pe^{1/2}$, including the prefactor, in which the influence of the solution composition and channel geometry is reflected. However, to compare the prefactor with the experiments, we need to prevent the oversaturated oil – especially the bulk droplets that nucleated out of it – from leaving the channel. To achieve this, a porous region consisting of circular pillars is put in the middle of the channel. For channels with such a porous region, the prefactor can collapse different groups of data onto one universal master curve for different channel geometries and mixtures. For channels without the porous region, the measured oil volume is smaller than the theoretical prediction because some of the oversaturated oil – including nucleated droplets in the bulk – leaves the channel.
2. Experimental procedure and methods
Solvent exchange is performed in a thin square channel with various height 
$h$, width 
$W$ and length 
$L$, see figure 1 for the definitions and table 1 for the parameters. The microfluidic channel is made of a glass wafer covered on a silicon wafer which is decorated with an inlet, an outlet and, for some of them, a porous region in the centre. The porous region is made of an array of circular pillars with pillar diameter 
$d$ ranging from 
$6.7\ \mathrm {\mu }\textrm {m}$ to 
$8.7\ \mathrm {\mu }\textrm {m}$ (see figure 1, inset), and the pillar spacings in the transverse and axial directions 
$a$ and 
$b$ are varied to change the porosity 
$\phi =1-{\rm \pi} d^{2}/(2ab)$ (see table 1 for parameters). The length between the inlet and the porous region is 
$L$. For ‘smooth’ channels, i.e. those without the porous region, 
$L$ is the entire length of the channel. The whole channel is made hydrophobic by an OTS (octadecyltrichlorosilane) coating: 1 M hydrochloric acid is first pumped through the chip at 
$50\ \mathrm {\mu }\textrm {l}\,\textrm {min}^{-1}$ for 20 min by using a syringe pump (Harvard, PHD 2000). The chip is then put into a vacuum chamber at 1.8 mbar overnight to dry. A solution of OTS dissolved in hexadecane (Sigma-Aldrich, 
$\geq 99\,\%$) at 0.4 v/v % is pumped through the chip at 
$50\ \mathrm {\mu }\textrm {l}\,\textrm {min}^{-1}$ for 20 min. The chip is then sequentially cleaned by chloroform, toluene and ethanol, and finally dried in vacuum for use.
Figure 1. Side view of the experimental setup. Solvent exchange is performed in a square microfluidic channel with or without the porous region, which is made of circular pillar arrays. Solution A (yellow) is decane (the oil) saturated aqueous ethanol solution and is injected into the channel first. Solution B (light blue) is decane saturated water and is injected later at a fixed pressure. In the mixture between the two solutions, oil is oversaturated so that droplets (red) nucleate out of it, both in the bulk and on the hydrophobic walls. A flowmeter is used to measure the flow rate of solution B. The channel height is 
$h$, the width is 
$W$ (in the 
$y$ direction, not shown) and the channel length before the porous region is 
$L$. Solvent exchange is observed from below by confocal microscopy. The inset shows the bottom view. Pillar diameter 
$d$ and pillar spacings 
$a$ and 
$b$ are varied to change the porosity 
$\phi =1-{\rm \pi} d^{2}/(2ab)$. See table 1 for the parameters.
Table 1. Geometrical parameters of the microfluidic chips. See figure 1 for a definition of 
$h$, 
$L$, 
$W$, 
$d$, 
$a$ and 
$b$.
Solution A, which is rich in oil, consists of decane (oil) saturated aqueous ethanol solution. To make a solution A with the desired concentration, its ethanol-to-water weight ratio 
$w_{e,A}/w_{w,A}$ is first determined. The solution is prepared by first mix 
$\sim 400\ \textrm {g}$ mixture of ethanol (Sigma-Aldrich, 
$\geq 99.8\,\%$) and water (Milli-Q) with this specific weight ratio 
$w_{e,A}/w_{w,A}$, then added to decane (Sigma-Aldrich, 
$\geq 95\,\%$) until phase separation is observed, so that the solution is saturated. The weight fractions of each species are calculated from the actual ethanol-to-water weight ratio 
$w_{e,A}/w_{w,A}$ in the solution. The oil weight fraction 
$w_{o,A}$ of solution A is increased by increasing its ethanol-to-water weight ratio 
$w_{e,A}/w_{w,A}$, since the solubility of oil (decane) in ethanol is higher than in water. Finally, solution A is labelled yellow by adding a small amount of perylene (Sigma-Aldrich, 
$\geq 99\,\%$) at 
$0.2\ \textrm {mg}\,\textrm {ml}^{-1}$.
Solution B, which is poor in oil, is made of decane saturated water and is labelled light blue by Rhodamine 6G (Sigma-Aldrich, 99 %) at 
$0.2\ \textrm {mg}\,\textrm {ml}^{-1}$. Its oil weight fraction is 
$w_{o,B}=5.2\times 10^{-8}$ at 
$25\,^{\circ }\textrm {C}$.
Solution A is first injected to fill the entire channel, then solution B is injected at a constant driving pressure to perform the solvent exchange: an oil oversaturation pulse is generated in the mixture of solutions A and B, which leads to oil droplet nucleation both in the bulk and on the hydrophobic surfaces in the channel. The contact angle of oil in water is 
$\theta =15\pm 3^{\circ }$ on the same treated silicon surface and 
$\theta =11\pm 2^{\circ }$ on the same treated glass (see supplementary material available at https://doi.org/10.1017/jfm.2020.1137 for more details). The flow rate 
$Q$ of solution B, measured by a flowmeter (ML120V21, Bronkhorst, Netherlands), is varied by changing the driving pressure. The Péclet number of the flow is calculated by 
$Pe_h=\bar {u}h/D$, where 
$\bar {u}=Q/(Wh)$ is the average velocity of solution B, and a typical diffusivity of water in ethanol 
$D=0.84\times 10^{-9}\ \textrm {m}^{2}\,\textrm {s}^{-1}$ is used.
After approximately 0.3 ml (more than 1000 times of the channel volume 
$\approx 0.25\ \mathrm {\mu }\textrm {l}$) of solution B is injected, the injection is stopped by closing the valve, and a three-dimensional (3-D) scan of the whole channel is recorded by confocal microscopy (Nikon A1, Nikon, Japan) from below.
3. Experimental results
Typical mid-plane snapshots of the upstream part of the channel with the porous region (chip no. 1, see table 1 for the geometrical parameters) are shown in figure 2(a). The channel without porous region (chip no. 6, see table 1) is shown in figure 2(b). Black signals pillars, light blue signals water and red signals oil (because ethanol must have been dissolved in and washed away by the excess amount of water). In general, more oil droplets (red) are found in the channel after the solvent exchange. Furthermore, oil droplets are observed before and inside the porous region (black), but not behind the porous region. However, for channels without the porous region, droplets are observed in the entire channel (see supplementary material for the snapshots of the full channel). A closer look at the porous region is shown in figure 2(c), where the oil droplets and the circular pillars are clearly shown. A typical 3-D scan of the same area is shown in figure 2(d). Note that, in this work, we only focus on the total oil volume 
$V$, not on the droplet morphology.
Figure 2. (a) Typical mid-plane snapshots of the upstream part of chip no. 1, see table 1 for the geometrical parameters. Red signals oil, light blue signals water and black the pillars. The flow rates are 
$1.8\ \mathrm {\mu }\textrm {l}\,\textrm {min}^{-1}$ and 
$10.8\ \mathrm {\mu }\textrm {l}\,\textrm {min}^{-1}$. The arrow indicates the direction of the flow. (b) Typical mid-plane snapshots of the upstream part of chip no. 6. The flow rates are 
$6.5\ \mathrm {\mu }\textrm {l}\,\textrm {min}^{-1}$ and 
$29.6\ \mathrm {\mu }\textrm {l}\,\textrm {min}^{-1}$, respectively. The scale bar indicates 0.5 mm. Oil droplets are found both before the porous region and inside the porous region. (c) A closer look at the porous region of chip no. 1, where the circular pillars are clearly seen. (d) Three-dimensional confocal image of the morphology of oil in the same area, resolution is 
$0.28\times 0.28 \times 0.85\ \mathrm {\mu }\textrm {m}^{3}\,\textrm {pixel}^{-1}$ in the 
$x$, 
$y$ and 
$z$ directions, respectively.
Solvent exchange is performed for all the 6 chips (see table 1 for the geometrical parameters) at different 
$w_{e,A}$ and flow rates 
$Q$. After the solvent exchange, 3-D scans of the whole channel, similar to that shown in figure 2(d) are performed. The total volume 
$V$ of these droplets is measured by counting the number of red pixels of the 3-D confocal image and then multiplied by the volume of one pixel, then they are plotted against 
$Pe_h$ in figure 4(a,b) in log–log scale. Results for channels with the porous region (chip nos. 1–5) are shown in figure 4(a), and results for channels without the porous region (chip no. 6) are shown in figure 4(b). It is found that, for all the chips, 
$V$ increases with the Péclet number as 
$\propto Pe_h^{\alpha }$, with 
$\alpha \approx 0.50$, see table 2 for details. Also, 
$V$ increases with 
$w_{o,A}$, i.e. the more oil is in solution A, the more oil is nucleated. In the next subsection we will develop a theoretical model to quantitatively account for these two observations.
Table 2. Control parameters of the solvent exchange, and the best fitted power laws: 
$V\sim Pe_h^{\alpha }$. Here, 
$w_{o,A}$ is the oil weight fraction of the good solvent, 
$w_{e,A}$ is the corresponding ethanol weight fraction and 
$\alpha$ is the fitted exponent of the power law. The 95 % confidence bounds of the fittings are also shown. The average exponent 
$\alpha$ is 
$0.50\pm 0.01$.
4. Theoretical model
Figure 3(a) shows the schematic of the initial oversaturation pulse (red shaded region), which is a mixture of the two solutions due to advection and diffusion. The initial oversaturation pulse takes a parabolic shape in laminar flow, which is the case here. It broadens because of (Taylor–Aris) diffusion of the three components: oil, ethanol and water. The concentration of the mixture changes continuously from that of the solution A to that of solution B (Ruschak & Miller Reference Ruschak and Miller1972). Figure 3(b) is the phase diagram of the oil–ethanol–water system, it shows the concentrations of the solutions and the mixture, with the ethanol weight fraction 
$w_{e}$ and the oil weight fraction 
$w_{o}$ being the 
$x$ and 
$y$ axes, respectively. Then the water weight fraction of the solution/mixture can be calculated by 
$w_{w}=1-w_{e}-w_{o}$. The blue curve shows the concentrations of (oil) saturated solutions, this is the so-called ‘binodal curve’. The curve is fitted from the measured data points (see appendix A for details). Then, the concentrations of solutions A and B are on the binodal curve, denoted by points A and B in figure 3(b). The concentrations of the mixture lie on a curve connecting points A and B, denoted by the red curve AB – this is the so-called ‘diffusion path’ (Ruschak & Miller Reference Ruschak and Miller1972). All the possible concentrations of the liquid – both of the solutions and the mixture – lie on this curve. The oil concentration 
$c_{o}$ in the mixture is higher than its saturated oil concentration 
$c_{o,s}$, thus the mixture is oversaturated with oil. The (absolute) oil oversaturation is then denoted as
\begin{equation} \Delta c_{o}\equiv c_{o}-c_{o,s}. \end{equation}
Figure 3. (a) Sketch of the oversaturation pulse (red region) during solvent exchange. Solution A is on the right-hand side, being displaced by solution B. Liquid parcel 
${\rm M}^{\prime }$ is on the oversaturation pulse. The concentrations of solutions A, B and liquid parcel 
${\rm M}^{\prime }$ are represented by A, B and M in (b), respectively. (b) Phase diagram of the decane–ethanol–water ternary system. The blue curve is the binodal curve (saturation curve of oil) fitted from the data points (blue diamond) measured as described in Tan et al. (Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016). Red line AB is the diffusion path, black line MS is the dilution curve, which is a straight line passing through coordinate 
$(0,1)$. (c) Oil oversaturation normalized by the oil density, 
$\Delta c_{o}/\rho _{o}$, as a function of the normalized ethanol concentration 
$w_{e}=c_{e}/\rho$, at different initial conditions: black solid line for 
$w_{o,A}=0.063$, and red dashed line for 
$w_{o,A}=0.017$. (d) Plug shaped oversaturation pulse because of Taylor–Aris dispersion which happens in a long and thin channel; 
$x_1=x-\bar {u}t$ is the distance to the centre of the plug, 
$\bar {u}$ is the average flow velocity and 
$L_{p}$ is the length of the plug. The porous region is represented by the vertical black pillars. (e) The normalized oil oversaturation 
$\Delta c_{o}/\rho _{o}$ as a function of normalized length 
$\kappa x_1$, for 
$w_{o,A}=0.063$ (black solid line), and 
$w_{o,A}=0.017$ (red dashed line), respectively. The normalized length 
$\kappa x_1$ is a function of time 
$t$ and space 
$x$. (f) The calculated prefactor 
$s$ as a function of 
$w_{e,A}$. The two different initial conditions with 
$w_{o,A}=0.063$ and 0.017 are indicated by the red dots.
We now use mass-per-volume concentrations 
$c$ in the equations for simplicity but keep using weight fractions 
$w$ elsewhere. Notice that the weight fractions 
$w$ are the concentrations 
$c$ normalized by the density of the liquid 
$\rho$: 
$w=c/\rho$.
For any oversaturated liquid parcel 
${\rm M}^{\prime }$ in the mixture, as represented by point M on the diffusion path, when the oversaturated oil nucleates, it is considered that only oil nucleates out of the mixture, and the ethanol-to-water ratio in the mixture is kept constant (Lu et al. Reference Lu, Peng and Zhang2016). The concentration of the mixture moves along the so-called ‘dilution curve’ (Lu et al. Reference Lu, Peng and Zhang2016) to point S on the binodal curve. The dilution curve MS is a straight line passing through point 
$(0,1)$ in figure 3(b) (this point means pure oil in the phase diagram, see supplementary material and appendix B for more details), as shown by the black solid line. Therefore, the (absolute) oil oversaturation of this liquid parcel 
${\rm M}^{\prime }$ is
\begin{equation} \Delta c_{o}({\rm M})=(c_{o}-c_{o,s})|_{\rm M}=c_{o}({\rm M})-c_{o}({\rm S}). \end{equation} The exact shape of the diffusion path is determined by the diffusion speeds of the three components. In a simplified case, we assume that the diffusion of each component only depends on its own concentration gradient. To get a first-order approximation of the problem, we further assume that the diffusion coefficients of oil and water are equal: 
$D_{o}=D_{w}\equiv D$. Then, the diffusion path AB becomes a straight line (see Ruschak & Miller (Reference Ruschak and Miller1972) and appendix B. These two assumptions are reasonable since the oil concentration is small in all the tested cases: 
$w_{o,A}\leq 0.063$, and the errors induced by them are small, see appendix C), and 
$c_{o}({\rm M})$ becomes linearly dependent on the ethanol concentration 
$c_{e}$ of the arbitrary liquid parcel 
${\rm M}^{\prime }$. Since point M is on the diffusion path AB, we have 
$0\leq c_{e}\leq c_{e,A}$, where the initial condition 
$c_{e,A}$ is the ethanol concentration of solution A; 
$c_{o}({\rm S})$ can also be expressed as a function of 
$c_{e}({\rm M})$ by finding the intersection of the dilution curve (passing through point M) and the binodal curve. Then the (absolute) oil oversaturation 
$\Delta c_{o}$ of any liquid parcel in the channel can be expressed as a function of its ethanol concentration 
$c_{e}$
\begin{equation} \Delta c_{o}=\Delta c_{o}(c_{e}), \quad 0\leq c_{e}\leq c_{e,A}. \end{equation} This function of course also depends on the initial condition 
$c_{e,A}$. Equation (4.3) is calculated numerically (see supplementary material for more details), and the oil oversaturation normalized by the oil density, 
$\Delta c_{o}/\rho _{o}$, is plotted as a function of the normalized ethanol concentration 
$w_{e}=c_{e}/\rho$ in figure 3(c), for different initial conditions 
$w_{o,A}=c_{o,A}/\rho$: 0.063 and 0.017, respectively. Here, 
$\rho$ is the density of the liquid parcel (see supplementary material for the calculation of 
$\rho$), and 
$\rho _{o}=730\ \textrm {kg}\,\textrm {m}^{-3}$ is the density of the oil (decane).
The oversaturation pulse shown in figure 3(a) evolves as a function of space 
$\boldsymbol {x}$ and time 
$t$, so that both the oil oversaturation 
$\Delta c_{o}$ and 
$c_{e}$ also depend on space and time: 
$\Delta c_{o}[c_{e}(\boldsymbol {x},t,c_{e,A})]$. As mentioned previously, the total oil volume 
$V$ is nucleated from all the oil oversaturation generated in the liquid. The total oil volume 
$V$ in the porous medium can then be calculated by integrating the oil oversaturation over the entire channel volume, at the time just before the oversaturation pulse reaches the porous region
\begin{equation} V(t,c_{e,A})= \frac{1}{\rho_{o}}\int_{\varOmega} \Delta c_{o}[c_{e}(\boldsymbol{x},t,c_{e,A})] \,\mathrm{d}\boldsymbol{x}. \end{equation} The spatial and temporal distribution of ethanol concentration 
$c_{e}(\boldsymbol {x},t,c_{e,A})$ of the parabolic shaped oversaturation pulse is non-trivial. However, the channel used here is long and thin. The time scale for the axial advection is 
$\tau _{a} \sim L/\bar {u} \geq 4\ \textrm {s}$, which is much larger than the time scale of the vertical diffusion 
$\tau _{v} \sim (h/2)^{2}/D=0.12\ \textrm {s}$. In this situation, the concentration gradient in the vertical direction is smoothed out because of Taylor–Aris dispersion (Taylor Reference Taylor1953; Aris Reference Aris1956, Reference Aris1959), and the oversaturation pulse becomes like a plug, as shown in figure 3(d).
The analytical solution of the ethanol concentration 
$c_{e}(\boldsymbol {x}, t,c_{e,A})$ for pure water entering a thin square channel which at time 
$t=0$ contains only an aqueous ethanol solution with (uniform) concentration 
$c_{e,A}$ (no oil present) in this plugged regime is (Taylor Reference Taylor1953)
\begin{equation} c_{e}(\boldsymbol{x},t,c_{e,A})=\tfrac{1}{2}[1+\text{erf}(\kappa x_1)]\cdot c_{e,A}, \end{equation}
where 
$x_1=x-\bar {u}t$ is the distance to the central position of the plug, as shown in figure 2(d). Here, 
$\kappa = ({1}/{2}) k^{-{1}/{2}}t^{-{1}/{2}}$ is the wavenumber (with units 
$\textrm {m}^{-1}$), where 
$k=h^{2}\bar {u}^{2}/(210D)$ for a thin square channel of thickness 
$h$ (Dorfman & Brenner Reference Dorfman and Brenner2001). The wavenumber 
$\kappa$ has the temporal dependence and its inverse scales as the length of the plug: 
$L_{p}\sim 1/\kappa$.
For the solvent exchange in our case, oil is present in all of the liquid, but (4.5) still holds here because of the previous assumption 
$D_{o}=D_{w}\equiv D$. With this assumption, any portion of the water can be replaced by oil without influencing the above equation. Substituting (4.5) into (4.3), we have
\begin{equation} \Delta c_{o}=\Delta c_{o}(\kappa x_1, c_{e,A}). \end{equation}
Equation (4.6) is solved numerically, and figure 3(e) shows how the normalized oil oversaturation 
$\Delta c_{o}/\rho _{o}$ changes as a function of 
$\kappa x_1$ at the three different initial conditions. It is worth noting that the oil oversaturation is not symmetric, but its peak position is to the right (
$\kappa x_1\approx 0.55$), or in other words, more to the downstream. This asymmetry originates from the asymmetry in 
$({\Delta c_{o}}/{\rho _{o}})(w_{e})$ shown in figure 3(c), whose asymmetry originates from the shape of the binodal curve shown in figure 3(b). Also, the oil oversaturation is higher everywhere in the channel for the case with higher initial oil concentration (larger 
$w_{o,A}$).
Since the oversaturation pulse becomes a plug-like shape, which means the vertical concentration gradient is negligible, (4.4) can be further simplified as the integration along the single dimension 
$x$. Substituting (4.6), we have
\begin{equation} V(t,c_{e,A})= \frac{1}{\rho_{o}}Wh\int_{-\infty}^{\infty} \Delta c_{o}(\kappa x_1, c_{e,A}) \,\text{d}x. \end{equation} Substituting 
$t=L/\bar {u}$ into 
$\kappa$ as the time when the oversaturation pulse is just entering the porous region then integrating with respect to 
$\kappa x_1$ leads to
\begin{equation} V(c_{e,A})=s(c_{e,A})\cdot W h^{{3}/{2}} L^{{1}/{2}}\cdot Pe_h^{{1}/{2}}, \end{equation}where
\begin{equation} s(c_{e,A})=\frac{2}{\sqrt{210}}\frac{1}{\rho_{o}}\int_{-\infty}^{\infty} \Delta c_{o}(\kappa x_1, c_{e,A})\,\text{d}(\kappa x_1) \end{equation}
is a dimensionless prefactor which is proportional to the area under the curve 
$\Delta c_{o}(\kappa x_1, c_{e,A})$, as shown in figure 3(e).
From (4.8) we see that the power law dependence of 
$V$ on 
$Pe_h$ is indeed predicted to be 
$1/2$. The prefactor 
$s$ is shown in figure 3(f) as a function of normalized initial condition 
$w_{e,A}=c_{e,A}/\rho$. The two red dots correspond to the initial conditions 
$w_{o,A}=0.063$ and 0.017, respectively. It is worth noting that 
$s$ changes sharply within this range.
Equation (4.8) with (4.9) is the main theoretical result of our paper. It is universal in the sense that it includes all the factors that could influence the total amount of nucleated oil: the flow rate, the channel geometry and the solution composition.
Note that, in the model, only Taylor dispersion in the post-free region is considered, because the porous structure immediately interrupts further development of the oversaturation pulse. Further dispersion/mixing in the porous region is negligible: the length of the oversaturation pulse, which covers 99 % of the oil oversaturation (see figure 3e), is calculated to be 
$L_{p}=3.56/\kappa =0.49\sqrt {Lh}\cdot \sqrt {Pe_h}>0.44\ \textrm {mm}$, much larger than the pore size 
$\approx 20\ \mathrm {\mu }\textrm {m}$. That means that 
$c_{e}$ does not really change on the length scale of the pore size, thus the porous medium does not induce further dispersion/mixing. On the other hand, although in the theoretical model the length of the post-free region 
$L$ is used to calculate the total oil volume, this does not mean that the oversaturation would accumulate in the pulse and then all of a sudden nucleate into droplets at some time point, for example, at 
$t=L/\bar {u}$ when the centre of the oversaturation pulse is at the entrance of the porous region. Instead, as the oversaturation pulse develops while moving downstream, the (newly generated) oversaturation nucleates into droplets along the way, so that some of the droplets are observed upstream of the porous region.
5. Comparison between universal theoretical result and experiments
To compare the universal theoretical result with the experimental results, the oil volumes 
$V$ measured for various chips and mixtures shown in figure 4(a,b) are normalized by 
$s\cdot Wh^{3/2}L^{1/2}$, and then plotted against 
$Pe_h$ in figure 4(c,d) in log–log scale. Results for channels with the porous region (chip nos. 1–5) are shown in figure 4(c), and results for channels without the porous region (chip no. 6) are shown in figure 4(d). It is found that, for chips with the porous region, the calculated prefactors (
$s=1.2\times 10^{-2}$ and 
$3\times 10^{-3}$) can indeed collapse all data onto one master curve, regardless of the channel geometry and composition of the mixture. This is in agreement with the calculated prefactor of our theoretical model. However, for channels without the porous region, the measured oil volume 
$V$ is smaller than the theoretical prediction – which is not true for channels with the porous region. This confirms that the porous region can indeed prevent the oversaturated oil from leaving the channel, by which the comparison between experiment and theory is enabled.
Figure 4. (a,b) Measured total volume of oil in the entire channel 
$V$ as a function of the Péclet number 
$Pe_h=h\bar {u}/D$, for different initial conditions – oil weight fraction of solution A represented by 
$w_{o,A}$, and different flow rates as represent by 
$Pe_h$, in different chips. Results for the channels with the porous region (chip nos. 1–5) are shown in (a), and results for the channel without the porous region (chip no. 6) are shown in (b). For all the results, the total oil volume increases as 
$V\sim Pe_h^{\alpha }$, with 
$\alpha \approx 0.50$. Dashed lines are the best fit for each group. See table 2 for the solution compositions of the solvent exchange and the best fitted power laws. (c,d) Value of 
$V$ non-dimensionalized by 
$s\cdot Wh^{3/2}L^{1/2}$ and plotted as a function of 
$Pe_h$. Results for the channels with the porous region (chip nos. 1–5) are shown in (c), and results for the channel without the porous region (chip no. 6) in (d). The black solid lines are calculated from (4.8). Each data point is the average of five measurements and error bars are the standard deviation.
6. Conclusions and outlook
We have performed solvent exchanges in thin square channels with and without a porous region in the middle, at different initial conditions 
$w_{e,A}$, flow rates 
$Q$ and channel geometries. The total volume of these oil droplets 
$V$ is measured by confocal microscopy, and is found to increase with Péclet number as 
$\propto Pe_h^{\alpha }$, with 
$\alpha \approx 0.50$. A theoretical model is then developed, based on the ternary phase diagram and Taylor–Aris dispersion, to predict the total amount of oil 
$V$ nucleated from solvent exchange. The theory indeed predicts a power law dependence 
$V\propto Pe_h^{1/2}$. In addition, the influence of the channel geometry and the initial condition 
$w_{o,A}$ can also be calculated and included in the model to give a complete prediction of the oil volume 
$V$ on one universal curve, thanks to the porous region which can prevent the oversaturated oil from leaving the channel. This model is found to be able to predict the total nucleated oil volume 
$V$, irrespective of the channel geometry and initial mixture.
The findings of this work contribute to a better understanding of the solvent exchange, and could guide further design and research on this topic. First, a porous medium may serve as a good tool to collect all the oversaturation, providing a macroscopic approach to study this bottom-up method.
The results of this paper encourage us to propose a three-step approach to study solvent exchange. (i) Find the concentration distribution by solving the advection–diffusion equations. The boundary conditions are defined by the flow geometries and the initial conditions are set by the initial solution concentrations. (ii) Based on the concentration distribution, the oversaturation distribution can be calculated by applying the knowledge of the phase diagram. (iii) Investigate the quantity that is of interest, such as the volume of the nucleated phase, the dynamical interaction between the nucleated phase and the oversaturation, etc.
This research can be considered as a demonstration of the above proposed approach, with the total nucleated oil volume 
$V$ being the subject of interest, and the aid of Taylor–Aris dispersion in a long and thin channel to obtain the analytical solution of the concentration distribution. For more general cases where the oversaturation pulse is not a plug, analytical solutions of the concentration distribution may be difficult to get, then numerical simulations should be employed to finish the first step. It is also easier to incorporate the ternary (or multicomponent) diffusion effect in numerical simulations. Moreover, the dynamical interaction between the nucleation and the oversaturation pulse can also be incorporated in numerical simulations.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2020.1137.
Acknowledgements
We thank X. Zhang and C. Sun for valuable discussions, and H. Le The for the scanning electron microscope (SEM) Images. We thank the Micro-Nano Fabrication Laboratory of Peking University for providing the chips.
Funding
We acknowledge support from the Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC), an NWO Gravitation programme funded by the Ministry of Education, Culture and Science of the government of the Netherlands, and D.L.'s ERC-Advanced Grant under project number 740479.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Fitted binodal from the data points
The 
$X$- and 
$Y$-coordinates in figure 3(b) are actually 
$w_{e}=c_{e}/\rho$ and 
$w_{o}=c_{o}/\rho$. The binodal 
$w_{o,s}(w_{e})$ is fitted by two piecewise cubic Hermite interpolating polynomials (PCHIPs) through the data points, which are measured by titration following the procedure as described by Tan et al. (Reference Tan, Diddens, Lv, Kuerten, Zhang and Lohse2016).
For the first part (first 10 points, lower half), the polynomial is
\begin{equation} w_{o,s}=w_{o,s}(w_{e})=\sum_{i=1}^{10}P_i(w_{e}), \end{equation}where
\begin{equation} P_i(w)=a_i(w-w_i)^{3}+b_i(w-w_i)^{2}+c_i(w-w_i)+d_i \end{equation}
is valid on the 9 intervals between the 10 data points. Here, 
$w_i$ is the ethanol weight fraction of the 
$i$th data point. For the second part (last 5 points, upper half), the polynomial is
\begin{equation} w_{o,s}=w_{o,s}(w_{e})=\sum_{i=11}^{15}P_i(w_{e}), \end{equation}
where 
$P_i(w)$ has the same form as (A2). Parameters are shown in table 3.
Table 3. Parameters of the PCHIP spline fitting of the first 8 points of the binodal. Here, 
$i$ is the number of the data points, 
$w_i$ is the ethanol weight fraction of the 
$i$th data point and 
$a_i$, 
$b_i$, 
$c_i$ and 
$d_i$ are the coefficients.
Appendix B. Transformation from a ternary phase diagram to the phase diagram in a Cartesian coordinate
Figure 5(a) shows a typical ternary phase diagram, using the decane–ethanol–water system as an example. The three vertices E, W, D stand for ethanol, water and decane, respectively. Figure 5(c) shows the same phase diagram in Cartesian coordinates, with the ethanol weight fraction 
$w_{e}$ being the 
$X$-coordinate and the oil weight fraction 
$w_{o}$ being the 
$Y$-coordinate. Here, we prove that any straight lines in the ternary phase diagram shown in figure 5(a) are still straight lines in figure 5(c).
Figure 5. Transformation from the ternary phase diagram to a right-triangle phase diagram in Cartesian coordinates. (a) The ternary phase diagram of the decane–ethanol–water system. The three apexes are E for ethanol, W for water and D for decane. Blue curve is the binodal curve fitted from the data points (blue diamonds) taken from Skrzecz et al. (Reference Skrzecz, Shaw, Maczynski and Skrzecz1999). Red line AB is the diffusion path when 
$w_{e,A}=0.837$, and black line DM is the dilution curve. (b) Basis (
$\boldsymbol {p}$, 
$\boldsymbol {q}$) of the ternary phase diagram shown in red dashed arrows and basis (
$\boldsymbol {i}$, 
$\boldsymbol {j}$) of the Cartesian coordinates shown in black solid arrows. (c) Right-triangle phase diagram of the same system. Diffusion path AB and dilution curve DM are still straight lines. Here, 
$w_{e}$ is replaced by 
$x$ and 
$w_{o}$ is replaced by 
$y$.
To prove this, we need to find the transformation matrix between these two vector spaces. Let 
$\boldsymbol {p}$ and 
$\boldsymbol {q}$ be the basis vectors of the vector space that the ternary phase diagram is in; 
$\boldsymbol {p}$ is parallel to line WE in figure 5(a), and 
$\boldsymbol {q}$ is parallel to line WD. Let 
$\boldsymbol {i}$ and 
$\boldsymbol {j}$ be the basis vectors of the Cartesian coordinates. These two bases are put together in figure 5(b), with (
$\boldsymbol {p}$, 
$\boldsymbol {q}$) shown in red dashed arrows and (
$\boldsymbol {i}$, 
$\boldsymbol {j}$) shown in black solid arrows. Then we have
\begin{equation} (\boldsymbol{p},\boldsymbol{q}) = (\boldsymbol{i},\boldsymbol{j})\left[ \begin{array}{@{}cc@{}} \dfrac{1}{2} & \dfrac{\sqrt{3}}{2} \\ 1 & 0 \end{array} \right]. \end{equation}This is a linear transformation, with
\begin{equation} \boldsymbol{\mathsf{T}} = \left[ \begin{array}{@{}cc@{}} \dfrac{1}{2} & \dfrac{\sqrt{3}}{2} \\ 1 & 0 \end{array} \right] \end{equation}being the transformation matrix.
With the assumption of no ternary diffusion and 
$D_{e}=D_{o}\equiv D$, Ruschak & Miller (Reference Ruschak and Miller1972) first proved that the diffusion path in the ternary phase diagram is a straight line. By definition, the dilution curve is also a straight line. Because linear transformation transforms straight lines to straight lines, diffusion path AB and dilution curve DM in figure 5(c) are still straight lines.
On the dilution curve DM in figure 5(a), the ratio of ethanol weight fraction to water weight fraction is kept constant, that is
\begin{equation} w_{e}/w_{w}=c, \end{equation}
where 
$c$ is a positive constant.
In Cartesian coordinates, 
$w_{e}$ is the 
$X$-coordinate and 
$w_{o}$ is the 
$Y$-coordinate. For easier notation, let 
$x=w_{e}$ and 
$y=w_{o}$, then 
$w_{w}=1-x-y$, and
\begin{equation} \left. \begin{gathered} 0\leq x \leq 1,\\ 0\leq y \leq 1,\\ 0\leq 1-x-y \leq 1. \end{gathered} \right\} \end{equation}Equation (B3) can be transformed to
\begin{equation} y=1-(1+c)x. \end{equation}This is a straight line passing through 
$(0,1)$ in figure 5(c). Its slope satisfies 
$-\infty\leq -(1+c) \leq -1$. This is why the phase diagram in figure 5(c) is a right triangle.
Appendix C. Influence of non-equal diffusivities
Normally, the diffusivities of water and oil in the mixture are not equal. Here, we briefly discuss the influence of non-equal diffusivities of water and oil 
$D_{w}\neq D_{o}$ on the prefactor 
$s$.
With the assumption that ternary diffusion effects are neglected, and in the meantime we ignore the advection term in the advection–diffusion equations, the transport equations for water and oil become (Ruschak & Miller Reference Ruschak and Miller1972)
\begin{gather} \frac{\partial w_{w}}{\partial t}=D_{w}\frac{\partial^{2} w_{w}}{\partial x^{2}}, \end{gather}
\begin{gather} \frac{\partial w_{o}}{\partial t}=D_{o}\frac{\partial^{2} w_{o}}{\partial x^{2}}. \end{gather}Solutions of (C1) and (C2) have the form
\begin{gather} w_{w}=\frac{w_{w,A}+w_{w,B}}{2}+\frac{w_{w,A}-w_{w,B}}{2} \mathrm{erf} \frac{x}{2tD_{w}}, \end{gather}
\begin{gather} w_{o}=\frac{w_{o,A}+w_{o,B}}{2}+\frac{w_{o,A}-w_{o,B}}{2} \mathrm{erf} \frac{x}{2tD_{o}}, \end{gather}
where 
$x$ and 
$t$ are the distance (as defined in figure 3d) and time, respectively. Here, 
$w_{w,A}$ and 
$w_{w,B}$ are the water weight fractions of solutions A and B, 
$w_{o,A}$ and 
$w_{o,B}$ are the oil weight fractions of solutions A and B. With 
$-\infty < x < +\infty$ and 
$w_{e}=1-w_{w}-w_{o}$, the diffusion path can be computed.
Using the same diffusivity of water 
$D_{w}=0.84\times 10^{-9}\ \textrm {m}^{2}\,\textrm {s}^{-1}$ and taking the oil diffusivity 
$D_{o}=0.43\times 10^{-9}\ \textrm {m}^{2}\,\textrm {s}^{-1}$ as calculated by Perkins & Geankoplis (Reference Perkins and Geankoplis1969), we obtain the diffusion path for initial conditions 
$w_{o,A}=0.063$, 
$w_{w,A}=0.119$, as shown in figure 6. The diffusion path is no longer a straight line, but curved. The value of 
$s$ is calculated to be 
$s=1.2\times 10^{-2}$ and 
$3\times 10^{-3}$ for the two initial conditions 
$w_{o,A}=0.063$ and 0.017, all less than 4 % change as compared to the case when assuming 
$D_{w}=D_{o}$.
Figure 6. Diffusion path (red curve) for the case when 
$D_{w}=0.84\times 10^{-9}\ \textrm {m}^{2}\,\textrm {s}^{-1}$ and 
$D_{o}=0.43\times 10^{-9}\ \textrm {m}^{2}\,\textrm {s}^{-1}$, with initial conditions 
$w_{o,A}=0.063$, 
$w_{w,A}=0.119$. Notice now the diffusion path (red curve AB) changes to a curve rather than a straight line as shown in figure 3(b). The blue diamonds are the measured data on the binodal, and the blue curve is the binodal. The black line is the dilution curve.
Open access
