The role of phase topology in hysteresis during fluid injection and withdrawal in porous media is not fully understood. We address this by providing experimental and theoretical evidence on three key findings. (1) The topological evolution of the nonwetting fluid is distinct from the capillary pressure and the specific interfacial area, as shown by experiments and a generalized model. (2) Saturation paths with identical capillary pressure and interfacial area show different topologies, revealing insights into energy dissipation and phase connectivity. (3) The topological evolution of the nonwetting phase follows predictable trajectories captured by a piecewise non-linear model. These findings offer practical implications for optimizing subsurface hydrogen and carbon dioxide storage systems and provide a novel approach to study complex systems with topological singularities.

Relative permeability plays an important role in the upscaling of multiphase flow in porous media from the pore scale to the Darcy scale. The entire concept of relative permeability is contingent on the existence of a representative elementary volume (REV). As we move to smaller samples to measure relative permeability, such as with digital core analysis, the concept of a classical REV has become increasingly unlikely when using the conventional approach to defining a representative volume. The “‘conventional”’ understanding of an REV is that a large enough volume must be considered such that spatial variability averages out. In digital rock methods, such as pore-scale simulations based on micro-computed tomography (CT) images, the domain size is typically 2 to 4 mm. This is approximately the length scale of a single-phase flow REV using the classic REV approach. However, the single-phase perspective does not consider the complex dynamics and fluctuations often observed in multiphase flow systems, even at centimeter-scale experiments and/or simulations. A fundamental question is, therefore, whether the domain size commonly used in digital rock simulations can provide a consistent energy budget such that the concept of relative permeability exists. Based on first principles, relative permeability accounts for the rate of energy dissipated in a stationary process. If the dynamics are fluctuating, the energy dissipated can vary but will average out over a long enough timescale. The key to determining the validity of the relative permeability is the timescale of the measurement, not the spatial scale. The conventional REV theory assumes that spatial, temporal, and ensemble averages are equivalent in an ergodic system, but it does not provide a way to test this assumption. We provide a formal way to identify the timescale where the relative permeability accurately captures energy dissipation as a way to validate relative permeability measurements and quantitatively assess their accuracy. This result will be tested for a practical SCAL test, determining how long a flow experiment needs to be run to accurately characterize the rate of energy dissipation by the flow. The outcome will be a best practice guide for the determination of relative permeability from core-scale experiments and/or digital core simulations that ensure the energy budget is fully accounted for in the relative permeability coefficient.

Direct numerical simulations are commonly applied to X-ray computed microtomography images of porous rocks to measure petrophysical properties, such as relative permeability. To complement and/or expand laboratory data with useful simulations, such simulations need to reflect the true wetting condition of reservoir rock. It is therefore important to investigate approaches for assigning surface wetting conditions in pore -scale multiphase flow simulations. Two-phase Lattice -Boltzmann simulations are conducted on a Ketton carbonate image using different effective wetting models and then compared to pore -scale experimental data on the same rock, which is published elsewhere in (Scanziani et al., 2018). Two different wetting models are established based on the aging of the Ketton sample with crude oil at connate water saturation. Model 1 provides uniform wetting on all oil treated surfaces. Model 2 attempts to capture the multiscale nature of carbonate by allowing for spatially varying wetting conditions on oil treated surfaces based on local microporosity. Assessment criteria are based on porescale morphological measures, which include oil surface area coverage, oil phase topology, fluid/fluid curvatures, and contact angles. Overall, Model 1 with an assigned contact angle of greater than +10 degrees above the geometrically measured contact angle from the experimental data provided the best results based on our assessment criteria.

The impact of wettability on the co-moving velocity of two-fluid flow in porous media is analyzed herein. The co-moving velocity, developed by Roy et al. (2022), is a novel representation of the flow behavior of two fluids through porous media. Our study aims to better understand the behavior of the co-moving velocity by analyzing simulation data under various wetting conditions. The simulations were conducted using the Lattice-Boltzmann color-fluid model and evaluated the relative permeability for different wetting conditions on the same rock. The analysis of the simulation data followed the methodology proposed by Roy et al. (2022) to reconstruct a constitutive equation for the co-moving velocity. Surprisingly, it was found that the coefficients of the constitutive equation were nearly the same for all wetting conditions. Based on these results, a simple approach was proposed to reconstruct the oil phase relative permeability using only the co-moving velocity relationship and water phase relative permeability. This proposed method provides new insights into the dependency of relative permeability curves, which has implications for the history matching of production data and solving the associated inverse problem. The research findings contribute to a better understanding of the impact of wettability on fluid flow in porous media and provide a practical approach for estimating relative permeability based on the co-moving velocity relationship, which has never been shown before.

Underground hydrogen storage (UHS) is receiving increasing attention to address the challenges in hydrogen storage. A crucial aspect of UHS is understanding the transport of hydrogen in subsurface porous media. In this work, hydrogen core flooding experiments were conducted in a sandstone sample and then the pore-scale distribution of hydrogen and brine was visualized using high-resolution X-ray micro-computed tomography (micro-CT). After CT image processing, we measured the surface contact angles (CAs) on rock surfaces and found that the measured CAs followed a log-normal distribution. To investigate the influence of rock surface wettability on the transport properties in the hydrogen-brine-sandstone system, the lattice Boltzmann (LB) method was used to simulate two-phase flows with different surface CA distributions; hydrogen-brine relative permeability and capillary pressure curves through the primary drainage, imbibition, and secondary drainage processes were obtained by LB simulations. X-ray CT scanning showed that hydrogen resided in both large pores and small pores and pore throats after the primary drainage process (i.e., hydrogen displacing brine), whereas after the imbibition process (i.e., brine displacing hydrogen) hydrogen stayed primarily in large pores where the capillary pressure barriers were low. The LB two-phase flow modeling illustrated that water’s relative permeability increased whereas hydrogen’s relative permeability decreased when the core flooding moved from primary drainage to imbibition; in contrast, when the core flooding moved from imbibition to secondary drainage, water’s relative permeability decreased whereas hydrogen’s relative permeability increased. Decreasing the surface CA increased the capillary pressure and reduced the hydrogen residual rate, which indicates high hydrogen retrievability and thus is favorable for UHS practice. The change in CA’s standard deviation did not cause noticeable changes in water’s relative permeability curves, whereas it resulted in noticeable changes in hydrogen’s relative permeability curves. This work is the first study that utilizes X-ray micro-CT scanning and pore-scale multiphase flow modeling to quantitatively and comprehensively investigate hydrogen-brine two-phase flow properties under different rock surface wettability distributions. The research findings from this study will advance the understanding of hydrogen transport in an underground storage system and provide essential data for large-scale field studies in UHS.

This study aims to bridge length scales in immiscible multiphase flow simulation by connecting two published governing equations at the pore-scale and continuum-scale through a novel validation framework. We employ Niessner and Hassnaizadeh's [“A model for two-phase flow in porous media including fluid-fluid interfacial area,” Water Resour. Res. 44(8), W08439 (2008)] continuum-scale model for multiphase flow in porous media, combined with the geometric equation of state of McClure et al. [“Modeling geometric state for fluids in porous media: Evolution of the Euler characteristic,” Transp. Porous Med. 133(2), 229–250 (2020)]. Pore-scale fluid configurations simulated with the lattice-Boltzmann method are used to validate the continuum-scale results. We propose a mapping from the continuum-scale to pore-scale utilizing a generalized additive model to predict non-wetting phase Euler characteristics during imbibition, effectively bridging the continuum-to-pore length scale gap. Continuum-scale simulated measures of specific interfacial area, saturation, and capillary pressure are directly compared to up-scaled pore-scale simulation results. This research develops a numerical framework capable of capturing multiscale flow equations establishing a connection between pore-scale and continuum-scale simulations.

Ion transport within saturated porous media is an intricate process in which efficient ion delivery is desired in many engineering problems. However, controlling the behavior of ion transport proves challenging, as ion transport is influenced by a variety of driving mechanisms, which requires a systematic understanding. Herein, we study a coupled advection–diffusion–electromigration system for controlled ion transport within porous media using the scaling analysis. Using the Lattice–Boltzmann–Poisson method, we establish a transport regime classification based on an Advection Diffusion Index (ADI) and a novel Electrodiffusivity Index (EDI) for a two-dimensional (2D) microchannel model under various electric potentials, pressure gradients, and concentration conditions. The resulting transport regimes can be well controlled by changing the applied electric potential, the pressure field, and the injected ions concentration. Furthermore, we conduct numerical simulations in a synthetic 2D porous media and an x-ray microcomputed tomography sandstone image to validate the prevailing transport regime. The simulation results highlight that the defined transport regime observed in our simple micromodel domain is also observed in the synthetic two- and three-dimensional domains, but the boundary between each transport regime differs depending on the variation of the pore size within a given domain. Consequently, the proposed ADI and EDI emerge as dimensionless indicators for controlled ion transport. Overall, our proof-of-concept for ion transport control in porous media is demonstrated under advection–diffusion–electromigration transport, demonstrating the richness of transport regimes that can develop and provide future research directions for subsurface engineering applications.

Niessner and Hassnaizadeh (2008) proposed a set of governing equations to model multiphase flow in porous media based on first principles that includes fluid/fluid specific interfacial area as a state variable. Herein, we shed light on the lesser-discussed aspects of their proposed model, referred to as the extended model. Firstly, we show how the extended model facilitates real-time tracking of specific interfacial area. Secondly, we compare numerical solutions of the extended model to that of the traditional multiphase flow model with capillary dispersion and the analytical Buckley–Leverett solution. Thirdly, we provide a sensitivity analysis for the additional unknown parameters of the extended model - specific interfacial permeability and specific interfacial area generation. Lastly, a stability analysis is performed for the numerical solutions. As an auxiliary outcome, we provide an approach to solve the non-linear partial differential equations of the extended model by using torch.autograd as the automatic differentiation engine of PyTorch. Our results demonstrate that the extended model can provide saturation profiles similar to those generated by traditional models, but with the added benefit of real-time tracking of specific interfacial area. The generated specific interfacial area profiles were comparable to expected trends from previous numerical and experimental studies and highlighted the importance of the capillary pressure relationship and specific interfacial area generation term. The results also demonstrate an important balance between the permeability of the porous media and the interfacial permeability that must be maintained; otherwise, generated profiles become physically unrealistic. Overall, the extended model explains the evolution of specific interfacial area during multiphase flow at the cost of additional non-linearity and unknown parameters.

Wettability is one of the critical parameters affecting multiphase flow in porous media. The wettability is determined by the affinity of fluids to the rock surface, which varies due to factors such as mineral heterogeneity, roughness, ageing, and pore-space geometry. It is well known that wettability varies spatially in natural rocks, and it is still generally considered a constant parameter in pore-scale simulation studies. The accuracy of pore-scale simulation of multiphase flow in porous media is undermined by such inadequate wettability models. The advent of in situ visualization techniques, e.g. X-ray imaging and microtomography, enables us to characterize the spatial distribution of wetting more accurately. There are several approaches for such characterization. Most include the construction of a meshed surface of the interface surfaces in a segmented X-ray image and are known to have significant errors arising from insufficient resolution and surface-smoothing algorithms. This work presents a novel approach for spatial determination of wetting properties using local lattice-Boltzmann simulations. The scheme is computationally efficient as the segmented X-ray image is divided into subdomains before conducting the lattice-Boltzmann simulations, enabling fast simulations. To test the proposed method, it was applied to two synthetic cases with known wettability and three datasets of imaged fluid distributions. The wettability map was obtained for all samples using local lattice-Boltzmann calculations on trapped ganglia and optimization on surface affinity parameters. The results were quantitatively compared with a previously developed geometrical contact angle determination method. The two synthetic cases were used to validate the results of the developed workflow, as well as to compare the wettability results with the geometrical analysis method. It is shown that the developed workflow accurately characterizes the wetting state in the synthetic porous media with an acceptable uncertainty and is better to capture extreme wetting conditions. For the three datasets of imaged fluid distributions, our results show that the obtained contact angle distributions are consistent with the geometrical method. However, the obtained contact angle distributions tend to have a narrower span and are considered more realistic compared to the geometrical method. Finally, our results show the potential of the proposed scheme to efficiently obtain wettability maps of porous media using X-ray images of multiphase fluid distributions. The developed workflow can help for more accurate characterization of the wettability map in the porous media using limited experimental data, and hence more accurate digital rock analysis of multiphase flow in porous media.

In this paper, the underlying problem with the color -gradient (CG) method in handling density-contrast fluids is explored. It is shown that the CG method is not fluid invariant. Based on nondimensionalizing the CG method, a phase-field interface-capturing model is proposed which tackles the difficulty of handling density-contrast fluids. The proposed formulation is developed for incompressible, immiscible two-fluid flows without phase-change phenomena, and a solver based on the lattice Boltzmann method is proposed. Coupled with an available robust hydrodynamic solver, a binary fluid flow package that handles fluid flows with high density and viscosity contrasts is presented. The macroscopic and lattice Boltzmann equivalents of the formulation, which make the physical interpretation of it easier, are presented. In contrast to existing color -gradient models where the interface-capturing equations are coupled with the hydrodynamic ones and include the surface tension forces, the proposed formulation is in the same spirit as the other phase-field models such as the Cahn-Hilliard and the Allen-Cahn equations and is solely employed to capture the interface advected due to a flow velocity. As such, similarly to other phase-field models, a so-called mobility parameter comes into play. In contrast, the mobility is not related to the density field but a constant coefficient. This leads to a formulation that avoids individual speed of sound for the different fluids. On the lattice Boltzmann solver side, two separate distribution functions are adopted to solve the formulation, and another one is employed to solve the Navier-Stokes equations, yielding a total of three equations. Two series of numerical tests are conducted to validate the accuracy and stability of the model, where we compare simulated results with available analytical and numerical solutions, and good agreement is observed. In the first set the interfacial evolution equations are assessed, while in the second set the hydrodynamic effects are taken into account.

The flocking of self-propelled particles in heterogeneous environments is relevant to both natural and artificial systems. The Vicsek model is a canonical choice to investigate such systems due to the minimal number of parameters required to define flocking. Prior research on the Vicsek model has investigated the effects of interaction rules, particle speed, and obstacle packing on the flocking behavior, but the effect of interaction radius remains an open question. Unlike obstacle-free domains, the locality of interactions not only affects how quickly the system can become polarized, but also how well the flocks can align or realign after colliding with obstacles. In this letter, we delve into this subtle relationship that exists in the scale of the perception of Vicsek particles in the presence of obstacles. We demonstrate that the presence of obstacles impacts group density, which provides the basis to identify distinct phases for collective behavior. This leads to the counter-intuitive result that obstacles, while generally confounding for macroscopic order, may enable global order even as noise in the system increases.

Flow fluctuations that are commonly associated with multiphase flow in porous media are studied using concepts from non-equilibrium thermodynamic and statistical mechanics. We investigate how the Green–Kubo formulation of the fluctuation dissipation theorem can be used to predict the transport coefficient from the two-phase extension of Darcy's law. Flow rate-time series data are recorded at the millisecond timescale using a novel experimental setup that allows for the determination of flow fluctuation statistics. By using Green–Kubo relations, a transport coefficient is predicted based on the integrated autocorrelation function. Notably, this coefficient aligned closely with the total effective phase mobility computed using Darcy's equation for multiphase flow, particularly in scenarios where a linear relationship between flow rate and pressure gradient was observed. Our results open a new field of coefficient explorations where microscale fluctuations during multiphase flow are directly linked to macroscale parameters.

Various researchers have studied fluctuations in pore-scale phase occupancy during multiphase flow in porous media using synchrotron-based X-ray microcomputed tomography (micro-CT). However, the impact of these fluctuations on the concept of a representative volume is not yet fully understood. In this study, we performed spatial and temporal averaging of multiphase flow experiments visualized with synchrotron-based micro-CT, focusing on oil saturation as the key parameter to determine a representative time-and-space average. Our findings revealed that a saturation value representative of both time and space was achieved during fractional flow experiments in drainage mode with fractional flows of 0.8, 0.5, and 0.3. Furthermore, we computed a range of relative permeabilities on the basis of whether momentaneous saturation or time-and-space averaged saturation was utilized for direct simulation. Our results highlighted the importance of time-and-space averaging in determining a representative relative permeability and indicated that the temporal and spatial scales covered in a typical micro-CT flow experiment were sufficient to obtain a representative saturation value for sandstone rock under intermittent flow conditions. This research study focuses on the analysis of fluid flow behaviors in porous media during multiphase flow experiments. We used synchrotron-based X-ray micro-computed tomography technology to investigate the observed fluctuations and intermittency in porous media flows. We conducted spatial and temporal averaging on experimental dynamic 3D images of fractional flow, specifically focusing on oil saturation as a key parameter. The results showed that for a steady-state fractional flow experiment conducted in drainage mode, saturation convergence was typically reached within the duration of a standard synchrotron-based micro-CT fractional flow experiment. However, when the fractional flow was half water and half oil, the system exhibited increased fluctuations and required a longer time to achieve convergence. We also simulated the relative permeability and found that using experimentally observed saturation provided a range of values for each fractional flow, highlighting the importance of time-and-space averaging for determining a representative relative permeability. Overall, we demonstrated that the temporal and spatial scales sampled during a typical micro-CT flow experiment are adequate to provide a representative saturation value for sandstone rock under intermittent flow conditions, and that time-and-space averaging provides a consistent way to define a representative saturation. Synchrotron X-ray micro-CT reveals flow intermittency during the co-injection of immiscible fluids in porous media Saturation fluctuations identified for different fractional flows averaged out over the time scale of a typical experiment Using only a spatial average saturation could lead to inaccuracies in relative permeability. Time-and-space averaging is recommended

Proton exchange membrane fuel cells, consuming hydrogen and oxygen to generate clean electricity and water, suffer acute liquid water challenges. Accurate liquid water modelling is inherently challenging due to the multi-phase, multi-component, reactive dynamics within multi-scale, multi-layered porous media. In addition, currently inadequate imaging and modelling capabilities are limiting simulations to small areas (<1 mm2) or simplified architectures. Herein, an advancement in water modelling is achieved using X-ray micro-computed tomography, deep learned super-resolution, multi-label segmentation, and direct multi-phase simulation. The resulting image is the most resolved domain (16 mm2 with 700 nm voxel resolution) and the largest direct multi-phase flow simulation of a fuel cell. This generalisable approach unveils multi-scale water clustering and transport mechanisms over large dry and flooded areas in the gas diffusion layer and flow fields, paving the way for next generation proton exchange membrane fuel cells with optimised structures and wettabilities.

In this paper, we investigate phase-field interface capturing equations for two-fluid systems to probe their accuracy and computational cost. Two different schemes are considered: In the first scheme, one of the two order parameters is numerically solved based on a phase-field equation, while the other order parameter is determined through the summation relation; the summation of order parameters equals unity. In the second scheme, the two order parameters are both obtained numerically by solving their respective phase-field equations. A phase-field model based on the color-gradient (CG) method is chosen, and available lattice Boltzmann models are employed for solving the interface-capturing equations together with the hydrodynamic equation. It is shown that for the first scheme, which includes the summation relation, numerical results become asymmetrical. Also, in some cases, it results in nonphysical interfaces. In terms of computational resources, this first scheme is about 11% faster with 25% less computational memory usage than the second scheme. It is shown that only for a zero velocity domain do the two schemes lead to equal results. Also, a theoretical analysis is conducted to highlight the differences between the two approaches.

Groundwater pollution poses a significant threat to environmental sustainability during urbanization. Existing remediation methods like pump-and-treat and electrokinetics have limited ion transport control. This study introduces a coupled advection-diffusion-electromigration system for controlled ion transport in the subsurface. Using the Lattice-Boltzmann-Poisson method, we simulate ion transport in various two- and three-dimensional porous media. We establish an ion transport regime classification based on the Peclet number (Pe) and a novel Electrodiffusivity index (EDI). By manipulating the electric potential, hydrostatic pressure, and ion concentration, we identify four transport regimes: large channeling, uniform flow, small channeling, and no flow. Large channeling occurs when advection dominates, while uniform flow arises when diffusion and electromigration are more prevalent. Small channeling happens when the advection opposes electromigration and diffusion, and no flow occurs when the advection or electromigration impedes ion transport via diffusion. Simulations in heterogeneous models confirm these transport regimes, highlighting the influence of pore size variation on transport regimes. Consequently, $Pe$ and $EDI$ must be tailored for optimal transport control. These findings enable better control over ion transport, optimizing processes such as heavy metal removal, bioremediation, and contaminant degradation in groundwater management.

Multiphase flows through reservoir rocks are a universal and complex phenomenon. Relative permeability is one of the primary determinants in reservoir performance calculations. Accurate estimation of the relative permeability is crucial for reservoir management and future production. In this paper, we propose inferring relative permeability curves from sparse saturation data with an ensemble Kalman method. We represent these curves through a series of positive increments of relative permeability at specified saturation values, which guarantees monotonicity within, and boundedness between, 0 and 1. The proposed method is validated by the inference performances in two synthetic benchmarks designed by SPE and a field-scale model developed by Equinor that includes certain real-field features. The results indicate that the relative permeability curves can be accurately estimated within the saturation intervals having available observations and appropriately extrapolated to the remaining saturations by virtue of the embedded constraints. The predicted well responses are comparable to the ground truths, even though they are not included as the observation. The study demonstrates the feasibility of using ensemble Kalman method to infer relative permeability curves from saturation data, which can aid in the predictions of multiphase flow and reservoir production.

For multi-phase flow through multi-scale heterogeneous porous media, such as the pore space of rocks, the interaction between multiple immiscible fluids and an intricate network of pores, creates a wide range of dynamic flow phenomena. At larger scales i.e. scales relevant for practical applications such as carbon sequestration, this interplay of dynamic phenomena is often referred to as "complexity". However, it is important to describe the persistent features of the flow in an adequate manner, to represent the "complexity" of the system. Dynamic Mode Decomposition (DMD) is a dimensionality reduction algorithm that computes a set of modes associated with fixed oscillatory behaviours. In this work, saturation data extracted from dynamic two-phase flow experiments are analysed with DMD. We show that DMD can reproduce the data. Furthermore, not all dynamic modes are required to reproduce key dynamic features; this highlights the important spatial and temporal scales for flow. We show that DMD was able to identify localized regions important to flow. Overall, DMD is proven as a useful diagnostic tool for complex 4D flow dynamics for multi-phase flow.

Electrokinetic in-situ recovery is an alternative to conventional mining, relying on the application of an electric potential to enhance the subsurface flow of ions. Understanding the pore-scale flow and ion transport under electric potential is essential for petrophysical properties estimation and flow behavior characterization. The governing physics of electrokinetic transport is electromigration and electroosmotic flow, which depend on the electric potential gradient, mineral occurrence, domain morphology (tortuosity and porosity, grain size and distribution, etc.), and electrolyte properties (local pH distribution and lixiviant type and concentration, etc.). Herein, mineral occurrence and its associated zeta potential are investigated for EK transport. The new Ek model which is designed to solve the EK flow in complex porous media in a highly parallelizable manner includes three coupled equations: (1) Poisson equation, (2) Nernst–Planck equation, and (3) Navier–Stokes equation. These equations were solved using the lattice Boltzmann method within X-ray computed microtomography images. The proposed model is validated against COMSOL multiphysics in a two-dimensional microchannel in terms of fluid flow behavior when the electrical double layer is both resolvable and unresolvable. A more complex chalcopyrite-silica system is then obtained by micro-CT scanning to evaluate the model performance. The effects of mineral occurrence, zeta potential, and electric potential on the three-dimensional chalcopyrite-silica system were evaluated. Although the positive zeta potential of chalcopyrite can induce a flow of ferric ion counter to the direction of electromigration, the net effect is dependent on the occurrence of chalcopyrite. However, the ion flux induced by electromigration was the dominant transport mechanism, whereas advection induced by electroosmosis made a lower contribution. Overall, a pore-scale EK model is proposed for direct simulation on pore-scale images. The proposed model can be coupled with other geochemical models for full physicochemical transport simulations. Meanwhile, electrokinetic transport shows promise as a human-controllable technique because the electromigration of ions and the applied electric potential can be easily controlled externally. Graphical abstract