Visual Analysis of Two-Phase Flow Displacement Processes in Porous Media

We present the visual analysis of our novel parameter study of porous media experiments, focusing on gaining a better understanding of drainage processes on the micro-scale. We analyze the temporal evolution of extracted characteristic values, and discuss how to directly compare experiments that exhibit processes at different temporal scales due to varying boundary and physical conditions. To enable spatio-temporal analysis, we introduce a new abstract visual representation showing which paths through the porous media were occupied to what extent, e.g., allowing for classification into viscous and capillary regimes. This joint work of porous media experts and visualization researchers yields new insights regarding immiscible two-phase flow on the micro-scale toward the overarching goal of characterizing flow based on boundary conditions and physical fluid properties.


Introduction
A thorough understanding of the underlying physical processes occurring during two-phase flow is of high importance, since it will enhance our abilities to effectively build models for prediction and simulation. Here, two-phase flow term refers to the interactive flow of two fluids in a confined space surrounded by solid walls. The two fluids interact via their common interfaces, and they have no or very limited mixing potential with each other. Based on this assumption, at least one of the fluids has to be a liquid so that immiscible flow can take place. The fluid which has the stronger affinity to the solid phase (increased spreading) is usually referred to as the wetting phase, while the fluid with the least affinity is referred to as the non-wetting phase. The processes involved, like drainage (where the non-wetting phase forcibly displaces the wetting phase) and imbibition (the inverse process which can also take place without external force), have been the topic of extensive work in porous media research due to the degree of degeneracy that is expressed via hysteresis in the capillary pressure-saturation relation [JS90,Hil06], and how this hysteresis can be lifted [HG93]. Such fundamental physical processes play a crucial role in the handling and optimization of industrial and environmental applications with a dominant socioeconomic impact, including enhanced oil recovery [AMW * 17], soil remediation [FHR04], printing [KLG10], etc. Given the fundamental significance of two-phase flow, numerical [JNHD10, PSW09, BB92] and experimental [WS13,CCR97] schemes have been developed to describe the underlying physical processes in the most efficient way given the complexity of the process and the physical constrains to acquire pore-scale experimental evidence. However, a better understanding of the complex dynamic interplay between viscous and capillary forces during a drainage event, and the impact of the pore space geometry, boundary conditions, and physical properties of the involved fluids is crucial for further advancements in the field.
The close collaboration between porous media researchers and visualization experts presented in this work yielded new insights toward the overarching goal of flow characterization based on the correlation of boundary conditions and physical properties of the fluids in a porous medium. We argue that this lays a crucial foundation for further advancements in the domain. For this, we conducted and analyzed eleven two-phase displacement experiments, employing an artificial, Poly-Di-Methyl-Siloxane (PDMS) micro-model under varied viscosity ratios M between the invading and the defending fluid under a drainage scenario, and capillary numbers Ca (which provides a measure of the competition of the viscous vs the capillary forces for the corresponding boundary conditions). In this work, we aim at linking the flow behaviour to Ca and M, with the ultimate objective being the inclusion of the local porous medium geometry as a determining measure. This is expected to enable the investigation of how and why flow characteristics can change with respect to the heterogeneity and the local geometrical properties of the pore space. The ability to visualize experiments in a compact way and to quantify values of interest allows us to identify similarities and differences between the processes, to categorize the experimental results and to further research the relation between boundary conditions, physical fluid properties and the resulting flow. This requires the extraction and expressive presentation of relevant quantities, which also need to be made directly comparable as-depending on the viscosity ratio M and capillary number Ca-the speed of flow progression is in different orders of magnitude (ranging between 330 ms and 593 s, cf. Fig. 1). In particular, we visually analyze the occupancy of the pore space by each phase, the breakup of the connected phases due to the competition between viscous and interfacial effects, and coalescence of disconnected blobs, either to form bigger blobs or to reconnect with the continuous phase, and the production of interfacial area and the mechanisms behind that. As the first work in the domain, we also study how these parameters correlate to the pore space and the local geometrical properties.
In the following, we first review related work in visualization as well as porous media research (Sec. 2), before describing our experiments and the extraction of phases (Sec. 3). This provides the basis for determining visualization data for comparable analysis, including the extraction of quantities of interest (e.g., interfaces, saturation, inlet connection), and accounting for different time scales across experiments via temporal normalization (Sec. 4). This then allows to conduct a quantitative comparative analysis of time-dependent processes (Sec. 5). For spatial analysis in the pore space, we employ an aggregate network representation of the porous structure that allows to investigate which and to what degree paths through the porous medium were taken by a fluid (Sec. 6). Finally, we summarize the impact of this work on porous media research and discuss limitations as well as directions for future work (Sec. 7). Positioning in Porous Media Research. Various macroscopical continuum models for immiscible two-phase flow in porous media have been proposed and widely applied. Van Genuchten [VG80] introduced a "Darcy-scale" model which proposed that the capillary pressure is a simplified non-linear relation of the wetting phase saturation. Based on thermodynamical principles and approaches, it has been shown that capillary pressure does not only depend on the fluids' saturation, but crucially also on its spatial distribution in the porous domain [SBR * 16]. Numerous attempts have been made to embed the pore-scale geometry and the spatial fluids' configuration in the pore space in continuum models. One of these is based on a multi-scale approach that includes the specific interfacial area between the two phases as a separate state variable in addition to phase saturation [HG93]. Such theories including additional state variables are commonly denoted as "extended theories". They can also include process variables dependent on higher-order morphological features [KS14] which can be expressed with the first three Minkowski functionals M0 -M2 [Web01]. One common underlying assumption is that none of the two phases gets disconnected during the process. However, this only holds for ideal conditions which are hardly ever observed in practice, even in very well-conditioned experiments for drainage and imbibition (like the one discussed in this paper).

Related Work in Visualization and
Micro-model experiments. Two-phase flow has been studied in experiments by means of natural or artificial porous media. As natural porous media mostly rock samples are used [BRB * 16, BOK * 13, ABB13], while for the latter there are many approaches to manufacture porous media with well-defined properties [Dar57,GCS88,Sha98,CCR97]. In this work, the artificial porous media used are PDMS micromodels. As first defined by Karadimitriou and Hassanizadeh [KH12], the term "micromododels"" is used to describe an artificial, transparent porous medium with similar average properties, like pore size distribution, porosity and permeability, to a natural porous medium, with the pores measuring on the micro-scale and the overall dimensions on the centimeter scale. Micro-models are commonly produced either by etching processes [MK61], or lithography [Tho83]. Soft-lithography specifically is widely employed due to its comparably low cost and the high precision with which copies of the same flow network can be made [XW98]. As this type of micro-models is transparent, images taken during flow and displacement can be acquired in high temporal and spatial resolution with optical means.
In the literature there is a significant volume of works which make use of micro-models and optical illumination in order to visualize the effects taking place, at a spatial resolution of a few microns per sensor pixel. With such a micromodel, Lenormand et al. [LTZ88] investigated the effect of the boundary conditions and the physical properties of the involved fluids, both numerically and experimentally. They performed drainage experiments for a combination of boundary conditions and physical fluid properties and created a wellknown map of the parameter space and the resulting displacement regimes. They categorized the flow regimes into capillary fingering, viscous fingering as well as stable front and were able to match these regimes both numerically and experimentally. Probably due to the scarcity of computational resources and image analysis tools they limited themselves to qualitative observations of the experimental images, without being able to draw any quantitative conclusions. Cheng et al.
[CPNNG04] performed drainage and imbibition experiments, in an attempt to investigate the role of specific interfacial area as a separate state variable in two-phase flow studies, under quasi-static flow conditions. This is the first experimental proofvia recording features on the micro-scale-that specific interfacial area can potentially be used as a state variable in two-phase flow studies. Karadimitriou et al. [KHJNK14] discuss the possibility of the inclusion of specific interfacial area as a state variable for twophase flow, but under dynamic conditions. They used a micro-model based on Poly-Di-Methyl-Siloxane (PDMS) and akin to Cheng et al.
[CPNNG04] they performed drainage and imbibition experiments, including scanning curves, for three different fluxes, ranging from a capillary to a highly viscous flow regime. They extracted information on phase saturation, local capillary pressure, local contact angle, location of the apex of the interfaces and constructed again the corresponding surface in the capillary pressure-saturation-specific interfacial area space. They concluded that specific interfacial area can be included as a state variable in two-phase flow studies under quasi-static but not under dynamic conditions. While they attributed the observed mismatch mostly to the disconnected phases induced during the process, they could not quantify the disconnected phases in terms of saturation and specific interfacial area (which we do in this work).

Experiments and Image Analysis
We conducted experiments to study the distribution of the two fluid phases depending on the flow boundary conditions and the change in the physical properties of the fluids in combination with the pore geometry (Sec. 3.1). Initially, the pore space is occupied by a wetting phase, which is then displaced by a non-wetting phase during the course of the experiment (this scenario is also known as primary drainage). The three involved phases-solid, wetting fluid, and non-wetting fluid-are extracted from the resulting time-dependent ensemble for further analysis (Sec. 3.2). While our experimental setup allows for the comparably simple distinction between phases in general, high-precision segmentation is challenging. In particular, trapped air and dust particles and model imperfections can lead to sizeable deviations that require manual correction to complement automatic processing.

Experiments
All experiments use an artificial porous medium made of PDMS via soft lithography [XW98, KMK * 13], replicating the C-type network used in the work of Sivanesapillai and Steeb [SS18]. The flow network has a mean pore size of 410 µm, and a porosity of 0.44 (as a comparative measure between the available space for flow and the bulk volume of the porous medium). A primary drainage scenario is considered in which the flow network that is initially filled with the wetting phase (Fluorinert FC-43) is invaded by the non-wetting phase (water dyed with ink), as mentioned earlier. To achieve higher viscosities for the invading phase, thus increasing the viscosity ratio, water is mixed with glycerol at given concentrations so as to achieve viscosities equal to and ten times higher than this of Fluorinert. A 1 mL glass syringe in combination with a mid-pressure CETONI neMESYS 100N [CET] generates the flow in the network via an 1/16 inch ID Poly-Tetra-Fluoro-Ethylene (PTFE) tube. Right before the inlet of the flow network, an MPS2 [ELV] flow-through pressure sensor (maximum pressure 1 bar) measures the inlet pressure continuously. Another identical sensor was also used right after the outlet of the flow network. With this, pressure drop across the whole network can be measured at a frequency of 10 Hz.
The transparent nature of PDMS allows capturing the processes taking place in the pore space in real time, by using transmitted light microscopy ( Fig. 2). For this purpose, a custom made microscope was developed, which is able to visualize samples with a resolution of 0.5 to 20 microns per camera pixel. In this case, flow can be captured at a resolution of 8 microns/pixel, at a frame rate varying from 1 to 15 fps, depending on the speed of the investigated process indicated by the capillary number. A total of eleven experiments was conducted; the flow boundary conditions included capillary numbers Ca ∈ {1e−5, 1e−4, 1e−3, 1e−2}, and the viscosity ratios were M ∈ {0.2, 1, 10}. For the combination of Ca = 1e−2 and M = 0.2, capturing the processes reliably was impossible due to (i) the very high flow velocity that (ii) is also sufficient to deform the pore space. Accordingly, this yields irreproducible results and naturally prohibits meaningful comparison to the other experiments.
A gray-scale, 8-bit image sequence is recorded for every experiment (e.g., Fig. 3) The wetting phase, PDMS, and some trapped air appear as light-gray to white pixels, while the non-wetting fluid, interfaces and dust particles appear as dark-gray to black pixels. The micro-model is flushed between experiments to remove residual liquids. This leads to translational and rotational displacements in the recorded image sections as well as slight perspective changes. To compensate for this, we capture one mask image for each experiment after saturating the micro-model with the wetting phase (here, the wetting phase appears in dark gray, cf. Fig. 4a), and employ the Enhanced Correlation Coefficient (ECC) for correction [EP08]. For this, we compute a homography between the mask of the experiment with (Ca = 1e−2 ,M = 1) and the masks of all other experiments. The mask is further used as reference for distinguishing the three phases in the input images below, representing the pore space.

Phase Segmentation
Phase segmentation is conceptually straight-forward: masks allow to distinguish between solid and fluid phases, and the experimental setup already relatively clearly distinguishes fluids as light and dark gray. However, automatic phase segmentation would erroneously recognize non-existent structures and the detected interfaces are faulty. Note that accurate masks are particularly important: besides segmentation, masks also crucially provide the basis to create the (e,f) The non-wetting phase is finally partitioned into parts with and without uninterrupted connection to the inlet (cf. Sec. 4.1).
structure of transport networks (cf. Sec. 6.1 (1)). In the recorded masks, the walls of the PDMS structure appear as ≈ 10 pixel wide contours with a dark-gray to black tone that becomes lighter towards the edge (cf. closeup in Fig. 4a). The contours exhibit noise in brightness and thickness, which is caused by the fact that the crosssectional shape of the solid edges in "depth" in 3D are not precisely rectangular, and the illumination is not perfectly homogeneous. Additional artifacts are caused by dust particles from the air, which deposit on the model or the optical setup. They appear as stains with different gray values in the solid or fluid phase and their contours might interfere with those of the PDMS surface. Accordingly, flow paths may not be recognized in areas where there is only a gap of a few microns in the porous media between the walls of the solid to begin with. Furthermore, the model geometry features dead-end pores at the top and bottom which are connected to the remaining structure by a single pore throat. These cavities trap air bubbles that cannot be flushed out, and the respective void space becomes inaccessible to the fluids. In the recorded images, these air bubbles appear similar to PDMS as translucent material with a dark-gray to black contour, and they are manually removed. Additionally, larger, non-persistent artifacts from the image that do not impact the actual experiment are accounted for by merging them with their surrounding (Fig. 4b). The images are further denoised by applying a median blur filter with a rectangle of kernel size 5 × 5, and thresholding subsequently converts it into a binary image. Remaining small artifacts (below the minimum pore throat size) are further reduced by applying three iterations of the morphological operations opening and closing with a rectangular kernel of size 3 × 3. Finally, contour extraction [S * 85] is applied by redrawing the binarized image and excluding all contours with an area < 150 pixels (Fig. 4c, cf. Fig. 5a for full image). Parameters in the automatic part of the pre-processing were determined empirically; they were adjusted to reduce artifacts as effectively as possible without removing fine geometry or corrupting the topology of the pore network.
Input images are pre-processed automatically akin to the masks (without manual adjustment). With this, each image can be segmented into the three investigated phases (Fig. 5c). The solid phase corresponds to all pixels that are on in the pre-processed, binarized mask. The wetting phase consists of all those pixels which are on in the binarized input images but off in the respective mask. All remaining pixels are accordingly assigned to the non-wetting phase.

Visualization Data For Comparative Analysis
We now extract visualization data from the segmentation results that allow us to investigate processes quantitatively (Sec. 4.1). Direct comparison between experiments further requires dealing with different temporal scales (Sec. 4.2).

Quantities of Interest
Saturation of Fluid Phases Strictly speaking, the input images are two-dimensional projections of three-dimensional processes. However, the loss in depth information of processes is negligible: the considered porous medium is orders of magnitude smaller in depth than in width and height, at 100 µm in comparison to 20 mm × 15 mm. Accordingly, we can simply approximate the saturation of the fluid phases by the relative fraction of pixels per phase in relation to non-solid pixels in total.
Interfacial Areas Given that depth can be neglected in our case, the interface between two phases is proportional to the summed length of the interface lines. Our interface extraction generates one binary mask per phase combination, which contains one-pixel-thin lines at the position of the respective interface. For this, the mask of the first phase is subtracted from a dilated version of itself (computed via a 3 × 3 rectangular kernel). Then, a bitwise AND is applied to the result and the second considered phase. The result then undergoes a series of filter operations that remove all those pixels that have both horizontal/vertical and diagonal neighbours (essentially removing staircase shapes along the interface line).
The interfacial area is then computed via the Euclidean distance between the centers of adjacent pixels (Fig. 6). First, we calculate for each pixel the sum of distances from its center to the border of adjacent pixels on the interface via 2D convolution (distance convolution). The entries of the employed 3 × 3 kernel accordingly exhibit the distances to the central element. Second, we apply the respective interface mask, and determine the total interface length via the sum of pixel values.
Disconnected Fluid and Inlet Connection A commonly employed assumption in the domain is that fluids do not get discon-nected (cf. Sec. 2). However, we consider the explicit investigation of connected components to be of crucial importance for the analysis for a number of reasons. The disconnection of the non-wetting phase during primary drainage is an effect that results from the synergy of the capillary number, the viscosity ratio, and the geometry of the porous medium. For low capillary numbers, there is no specific reason for the disconnection of the non-wetting phase. However, it is expected for high capillary numbers that the non-wetting phase will get disconnected due to the developed shear stresses, and their competition against interfacial forces, as the only means to absorb the energy fed into the system. However, the location of these disconnections also depends on the local geometry of the porous medium, and how they can locally enhance the shear stresses due to local constrictions. We expect that considering number and location of the disconnected non-wetting parts helps in quantifying the effect of the pore geometry on the evolution of flow for various boundary and physical conditions. In addition to this, the extended theories of two-phase flow can be better evaluated, at later times, since the discrimination of the type of the newly formed interfaces with respect to the phases involved and their connectivity to the corresponding phase is an intrinsic characteristic of them.
To this end, we split the fluid phases into parts that are connected or disconnected from their respective inlet (e.g., on the right for the non-wetting phase, Fig. 5e and f). We check which parts of fluid are connected by determining their distance to the respective image border where their inlet is located (a distance of 40 pixels in the 2448 pixels wide frame yielded reliable detection of both fluid parts across all experiments). Disconnected components in the non-wetting phase correspond to blobs and ganglia, whereas we do not consider small droplets below 50 pixels (blobs occupy a single pore, while ganglia stretch across pores). We do not distinguish between blobs and ganglia here, but this is planned for future work.

Temporal Normalization
Our analysis particularly focuses on the comparison of the experiments. However, the experimental parameters have a strong influence on the speed at which the displacement processes take place, and accordingly experiments are recorded at different frame rates. Since the image recordings are also started and stopped manually, there are no directly given reference points. To make processes at different temporal scales directly comparable, we derive three reference points in time based on the flow of the non-wetting phase: entry t i of a phase into the domain, breakthrough t b when an uninterrupted connection between the left and right boundary is formed (Fig. 3), and finally end te when the experiment is stopped (this is done manually when all relevant processes have been captured). Especially the universal and dimensionless time stamp for breakthrough is critical for further analysis: depending on the boundary conditions, an event, like the continuation of saturation change, may or may not happen after breakthrough. In combination with the actual physical time, this provides the basis for categorizing experimental data. Based on these reference points, we employ different temporal normalization schemes: relative timet = t−ti te−ti (t ∈ [0, 1]), breakthrough normalizationt = t−ti t b −ti (t ∈ [0, 1]), and breakthrough scalingt = t−ti t b −ti (t ∈ [0, te−ti t b −ti ], with fixed point t b = 1). Relative time has the advantage that it gives an impression of the speed with which the displacement processes take place. The behaviour of the fluid before and after breakthrough is of particular interest in the analysis. Breakthrough normalization and breakthrough scaling enables direct comparability of the experiments in this respect.

Temporal Evolution of Characteristic Values
The extracted visualization data and different means to deal with temporal heterogeneity now allow us to quantitatively compare occurring processes across different experiment runs (Fig. 7).
Interface between wetting and non-wetting phase (Fig. 7a) Four groups of two to three curves can be identified, whereas the interface length starts increasing at similar points in time: Interestingly, this means that a configuration of Ca and M yields similar behavior in this regard to a configuration with Ca being one order lower and M one order higher for the considered parameter range. It is also indicated that when we are in a capillary regime already (with (Ca 1e−5, M 1) and (Ca 1e−4, M 10)), the next "step" of (Ca 1e−5, M 10) will still be capillary dominated, which means that the phase topology does not change and the resulting impact on the interface length is comparably small. The investigation of respective transport networks that exhibit similar flow paths in Fig. 1 complements this observation. Generally, it can be seen that higher capillary numbers yield a larger interfacial area, which can largely be attributed to increasing phase disconnection (cf. Fig. 7c). This effect is also one of the main reasons why we need to differentiate between connected and disconnected phases, since the extended theories and their validation only account for connected phases.
Saturation of the Connected Wetting Phase (Fig. 7b) Two groups of curves can clearly be identified in the diagram: one where the saturation settles between 0.3 and 0.5 ({(1e−5, 1), (1e−5, 10), (1e−5, 0.2), (1e−4, 1), (1e−4, 10)}), and one where the saturation after the breakthrough approaches zero (the rest). This graph clearly reflects that for increasing capillary number, the breakthrough saturation also increases. In addition to this, for increasing viscosity ratios, the saturation change after breakthrough evolves smoothly and nearly undisturbed, resulting in a nearly total recovery of the wetting phase from the pore space (the saturation of the connected wetting phase approaches zero). These findings indicate the following behavior. For low capillary numbers, where the flow regime is capillarity-dominated and capillary fingers make their appearance, as soon as breakthrough takes place, there is limited possibility for the saturation to change, since the connectivity of the wetting phase is also limited due to the two-dimensional nature of the porous medium, and topological changes are not favored. As the capillary number increases, we gradually switch from a capillary regime to a mixed capillary and viscous regime. The contribution of the viscous forces enables the re-mobilization of the the disconnected wetting phase, which eventually leads to an increased final saturation than this of breakthrough. When we increase the capillary number even further (for viscosities M = 0.2  and M = 1), then blobs of disconnected non-wetting phase-caused by the competition between shear stresses and capillary forces-are rushed into the network (e.g., (Ca 1e−5,M 02) in Fig. 3a, also see below in Fig. 7c). This further increases the saturation, but in a erratic way due to the synergy with viscous fingering. As the viscosity of the non-wetting phase also increases (M = 10), then we switch to a piston-like movement of the non-wetting phase, which invades a large part of the pore space with fingers (i.e., elongated tail structures, e.g., (Ca 1e−2,M 10) in Fig. 3k), establishing a steady state of increased saturation after breakthrough.
Disconnected Components (Fig. 7c) We now consider the number of disconnected components (i.e., blobs and ganglia) between entry of the non-wetting phase and breakthrough. This gives us a complementary view of the behaviour described my means of Fig. 7a. In this diagram, the experiments (Ca 1e−2, M 1) and (Ca = 1e−3, M = 0.2) stand out, whose number is more than ×3 above the values of the other experiments. Accordingly, it shows that the combination of high capillary number and low viscosity ratio creates many more isolated and disconnected blobs of the non-wetting phase than the other configurations. From these results, we would expect that for the configuration (Ca = 1e−2, M = 0.2) we would see that the number of blobs would be even higher, since the high capillary number and the low viscosity ratio would favour such a behaviour. This would have to to be confirmed with additional experimental measurements. It is also clearly reflected that as the viscosity ratio increases for high capillary numbers, the system switches from creating blobs to creating thin fingers. This leads to a corresponding increase of the interfacial area (as seen in Fig. 7a), but to a lesser degree since there is also a different mechanism creating this interfacial area.

Space-Time Analysis via Transport Networks
The analysis in Sec. 5 is based on values of interest to reflect the state of the experiment at a given time. While these values quantitatively show how Ca and M affect the displacement processes, they do not convey how this impacts the resulting flow patterns and how they differ spatially between experiments. For this, we introduce a compact representation based on node-link diagrams of so-called transport networks, which aggregate the data over time while maintaining spatial flow information. The porous medium is modeled as a network of pore bodies (nodes) connected by pore throats (links). This concept is widely used in the estimation of material properties and the simulation of different processes in porous media research; however, no precise definition of pore bodies and throats exists [Blu01,XBJ16]. We developed an approach based on the triangulation of the whole domain to not only extract topological information but also capture covered pore and throat areas.

Transport Network Visualization
We first create a graph capturing the pore network structure w.r.t. the solid (1). In doing this, we associate each node and edge (i.e., pore and pore throat) with a certain region in the original data. This then allows to determine the fluid coverage for each edge and node across experiments and time steps (2).
(1) Extraction of porous structure network graph. Transport networks and their topology are based on a triangulation of the void space, i.e., the part of the porous medium that is available for fluid transport (cf. Fig. 8a). For this, we consider the extracted mask interfaces (cf. Fig. 4.1). Each contour is described by a clockwise ordered list of coordinates of the associated pixels which are sampled by selecting indices with a fixed step size (blue points in Fig. 8a). In our application, a step size of 20 pixels has proven to be sufficient to adequately resolve the domain. From this set of points, an unstructured triangle grid is generated by computing a Delaunay triangulation.
For each triangle edge, we then consider its midpoint and check whether it lies in void space. For this we use the preprocessed, binary mask image of the micro-model, as described in Sec. 3.2. To compensate for rounding errors due to converting the calculated floating point positions of the midpoints into integer pixel coordinates, the midpoints are checked against a version of the mask image where the solid phase is dilated by one pixel. If the midpoint lies in void space it becomes an intermediate node (yellow circle in Fig. 8a) and is added to our transport network. We then consider the triangles again and check for the numbers of intermediate nodes at their edges. In case there are two intermediate nodes, we add an intermediate edge connecting them to our transport network (yellow line in Fig. 8a). If there are three intermediate nodes, a new pore body node is inserted at the barycenter of the triangle (red circle in Fig. 8a). In addition, connections of the three intermediate nodes to the pore body node are inserted as well. Thus, a pore throat connecting two pore bodies is formed by at least two intermediate edges.  (2) Checking for fluid progression. Each pore body node and intermediate node is associated with one triangle each which again represents a portion of the porous medium that is available for fluid transport. Each of these elements is indexed via n below. A high coverage of the area of a triangle by the transported fluid indicates that the corresponding element of the network is used for fluid transport. The normalized coverage Cn,t = An,t A∆ n is computed as the ratio between triangle area An,t and the area of the full triangle A ∆n Based on this, we now associate a usage rate U with each intermediate edge and pore body node of the network. We define this measure Un as the number of time steps t at which the normalized coverage of the triangle of element n is larger than a fixed threshold k, divided by the total number of time steps of the investigated time period T : In our analysis we set a comparably low threshold of k = 0.1 to exclude small dust particles in the micro-model from being detected as transported fluid. However, the sensitivity w.r.t. k is rather low, and our measurements showed only little impact for k ∈ [0.1, 0.7].
(3) Visual representation. For visualizing the transport network, intermediate line segments between two pore bodies are contracted to form a single edge connecting pore bodies (red line in Fig. 8a). Their thickness is scaled linearly with the minimum width of the represented pore throat. We approximate this width by the minimum length of corresponding triangle edges (orange in Fig. 8a). The usage rate U of each segment is mapped onto the contracted edge by orthogonally projecting its intermediate nodes onto the single edge and assigning the value to this segment (projection indicated by purple arrows). Finally, the contracted edge is drawn segment by segment, mapping the usage rate to color using the Inferno color map [SvdWF15]). Complementing U, an important characteristic of interest for the analysis is the length of the considered time interval. To incorporate this information, we use unoccupied pores and throats with U < 0.01 (i.e, that would otherwise be mapped to black in our case). For this, we use a secondary isoluminant color map based on CET-I1 by Kovesi [Kov15] (which we consider to be rather low-key visually compared to Inferno).
An aspect of particular interest is how the flow paths change after the breakthrough, but to clearly show that the transport network in Fig. 9 is created in a slightly modified way. In the first step, the usage rate is calculated for the complete time interval after breakthrough until the end of the recording. In the second step, these results are compared to the usage rate values before breakthrough and mapped to black for any element that was used for fluid transport before breakthrough (U > 1%).
The generated void space triangle mesh, nodes and connections can be applied to all experiments that have been conducted with the same micro-model, after they have been aligned as described in Sec. 4. This allows to analyze and compare the experiments with a uniform representation, clearly conveying differences in flow processes.

Transport Network Matrix
For investigating the effects of Ca and M on flow behaviour, we create a small multiples visualization that arranges transport networks such that it reflects their position in the parameter space. Showing all available combinations of capillary numbers and viscosity ratios allows to directly compare their impact on the spatio-temporal distribution of phases in a static visualization. In our analysis, we investigate two kinds of transport matrices, differing in the considered time span: before (Fig. 1) and after (Fig. 9) breakthrough. Their combination conveys a comprehensive qualitative and quantitative summary of the evolution of flow, and also allows to identify isolated effects and evaluate their contribution to phase distribution.
First, we consider the time span from entry to breakthrough in Fig. 1. Generally, the progression of the flow through the porous medium can clearly be seen, especially also indicating for how long certain structures were occupied. It is of particular interest to see that late in the experiment new flow paths are covered close to the inlet of the model (as indicated via low usage rate U), with a residence time similar to the occupied areas close to the outlet. This means that these effects occur later than could be expected, which happens to some extent in all experiments, but is particularly significant for larger Ca and M. This provides an indication of the importance of the phase occupancy at the entry pores of the flow network in these cases especially, and how they contribute to the breakthrough saturation. We further clearly see that for low capillary numbers (Ca = 1e−5), the viscosity ratio M does only have little impact on the distribution of phases before breakthrough, as they show approximately the same spatial distribution.
From Fig. 9, it becomes apparent in which experiments the flow paths still change significantly after breakthrough. This behavior indicates a flow dominated by viscous rather than capillary forces, and allows us to partition the parameter space into different regimes. For low capillary numbers Ca, we see that after breakthrough the usage rate does not change significantly since the capillary forces dominate the process, and the geometry of the porous medium is in charge of the process. The minor changes which take place can mostly be attributed to so-called capillary end effects, which would not be of essence in a bigger porous medium (capillary end effects arise from the discontinuity of capillarity in the wetting phase at the outlet). For intermediate capillary numbers, we observe some similarities with the previous case, in the sense that the distribution of phases changes so as to increase the breakthrough saturation (but not drastically). However, there is a stronger tendency toward newly visited pores close to the inlet. This is the outcome of the effect of the more dynamic boundary conditions in the network of channels outside of the field of view, which were used to feed the porous medium with the non-wetting phase. We also discover significant mobilization of the non-wetting phase at later times which increased the final saturation. We attribute this effect to the fact that the non-wetting phase has reached the outlet of the model and re-mobilization was inevitable due to extensive viscous effects. As the capillary number is increased even further to Ca = 1e−2, we see the breakthrough saturation to increase accordingly, with extensive re-mobilization of the non-wetting phase due to the increased viscous effects. Crucially, the re-mobilization of the non-wetting phase covers (nearly) the whole pore space, which is an important property for numerous applications. Another interesting aspect we learn from the visualization is the effect of the phase occupancy at the boundary pores. The degree of occupancy mostly depends on the boundary conditions, but the visualization shows that for high capillary numbers-i.e. in a viscous flow regime-the breakthrough saturation and the phase distribution are both highly affected by the fact that the occupancy of boundary pores varies significantly. This can impact the comparability of the experiments. Resolving this issue by adapting the experimental setup is almost impossible in practice since for high capillary numbers even a small deviation in terms of surface roughness or pump fluctuations can have a significant effect. This emphasizes the importance of showing such deviations as a crucial part of a comprehensive analysis.

Pore Throat Histogram
Abstracting the porous medium as a network allows us to characterize the displacement process by discrete paths through which the non-wetting phase flows. It further crucially allows to create histograms that make a connection between pore throat width and fluid transport. In the presentation of the histograms, we differentiate between pore throats that already transported fluid before the breakthrough and those that only do so after the breakthrough (akin to the analysis in Sec. 6.2). The two types are visualized as stacked bars with blue and orange color, respectively. The remaining pore throats, which do not transport fluid, are added on top as light-gray bars to convey the total distribution of pore throat widths in the porous medium.
These histograms in Fig. 10 complement the transport network visualization in Fig. 1 and Fig. 9, covering both time periods (before/after breakthrough) in one view and allowing for a more quantitative comparison of throat occupation. Different aspects already discovered in the transport network visualizations can be investigated in more detail. Generally, larger Ca lead to a higher saturation, with Ca ∈ {1e−3, 1e−2} yielding (almost) full saturation of the porous medium eventually. This complies with the general understanding that high capillary numbers in combination with low viscosity ratios lead to viscous fingering, low capillary numbers with practically any viscosity ratio lead to capillary fingering, and high capillary numbers in combination with high viscosity ratios lead to a stable front. Here, it can be seen that the viscosity parameter M crucially impacts when parts of the medium are saturated: with lower M large portions of the domain are only saturated after breakthrough. It is also confirmed that-across all experiments-a lower coverage is achieved with smaller throat widths, due to the increased entry capillary pressure of the corresponding throats in a drainage process.

Discussion and Future Work
In this work, we conduct an in-depth visual analysis of the characteristics of two-phase fluid flow and its preferred transport paths based on new experimental data. In particular, we study the characteristic flow patterns resulting from different combinations of viscous and capillary forces, and the impact of the porous medium structure. For this, we extract quantities of interest from the experiments with the high accuracy required, and employ breakthrough as a reference point to achieve comparability between processes running at different time scales across experiments. This enables quantitative analysis, but abstracts spatial effects. To additionally consider this aspect, we introduce a new spatio-temporal visual representation to depict which paths through the porous medium were occupied to what extent. Apart from visualizing it directly, it can also be useful to extract abstract information regarding spatial aspects, namely the occupancy of throats depending on their width.
Identifying and classifying experiments with similar behavior on this basis allows us to partition the experiment parameter space into regimes and relate it to the dominant driving forces. An important insight gained from the visual analysis is that we can now see that the reason for the creation of the interfacial area is completely different in two characteristic cases. One is due to the breakup of the non-wetting phase due to the competition between the viscous and capillary forces, while the other one is due to the creation of thin fingers. Based on our conducted experiments, we can now not only quantify the two different effects, but we can also spatially correlate them, and map them to the pore space using our spatio-temporal representation. There are significant implications driven by this ability to differentiate between degenerative identical states in a robust, compact and efficient way, since these findings can be attributed to different effects, which correspond to different approaches regarding the physical context.
In general, our experiments exhibit physical effects that are challenging to model due to the increased physical complexity of the processes involved. While micro-scale approaches exist, they cannot adequately model all observed effects, whereas macro-scale models do not consider some effects at all (e.g., the reduction of effective interfacial area between solid and non-wetting fluid not appropriately reflected). Based on the extracted data and findings of this work, we plan to make a first attempt to effectively model the respective processes. We aim to investigate the role of the solid phase geometry more closely, especially in the interplay with the boundary and physical conditions of the fluids involved. We also intend to extend our research to scenarios further than primary drainage, and evaluate the role of the pore bodies and throats of the porous structure as the dominant geometrical measure of each displacement process. This is directly linked to the Young-Laplace equation, where capillary pressure is dominantly affected by the local geometrical properties of the porous medium. Pore bodies are the determining feature for drainage, due to their increased size and lower entry capillary pressure for an invasion process, while pore throats are dominant during imbibition due to their reduced size and the spontaneous response of the system to the increased capillarity. Additionally, the role of the disconnected phase will be investigated more closely, both in terms of their effective influence in the evolution of flow as a pressure barrier, but also in terms of the local, boundary, and physical conditions under which a phase can become disconnected.
These efforts will require a much more densely sampled experiment parameter space (via further experiments or conducted simulations), posing additional challenges for the analysis as well. In particular, our visualization approaches should also be made more scalable visually, and we aim to consider different directions in this regard. First of all, while our transport networks proved to be crucial in our analysis, they cannot directly scale visually to larger experiment and/or simulation ensembles. However, they already provide a simplified representation of the pore structure, that could further be scaled towards more expressive smaller representations (like glyphs) by reducing the underlying graph accordingly. In addition, we aim to conduct further extraction of other quantities that could be relevant (e.g., distinguishing between ganglia and blobs and quantifying their shape via descriptors). This accordingly results in a high-dimensional feature vector describing each experiment, that could further be the basis for further analysis techniques based on projection and/or clustering.
Furthermore, in this work, we employ the same porous media model for the sake of direct comparability. A larger, more comprehensive study, however, could include different models for a more targeted investigation of specific effects. This would be supported by our full pipeline: the porous model for the experiments was generated using AutoCAD (and could easily be adapted) and all other steps, including the transport network visualization, could be done in exactly the same way. However, naturally, direct spatial comparison across different pore structures is not possible directly anymore, but needs to be substituted by more abstract alternatives. Perspectively, 3D experiments (and simulations) could also be conducted for further investigation. The extraction of quantities generally works akin to the 2D case, and the transport networks can be generalized to 3D as well: based on the discretization of the porous medium by a tetrahedral mesh, triangle edge midpoints are simply replaced with the barycenters of the triangular tetrahedron faces, and midpoint checks can also be adapted accordingly in a straight-forward fashion.