Examination of Haines Jump in Microfluidic Experiments via Evolution Graphs and Interface Tracking

Abstract

This work examines a type of rapid pore-filling event in multiphase flow through permeable media that is better known as Haines Jump. While existing microfluidic experiments on Haines Jump mostly seek to maintain quasi-steady states through very low bulk flow rates over long periods of time, this work explores the combined use of a highly structured microscale transport network, high-speed fluorescent microscopy, displacement front segmentation algorithms, and a tracking algorithm to build evolution graphs that track displacement fronts as they evolve through high-speed video recording. The resulting evolution graph allows the segmentation of a high-speed recording in both space and time, potentially facilitating topology-cognitive computation on the transport network. Occurrences of Haines Jump are identified in the microfluidic displacement experiments and their significance in bulk flow rates is qualitatively analyzed. The bulk flow rate has little effect on the significance of Haines Jump during merging and splitting, but large bulk flow rates may obscure small bursts at the narrowest part of the throat.

1. Introduction

Multiphase flow through permeable media occurs in natural and artificial processes where fluid flows in void spaces that are interspersed in a solid matrix. Examples include blood flow through capillaries [1,2] and transport through membranes [3]. This work is motivated by geological applications such as hydrocarbon recovery from subsurface formations [4,5] and geological storage of greenhouse gases [6,7]. Void spaces in subsurface formations often have complicated geometries and are formally described as a transport network consisting of spacious intersections (pores) of narrow bridges (throats) [8]. Accurate characterization of fluid dynamics at the same spatial scale as pores and throats is required to understand multiphase flow through complex networks. With pores and throats at the micrometer scale, capillary forces dominate flow dynamics such that the dynamics are contingent on void space geometry and matrix/fluid wettability [9,10]. In immiscible displacement processes, the injected fluid may form a finger-like pattern as it displaces the resident fluid [11]. The exact mechanism by which microscale interfacial phenomena translate to observable macroscale phenomena is an active area of research, a prominent example being fingering instabilities [12,13]. The complexity of multiphase flow in permeable media is often attributed in the literature to two types of interface-related events:

• Haines Jump [14], originally conceived to address rapid pore-filling events during drainage of a wetting phase where it is rapidly displaced from a pore by a non-wetting phase from the most available throat, and
• snap-off events [12,15,16,17], where interfaces enter a throat from both sides then coalesce to ensnare the wetting phase (during drainage) or the non-wetting phase (during imbibition);

Haines Jump is regarded as a promoting factor for flow instabilities [18,19,20], and it is conjectured that Haines Jump velocity is independent of bulk flow rate, defined as volumetric flow rate across the entire medium [21]. This work addresses a gap in existing literature regarding the occurrence and significance of Haines Jump during high bulk flow rates; refer to

Table 1. Summary of existing literature reporting experiments similar to this work.

Article	Sample	Geometry	Fluids	Drainage	Imaging
Abdulla et al. [22]	oil-wet rock sample	heterogeneous; wide pore size distribution of 3.5–120 μm	ⓝ brine (20% wt. KI); ⓦ oil mixture (85% n-decane, 15% i-iododecane)	0.15 μL/min	synchrotron X-ray µcT @ 70 s interval
Armstrong et al. [19]	microfluidic device	hex-tiling; special 60 μm pore body with 13 μm-wide, 20 μm-long, 5 μm-deep necks	ⓝ decane; ⓦ Millipore water with 50 mg/L fluorescent particles;	2, 20 and 200 nL/min	Olympus iSpeed high-speed video camera @ 2 kfps
Sun et al. [27]	PDMS microfluidic device	sinusoid pore with 0.5 mm radius and 0.2 mm throat; with/without gas chamber	ⓝ air and oil; ⓦ water	100 μL/h max	Zeiss Stemi 2000-CS, unknown fps
Moebius et al. [21]	glass model	monolayer 4 mm-radius glass beads	ⓝ air; ⓦ DI water or ethanol	10–80 mL/min	Marathon Ultra @ 1.2 kfps
Singh et al. [12]	glass model	randomly packed glass beads	ⓝ filtered dodecane; ⓦ water	61 nL/s	µcT
Edery et al. [28]	polymer microfluidic device	$5\times 5$ array of 50 μm3 square posts 75 μm apart	ⓝ water; ⓦ oil	1, 10, 100 μL/h	confocal microscope with camera @ 10 fps
This Work	borosilicate glass microfluidic device [29]	see Section 2	ⓝ DI water with 0.5% vol. microspheres; ⓦ decane	40 μL/min with 1500 mbar drop	Phantom VEO 440L @ 1 kfps

Note: ⓝ stands for “non-wetting” phase, ⓦ stands for “wetting” phase, as defined in [9] p. 7, § 1.2.1.

for comparison. In the work of Moebius et al. [21], the highest bulk flow rates are approximately three orders of magnitude larger than those used in this work. Aside from bulk flow rates, the following two aspects differentiate this work from [21]: (i) the use of a pore network that is approximately 100 times smaller, and (ii) the use of an evolution graph in the spatiotemporal indexing of potential rapid pore displacement events.

1.1. Experiment Framework

Two major technologies are typically used to observe the distribution of fluids as they flow through permeable media: X-ray tomography and microfluidic device imaging. X-ray tomography is capable of taking high-resolution 2D and 3D X-ray scans of a permeable media with two subsequent scans separated by time intervals on the order of seconds or minutes [22,23]. In 3D scans, injection schedules may need to maintain quasi-steady-state (QSS) flows to be able to observe the flow dynamics in successive scans. In the case of 2D scans, the medium is scanned at prescribed locations along the medium, which means there is a time difference between individual scans in each sequence. A microfluidic device, on the other hand, allows quasi two-dimensional fluid flow to be observed in a complex pore network continuously [24]. Glass microfluidic devices provide a platform for reproducible flow experiments because it may capture real or artificial geometries and is made of transparent materials that are chemically-inert to the injectant [25]. The fabrication procedure for such devices has matured over the past several decades, enabling a trove of original quantitative research, especially with the help of high speed cameras that can capture fluid distributions with high frame rates (hundreds, thousands, or more frames per second). Microfluidic flow experiments, therefore, allow examination of high injection rates that induce non-steady flow states. Since a large number of images are captured (hundreds to thousands)in these experiments, proper analysis lies well beyond manual inspection and instead calls for digital image processing workflows.

Existing microfluidic experiments on Haines Jump [19] often attempt to maintain QSS by injecting the non-wetting phase at very low bulk flow rates for prolonged periods of time (e.g., hours), and postulating that the velocity of Haines Jump is independent of the bulk flow rate. However, as the authors confirmed in a previous single-phase microfluidic investigation [26], a transient flow state of sustained duration (in seconds) occurs even when the medium experiences small perturbation of the inlet pressure condition. It is not a stretch to hypothesize that the transient state for immiscible flow is more sophisticated and prolonged due to the presence of evolving fluid–fluid interfaces and loosely trapped fluids that may be drained sometime during the displacement pressure. As part of a series of ongoing experimental efforts that explores the transient state of flow in permeable media, this paper addresses immiscible flow through the automatic identification of rapid pore filling events from high-speed recordings. A summary of relevant literature is shown in in Table 1.

1.2. Interface Tracking

The identification of fluid–fluid interfaces from digital images lies at the core of this study. Different imaging setups produce images with different features that require unique segmentation techniques. Images obtained within the visible spectrum accurately record the shape and location of fluid interfaces as thick darker curves. In principle, edge detection algorithms can then segment the fluid interfaces so that changes to the interface may be identified through frame-by-frame comparison procedures that leverage level set methods. However, all work requires manual analysis and complete automation is still an open problem.

In this work, fluorescent microscopy is conducted with the premise that fluorescent particle concentration (seeding density) is different among different fluid phases. In this framework, interfaces are not explicitly recorded but can be implicitly reconstructed as imaginary curves that separate fluids into different phases by their particle seeding density [30]. More concretely, a Delaunay triangulation [31] is constructed using all segmented particles and is then thresholded and filtered using the distribution of edge lengths to obtain the phase mesh of the non-wetting fluid. The boundary of the phase mesh is the fluid interface. Individual displacement fronts are obtained from the fluid interface by removing the boundaries of the solid matrix from the boundary of the phase mesh. While a subtraction using two consecutive frames is able to remove sensor-attributed background noise that varies slowly, it is unfaithful to the location of particles in either frame when large clusters of particles stay completely still, especially for displacement fronts that stay put for multiple frames. This paper later introduces a classification-based removal of stationary particles so that temporarily stationary particles are not excluded from the fluid mesh.

Tracking graphs provide a direct tool to string the displacement fronts together based on their spatiotemporal proximity [32,33]. As a directed graph, vertices in tracking graphs that have multiple parents or leafs (children) correspond to merging and splitting events. Merging and splitting fronts are crucial to this study because merging signifies the start of a pore filling event while splitting signifies the the end of another pore filling event.

1.3. Quantitative Analysis

Various quantitative descriptions of multiphase flow through permeable media are of academic interest. Quantities related to dynamics, (e.g., capillary pressure, contact angle, curvature, velocity) and quantities regarding flow state (e.g., transient state duration, volume of displacement) are all subject to quantitative analysis that serves various purposes [8,11]. Instead of using the complete video, computation may be performed within a much smaller spatiotemporal window that encloses a phenomenon of interest, providing speed and accuracy boosts to algorithms such as those used in particle image velocimetry (PIV), e.g., [34]. Local computation is made possible by using tracking graphs in combination with traditional quantitative analysis.

2. Materials and Methods

2.1. Experiment

The microfluidic experiments in this work use artificial geometries (Figure 1). This device is a 1′′-by-2′′ (2.54 cm × 5.08 cm) borosilicate glass model that encases a 0.5′′-square (1.27 cm × 1.27 cm) flow field connected to a pair of injection/outlet ports through distribution channels that obey Murray’s law [35]. Its fabrication uses photolithography and wet etching [36,37]. To create observable interfaces, decane serves as the displaced phase while a water-based fluorescent particle solution serves as the displacing phase. The fluorescent particles used in this experiment are yellow-green Carboxylate-Modified FluoSpheres^® with an excitation/emission maxima of 505/515 nm and a diameter of 1 μm. These microspheres are supplied by ThermoFischer Scientific and are diluted to 0.5% solids from the supplier-provided concentration of 2% solids.

The experiment setup, shown in Figure 2, uses a syringe pump to inject decane and uses FLUIGENT pumps to inject the fluorescent particle solution and DI water. In preparation of each flow experiment, the microfluidic device is saturated with the fluorescent particle solution followed by refocusing of the microscope. Next, DI water is injected into the device to displace the fluorescent particles. Background particles deteriorate the performance of the proposed interface tracking algorithm, so it is important that the medium is as clean as reasonably possible. As the final preparation step, decane is injected until the medium is fully saturated. To run the experiment, the fluorescent particle solution is injected at 150, 750 and 1500 mbar and then video recording commences with a resolution of at 2048 × 1536 px at 1000 fps.

2.2. Image Segmentation

Proper image enhancement and segmentation are required to accurately obtain particle centroids in preparation for interface tracking. Different video data require different treatment, but all involve solid mask generation, particle segmentation, and stationary particle removal. This subsection describes various algorithms scripted for this research that make use of the image processing package by MATLAB (MATLAB and Statistics Toolbox Release 2020b, The MathWorks, Inc., Natick, MA, USA). All manuscript figures apart from Figure 1, Figure 2, and 5 are generated in MATLAB.

The maximum image, minimum image, range image (Figure 3a), and a binary solid mask are required for each experiment, especially if the field of view is not constant from one experiment to the next. The maximum/minimum image consists of the maximal/minimal pixel value for each pixel throughout the video. The range image is their difference, and the solid binary mask is obtained from its binarization. Global thresholding can be used when the lighting in the field of view is uniform, otherwise lighting-aware adaptive binarization should be used instead [38,39,40]. Both binarizations produce a preliminary solid mask that requires fine-tuning (refer to Figure 3b), see Listing 1. Conversion from the artificial integer coordinates $(u,v)$ to centroid coordinates $(x,y)$ requires an affine transformation with six degrees of freedom through homogeneous coordinates

$$\left[\begin{array}{ccc}a& b& c\\ d& e& f\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$

(1)

which can be solved or fitted as long as three or more non-collinear solid grains are selected during step (f) of Listing 1. In addition, steps (c) through (e) remove small imperfections from the preliminary solid mask generated by step (b), which may be present in either the void space or the matrix space.

For particle segmentation, the minimum image is first subtracted from each frame to remove sensor-attributed static background noise. Because the illumination is mostly even throughout the video while the random background noise varies slowly throughout the video, a simple procedure is used for noise removal and particle centroid computation, (Listing 2). Note that the centroid of each connected component in the particle mask is recognized as the position of a particle only when weighted by pixels in the clean frame obtained from step (a).

Listing 1. Computer -aided procedure to generate a fine-tuned solid mask from the range image.

(a)  Adaptively binarize the range image with dark foreground polarity and a default 50% sensitivity; (b)  Isolate the largest connected component, remove all other pixels, then invert the binary mask; (c)  Fill all holes, remove connected components smaller than 4000 px in area (approximately 0.127% of the image size), and dilate by 1 pixel; (d)  Invert the binary mask, remove connected components smaller than 4000 px in area again, and dilate by 1 pixel; (e)  Invert the binary mask to obtain a half-tuned solid mask. (f)  Hand-select solid grains that are round, not on the boundary, and not involved in snap-off to tally their mean radius, centroid coordinates, then assign artificial integer coordinates to them (refer to Figure 3c). (g)  Use least squares analysis (refer to Equation (1)) to recompute all solid grain centroids based on their artificial integer coordinates and replace the half-tuned solid mask with a tiling of the circular solid grain with mean radius at the recomputed centroids (refer to Figure 3d).

Listing 1. Computer -aided procedure to generate a fine-tuned solid mask from the range image.

(a)  Adaptively binarize the range image with dark foreground polarity and a default 50% sensitivity; (b)  Isolate the largest connected component, remove all other pixels, then invert the binary mask; (c)  Fill all holes, remove connected components smaller than 4000 px in area (approximately 0.127% of the image size), and dilate by 1 pixel; (d)  Invert the binary mask, remove connected components smaller than 4000 px in area again, and dilate by 1 pixel; (e)  Invert the binary mask to obtain a half-tuned solid mask. (f)  Hand-select solid grains that are round, not on the boundary, and not involved in snap-off to tally their mean radius, centroid coordinates, then assign artificial integer coordinates to them (refer to Figure 3c). (g)  Use least squares analysis (refer to Equation (1)) to recompute all solid grain centroids based on their artificial integer coordinates and replace the half-tuned solid mask with a tiling of the circular solid grain with mean radius at the recomputed centroids (refer to Figure 3d).

Listing 2. Automatic procedure for noise removal and particle segmentation.

(a)  Subtract the minimum image from each frame to obtain a clean frame; (b)  Globally threshold the clean frame to obtain a binary particle mask; (c)  Apply a median filter to the particle mask, fill holes, then remove connected components smaller than 8 px in area; (d)  Compute the centroids of bright clusters in the clean frame using the particle mask.

Listing 2. Automatic procedure for noise removal and particle segmentation.

(a)  Subtract the minimum image from each frame to obtain a clean frame; (b)  Globally threshold the clean frame to obtain a binary particle mask; (c)  Apply a median filter to the particle mask, fill holes, then remove connected components smaller than 8 px in area; (d)  Compute the centroids of bright clusters in the clean frame using the particle mask.

Removal of stationary particles is conducted as a three-part process. The prerequisite is to ensure that an extra number of frames (e.g., 20 frames) are recorded before the flow enters the field of view. The first part is to cluster all particle centroids in these extra frames with a distance criterion of 1 px to generate the list of stationary particles. The second part computes the stationary lifespan of each stationary particle before which they start moving. It requires checking the presence of each stationary particle throughout the video as judged by the same spatial proximity of 1 px. Note that the presence of each particle throughout a video frame is effectively a binary array and is subject to falling edge detection through logistic regression. The final part is to loop through all frames and remove any particle that is within 1 px of a stationary particle still in its stationary lifespan.

2.3. Interface Segmentation

The remaining particles lie in the void space and the following computer-aided procedure extracts from it the displacement fronts; see Figure 4.

To begin, a frame is chosen where the displacement is roughly halfway through the field of view. This frame can be chosen manually or automatically based on a plot of the average pixel value versus the frame number. This choice is made to guarantee that the two fluids of different seeding densities are both well-represented in the frame. A Delaunay triangulation D is then constructed from all remaining particles in this frame and the length of all edges in D are calculated. Next, given a fixed edge length threshold, remove from D all triangles whose side lengths all exceed the threshold, then remove all hanging edges and vertices. After removal, a phase mesh M is obtained for the phase with higher seeding density. Different edge length thresholds are tested en masse to arrive at an optimal value L where erroneous faces are mostly removed and correct faces are mostly kept in the chosen frame.

Next, use L to generate a phase mesh for each frame. The segmentation of displacement fronts from a phase mesh M is grounded on the idea that displacement fronts are the parts of the phase mesh boundary $\partial M$ that do not trace along the boundary of the solid grains. Note that each frame may contain multiple phase mesh boundaries; therefore, those consisting of 20 or fewer particles are ignored. The automatic per-frame procedure described in Listing 3 describes the extraction of the phase boundary into individual displacement fronts.

Listing 3. Automatic procedure for segmenting a phase mesh boundary into displacement fronts.

(a)  Find all connected components of the reconstructed solid mask; (b)  For each vertex and each edge in the phase mesh boundary, compute its closest connected component and the distance; for edges, use its centroid; (c)  Remove all edges with a distance not exceeding 20 px; (d)  Partition the remaining edges into displacement fronts by vertices with a distance smaller than 10 px; (e)  Remove displacement fronts that fail to straddle two different solid grains.

Listing 3. Automatic procedure for segmenting a phase mesh boundary into displacement fronts.

(a)  Find all connected components of the reconstructed solid mask; (b)  For each vertex and each edge in the phase mesh boundary, compute its closest connected component and the distance; for edges, use its centroid; (c)  Remove all edges with a distance not exceeding 20 px; (d)  Partition the remaining edges into displacement fronts by vertices with a distance smaller than 10 px; (e)  Remove displacement fronts that fail to straddle two different solid grains.

2.4. Evolution Analysis

With all displacement fronts extracted, a tracking graph G can be built by exploiting the spatiotemporal relations between fronts segmented throughout the video. A tracking graph for displacement fronts is a directed graph that tracks their movement as they move or experience topology-changing events (emerge or disappear, enter or exit, merge or split). Instead of building the full graph G, its contraction $G^{\prime }$ called the evolution graph is built.

The vertices of the evolution graph $G^{\prime }$ are branches of evolution, which are collections of moving displacement fronts from G after (inclusive) and before (exclusive) topology-changing events. A per-frame incremental construction of $G^{\prime }$ is possible by gauging the spatiotemporal proximity of two displacement fronts through their bounding box. The first displacement front is the end of a branch of evolution that is still actively collecting displacement fronts as of frame k, whereas the second displacement front are those yet to be collected in frame $k+1$. The bounding boxes used could be level curves surrounding the displacement fronts by a certain distance, but this work uses circumcircles [41].

The collisions between the bounding boxes of the final displacement front of each active evolution branch in frame k and each displacement front in frame $k+1$ are used to build the so-called adjacency matrix (see Figure 5). The adjacency matrix dictates whether the active evolution branches should be terminated or allowed to collect new displacement fronts, and also create directed edges when proper.

It is important to resolve the so-called “unresolved events” where multiple branches are involved in simultaneous merging and splitting. An unresolved event is indicated in the adjacency matrix as the existence of a row and a column, each containing at least two non-zero entries that share a common non-zero entry. To resolve these events, binary integer programming is applied to the unresolved submatrix such that merging is forbidden (if the submatrix has more rows than columns) or splitting is forbidden (the other way around). The objective is to minimize the sum of distances between the resolved active evolution branches and their corresponding resolved displacement fronts in the upcoming frame. The solution of the binary programming is used to replace the unresolved submatrix entry-by-entry.

Finally, an interactive viewer is available as part of the software package for users to visually browse through the evolution graphs (see Figure 6 and Data Availability Statement). Manual branch merging, interface splitting, pruning, and automatic updates to parent–child relations can improve imperfections in the automatically generated evolution graph. Visual inspection allows the correspondence between throats and evolution branches to be established for quantitative analysis of desired phenomena.

3. Results

The quantitative analysis of rapid pore filling events lies on a special type of evolution branch with exactly two parent branches and two child branches. Such branches are referred to as “nuclear branches” as each involves the merging of two parent branches and the splitting into two children branches. They form a basic study unit within the evolution graph. For the 1500 $\mathrm{m}$$\mathrm{bar}$ experiment, Figure 7 shows the correspondence between throats in the field of view with branches in the evolution chart. Three nuclear branches, numbers 7, 16, and 31 are selected. Note from Figure 7 that the parents of branch 7 are branches 4 and 5, and the children of branch 7 are branches 11 and 12.

To identify the occurrence of Haines jump, the volume of displacement is tracked for each nuclear branch starting from slightly before the merging of the parent branches and ending slightly after the splitting into the child branches. As the recording is two-dimensional, it is sufficient to keep track of the filled pixels and the pixel filling speed, see Figure 8. A derivative plot called the filling speed plot is computed for each branch. Several notable features are visible from Figure 8. For all nuclear branches, there exist noticeable rapid pore filling events the moment when displacement fronts from parent branches merge together. Noticeable rapid pore filling events also occur when the displacement front in a throat splits.

Figure 9 shows distribution of the two immiscible phases before and after rapid pore filling events.

4. Discussion

The behavior of the curves between the dashed merge and split lines in Figure 8 varies dependent on the bulk flow rates.

• For the higher bulk flow rate (Figure 8a–c) of 24.7 μL/min, the filling speed plot is concave upwards and almost smooth. This indicates that, as the merged displacement front travels through the throat, the filling rate noticeably decreases to the lowest value at around the narrowest part of the throat and then noticeably increases as the merged displacement front enters the other pore of the throat until it splits into two. The number of frames between the lowest filling rate and the split lasts 15, 14, 12 ms, respectively.
• For the lower bulk flow rate (Figure 8g–i) of 4.2 μL/min, the filling speed within the throat is mostly slow, but there is a single burst of rapid filling around the narrowest part of the throat.
• For the intermediate bulk flow rate (Figure 8d–f) of 13.6 μL/min, single bursts occur around the narrowest part of the throat but are significantly less noticeable due to a meaningful bulk flow rate. The single burst typically follows a stationary waiting period.

Through observation, it is clear that Haines jump is most significant when two displacement fronts merge together. It is somewhat significant as a displacement front splits and enters into two throats. When the bulk flow rate is low, Haines jump is observable as a single burst near the narrowest part of the throat, which is less or completely not observable when the bulk flow rate is high.

The evolution graph provides an opportunity to design and study other metrics of the fluid flow. Another example metric is the “nuclear rate” defined as:

$$\mathrm{Nuclear}\mathrm{Rate}=\frac{\mathrm{Number}\mathrm{of}\mathrm{Nuclear}\mathrm{Branches}}{\mathrm{Number}\mathrm{of}\mathrm{Good}\mathrm{Branches}},$$

(2)

where a “good branch” is an evolution branch with at least one parent and at least one child. While the 150 $\mathrm{m}$$\mathrm{bar}$, 750 $\mathrm{m}$$\mathrm{bar}$, and 1500 $\mathrm{m}$$\mathrm{bar}$ experiments each have 38, 40, and 37 branches, only 28, 28, and 29 good branches are filtered of which only 5, 8, and 9 are nuclear, resulting in nuclear rates of 18%, 29%, and 31%, respectively. Because higher nuclear rate suggests that the flow pattern is more regular, one may compare the transport network evolution graphs (refer to SI) to conclude that the disruptive effect that Haines Jump has on flow patterns is remedied by higher bulk flow rates.

Lastly, high-speed videos (see Data Availability Statement) contain the nine branches depicted in Figure 8. These videos are cropped using spatiotemporal information stored in the evolution graph, e.g., bounding boxes of the displacement fronts, the first-seen frame, and the last-seen frame. Such videos, produced from a spatiotemporal segmentation of high-speed recordings, could be stored for localized PIV analysis, e.g., if high-density PIV is used for the entire video, low-density PIV could be used for particles that form the displacement front. Potential future work includes the joint application of PIV and evolution graphs.

5. Conclusions

This paper introduces an evolution chart of displacement fronts for multiphase flow in permeable media through image segmentation and tracking graphs. High-speed video recordings are obtained through microfluidic flow experiments and are analyzed for merging and splitting events. Displacement through throats is isolated from the evolution graph to study the volumetric rate of displacement through the throats. Occurrences of Haines Jump are identified from the pixels-filled-each-frame plot at the onset of interface merging, interface splitting, and at around the narrowest part of the throat. Specifically, Haines jump is signified by rapid shifts in the pore filling process, i.e., total pixels filled, and associated sudden changes in velocity, i.e., pixels filled per frame. The bulk flow rate does not affect the significance of Haines Jump at merging and splitting, but the single burst at the narrowest part of the throat may be eclipsed by high bulk flow rates.

While the original goal of this research is the fully automated generation of evolution graphs from high-speed microfluidic experiment recordings, certain aspects of this workflow fall short of this expectation, and are either manual or, at best, computer aided. For instance, all thresholds used in this work are globally thresholded from representative frames and identical thresholds are applied to all frames. The fact that such a workflow functions correctly is mostly attributed to the highly constructed nature of the pore network and grain shape. Complex networks, especially multiscale networks, could benefit from the use of adaptive algorithms, neural networks, or a combination thereof. Rule-based pruning of the evolution graph, including the merging of branches (in parallel or in series) and their deletion, could also be formulated to improve the level of automation.