## 2 Materials and Methods
### Governing Equations  
For the simulations, we used a single-phase laminar fluid flow interface to model the diffusive/convective transport properties of the solutions in the mixer microchannel. All fluids were considered to be incompressible, and steady-state conditions were applied. The model solves the continuity equation:  
\[\rho\
abla\cdot\mathbf{u}\ =\ 0 \tag{1}\]  
Figure 1: Schematic of the stable operating conditions for our mix-and inject microfluidic devices [(40; 41)] based on optical experiments. Uniform mixing was observed to occur over a flow rate range of 80–220 \(\upmu\)L/min. The green area indicates the optimal conditions where both uniform mixing and stable jetting were observed. Figure 2: The mixing component of the microfluidic mix-and-inject device, (**a**) optical microscope image of the serpentine mixer fabricated using SU8 on glass; (**b**) schematic of the simulated geometry employed for the present study. The points at which the device characteristics (e.g., concentration) were probed are shown by the light blue squares in the cross-section of the mixer channel. created by the microchannel walls because of the serpentine geometry. The value of Pe also peaks near the outlet because of the sudden release of the flow due to the lack of microchannel wall resistance near the outlet. The simulated concentration profile is shown in Figure 3c; this agrees well with the corresponding experimental optical fluorescence data, an example of which is shown in Figure 3d, providing additional confirmation of the accuracy and interpretation of our FEM model. In both simulated and experimental concentration profiles, there are light and dark bands visible along the serpentine channel due to the effects of chaotic advection. However, we note that there are some variations between the simulated and experimental data which we attribute to several factors. These include the fact that out-of-plane, the optical concentration profile data are measured within a narrow, but finite, region within the channel defined by the depth-of-focus of the microscope, whereas the simulated concentration profile is precisely defined at the mid-height (with respect to the z-axis) of the microfluidic channel. Another point to consider is that any errors in the fabrication, though confined to the (sub-) micron scale, may have an effect and this will also impact the concentration profile. In order to further investigate the mixing characteristics, the slice concentration distributions and the cross-section iso-concentration contours for 12 probe positions at a flow rate of 80 \(\upmu\)L/min are shown in Figure 3e. The velocity shift towards the centre affects  
Figure 3: Distribution plots in the x–y plane at a flow rate of 80 \(\upmu\)L/min for (**a**) velocity, (**b**) Peclet number, and (**c**) concentration. The intensity of the Peclet number (Pe) distribution in (**b**) directly corresponds to the intensity of the chaotic advection effect (chaotic mixing is higher in darker areas). The chaotic mixing of the two solutions is also apparent from the concentration distribution shown in (**c**). (**d**) An experimental optical fluorescence image of mixing in the serpentine, also at 80 \(\upmu\)L/min, shows the same chaotic effect as was observed in the serpentine microchannel during simulations. Note that the outline of the channel is indicated by the dotted white line in (**c**,**d**). (**e**) shows the cross-sectional concentration distributions and iso-concentration contours as a function of distance travelled along the serpentine (at the 12 probe points which are indicated in Figure 2). the concentration profile by disrupting the interface between solution #1 and solution #2, which otherwise would be dominated by diffusion. Example concentration profiles for six probes along the serpentine micromixer at a flow rate of 80 \(\upmu\)L/min are plotted in Figure 4. The solutions being mixed can be distinguished by the sharp concentration gradient between them near the Y-junction.