The Effect of Junction Gutters for the Upscaling of Droplet Generation in a Microfluidic T-Junction

The influence of drop formation due to micro rib-like structures, viz., the Junction Gutters (JGs) within a standard microfluidic T-junction, is numerically investigated. Hydrodynamic conditions that lead to various flow regimes are identified characterized by the Capillary number (Ca) and velocity ratios of the dispersed and continuous phases (q) within a standard T-junction. Subsequently, under such conditions, a range of gutter configurations is introduced in the standard channel. The results predict that the introduction of JGs can favourably alter the formation frequency and morphology of drops and, consequently, promote upscaling significantly for the hydrodynamic conditions associated with low Ca. Detailed flow maps are presented that reveal a plethora of transitions during the formation of droplets with higher Ca and q that would otherwise signify a dripping or a jetting regime in a standard junction. However, specific gutter configurations are identified where JGs are unfavourable for generating monodisperse droplets.


Introduction
Droplet-based microfluidics is continuously evolving with applications associated with several aspects of science and engineering due to the reliable manipulation of drops (Zhao and Middelberg 2011;Stone et al. 2004;Link et al. 2004). These include chemical and medical applications wherein highly uniform droplet formation is a constant requirement, such as compartmentalization in biological assays (Scheler et al. 2019), medical diagnostics (Rivet et al. 2011), cellencapsulation (Köster et al. 2008), DNA-sequencing (Lan et al. 2017) and drug release for which microfluidic Labon-Chip devices are employed (Cui and Wang 2019). More recently, microfluidic devices have begun to evolve as excellent platforms for detecting RNA viruses (Basiri et al. 2020;Dolan et al. 2018). However, despite such substantial breakthroughs achieved in the society using microfluidic technologies, there are significant underlying challenges wherein a) droplet production rates using microfluidic devices owing to exorbitant handling costs and requirements on b) highest precision, control, and stringent quality standards during fabrication of microfluidic devices, bearing in mind their life-saving applications. To complicate matters further, it becomes essential to ensure the desired chemical and biological transformations during droplet fragmentation are intrinsically safe and environmentally friendly. Typically, in two-phase, liquid-liquid microfluidic systems within the scope of applications mentioned above, to precisely control, manipulate and enhance the droplet throughput, conditions such as the flow behaviour, droplet size, conditions of wettability and the geometry of the microfluidic device become critical (Sattari et al. 2020). Therefore, various microfluidic channels exist, such as but are not restricted to i) cross-flow, ii) co-flow, and c) flowfocusing and several geometric variations within them are possible. The T-junction is a form of a cross-flow configuration where the dispersed and the continuous phase fluids are fed orthogonally to generate droplets. Since its inception (Thorsen et al. 2001), the T-junction has significantly gained popularity owing to its simplicity and ability to produce monodisperse droplets with a coefficient of variation of < 2% (Xu et al. 2006). Furthermore, considering its capability to robustly upscale through minimal modifications through methods such as integrating several parallel-T junctions (Nisisako and Torri 2008), split the primary and secondary droplets (Sun et al. 2018). More recently, their potential to be configured within an Interactive Learning Control (ILC) (Huang et al. 2020) makes T-junctions favourable for mass production of droplets with high break up rates (Zhu and Wang 2017).
Nevertheless, the complexities involved with the dynamics of the fluids and their interactions with the T-junction configuration and a constant need for device miniaturization, scaling and upscaling continue to pose challenges associated with achieving the control and breakup of droplets (Chiu et al. 2017). Therefore, to overcome such challenges, several passive and active methods have been proposed where the former does not require external actuation. In contrast, the latter typically makes use of additional energy through electrical (Singh et al. 2020), thermal (Sohel Murshed et al. 2008), magnetic (Tan and Nguyen 2011), and mechanical actuation (Churski et al. 2010) with which the droplets are generated within the framework of a T-junction configuration. Although most active methods yield an excellent coefficient of variation and present a range of possibilities for drop generation, the challenges for parallelization, an additional level of complexity in handing the external input, and cost-based constraints may persist depending on the nature of the external input employed (Chong et al. 2016). In passive methods, the hydrodynamic conditions such as capillary number (Ca), the flow rate or velocity ratios of the dispersed and continuous phases, and at times, trivial geometric changes to the standard cross-flow or a co-flow configuration can be harnessed to generate droplets that emanate from flow regimes in a T-junction such as squeezing, dripping, and jetting (De Menech et al. 2008;Li et al. 2012).
Several research works have successfully modified the standard T-junction configuration both experimentally and numerically to enhance the process of drop formation through passive methods. Abate and Weitz (2011) experimentally proposed a modified T-junction that consisted of three junctions; a jetting junction, a bubbling junction where air bubbles were formed and a triggering junction in which the air bubbles deformed the jet to form droplets due to Rayleigh-Plateau instability. Shui et al. (2009) developed a head-on T-junction configuration in which two identical inlet channels, a constriction channel a wide outlet channel to examine the drop formation at different flow regimes. Their results suggested that the head-on devices have wider applicability to generate a broad range of droplet sizes in regimes driven by capillary instability, squeezing and shearing. Various numerical studies on the modified T-junctions have evolved in the recent past that investigated the head-on T junction in the form of an orthogonal double junction (Ngo et al. 2015;Han and Chen 2019;Raja et al. 2021) to generate alternating droplets and investigated the resulting drop sizes due to channel width and with different junction injection angles in standard T-junctions (Jamalabadi et al. 2017) and for double junctions ranging between 30-90 degrees (Ngo et al. 2016) that suggested the drop formation in an alternating pattern increases with injection angles. Consequently, the studies described above suggest that topological changes to the standard T-junction can be utilized to significantly promote droplet/ bubble breakup (Arias 2020), as detailed in the review of Cerderia et al. (2020).
Recently, a step like modification and a capillary device in the standard T-junction has shown to have substantial potential to produce monodisperse droplets under jetting regimes wherein polydisperse droplets are most often realized due to the uncontrollable Rayleigh instability (Li et al. 2015). However, with the step-like modification to the standard T-junction in place, the experiments of Li et al. (2015) demonstrated that monodisperse microdrops are accomplished owing to the narrowing jetting flow that can be realized. It was further demonstrated that monodisperse drops as small as 15 µm could be achieved under stable jetting conditions agreeing with the theoretical scaling. (Li et al. 2016). The most recent work of Zhang et al. (2022) that implemented the deep learning technique to enhance the ease of measuring the microdrops formed with the step arrangement in the T-junction further endorses the efficacy of such a step-like arrangement for forming monodisperse drops passively. The use of rib-like structures within the T-junction channel similar to the step arrangement mentioned above was numerically investigated comprehensively by Li et al. 2019. Four different rib structures viz., two triangular structures that are of the same width and height but with different orientations in reference to upstream and downstream of the channel, a streamlined structure formed by a semicircle complemented by a quarter of a circle, and a rectangular structure were embedded in the junction. The superiority of the rectangular rib was adequately demonstrated in their work which suggested that in the flow regime phase map, the jetting regime is curtailed by the rib, which resulted in favourably expanding the phase space of the dripping regime. New scaling laws for squeezing and dripping regimes were derived that indicated that the droplet diameter from the T-junction with the rectangular rib decreases linearly with the micro rib intrusion into the channel. The work of Li et al. (2019) opens further avenues for exploration of the rib-like structure closer to the junction, which acts as a droplet gutter by effectively promoting passive droplet generation in T-Junctions.
The present numerical work developed in this paper derives motivations from the works of Li et al. (2015Li et al. ( , 2016Li et al. ( , 2019. The current work furthers the investigation by exploring i) the behaviour of a wide range of rectangular/square junction gutters when embedded onto the standard T-junction subjected to various flow regimes, ii) transitions during drop formation that occur under the same hydrodynamic conditions within various flow regimes by purely varying the gutter topologies, and iii) the flow regime maps of the gutter phase-space topologies that promote and deteriorate droplet upscaling are identified.

System Description
The schematic of the microfluidic T-junction together with the Junction Gutter (JG) of length ( a ) and depth ( b ) is described in Fig. 1. A two-dimensional modelling approach is chosen for this study wherein the continuous phase inlet flows through the main channel and interacts with the dispersed phase injected through the side channel. The width of the continuous phase inlet and the main channel ( W c ), the width of the dispersed phase inlet ( W D ), the fluid properties of the continuous and dispersed phases are maintained as identical to that presented by Glawdel et al. (2012) since the current numerical work is validated against their experiments as demonstrated in System Description.
To ensure a fully developed laminar flow, the entrance lengths of the main channel ( L ec ) and the side channels ( L ed ) are chosen in accordance with Eq. 1. (Nekouei and Vanapalli 2017) that is used previously by Li et al. (2019): where W c corresponds to the hydraulic diameter, which is essentially the width of the continuous phase channel, and Re corresponds the Reynolds number. The Eq. 1. stated herein is for the continuous phase entrance length; however, it takes the same form with the corresponding parameters to obtain the dispersed phase entrance. The chosen entrance lengths correspond to the largest Re that is investigated in this study for both the continuous and dispersed phase channels, respectively, which is in line with the work of Li et al. 2019. The dimensions of the microchannel and the properties of the fluids used in the current work is presented in Table 1.

Numerical Description and Boundary Conditions
A range of numerical approaches such as the level set method (LSM) (Bashir et al. 2011;Wong et al. 2017), the volume of fluid (VOF) (Kashid et al. 2010;Ngo et al. 2015;Mastiani et al. 2017;Zhang et al. 2015), coupled-level set-and volume of fluid (CLSVOF) (Chakraborty et al. 2019) lattice Boltzmann (LBM) (Liu and Zhang 2009;Chen and Deng 2017) have been successfully employed in describing various regimes of drop formation in a microfluidic T-junction within the 2D framework. An overview of the numerical techniques, together (1) L ec = W c 0.6 1 + 0.035Re + 0.056Re , with their advantages and disadvantages, are comprehensively described in the work of van Sit Annaland et al. (2005). The VOF is a free surface reconstruction method that offers a simpler but robust treatment of the topological changes of the interface (Viswanathan 2019(Viswanathan , 2020 and can be applied for effectively describing droplet breakup and coalescence. However, the CLSVOF is a hybrid approach that couples the level set function to the VOF to estimate the curvature of the interface more adequately. In the present work, both the VOF and the CLSVOF are assessed within a standard T-junction (presented in Appendix A), and their suitability is examined for their applicability into T-junctions embedded with JGs. The equations that govern the flow in the system are described as follows: The continuity equation is given as and the momentum equation is described by where � ⃗ v is velocity, is the density, p is the pressure, is the viscosity, and t is the time. �� ⃗ F is the continuum surface tension force on the interface (Brackbill et al. 1992) of the volume fraction field , given by where k is the local curvature on the interface and is computed as In terms of describing the interface between two immiscible fluids, namely, the continuous and dispersed phases, and providing the volume fraction conservation throughout The volume fraction gives the portion of the cell which is filled with either phase, where The density and viscosity can be expressed as linear contributions from the continuous and dispersed phases indicated by the subscripts d and c as follows: For the CLSVOF method, firstly the level set function is defined ∅ is defined by The liquid phase properties that are interpolated across the interface are given by where the Heaviside function H(∅) can be written as follows (Sussman et al. 1994) where 2 is the finite interface thickness over which the fluid properties are smoothed. The value of is typically chosen between one to four times the length of the smallest computational cell so that numerical instability owing to parasitic currents are prevented. The surface tension force in the case of CLSVOF is given by, where the interface curvature is determined by the cell is f illed with the dispersed phase f luid In both cases, the wall adhesion is taken into account through a contact angle at the channel wall given by where n is the surface normal and n w , t w are, unit vectors normal and tangential to the wall, respectively. In the present simulations, it is assumed that the continuous phase perfectly wets the wall of the channel walls. In all the simulations in this study, the value of is fixed to be at 140º. As shown in Fig. 1, a uniform velocity is provided at the continuous and dispersed phase inlets, and at the channel boundaries C b1 and C b2 , a no-slip condition is employed. A zeropressure condition is prescribed at the outlet boundary. The dimensionless parameters that characterize this system are the Capillary number Ca = c U c , the channel Reynolds num- where U d and U c correspond to the maximum velocity of the dispersed and continuous phases. Considering that the maximum Reynolds number Re ≪ 1 , and Ψ , are fixed as given in Table 1, therefore, the two governing parameters are Ca and q for a standard microchannel T-junction. However, the presence of JGs in this system introduces two additional dimensionless parameters viz., the dimensionless gutter length a * = a W c and depth b * = b W c respectively.

Numerical Solution Procedure
The governing equations described in Numerical Solution Procedure are solved using the commercially available finite-volume based commercial software, ANSYS Fluent, Version 2020 R2. For the VOF method, the pressurevelocity scheme used was PISO that splits the solution into predictor and corrector steps together with the Non-Iterative Algorithm (NITA). However, in the case of the CLSVOF method, the PISO scheme was used in conjunction with the Iterative Algorithm (ITA) that required at least 45 inner iterations to ensure that all residuals of the CLSVOF method meet the sufficient convergence tolerance as demonstrated by the VOF, i.e., < 10 -6 . In both cases, the momentum equation was discretized using the QUICK scheme, and the gradients of the scalars were computed by using the Least-Squares cell-based method. The Least-Squares cell-based was chosen since it is directly comparable to node-based gradient methods, is much more accurate than cell-based methods, and is computationally less (16) n =n w cos +t w sin intensive. The "PRESTO!" (PREssure STaggering Option) scheme, although computationally more expensive, was used to interpolate the pressure term as it directly calculates the pressures at cell faces and avoids interpolation errors. The interface was determined by the Geo-Reconstruct Algorithm (Youngs 1982) that uses a piecewise-linear approach to determine the interface between fluids. To calculate the interfacial forces, the Continnum Surface Force (CSF) model was used (Brackbill et al. 1992). In the case of CLSVOF, the level set function described by Eq. 10 was discretized using the QUICK scheme. In both cases, the time-step is chosen to ensure that the Courant number is lesser than 0.25. For more information on the VOF solution process, the readers are directed to references (Viswanathan 2019(Viswanathan , 2020.

Grid Verification, Choice of Methods and Validation
The selection of a grid and the choice of methods was a multi-fold process. Firstly, the grid verification test is conducted by evaluating meshes of sizes 5 µm (coarse), 4 µm (medium), and 3 µm (fine), as shown in Fig. 2.
The spatially evolved droplet profiles superimposed with the grid details, just before pinch-off, are shown for the three mesh cases in Fig. 2a). The curvature associated with the evolution of the neck for the 3 µm (fine) is more prominent compared to the other cases. The value associated with t * denotes the time just before the droplet is formed with 3 µm (fine) case and it is observed that time for formation of the droplet relatively increases with the increased mesh  Table 1 for details of L ec ) for different grids evaluated against the theoretical description of the fully developed flow size. Further, the theoretical velocity profile u c (y) for a fully developed laminar flow based on entrance length estimate from Eq. 1 is given by.
A comparison of the numerically predicted velocity profiles against the theory (Eq. 17) for the three mesh cases at the same location, is shown in Fig. 2b). The differences between the fine and medium grids appear to be much lesser compared to the coarse grid, which is consistent with the observation in Fig. 2a). Although a good agreement is seen between the theory and numerical predictions for all the grids, the fine grid appears to show the least difference with the overall theoretical profile, and therefore, the 3 µm (fine) grid was employed for all the cases examined in this study. In the present work, a comparison between two interface capturing methods viz., the VOF and CLSVOF is presented in Appendix A to justify the choice of methods. As shown in Table 2, despite a close numerical prediction, the CLSVOF method, together with the iterative time advancement (ITA) scheme, requires ~ 3.6 times higher time to compute than compared to its VOF counterpart using the non-iterative time advancement (NITA) scheme. Therefore, the VOF method, together with the 3 µm (fine) grid alongside the numerical procedure described in Numerical Solution Procedure, was employed throughout the rest of the analysis.
The prediction from the VOF model is validated by comparing against the experimental data of Glawdel et al. (2012). Both the spatial and temporal periods of the drop formation and pinch-off such as the i) lag stage where the interface recedes to a small distance back into the dispersed phase inlet before it enters into the main channel, ii) filling stage wherein the interface penetrates and fills the main channel by proceeding towards the channel boundary C b1 , iii) transitioning into necking and iv) detachment are directly compared against the experimental data as shown in Fig. 3a).
The results show that the numerical predictions agree well with the experiments at various stages of drop formation and the overall droplet dimensions; an additional validation reinforces this with a different dispersed phase inlet dimension and hydrodynamic conditions presented in Appendix B. Nevertheless, the necking behaviour observed in the experimental data appears more pronounced than the numerical results suggesting that the stage between necking and pinchoff is quicker with the model. The apparent differences between model and experimental data during the necking stage before pinch-off due to the enhanced influence of wall boundary (C b1 ) could be attributed to a) the presence of dynamic wetting conditions in the experiments, whereas the simulations employ a static contact angle and b) the formation of thin liquid film between the droplet and wall in the experiment that is not resolved in the current simulations. However, the drop formation frequency predicted by the model shows a difference of ~ 1% w.r.t to the experimental data of Glawdel et al. (2012), as shown in Table 2, suggesting a good agreement between the model and experiment. Figure 3b) shows the intricate features of the overall pressure and continuous phase velocity fields during the formation of the droplet. Features such as the a) drop emergence into the channel at 5 ms shows an increase in the velocity field due to distortion of the fully developed continuous fluid, b) during the filling stage at 15 ms, the region between the boundary C b1 and the interface experiences significantly higher velocity, c) diminishing of the velocity vectors at 30 ms just prior to breakup, a significant volume of fluid fills the channel, and finally d) the high flow is directed into a region where the droplet pinches off. These critical features are consistent with that observed by the experimental demonstrations of van Stejin et al. (2007) and the 3D numerical work of Soh et al. (2016). Both spatially and temporally, resolving the formation and migration of droplets in microchannel T-junctions using the 3D numerical simulations is  (2012)) ---29.7 -highly challenging and computationally expensive. Despite being computationally prohibitive, well-resolved 3D simulations have shown the ability to accurately resolve the lubricating film's dynamics that signify the droplet-wall interactions (Ling et al. (2016)) for a wide range of Ca values. Unlike the 3D direct numerical simulations, the 2D numerical framework adopted in this study is inherently limited to predicting the leakage and corner flow characteristics revealed in the experiments (Korczyk et al. (2019)). Nevertheless, the comparisons illustrate that the 2D approach adopted in this work is able to replicate the essential features of drop formation, viz., the drop length and the formation frequency of  Figure b) shows the pressure distribution in the channel superimposed with the velocity vectors of the continuous phase the droplet, which are crucial parameters of interest for the present study. In the following sections, firstly, the hydrodynamic conditions that lead to various flow regimes within the microfluidic T-junction are identified. As a next step, the effect of upscaling is analysed in the context of the droplet morphologies, formation frequencies, transitions, phase maps, and scaling laws by embedding a wide range of JGs within the microfluidic T-junction that is subjected to the hydrodynamic conditions corresponding to the identified flow regimes.

Microfluidic-T Junction
As described in Numerical Description and Boundary Conditions, the variation of two governing parameters viz., the Ca and q leads to several flow regimes in the present system. Earlier studies have shown that the flow regimes that are present within a microfluidic T-junction can be classified as squeezing, dripping, jetting, and parallel flow where no droplets are formed (Li et al. 2019(Li et al. , 2012Santos and Kwaji 2010). In addition, the transition between each of these regimes exists that can potentially result in a wide range of flow regimes. Furthermore, the experimental results of Santos and Kwaji (2010) illustrated a snapping regime at very low Ca and with low dispersed phase velocities. However, based on the competition between the forces that are involved within the flow regimes, viz., the surface tension force, shear force on the interface, and the hydrostatic pressure differences on the sides of the emerging droplet, the nature of slugs that are formed differ (Gastecki et al. 2006;Li et al. 2012;Christopher et al. 2008). Figure 4 shows the blocking (B s ), squeezing (S s ), dripping (D s ), jetting (J s ), and parallel flow (PF s ) regimes that are identified with the Ca and q values respectively. The subscript 's' denotes that the regimes correspond to the hydrodynamic conditions associated with the standard microfluidic T-junction. In each case, the transient pressure at P j is plotted since it provides adequate information on the inherent droplet breakup characteristics as adopted previously (Wong et al. 2017;Li et al. 2012). At lower capillary numbers ( Ca ≪ 0.032 ), two regimes viz., the B s and the S s are identified. In Fig. 4a) within the images show by a-e, the B s regime is identified with Ca = 0.0015 and q = 2.001 , the dispersed phase enters the main channel that encounters the junction point leading to a pressure increase as shown by the image a. Soon after blocking the channel shown by c, a sluglike feature evolves with a localized neck that appears closer to the junction, as shown by d, leading to pressure buildup at the junction from the point c. However, this slug-like feature continues to grow wherein the neck moves away from the junction and shows significant resistance to breaking by growing and therefore, the entire channel cross-section is fully blocked as shown in e. The B s regime described in this study inherits the some characteristics of the snapping regime experimentally described by Santos and Kwaji (2010), except that in the present study except that the pinch-off did not occur at least until t*(s) ~ 0.13 s. Experimentally, the work of Arias and Montlaur (2020) presented the B s like regime for small Ca values (as small as 0.6 × 10 -3 ) when analysing the bubble breakup in a microfluidic T-junction. In addition, the 3D numerical results of Li et al. (2014) predicted a longslug formation that exhibits the same characteristics as that identified in the present study. In their results, Li et al. (2014) observe that such long slugs break with higher wall adhesion forces due to longer wall contact times. The squeezing regime S s , as shown by the images a-f within the Fig. 4b), exhibits the features such as filling (show by c) and necking (shown by the image d) that are similar to the B s regime but eventually the breakup (shown by images e and f)) is evidenced through the higher pressure in the image e relative to that predicted in the image d. Although the S s regime undergoes same sequence in terms of emerging, filling and blocking the channel, it is quite well known that the pressure difference across the slug causes breakup at low Ca such as the squeezing regime (Gastecki et al. 2006;Li et al. 2012Li et al. , 2019Christopher et al. 2008;De Menech et al. 2008). However, interestingly, this characteristic increase in transient pressure during breakup witnessed in the S s regime is not evidenced with the B s regime suggesting that in the B s regime, the transient pressure that has evolved at the junction was not sufficient to induce a breakup at least until the simulated time t * (s)~0.13 s. The dripping regime (D s ), as shown in Fig. 4c) with Ca = 0.0322 and q = 0.1458 , is shear-dominated, where the breakup occurs when the shear forces overcome the interfacial force.
Unlike the B s and the S s regimes, the droplet does not block the main channel entirely in the case of the D s regime where the droplet pinch-off is confined to the lower boundary C b2 (see Fig. 1); therefore, the droplet breakup is accompanied by the continuous fluid that emerges through the gap between the droplet and the upper boundary C b1 as seen in the image d in Fig. 4c). With Ca = 0.049 and q = 0.3296 eventually leads to the jetting regime J s where the dispersed phase fluid shows a thread-like structure. In contrast to the B s and the D s regimes where the breakup occurs at the junction, the droplets break up in the J s regime occurs downstream of the channel as shown by the images d, g in Fig. 4d). Finally, the parallel flow regime PF s is realized with Ca = 0.0725 and q = 0.2777 where the dispersed phase eventually flows parallel to the channel boundary C b2 (shown in image e within Fig. 4e)), where no droplets are generated. In addition to the validation shown in Fig. 3 and that presented in Appendix B, the parameters , , q, and Ca established in this study to characterize the S s , D s , and the J s regimes, which are in very good agreement with conditions identified for squeezing, dripping, and jetting in a standard T-junction as shown in other previous works (De Menech et al. 2008;Li et al. 2019;Liu and Zhang 2011). In the context of a gutter that is positioned closer to the junction, as shown in Fig. 1, the J s is representative of the PF s considering the similarities between the flow behaviour closer to the standard junction as shown by the images e and f within Fig. 4d). Therefore, as a next step, a wide range of junction gutters are introduced into the channel that is subjected to the hydrodynamic conditions pertaining to B s , S s , D s , and the J s but excluding the PF s . Figure 5 shows the effect of drop formation due to junction gutters ranging between a * =0.05 to 1.00 and b * =0.10-0.70 when the microchannel T-junction is subjected to hydrodynamic conditions pertaining to the B s regime as described previously in Fig. 4a). The introduction of a gutter in the B s regime promotes the droplet breakup in the junction, In each case (a-e), the transient pressure at the junction probe (P j ) is shown alongside the sequence of events described by (a-f). In each case, t* = 0 s corresponds to the time that the dispersed phase volume fraction first emerges into the continuous phase channel unlike that witnessed in the standard T-junction. The slug size appears to be largest for the case with the lowest gutter depth b * =0.10 together with a greater gutter length a * =1.00 analyzed in this study. As shown in Fig. 5, although both a * and b * influence the length of the droplet L d , however, the gutter depth b * appears to predominantly influence the droplet size compared to the gutter length a * with the hydrodynamic conditions pertaining to the B s regime. In addition to the droplet size, the variation of gutter dimensions has a broader consequence in terms of the morphological characteristics pertaining to the breakup, as demonstrated in Figs. 6 and 7.

Influence of Junction Gutters in the Hydrodynamic Conditions Pertaining to the Bs and the Ss Regimes
The influence of droplet breakup for a fixed value of b * =0.10 but with varying a * =0.10-1.00 when subjected to the hydrodynamic conditions of B s is shown in Fig. 6. For each gutter configuration, the sub-figure a) shows the necking just prior to slug pinch-off whilst the figure b) shows the incipience of the slug breakup.
For a * <1.00, the necking process is at the junction, and the images (at the bottom) show that the interface recedes into the main channel after the slug is formed. Interestingly, the necking and breakup sequences described in all cases with a * <1.00, and b * =0.10 are reminiscent of the S s regime, although the channel is subjected to hydrodynamic conditions pertaining to the B s regime. However, with a * =1.00 and b * =0.10, the slug grows into the main channel with a thin liquid thread but eventually snaps at a distance from the junction, unlike in the other cases. Eventually, the thin liquid tail tends to retract into the T-Junction as shown by the dotted square box in the sub-figure b) a * =1.00 and b * =0.10. It can be observed that the process of (i) tail thinning away from the junction point and (ii) the breakup of the large slug closely resembles the snapping slug described by Santos and Kwaji (2010). The effect of variation in gutter depth ( b * ) for fixed gutter length ( a * ) is shown in Fig. 7. The post-breakup images shown by the cases with b * =0.50 and 0.70 suggest that an increase in the gutter depth significantly reduces the slug size, and the slug evolves eventually by squeezing in between the gutter and channel walls. It is interesting to note that, for every gutter topology besides the small gutter case ( a * =0.10 and b * =0.10) shown in both  to left), the necking, droplet detachment, and the location of the droplet in the channel extracted at t * =0.03125 s; where t * =0 is identical to that presented in Fig. 4b). Figure b) shows the transient pressure evolution at the junction probe (P j ) for two cases where b * =0.3 and a * =0.05 and 0.5. The droplet formation times are indicated by green circles (for a * =0.05) and blue squares (for a * =0. Figs. 6 and 7, the dispersed phase channel experiences lower pressure distribution during post-breakup than compared to the necking stage. Figure 8a) presents the variation of gutter length ( a * =0.05-1.00) with fixed gutter depth ( b * =0.30) on drop formation when the microfluidic T-junction is subjected to the hydrodynamic conditions pertaining to the S s Regime (see Fig. 4b)). In each case, the images show necking (right column), breakup (middle column) and location of the droplet in the channel (left column) extracted at the same instance of time ( t * =0.03215) for all cases to assess the speed of droplet traverse in the main channel. With small values of a * =0.05 and 0.10, the necking and pinch-off responses shown in Fig. 8a) no longer exhibit the characteristics of the S s regime, i.e., during necking or pinch-off, the droplet does not block the channel by adhering to the channel boundary (C b1 ) instead, characteristics like that of a dripping regime D s (see Fig. 4c)) is observed. However, for cases with a * ≥0.30, the necking and pinch-off characteristics of S s re-appear and are well preserved. To distinguish the nature of breakup shown by a * <0.30 against that observed with a * ≥0.30, the transient pressure evolution at the junction point (P j ) for two different cases is presented in Fig. 8b). For a * =0.05, the incipience of pinch-off indicated by blue squares shows high-pressure peaks, which is a characteristic of the D s regime observed in Fig. 4c) (shown by the images d) and e)). With a * =0.50, although the pressure at the junction increases during pinch-off (shown by green circles), the peak observed at the incipience of droplet pinch-off is markedly different to a * =0.05 but shows similarities to the pressure profile of the S s regime. This suggests that the choice of gutter configurations in the microchannel can influence transitions at least during the necking phase and drop detachment at low Ca . For all cases shown in Fig. 8a) (left column), the results predict that the speed of the evolved droplet in the main channel appears to increase with incremental values of a * .
The effect of gutter dimensions on the size of the droplet ) and the droplet formation frequency (f) is summarized in Fig. 9 when the channel is subjected to conditions corresponding to the B s (show by symbols with solid lines) and the S s (shown by symbols with dashed lines). For cases pertaining to the B s , except for b * =0.10, both L D * and f show no noticeable variation w.r.t to the gutter length a*. However, the variation of a * has a significant impact when b * is small, viz., with 0.10, which is supported by the images of drop evolution shown in Fig. 6). The behavior is similar with the conditions of S s where, for most cases, the change in a * has negligible influence on L D * and f. However, the presence of a transition shown for a * <0.30 to a * ≥0.30 suggests that the choice of a * can influence a transition in the system, as shown in Fig. 8. For all cases investigated, the predictions suggest that increasing the gutter depth b * for all values of a * significantly reduces the size together with increasing the frequency of the droplet formation when the channel is subjected to the B s as well as the S s conditions.  Fig. 4a) with Ca=0.0015, q=2.001, and (S s ) regime (See Fig. 4b)) with Ca=0.0043, q=0.6487 respectively Several researchers have proposed scaling laws for predicting the size of the slug ( L D * ) for over a range of Ca values within the (i) squeezing (Garstecki et al. (2006)) and (ii) transition from squeezing to dripping regimes (Xu et al. (2008)) for standard T-junctions. More recently, the work of Li et al. (2019) proposed a scaling relationship within the squeezing regime for T-junction microchannels with rectangular ribs, which suggested that the resulting agreement between the numerical data and derived scaling law could be as close as ~ 15%. However, their study was limited to rib depth of up to 50%. In the current work, the general equation proposed by Li et al. (2019) is adopted in Eq. (18), but as a function of velocity ratios ( q ) of the continuous and dispersed phase fluids, as follows: In Eq. 18, and are fitting parameters. To verify the predictions of the numerical model, the resulting slug size is compared against theory (Eq. 18) for a * =0.30 but tested for a (18) wide range of gutter depths ( b * ). The estimated using Eq. 18 for both B s and S s conditions resulted in Eqs. 19 and 20.
When the channel is subjected to the Bs conditions, the Ca =0.0015 (small) and therefore, the overall agreement with the theory (Eq. 18) is < ~ 4% for all cases except for b * =0.70 where the deviation is up to ~ 16% as shown in Fig. 10. With hydrodynamic conditions pertaining to S s , the resulting increases to Ca =0.0043, and consequently, the nonlinearities associated with the breakup process significantly increases, thereby resulting in deviations of ~ 16% with b * =0.50 and as high as 88% for b * =0.70. Nevertheless, the scaling law described in Eq. 18 agrees well for b * ≤ 0.50 and b * <0.5 for conditions pertaining to B s and S s , respectively. To further analyse the consequence of nonlinearities associated with the breakup process, the channel is subjected to D s and J s conditions with several junction gutters.  Figure 11a) presents a phase space over a wide range of gutter depth ( b * ) and gutter length ( a * ) that occur at Ca=0.0322, q =0.1458, which correspond to the conditions of the D s regime. Three distinct droplet morphologies in terms of adherence to channel walls of the T-junction are predicted for over a wide range of ( a * ,b * ) in the flow map as shown in Fig. 11b). The limits marked by the solid red and blue lines in Fig. 11a) present a clear transition between droplets' adherence to the channel boundaries. For most cases when b * ≤0.30, the droplet formation and evolution processes are much like the dripping regime wherein the drops adhere to both the upper (C b2 ) and lower (C b1 ) channels, indicated by circles in the map. However, with b * =0.30 and for a * >0.60, a transition appears where the evolved drops are much smaller and are unbounded to either of the channel boundaries indicated by open squares. When b * ≥0.5, another transition appears where the drops adhere to the boundary C b2 as shown by open triangles. This regime (shown by open triangles) diminishes with the increase in b * and for increasing values of a * where the unbounded drop regime starts to predominate in the map. Like the earlier conditions of B s and S s , for any given configuration of the gutter, the droplet spacing ( ), (see Fig. 1), predicted by the numerical model between any two droplets with the D s conditions,  Fig. 4(c) are identical. Consequently, all droplet shapes and drop formation frequencies predicted are identical for every subsequent droplet formed after the first drop, which reinforces that the gutters exhibit strong potential to produce monodisperse drops for conditions pertaining to D s . Figure 12a) and b) complement the information from the regime map  Fig. 13 a) flow regimes in ( a * , b * ) with droplets generated through three characteristic modes predicted viz., the uniform dripping, the non-uniform dripping that proceeds to parallel flow, and non-uniform dripping. The solid red lines indicate the boundary of the uniform dripping region. b) shows the transient pressure evolution at the junc-tion probe (P j ) with different gutter configurations that describe the three regimes and morphologies of the evolved droplets. In all cases, the hydrodynamic conditions are fixed with Ca=0.049, q=0.3296 which correspond to Js regime as shown in Fig. 4d) by showing the variation in L D * and f due to a * . Unlike the S s and the B s conditions, as shown in Fig. 9, although negligible, the size of drops shows some noticeable variations with change in a * due to a plethora of morphological transitions that are evidenced in the ( a * , b * ) phase map shown in Fig. 11a). With the increase in b * , that the drop formation frequency can significantly increase, as shown in Fig. 12b), thereby suggesting that the inclusion of gutters can substantially promote upscaling correspondingly in the D s regime.

Influence of Junction Gutters in the Hydrodynamic Conditions Pertaining to the Ds and the Js Regimes
When the microchannel T-junction is subjected to the J s regime with Ca=0.049, q=0.3296, the gutters alter the flow pattern more drastically, as shown by the ( a * , b * ) flow regime map in Fig. 12a). The attributes of the J s regime as shown in Fig. 4d) are no longer preserved; instead, with b * <0.40 and a * <0.4, a uniform dripping regime emerges. The uniform dripping regime inherits the features of the D s regime and is shown by filled circles in the phase map; its boundaries are defined by the solid red lines. One such formation of droplets in the uniform dripping regime is evidenced in Fig. 12b) with gutter dimensions a * =0.30 and b * =0.40; together, the transient pressure at the junction (P j ) shown by the solid green line suggests that the once the junction pressure stabilizes, the droplets continue to evolve with uniform droplet spacing ( ), the droplets are monodisperse, and the frequency of formation (f) of the droplets are uniform and stable (Fig. 13).
Crossing the threshold of the uniform dripping leads to the onset of non-uniform jetting (marked by filled squares in Fig. 12a). One instance of this regime is shown by a * =1.00, b * =0.60 in Fig. 12b); the corresponding transient pressure evolution at the junction is shown by the solid black line. This regime shows attributes of dripping at the initial stages, wherein some drops begin to emerge but with non-uniform spacing and size. However, eventually, the dispersed phase liquid continues to evolve with characteristics similar to a parallel flow regime (PF s ) as previously shown in Fig. 4e). For values of b * ≥0.70, the parallel flow characteristics are inhibited, but the non-uniform drips continue to evolve in the channel; this regime is shown by filled triangles in regime map in Fig. 12a) and is shown by a * =0.50 and b * =0.70 with the corresponding transient pressure at the junction shown by the solid red line in Fig. 12b). Unlike the uniform dripping regime, both the non-uniform jetting and the nonuniform dripping tend to deteriorate the monodispersity in drop formation significantly and tend to become unfavorable.

Conclusions
In this study, the effect of junction gutters on drop formation in a microchannel T-junction was numerically investigated. The numerical model was comprehensively assessed in the form of a grid verification test, evaluating the choice of interface capturing methods such as VOF and CLSVOF and identifying the rationale behind choosing the VOF method and validated with the experimental data of Glawdel et al. (2012) for a standard microchannel T-Junction. The hydrodynamic conditions leading to flow regimes such as the blocking (B s ), squeezing (S s ), dripping (D s ), jetting (J s ) and parallel flow (PF s ) that are characterized by the capillary number ( Ca ) and velocity ratio ( q ) were identified for the standard channel. An extensive range of Junction Gutters (JGs) was then embedded onto the junction of the standard microchannel, with gutter lengths a * , and depth b * , ranging from 0.05-1.00 and 0.10 and 0.80, respectively. With the introduction of the JGs under the same hydrodynamic conditions pertaining to regimes mentioned above, the following findings are summarized.
(i) In the hydrodynamic conditions of Bs and Ss where the Ca and flow velocities are small, gutter depth b * significantly influences the drop size and frequency of formation of droplets than the gutter length, a * . However, a * tends to promote transition by invoking changes to the morphology of the breakup of drops for small values of a * . In these regimes, the theoretical scaling law for predicting the size of the droplet with the gutters appears to strongly depend on b * and matches reasonably well with the numerical predictions for b * ≤0.50 for both the regimes.
(ii) With the presence of JGs in the Ds regime, the results suggest the presence of three distinct droplet morphologies of drops which are detailed by the ( a * , b * ) phase map in terms of their nature of adherence to the channel boundary of the T-junction. Unlike in the Bs and the Ss conditions, the size and formation frequency of the drops show noticeable variations with a * when the channel is subjected to Ds conditions. However, when the channel is subjected to Js conditions, flow transitions such as uniform dripping, nonuniform jetting, and non-uniform jetting occur that are presented through the ( a * , b * ) flow map, which details both the favourable and unfavourable topologies of gutters.
(iii) For the range of flow regimes identified through Ca and q within the standard T-junction, the JGs tend to influence the drop generation rates by promoting the upscaling favourably. Nevertheless, a careful selection of JGs in the hydrodynamic conditions of the Js regime is vital to foster monodisperse drop generation by transitioning into a uniform dripping regime.
Further numerical and experimental investigations are necessary to underpin (i) the optimal shape of JGs for drop formation, (ii) the influence of wall wettability between JGs and channel boundaries, (iii) critical conditions that can alter flow transition, and (iv) modified scaling laws with gutters that can predict the size of droplets during transitions.

Appendix A
To justify the choice of methods in capturing the interface of the two phases, a comparison between VOF and CLSVOF was undertaken as shown below in Fig. 14 with the fine mesh with the same hydrodynamic conditions and channel dimensions as in Fig. 2.
The spatial and temporal evolution of the droplets' interface for both the CLSVOF and VOF methods appear to be identical. However, at the incipience of the breakup at 30.25 ms, subtle differences exist between these methods on the interface curvature, as seen in Fig. 14c). The pressure difference P diff , which is the difference in pressure predicted by the CLSVOF (P CLSVOF ) and the VOF (P VOF ), suggests that the CLSVOF predicts a marginally higher pressure at the interface just after the breakup. Nevertheless, such differences are negligible considering that the two important parameters for the current study, viz., a) the final drop shape and b) the frequency of drop formation predicted by the two methods, show differences of ~ 0.17% and ~ 0.02%, respectively suggesting excellent agreements between CLSVOF and VOF methods (Table 2).

Appendix B
An additional comparison is presented in Fig. 15 for different operating conditions to fortify the validation of the VOF model against the experimental data of Glawdel et al. (2012) using a fine grid size of 3 µm.  b) shows the differences in the pressure distribution (P diff = P CLSVOF -P VOF ) in the channel between CLSVOF and VOF in the channel at the same instances In this case, q = 0.3514, parameters such as the width of the dispersed phase inlet ( W d ) is 45 µm, and Ca = 0.0027 were maintained the same as the previously published experimental result of Glawdel et al. (2012). The remaining parameters were maintained the same as provided in Table 1 of the manuscript. As shown in Fig. 15, the drop formation frequency predicted by the numerical result is 45.454 Hz, whereas the experimental drop formation frequency reported was 45.8 Hz resulting in a difference of ~ 0.75%. The numerical results agree well with the experimental measurements for various stages of the drop formation process.