Relating Pore‐Scale Observations to Continuum‐Scale Models: Impact of Ganglion Dynamics on Flow Transport Kinetics

Continuum‐scale models used to model and predict two‐phase flow in the subsurface are often based on averaged flow parameters and do not consider pore‐scale fluid flow phenomena, for example, ganglion dynamics and thin‐film flow. As such, a major challenge in upscaling two‐phase flow for groundwater engineering applications is understanding the impact of disconnected flow and ganglion dynamics on continuum‐scale flow functions such as relative permeability‐saturation and capillary pressure‐saturation curves. In this study, we explored how changes in wettability and fluid velocity affect ganglion dynamics. We conducted pore‐scale numerical simulations with OpenFOAM to investigate the displacement of decane by water. Additionally, we examined how displaced phase saturation (a continuum‐scale flow function) responds to changes in dynamic fluid connectivity. We identified three different fluid flow regimes, that is, the connected pathway flow regime, ganglion dynamics (GD) flow regime, and droplet traffic flow regime, and studied the effects of changes in the wettability of the porous medium and the velocity of the invading fluid on the transitions between these different regimes. Our research showed that transitions between connected and disconnected pore‐scale flow regimes, which are induced by changes in fluid velocity and wettability, have a significant impact on both fluid displacement efficiency and average fluid flow transport kinetics.


Introduction
The process of immiscible displacement, where one fluid phase is replaced by another in a porous medium, is a key feature of many subsurface engineering processes.This includes activities such as enhanced oil recovery, groundwater remediation, and geological storage of carbon dioxide (CO 2 ) (Singh, Scholl, et al., 2017;Wolf et al., 2020).During these processes, fluid ganglia, which are small or large clusters of individual fluid, become disconnected from the bulk fluid.This leads to the trapping of fluid phases (Singh, Menke, et al., 2017).
Since the work by Richards (1931), it has long been believed that immiscible fluid flow occurs only through connected pathways and that disconnected fluid ganglia are stranded and immobile (Payatakes et al., 1996).However, Avraam and Payatakes (1995) conducted experiments in two-dimensional microfluidic glass micromodels and identified three different pore-scale fluid flow regimes: the connected pathway flow (CPF) regime, the ganglion dynamics (GD) flow regime, and the droplet traffic flow (DTF) regime (Armstrong et al., 2017;Avraam & Payatakes, 1995).
In the CPF regime, as the name suggests, the fluid flows via connected pathways, while the GD flow regime involves periodic disconnection, reconnection, and mobilization of large, disconnected fluid clusters.These events are collectively termed "ganglion dynamics."Similarly, the DTF regime is also a disconnected flow regime; however, the size of the mobile fluid clusters is much smaller than that in the GD regime (minute fluid droplets instead of large clusters) (Picchi & Battiato, 2018).These flow regimes have been observed not only during flow in custom-made microfluidic pore networks but also in real porous media (Valavanides, 2018).In recent years, they have been observed and studied in rock cores using microcomputed tomography imaging techniques (Armstrong et al., 2016;Rücker et al., 2015;Schlüter et al., 2016;Spurin et al., 2021) and numerical simulation methods (Armstrong et al., 2016;Nhunduru et al., 2022;Zhang et al., 2022).
There are still uncertainties concerning the conditions under which a given flow regime (CPF, GD, or DTF) prevails (Alhammadi et al., 2017) and how changes in the balance of forces influence transitions between regimes (Spurin et al., 2019).Understanding disconnected flow, ganglion dynamics and other complex fluid displacement phenomena in subsurface geological formations such as aquifers is of paramount importance.Such an understanding has significant implications for underground water resources and subsurface engineering operations, for example, remediation of nonaqueous liquid contaminants and pollutants in aquifers and the injection and storage of CO 2 in deep saline aquifers.
Most of the world's groundwater is contained and stored in underground aquifers (NASA Jet Propulsion Laboratory, 2023).Thus, it is critical to understand complex fluid flow phenomena such as fluid fragmentation and ganglion dynamics to ensure effective management of water resources.In the remediation of nonaqueous phase liquids (NAPLs) and oil spills from groundwater, disconnection and fragmentation of fluid ganglia can help to mobilize the trapped nonaqueous phase and reduce NAPL residual saturation (Pak et al., 2015(Pak et al., , 2020)).Furthermore, fluid fragmentation increases the fluid fluid contact area, which in turn enhances the effectiveness of added surfactants and accelerates the rate at which microbial agents and/or inorganic reagents breakdown the oil/ NAPL (Pak et al., 2015(Pak et al., , 2020)).
The dominance of the CPF, GD, or DTF regime and fluid fragmentation at the pore scale depends on several factors, including the physical and structural properties of the porous medium, such as the pore geometry, pore connectivity, and wettability of the porous medium.Fluid properties such as flow velocity and fluid viscosity contrast (Jahanbakhsh et al., 2021) also have an impact on fluid fragmentation and pore-scale fluid flow regimes.Avraam and Payatakes (1999) studied steady-state two-phase flow in a glass model pore network and observed different pore-level mechanisms.Their findings revealed that strong wettability leads to increased fluid fragmentation and ganglion dynamics at the pore scale.These researchers also revealed that large viscosity ratios promote the formation of smaller ganglia and oil droplets.In another study, Spurin et al. (2019) investigated porescale flow regimes in a water-wet Estaillades rock sample.These scholars studied a range of capillary numbers and fluid mobility ratios and showed that fluid intermittency depends on the capillary number (Ca) and mobility ratio (M).The mobility ratio (M) is the ratio of the dynamic viscosity of the invading fluid to that of the defending fluid.They reported fluid intermittency for M < 1 at low capillary numbers and CPF for M > 1.
In Valavanides (2018) conducted simulations to investigate oil fragmentation under intermediate-wet conditions by testing a range of flow velocities and mobility ratios.Their findings suggest that the flow of the nonwetting phase in disconnected form is significant, and in certain cases of flow conditions, this is the dominant flow mode.These investigators also identified systematic flow structure mutations that occur with changing flow conditions and reported that some of these mutations can surface on a macroscopic scale and be measured, for example, through a reduced pressure gradient, while others remain latent within the interstitial flow structure.Deeper within the disconnected oil flow, mutations between ganglion dynamics and drop traffic flow prevail.The flow mutations become more pronounced with increasing fluid viscosity contrasts.
It is commonly understood that fluid displacement in porous media can undergo a transition from connected to disconnected flow.However, there is limited knowledge regarding the sensitivity of this transient process to changes in pore-scale flow and transport properties such as wettability and fluid viscosity contrasts.This lack of understanding also extends to the subsequent effects on average flow functions such as relative permeability, saturation function, capillary pressure, etc. (Datta et al., 1994;McClure et al., 2022).
Relative permeability-saturation curves are one of the primary functions used to describe the dynamic behavior of multiphase flow systems at the continuum scale.In such instances, relative permeability is tabulated as a nonlinear function of phase saturation through basic parameterization.Relative permeability-saturation curves provide an average representation of immiscible fluid flow behavior.These curves are widely used for predicting continuum-scale flow behavior for subsurface applications such as improved oil recovery, carbon dioxide storage, and contaminant transport (Gharbi & Blunt, 2012).Most continuum-scale relative permeability models assume that during immiscible fluid displacement, each fluid flows through a single connected pathway in the porous medium.However, as the fluid displacement process progresses and fluid volume fractions change, connected pathways break apart, and fluid mobilization occurs.This results in fluctuations in fluid topology and meniscus configuration.Dynamic changes in fluid connectivity can cause fluctuations in pressure and energy within a given fluid system during fluid displacement processes at steady state.Consequently, this can alter the fluid relative permeability, as reported by several authors (Armstrong et al., 2016;Khanamiri & Torsaeter, 2018;Ott et al., 2020;Picchi & Battiato, 2018, 2019;Rücker et al., 2019).
Authors, Rücker et al. (2019) conducted a study to image the pore-scale flow regime in a carbonate that was altered to a mixed-wet condition by aging with crude oil to represent the natural configuration in an oil reservoir.These authors used fast synchrotron-based X-ray computed tomography to directly image the pore-scale flow regimes and observed ganglion dynamics, wherein the pore space is intermittently filled with oil and brine.They found that the frequency and size of these intermittent fluctuations are greater than those in water-wet rock and that the impact of ganglion dynamics on the overall flow and relative permeability cannot be neglected in modeling approaches.
In another study, Ott et al. (2020) studied the displacement of crude oil by alkaline water injection in a microfluidic environment and revealed that changes in injection-water chemistry (wetting properties) and fluid connectivity have a considerable impact on the relative permeability.The results of their study demonstrate that statistical and topological information can be used to analyze displacement physics.Other authors have also used topological information to analyze displacement physics (Armstrong et al., 2016;Khanamiri & Torsaeter, 2018).Moreover, in multiphase flow studies (Armstrong et al., 2016(Armstrong et al., , 2019;;Herring et al., 2013Herring et al., , 2015;;McClure et al., 2019McClure et al., , 2020)), the Euler characteristic, χ, is commonly used to describe fluid topology and configuration.The Euler characteristic is defined as follows: where β 0 is the number of isolated clusters (zeroth Betti number), β 1 represents the number of loops (first Betti number), and β 2 is the number of cavities enclosed within a given system (second Betti number).The Euler characteristic, χ, defines topological changes by quantifying the evolution of the number of isolated clusters, loops and cavities as the displacement process progresses.
Cluster disconnection events result in an increase in the number of isolated clusters (β 0 ) and thus an increase in the Euler characteristic χ.Cluster coalescence leads to a decrease in the number of isolated clusters (β 0 ) and hence a decrease in the Euler characteristic χ.Cluster disconnection and reconnection events culminate in changes in the connectivity of the fluid system.
In Khanamiri and Torsaeter (2018) conducted drainage and imbibition flow experiments on water-wet Berea sandstone.These researchers studied and imaged the pore-scale fluid topology in a Berea sandstone rock core using synchrotron X-ray microtomography.During their research, they observed both connected and disconnected flows, as well as dynamic ganglion behavior.The results of their study showed that the Euler characteristic of disconnected ganglia combined with the connected pathway versus saturation exhibited hysteresis in imbibition and drainage.In contrast, the Euler characteristic of only the connected pathway versus saturation did not show any significant hysteresis in imbibition and drainage.Furthermore, it was reported that disconnected ganglia contributed to the flow by creating internal redundant connections, volume exchange with the connected pathway, and saturation changes.Their results provide evidence that the contribution of disconnected flow and ganglion dynamics to average (continuum-scale) flow functions should not be overlooked.
In another study, Armstrong et al. (2016) used fast X-ray microtomography to observe pore-scale two-fluid configurations during immiscible flow and initialize lattice Boltzmann simulations.According to their simulations, the mobilization of disconnected nonwetting phase clusters can account for a significant fraction of the total flux.These researchers also reported that fluid topology can undergo substantial changes during flow at constant saturation, causing hysteretic behavior.They concluded that traditional assumptions regarding fluid configurations are oversimplified and that the role of fluid connectivity cannot be ignored for multiphase flow.At the continuum scale, fluid topology and phase connectivity are accounted for by interfacial area and Euler characteristics.These are parameters that are currently missing from our current models.The results of their work suggest that average topological properties may be sufficient for predicting relative permeability, even for cases where fluid configurations are continuously changing at the microscale.
In Picchi and Battiato (2018) presented analytical expressions for the relative permeability of wetting and nonwetting phases in different flow regimes.These scholars demonstrated that changes in topology and dynamic connectivity of the flowing phases have a significant impact on relative permeabilities.Their research revealed that relative permeabilities for the CPF regime depend only on saturation.The study also highlighted that twophase model predictive capabilities are limited when flow regimes are not accounted for.Picchi and Battiato (2018) proposed a new homogenization framework that enables the derivation of regimespecific upscaled equations.The key advantage of this approach is that it incorporates the interactions between the flowing phases and the solid while maintaining a simple final upscaled formulation.
Specifically, the complexity of the problem is reduced by proposing an analogy between the flow regimes in real porous media and in a capillary tube, where the effective parameters of the macroscopic equations are determined by postulating the pore-scale flow regimes (Picchi & Battiato, 2019).Using this approach, the geometry of the flowing phases is simplified, however, the pore-scale physics of the interaction between different phases is retained, and the relative permeabilities of both the wetting and nonwetting phases can be analytically derived for the CPF and small and large ganglion dynamics flow regimes (Picchi & Battiato, 2019).Although this approach has proven successful in quantifying the impact of flow regimes on relative permeabilities, its use is limited to highly idealized configurations, and its ability to predict relative permeabilities in complex porous media is questionable.Importantly, the question of whether any type of correction needs to be introduced in the context of realistic porous media remains open (Picchi & Battiato, 2019).
In this work, the combined effects of wettability and fluid velocity on flow transport kinetics, the transition between pore-scale flow regimes and fluid topological behavior are studied and evaluated.To investigate and evaluate the sensitivity of the transition between connected and disconnected flow regimes to changes in wettability and invading fluid velocity, pore-scale numerical simulations are performed in OpenFOAM.The connectivity index (λ) is used to monitor and track temporal changes in fluid topology and dynamic fluid connectivity.The dynamic behavioral response of the saturation function (a continuum-scale function) to wettability and flow rate-induced changes in dynamic fluid connectivity are examined.In addition, the observed dynamic response is related to existing relative permeability-saturation models.

Numerical Simulations
Pore-scale simulations were conducted using the open-source computational fluid dynamics software Open-FOAM version 6 (OpenFOAM 6, 2018) to study the displacement of a decane oil phase by water in a heterogeneous porous matrix with a porosity of 32%.The porous matrix, depicted in Figure 1, is a small-scale representation of the microstructure of a Berea sandstone rock sample.
The microstructure of Berea sandstone was created by replicating an X-ray microtomography image of a Berea sandstone rock, as detailed by Boek & Venturoli (2010) and Wlodarczyk et al. (2019).The microstructure was engineered at Schlumberger Cambridge Research using a thin section of a 3D Berea sandstone rock sample that was imaged using X-ray microtomography.
The matrix was discretized on a lattice and converted into a bitmap format, where 0 represented pores and 1 represented solid grains or obstacles.This bitmap was then transformed into the stereolithography (STL) format, which is used for flow simulations in OpenFOAM software.The dimensions of the digital simulation model were scaled down by a factor of 10 from the physical dimensions documented by Wlodarczyk et al. (2019).The 2 cm × 0.91 cm patterned area of the microfluidic model was scaled down to 0.2 cm × 0.091 cm.This resulted in a complete miniature image of the original model.
The pore-scale numerical simulations were conducted using the volume of fluid (VOF) method in OpenFOAM.The transient, two-phase flow solver InterFoam for incompressible, isothermal fluids (Datta et al., 1994;Gharbi & Blunt, 2012;McClure et al., 2022) was used to discretize the numerical solution.The InterFoam solver was used because it is built into the OpenFOAM package, and its ability to simulate two-phase flow displacement problems has been proven in many previous studies (Deshpande et al., 2012;Friedemann et al., 2021).Details of the numerical method and InterFoam solver are discussed in previous work by Nhunduru et al. (2022).For the mesh used for the simulations, the porous microstructure was composed of 104,587 cells and 224,302 points.The mesh was generated using the blockMesh and snappyHexMesh utilities in OpenFOAM.

Mesh Validation
To ensure that all observed phenomena were a result of the simulation conditions under study and not due to spatial discretization errors, a mesh independence study was conducted.The chosen mesh density of 104,587 cells and 224,302 cells was refined by factors of 0.75 and 1.5, 2, and 4 respectively.The changes in the residual values for the four levels of mesh refinement are shown in Figure 2.
The results showed that the chosen mesh density (refinement factor of 1) was optimal and appropriate for the analysis.This was confirmed by the residual plot (Figure 2), which revealed that the solution converged for all the   75, 1, 1.5, 2, and 4.

Water Resources Research
10.1029/2023WR035624 refinement factors considered.As shown in Figure 2, increasing the mesh density up to a factor of 4 had no significant effect on the numerical solution.This indicates the independence of the solution values from mesh resolution.

Simulation Validation
To validate the ability of the numerical method to accurately model and simulate natural fluid flow phenomena in porous media, an experiment involving the displacement of decane by water was conducted in the laser-fabricated microfluidic device (micromodel) shown in Figure 1b.
In the experiment, the micromodel was saturated with decane, followed by subsequent water injection.The micromodel was oil-wet (contact angle of 135°), and the experiment was conducted with a fluid inlet velocity of 0.03 m/s at ambient temperature and pressure.The microfluidic experiment was then modeled using the VOF method in OpenFOAM for comparison.The simulation results were compared to and validated against the experimental results (Figure 3).
Figure 3 shows that the experimental and simulation results exhibit a good dynamic flow match.This proves that the numerical method used is both valid and accurate for modeling fluid flow in porous media.
There is, however, some variability observed between the times at which the fluid displacement patterns in the microfluidic model and the simulation match, as shown in Figure 3.This is because the physical model has surface heterogeneities, such as surface roughness and random variations in channel depth, due to the nature of the laser fabrication process, as described by Wlodarczyk et al. (2019).These surface heterogeneities did not exist in the simulation model, which had ideal smooth surfaces.

Investigating the Effect of Wettability and Invading Fluid Flow Rate on Pore-Scale Fluid Displacement
Several simulations were performed to study how the wettability and the flow rate of invading fluid affect fluid displacement at the pore scale.The simulations were conducted in the digital Berea sandstone model, shown in Figure 1a.In the simulations, water was used to displace a decane oil phase from the porous matrix.At the beginning of the simulations, the Berea sandstone matrix was fully saturated with decane, and water was injected from left to right to displace the decane.The properties of the water and decane phases used in the simulations were obtained from the literature (Dymond & Oye, 1994;Sayed et al., 2019;The Engineering Toolbox, 2004) and are summarized in Table 1.Viscous fingering is a fluid flow phenomenon that occurs when the dynamic viscosity of the invading fluid is lower than that of the defending fluid, that is, μ 1 /μ 2 < 1.When viscous fingering occurs, the invading fluid preferentially flows through the pore network (Zheng et al., 2017).In this work, the value of μ 1 /μ 2, that is, the viscosity ratio (M), was greater than 1; therefore, no viscous fingering occurred in any of the simulations.
Simulations were performed to investigate the effects of three different wetting conditions (contact angles) on pore-scale fluid displacement.These conditions were moderately water-wet (contact angle of 45°), neutrally wet (contact angle of 90°), and strongly oil-wet (contact angle of 150°).Additionally, four different water velocities were tested, and their values listed in Table 2.In total, 12 scenarios were simulated.
The capillary number is defined by the expression in Equation 2 as follows: where ρ inv is the invading fluid density, d H is the hydraulic diameter at the inlet, v inv is the invading phase velocity at the inlet, μ inv is the dynamic viscosity of the invading phase, and σ is the interfacial tension.
In the applicational domain, the chosen capillary numbers represent near-wellbore scenarios where the capillary numbers and fluid velocities are high.In all the simulations, the same boundary conditions were used.The walls were set to a no-slip condition.At the outlet, a zero-gradient velocity condition was applied, which means that the velocity of the fluid at the boundary is the same as the velocity of the fluid inside the simulation domain.A zeropressure boundary condition was also used at the outlet.The only differences between the simulations were the wettability (contact angle at the fluid-solid interface) and the inlet fluid velocity.

Measurement of Fluid Topology and Dynamic Fluid Connectivity
The Euler characteristic, χ, defined previously in Equation 1, is often used to describe fluid topology and connectivity.One commonly encountered challenge when using this characteristic is that due to edge effects, the Euler characteristic of a given porous specimen provides a biased estimate of the Euler characteristic of the region from which the specimen was taken (Odgaard & Gundersen, 1992).As such, in order to measure the changes in connectivity within the Berea sandstone microstructure during fluid displacement processes, a special connectivity index λ, (Domander et al., 2021) was used.This index was derived from the Euler characteristic χ as described in Equation 3. The measurements were obtained by using the "Bone J" plugin (Doube, 2015(Doube, , 2020) ) in ImageJ for a sequence of images extracted from the simulations.
The Bone J plugin employs the "purify" method for connectivity analysis.In this method, it is assumed that there is only one foreground particle and no cavities to remove small foreground and background particles.Bone J then iterates through the stack of images, calculating the Euler characteristic for Furthermore, Bone J calculates the contribution of disconnected elements to the Euler characteristic of the structure it was originally connected to (Δχ) by checking the intersections of voxels and stack edges.The connectivity is then calculated as follows: The Bone J plugin also determines the connectivity density (Conn.D), which can be interpreted as the number of trabecule per unit volume in mm (Doube, 2015(Doube, , 2020;;Singh, Menke, et al., 2017): The term ∆χ (Domander et al., 2021) corrects for the change in the topology of a system when it is broken down into smaller elements.

Pore-Scale Flow Regimes
The simulations revealed three distinct fluid flow regimes, namely, DTF, ganglion dynamics (GD), and CPF, as illustrated in Figure 4.The DTF and GD regimes were identified as disconnected flow regimes, while the CPF regime was considered to be a connected flow regime.
The flow regimes were classified based on the connectivity of the invading phase.The DTF regime refers to the situation where the invading fluid (shown in dark blue in Figure 4) flows as multiple disconnected fluid elements, each with a size comparable to the diameter of the pore throat.In this regime, the fluid elements are densely distributed within the porous network and periodically connect and disconnect during the fluid displacement process.In the ganglion dynamics (GD) flow regime, the invading fluid also flows as multiple disconnected fluid elements; however, these fluid elements are larger than the diameters of the pore throats and are sparsely distributed within the pore network.As in the DTF regime, the fluid elements periodically connect and disconnect during the fluid displacement process.
The DTF regime is characterized by more disconnected fluid elements than the GD regime.However, in the DTF regime, these fluid elements are closely packed together, resulting in highly cohesive flow.Consequently, the DTF regime has a high fluid displacement efficiency of approximately 60%-80%, similar to that of the CPF regime (see Figures 4a and 4c).In the GD regime, the invading fluid flows as sparsely distributed, disconnected fluid elements, as shown in Figure 4b.This causes the flow of the invading fluid to lack cohesion, leading to incomplete and inefficient displacement of the defending fluid from the pores.In this scenario, the displacement efficiency is much lower (40%-45%).
The CPF regime, on the other hand, is characterized by invading fluid flowing as a single pathway or fluid element that remains connected and does not disconnect during the fluid displacement process.Consequently, the CPF regime exhibits high cohesion and high displacement efficiency (over 70%).

Impact of Wettability and Invading Fluid Velocity on Pore-Scale Flow Regimes and Fluid Invasion Patterns
To understand the impact of wettability and fluid velocity on fluid invasion patterns, a comparison between the fluid displacement patterns for the 12 simulations was made, as shown in Figure 5.The comparison was carried out at a similar dimensionless time of τ = 0.8.The dimensionless time, τ, is defined by the following expression:   , 1995;Picchi & Battiato, 2018;Valavanides, 2018).The ganglion behavior varied significantly between the θ = 45°and 150°c ases.This difference was due to variations in the wetting properties.
In the 150°simulations, the porous medium was oil-wet, and the walls of the porous medium repelled the invading water phase.Consequently, the invading water phase did not wet the surface of the porous medium completely.In the 45°simulation cases, on the other hand, the porous medium was water-wet, and the walls of the porous medium strongly attracted the invading water phase.This resulted in complete wetting of the surface of the porous medium during the fluid displacement process by the invading water phase.Prior to fluid propagation via the GD regime, fluid displacement in the log Ca = 4 and log Ca = 3 simulation cases for θ = 45°and θ = 150°o ccurred via the connected flow pathway (CPF) regime in the initial stages of the displacement processes.The onset of fluid disconnection and the first GD event triggered an avalanche of similar events and a sudden acceleration in the flow transport kinetics.

Impact of Flow Transitions on Fluid Dynamic Connectivity and Flow Transport Kinetics
For simulation cases in which a transition from the CPF regime to the GD regime occurred, the onset points of fluid disconnection and GD events are marked "GD" in Figure 6. Figure 6 shows the temporal evolution of saturation and connectivity (λ) for the investigated cases.Figures 6f, 6g, 6i, and 6l indicate that the occurrence of fluid disconnection and dynamic ganglion (GD) events leads to a sudden acceleration in the flow kinetics.This results in a sudden increase in the gradients of the saturation and connectivity (λ) curves.The saturation curves are shown in gray for all simulation cases in Figure 6, and the connectivity (λ) time evolution curves are shown in blue for θ = 45°, green for θ = 90°, and red for θ = 150°.In all the simulation cases investigated, the microstructure of Berea sandstone was initially fully saturated with decane (initial decane saturation = 100% and decane connectivity index = 34.5).As the water phase invaded the pores of the porous matrix, the decane saturation decreased, and the defending decane phase became more disconnected (i.e., λ decreased).
Disconnection and break-up of the defending phase resulted in more isolated clusters, that is, an increase in β 0 in Equation 1, leading to an increase in the value of the Euler characteristic (χ).According to Equation 3, an increase in the value of χ causes the connectivity index (λ) to decrease with time.
The points at which breakthrough was achieved in each of the scenarios are labeled "BT" in Figure 6.Tailing-off of the saturation curves is observed, as no significant changes in fluid displacement or changes in the defending phase saturation occur after fluid breakthrough.

Non-Neutrally Wet, θ = 45°and θ = 150°Simulation Cases
At low capillary numbers (i.e., log Ca = 4 and log Ca = 3), the connectivity index remained constant after breakthrough for the non-neutrally wet θ = 45°and 150°simulation cases, and the connectivity index oscillated about a mean value (Figures 6g,6i,6j,and 6l).This is because GD events continued to occur after breakthrough in the θ = 45°and θ = 90°cases; thus, the system remained dynamic with periodic disconnection and reconnection of flow pathways, causing fluctuations in the overall fluid topology and connectivity.The dynamic connectivity (λ) curves for the non-neutrally wet θ = 45°and 150°cases are shown in blue and red, respectively, in Figure 6.
Significant fluctuations in connectivity are observed post-breakthrough (BT) for the log Ca = 4, 3, and 2 for the θ = 45°and 150°simulation cases.These post-breakthrough fluctuations in connectivity are seen to weaken as the capillary number increases from log Ca = 4 to log Ca = 1 in Figure 6.This indicates that for these cases, the system becomes more dynamically stable with increasing capillary number.This is because at low capillary numbers (i.e., log Ca = 4 and log Ca = 3) for the non-neutrally wet cases where the contact angle is θ = 45°and θ = 150°, the GD regime was dominant.However, at high capillary numbers (e.g., log Ca = 1), the DTF regime dominated.The DTF regime is more dynamically stable and cohesive than the GD regime.As a result, the DTF regime has fluid invasion and fluid displacement characteristics similar to those of the CPF regime.

Neutrally Wet, θ = 90°Simulation Cases
Since the CPF regime dominated throughout the entire fluid displacement process for the neutrally wet θ = 90°s imulation cases at log Ca = 3 and 4, no fluctuations in connectivity were observed after breakthrough (BT) in Figures 6h and 6k.
The green curves in Figure 6 represent the dynamic connectivity curves for the 90°contact angle scenarios.By comparing the connectivity fluctuations for contact angles of 45°, 90°, and 150°at log Ca values of 4, 3, 2, and 1 in Figure 6, we can observe that the systems with θ = 90°were generally more stable after the breakthrough when compared to those with contact angles of 45°and 150°.
This can be attributed to the dominance of highly cohesive flow regimes that is, the CPF and DTF regimes.At low capillary numbers (log Ca = 4 and log Ca = 3) for θ = 90°, there was dominance of highly cohesive flow regimes, namely the CPF and DTF regimes.However, as the capillary number increased, there was a transition in dominance from the CPF to the DTF regime.This transition led to a slight increase in fluctuations in the connectivity curves for the θ = 90 simulation cases.The reason for this is that the DTF regime is less dynamically stable than the CPF regime.
For θ = 90°, there was a shift from CPF dominance at low capillary numbers to DTF dominance at high capillary numbers.However, this transition did not have as much impact on the fluctuations in post-breakthrough dynamic connectivity, compared to the cases where θ was 45°and 150°.In those cases, an increase in capillary number resulted in a shift from the GD regime to the DTF regime, which caused a significant change in displacement characteristics.This is because the difference in displacement characteristics between the DTF and CPF regimes is not as significant as the difference between the DTF and GD regimes, as shown in Figure 4.

Impact of Wettability and Invading Fluid Velocity on Ganglion Dynamics Onset Times, Displacement Efficiency, and Fluid Breakthrough
Figure 7a shows a plot of fluid disconnection onset times for cases where a transition from connected to disconnected flow occurred in the simulations performed.
The graph in Figure 7a demonstrates that as the capillary number increases, the onset time for fluid disconnection decreases.This means that fluid disconnection happens earlier.The reason for this is that the throats in the structure of the Berea sandstone being studied are much narrower than the pore bodies, with an average pore-tothroat size aspect ratio of 2.98.At high capillary numbers or fluid velocities, the invading fluid is pushed through the narrow throat constrictions at a high speed, which causes it to break up while passing through the throats.As the capillary number increases, the fluid velocity also increases, leading to earlier disconnection of invading fluid during the displacement process.
As the fluid velocity increases, the onset times for fluid disconnection become more similar, regardless of the wetting state at high capillary numbers.This is evident from Figure 7a, where the onset times converge toward the DTF regime.This indicates that the differences in wettability become less significant with increasing fluid velocity.The tendency of the flow toward the DTF regime with increasing capillary number for all wetting states is also visible in Figure 7b, which shows the log capillary number versus the log breakthrough plot.In Figure 7b, increasing capillary number leads to earlier breakthrough and more negative log BT for all wetting states.At low capillary numbers, a transition from the CPF to GD regime occurred in the θ = 45°and 150°cases, resulting in earlier breakthrough times.As the capillary number increased, the log BT values converged for all wetting conditions.Finally, we studied the changes in displacement efficiency with increasing capillary number.The displacement efficiency refers to the amount of the defending phase displaced and is expressed as a percentage of the original amount of the defending phase as follows: Displacement Efficiency = Amount of the Displaced Defending Phase Original Amount of the Defending Phase × 100 % (7) Table 3 presents the displacement efficiencies for the 12 simulated cases.
The fluid displacement efficiency increases with increasing invading fluid velocity (capillary number) at low capillary numbers, and the displacement efficiency is greater for the θ = 90°simulation case than for the θ = 45°a nd θ = 150°simulation cases.
Figure 7c illustrates the relationship between log displacement efficiency and log Ca for the simulated cases.At low capillary numbers, there is a significant contrast in the displacement efficiency between the neutrally wet state at θ = 90°, wherein the CPF regime dominates, and the non-neutrally wet states at θ = 45°and 150°, which are dominated by the GD regime.
The CPF regime results in slow but complete and efficient fluid displacement, while the GD regime leads to fast but incomplete and inefficient displacement of fluid from the pore spaces, as shown in Figure 4.
Although there are more disconnected individual fluid elements in the DTF regime than in the GD regime, the DTF regime presents better coherence of the invading fluid phase during the fluid invasion process.This happens because the small and disconnected individual fluid elements flow in close proximity to each other in the DTF regime, as shown in Figure 4.
In the DTF regime, these fluid elements periodically connect and disconnect, forming continuous flow paths at a high frequency.This, in turn, causes the fluid elements to remain connected for a considerable amount of time.
The intermittent flow behavior observed in the DTF regime causes high cohesion of the invading phase, similar to the CPF regime.As a result, there is very little difference in the displacement efficiency between the CPF regime and the droplet traffic (DTF) regime.
In contrast, the difference in fluid displacement efficiency between the ganglion dynamics (GD) and droplet traffic (DTF) regimes (see the blue 45°and red 150°series in Figure 7c) is very large due to the differences in the sparsity and cohesion of the flowing fluid elements, as previously shown in Figure 4.This is also true for the difference in displacement efficiency between the ganglion dynamics (GD) and CPF regimes (Figure 7c).As the capillary number increased, the displacement efficiencies for all the wetting states increased and converged as the flow tended toward the DTF regime regardless of the wettability.
By studying the temporal evolution of ganglia for different simulation cases, the differences in fluid disconnection onset times, dynamic fluid connectivity, and kinetic behavior can be identified.Figure 8 shows the temporal evolution of the number and average size of the defending phase ganglia in the investigated cases.
Figure 8 shows that at low capillary numbers, such as log Ca = 3 and log Ca = 4, there are noticeable variations in the rates of increase in the number of ganglia and the rates of decrease in average ganglion size for the 45°, 90°, and 150°wetting states.This is due to differences in the dominant pore-scale flow regimes for the three wetting states (Figures 8e-8h).A transition from the CPF regime to the GD regime occurred for θ = 45°and θ = 150°, which is apparent in the red and blue trends in Figures 8e and 8g, where the defending phase ganglia increase at a very slow rate during the initial stages of fluid displacement, where the CPF regime dominates.The occurrence of ganglion triggers an avalanche of GD events, resulting in a drastic increase in the number of defending phase ganglia.In Figures 8e and 8g, when the angle (θ) is 90°(green trend line), the rate of change in the number of ganglia is constant and much slower at low capillary numbers than at the blue θ = 45°and red θ = 150°trend lines.This occurs because the CPF regime is dominant, and there is no transition to disconnected flow.As the capillary number increases, differences in the rates of increase in the number of fluid ganglia and average ganglion size become less significant and less apparent for the three wetting states in Figure 8 as they tend toward the DTF regime.
Although there is a notable difference in the rates of decrease of the average ganglion size for θ = 45°, 90°, and 150°at log Ca = 4 (Figure 6h) at log Ca = 1, the trendlines for θ = 45°, 90°, and 150°overlap (Figure 8b).This shows that the flow tends toward the DTF regime with increasing capillary number and that flow patterns converge with increasing invading fluid velocity for all wetting states.This is also seen in the plots showing the temporal evolution of the number of ganglia (Figures 8b, 8d, 8f, and 8h).

Impact of Flow Transitions on the Shapes of the Saturation and Connectivity Curves
Apart from influencing dynamic fluid connectivity and post-breakthrough fluctuations in connectivity (λ), transitions from connected and disconnected to disconnected pathway flow regimes also impact the shape of the saturation and connectivity curves.Table 4 summarizes the flow regime transitions and shapes of the saturation and connectivity curves ("S"-shape or "L"-shape) in Figure 6.
When a transition occurs from the CPF regime (connected flow) to the GD regime (disconnected flow), for example, for θ = 45°and θ = 150°at a low log Ca of 4 or 3, the rate of decline in the defending phase is initially slow in the period where the CPF regime dominates.Once the first GD event occurs, it triggers system instability and a series of similar events.The onset of ganglion dynamics (denoted by "GD" in Figure 6 for each relevant simulation case) rapidly accelerates the kinetics of the fluid displacement process.
The first occurrence of a GD (gas displacement) event causes instability in the fluid system that increases the frequency of similar events, leading to a cascading effect.This results in a rapid increase in flow transport kinetics and abrupt changes in the gradient at the onset of GD in the saturation and connectivity (λ) curves at low capillary numbers (log Ca = 4 and log Ca = 3) for the non-neutrally wet θ = 45°and θ = 150°simulation cases (as shown in Figures 6g, 6i, 6j, and 6l).This is why "S-shaped" curves are observed in these cases.
As shown in Figures 6a-6d, 6h, and 6k, "L-shaped" saturation and connectivity curves occur when no transition transpires between flow regimes and only one flow regime dominates during the entire fluid displacement process.Examples can be seen in the log Ca = 4 and log Ca = 3 simulation cases for θ = 90°, where the CPF regime dominated during the entire fluid displacement process, and in all the log Ca = 1 simulations, where the DTF regime dominated for the entire fluid displacement process.
In cases where a transition occurred from the CPF regime to the DTF regime, for example, in the log Ca = 2 simulation case for θ = 90°, "L-shaped" curves are also observed in Figure 6e, as in cases where only one flow regime dominated throughout the entire fluid displacement process.The reason is that although this is a transition from connected to disconnected flow, it does not have much of an impact on flow transport kinetics or the shape of the saturation curve.This is because the DTF and CPF regimes are very similar in terms of fluid displacement patterns, flow cohesion and displacement efficiencies, as illustrated previously in Figure 4.

Implication of the Results
Two relative-saturation models/approximations are commonly used to represent two-phase flows at the continuum scale.These are the Corey approximation (Brooks & Corey, 1964;Corey, 1954) and the L-E-T approximation (Lomeland et al., 2005).The Corey approximation follows the power law in the saturation function and shows L-type behavior in the shape of relative permeability-saturation curves.The L-E-T approximation, on the other hand, is defined by a logistic function and shows S-type behavior in the shape of relative permeability-saturation curves.
The Corey correlations for the relative permeability of oil and water are defined as follows: where K row is the oil relative permeability, K rw is the water relative permeability, and K rot and K wr are endpoint parameters.N o and N w are empirical parameters and are referred to as curve shape parameters.N o and N w are obtained from measured data through analytical interpretation or optimization using numerical simulations to match experimental data (history matching).S w is the fraction of the pore volume occupied by water.Sw can be expressed as follows:  where S o is the fraction of the pore volume occupied by the oil phase.S wn is the normalized water saturation, which is defined as follows: S wn (S w ) = S w S wir 1 S wir S orw (11) Here, S wir is the residual water saturation, and S orw is the residual oil saturation after water flooding.The Corey model (Corey, 1954) is one of the oldest relative permeability-saturation approximation models and has been widely used since 1954.The L-E-T model was developed in 2005 by Lomeland et al. (2005) after the limitations of the Corey model became evident.Specifically, modeling of the upper and lower saturation endpoints, where the rates of fluid displacement and change in saturation can be slow for some fluid displacement processes, may not be accurate.
The L-E-T-type approximation is described by 3 parameters: L, E, and T. The L-E-T correlations for the relative permeability of oil and water are defined as follows: The L parameter defines the lower part of the curve, and it is comparable to the appropriate Corey parameter based on similarity and experience.The T parameter, on the other hand, describes the upper part of the curve in a similar manner as the L parameter describes the lower part of the curve.Finally, the E parameter defines the position of the slope of the curve.
The Corey model possesses only one degree of freedom to accommodate the shape of each relative permeabilitysaturation curve.The L-E-T model (Lomeland et al., 2005) adds more degrees of freedom to better match the Sshape of the relative permeability-saturation curves observed in the special core analysis (SCAL) experiments and in 3D reservoir models.Moghadasi et al. (2015) evaluated the Corey, Chierici, and L-E-T approximations for oil/water relative permeability.These researchers used a sophisticated method that accounted for the number of uncertain model parameters.The study revealed that the L-E-T approximation, which has the largest number of uncertain parameters (three), demonstrated the strongest correlation with both the oil and water relative permeabilities.Sakhaei et al. (2016) assessed 10 common and widely used relative permeability approximations for gas/oil and gas/ condensate systems.These scholars discovered that the L-E-T approximation exhibited the best correlation with the experimental values for the gas and oil/condensate relative permeabilities.
Wettability impacts the relative movement and distribution of fluids in a porous medium and alters the dynamics of relative permeabilities at the continuum scale.To better understand the impact of fluid interactions in porous media, Kassa et al. (2021) developed a framework to model the impact of wettability changes on relative permeabilities.In their study, they included a time-dependent wettability change in the relative permeabilitysaturation relation by modifying the existing relative permeability function.These investigators found that the L-E-T model is ideal for eliminating the hysteresis in relative permeability induced by contact angle changes caused by exposure to a wetting alteration agent during drainage and/or imbibition displacements.Furthermore, they reported that the dynamic L-E-T model has the potential to capture the underlying wettability process at the pore scale and suggested that further investigation may be needed to re-evaluate the potential of capturing the wettability alteration process in the relative permeabilities based on the interpolation approach.
In this work, the behavioral response of the saturation function to wettability and flow rate-induced changes in dynamic fluid connectivity at the pore scale was investigated.We show that shifts in pore-scale flow regimes and changes in dynamic fluid connectivity have a significant impact on flow transport kinetics, fluid displacement NHUNDURU ET AL.
rates and, consequently, the shape of the saturation curve/function.Two distinct shapes for the saturation curves are observed: the S-type curve and the L-type curve.
Since relative permeability-saturation functions are strong functions of the saturation distribution, changes in pore-scale flow regimes and dynamic fluid connectivity subsequently impact the relative permeability-saturation function.The results of this study suggest that the Corey model is suitable for describing displacement processes where a single flow regime dominates the fluid displacement process and processes where there is a transition from the CPF to the DTF regime.Such processes produce "L-Type" curves.
On the other hand, the L-E-T model is ideal for describing fluid displacement processes in which there is a transition from connected to disconnected flow, that is, a transition from the CPF regime to the ganglion dynamics flow regime (CPF to GD).This occurs in the early stages of fluid displacement at low capillary numbers (log Ca ≤ 3) in non-neutrally wet cases.When the CPF regime dominates, slow displacement occurs, and the onset of the GD regime accelerates the rate of The L-E-T model has more degrees of freedom than the Corey model, which allows it to capture the slow displacement in the upper saturation end points where the transition occurs.
The results of this work suggest that the L-E-T model is a true-to-mechanism model that can capture the impact of pore-scale changes in flow regimes on average flow transport kinetics.It is recommended that numerical upscaling studies be performed to link pore-scale observations with continuum-scale models.
Authors Kassa et al. (2021), andPicchi andBattiato (2019) have begun paving the way for such studies.Kassa et al. (2021) quantified the link between pore-scale model parameters and core-scale model parameters.Their simulation results show that a very simple scaling can relate a pore-scale process to the core scale.Picchi and Battiato (2018) proposed a set of upscaled equations based on pore-scale flow regimes, that is, the topology of flowing phases.In their study, the incompressible Navier-Stokes equation was upscaled using multiple-scale expansions, and its closures were derived from the mechanical energy balance for different flow regimes at the pore scale.These researchers provided analytical expressions for the relative permeability of wetting and nonwetting phases in different flow regimes and demonstrated that the effect of the flowing-phases topology on the relative permeabilities is significant.In addition, they show that the classical two-phase Darcy law is recovered for a limited range of operative conditions and that specific terms accounting for interfacial and wall interactions should be incorporated to accurately model ganglia or drop traffic flow.

Conclusion
This work explores and provides insight into the pore-scale transient process from connected to disconnected flow for subsurface flow applications.We investigate the sensitivity of this process to changes in flow and transport properties (fluid velocity and wettability).We demonstrate that the effect of wettability is significant at low capillary numbers (log Ca ≤ 3), where the dominance of the ganglion dynamics (GD) flow regime under nonneutrally wet conditions (θ ≠ 90°) reduces the fluid displacement efficiency by up to 21%.
We demonstrate that flow transport kinetics can, in fact, be strongly nonlinear functions of the invading fluid flow rate in cases where dynamic ganglion behavior occurs.At fixed saturation, fluid flow can occur through a series of periodic ganglion breakup, mobilization, and coalescence events (ganglion dynamics), which constantly rearrange the fluid topology in a network of connected pores.
We show that ganglion dynamics have a significant impact on the saturation function and flow transport kinetics.
The rate of change of the displaced phase saturation is sensitive to fluid velocity and wettability-induced transitions between connected and disconnected flow regimes, with sigmoidal or linear saturation curves being observed that depend on whether a single flow regime dominates the displacement process or if a transition between flow regimes occurs during the fluid displacement process.
The results of this work suggest that the L-E-T relative permeability model is a true-to-mechanism model that accurately captures the impact of changes in pore-scale flow regimes on average flow transport kinetics.Our current approach to continuum-scale flow modeling needs to be revised to account for the contribution of porescale flow phenomena such as disconnected flow and ganglion dynamics.This will enable the development of more precise subsurface flow models and resolve the discrepancies that are often encountered between pore-scale and field-scale predictive flow models.

Figure 1 .
Figure 1.(a) A digital model of the Berea sandstone microstructure used for two-phase displacement simulations (waterdecane displacement) in OpenFOAM.(b) Physical model (microfluidic model) of the Berea sandstone microstructure described previously by Wlodarczyk et al. (2019).The dimensions of the digital model are scaled down 10 times from those of the microfluidic model.

Figure 2 .
Figure 2. The effect of mesh refinement on the numerical solution.Mesh sensitivity analysis of the Berea sandstone microstructure for refinement factors of 0.75, 1, 1.5, 2, and 4.

Figure 3 .
Figure 3. Validation of the numerical method.Comparison of fluid invasion patterns in microfluidic and digital models of the Berea sandstone microstructure.Regions of the porous medium occupied by the invading water phase appear in blue, and regions occupied by the defending decane phase appear in red.
the time from the start of the simulation and breakthrough time (BT) is the time at which fluid breakthrough occurred.

Figure 5
Figure5shows that at high capillary numbers, for example, log Ca = 1, the fluid displacement patterns are uniform and quite similar under all wetting conditions.For neutrally wet conditions (contact angle = 90°) at low capillary numbers, that is, log Ca = 4 and log Ca = 3, the invading fluid (in navy blue) propagates via a single flow pathway (Figures5h and 5k), namely, the CPF regime.However, for non-neutrally wet conditions (θ ≠ 90°), at low capillary numbers, log Ca = 4 and log Ca = 3, for θ = 45°and 150°, the invading fluid (in navy blue) propagates in the porous structure in the form of large, disconnected fluid clusters, as shown in Figures5g, 5i, 5j, and 5l.These disconnected fluid clusters are mobile and periodically disconnect and reconnect during the fluid displacement process.This phenomenon is known as ganglion dynamics(Avraam & Payatakes, 1995;Picchi & Battiato, 2018;Valavanides, 2018).The ganglion behavior varied significantly between the θ = 45°and 150°c ases.This difference was due to variations in the wetting properties.

Figure 4 .
Figure 4. Flow regimes observed at breakthrough time: (a) droplet traffic flow regime observed in the Log Ca = 2; θ = 150°s imulation, (b) ganglion dynamics (GD) flow regime observed in the Log Ca = 4; θ = 150°simulation and (c) connected pathway flow regime observed in the Log Ca = 4; θ = 90°simulation.Here, χ is the Euler characteristic, and λ is the connectivity index.

Figure 7 .
Figure 7. Plots showing (a) fluid disconnection onset time versus log capillary number (log Ca), (b) log breakthrough time versus log capillary number (log Ca), and (c) log displacement efficiency versus log capillary number (log Ca) for contact angles of 45°, 90°, and 150°.DTF-droplet traffic flow, CPF-connected pathway flow, GD-ganglion dynamics.

Table 1
Properties of the Fluid Phases Used in the Two-Phase Flow Simulations voxel (δχ) and summing to give a Euler characteristic of the structure (χ).The plugin calculates the Euler characteristic of the structure as though it is floating in space (χ = ∑δχ). each