4D Imaging of Two‐Phase Flow in Porous Media Using Laboratory‐Based Micro‐Computed Tomography

Dynamic three‐dimensional computed tomography (CT) imaging of liquid transport in porous media has primarily been conducted at high brilliance synchrotrons thus allowing fast, sometimes sub‐second, temporal resolution to be obtained. University laboratory CT instruments lack the photon flux available at synchrotrons, limiting the obtainable spatiotemporal resolution. Here, we discuss our experiences with instrumentation and software methods to conduct time‐resolved micro‐computed tomography (4D‐CT) experiments of flow in porous media, based on a conventional CT instrument operated with a highly undersampled number of projections. An experimental stage outfitted with syringe pumps placed on a slip ring allowed two‐phase flow experiments to be carried out with continuous unidirectional rotation and without obstruction of the liquid supply lines. An iterative reconstruction algorithm based on a priori information was used to provide high image quality and ∼30 s time resolution despite the few and low‐exposed projections compared to standard protocols. The experimental technique was demonstrated with imbibition and drainage in glass bead‐pack and Bentheimer sandstone samples with sub‐minute temporal resolution, allowing the liquid configurations just before and after fast dynamic phenomena such as cooperative pore‐filling events and Haines jumps to be captured. Power law scaling exponents for burst volumes associated with imbibition and drainage were estimated and compared with the literature. That 4D‐CT experiments can be carried out using conventional CT instruments to challenge contemporary permeability models is of high importance for many geo‐, bio‐ and environmental physics challenges.


Introduction
Over the last decade, the scientific community has witnessed significant improvements in the field of temporal CT imaging, enabling dynamic experiments to be carried out with increasingly good time-resolution.For synchrotron experiments, which offers both high brilliance sources, bespoke experimental setups, and rapid readout X-ray detector systems, sub-second time-resolution can be obtained for certain CT experiments (García-Moreno et al., 2021;Piovesan et al., 2020;Singh, Scholl, et al., 2017).Limiting the scope to multi-phase flow in porous media, time-resolved "4D CT" (i.e., 3 spatial dimensions and time) allows fluid phases and their intricate dynamics to be tracked at the pore-scale.For instance, Berg et al. (2013) carried out a three-dimensional study of liquid dynamics at the pore scale with a time-resolution of 16.8 s, demonstrating the feasibility of conducting dynamic two-phase flow imaging experiments with time-resolved CT.Since then, many dynamic studies of porous media multi-phase flow have been conducted at synchrotrons, including flow with transient, steady-state, and reactive flow (Piovesan et al., 2020;Schlüter et al., 2016;Singh, Menke, et al., 2017).Most recently, we studied the interfacial dynamics of Haines jumps in a glass shard sample with sub-millisecond temporal resolution using a stroboscopic acquisition scheme, providing new ways of studying this fascinating phenomenon (Tekseth et al., 2024).However, access to synchrotrons is a scarce and costly resource, and cannot substitute the regular research work that can be done in home-laboratories.Albeit being a tool that is coveted by many academic and industrial researchers, 4D-CT techniques for dynamic studies of fluid flow in three-dimensional porous media are still challenging and laborious to establish and execute.
Dynamic studies of two-phase flow experiments with home-laboratory μCT instruments are still in an early stage of development, primarily because these instruments lack the high photon flux of synchrotron sources.The primary parameters defining CT scanning times are the exposure time and number of projections needed for faithful reconstruction, and the overhead time for motor movements.Conventional reconstruction schemes such as the Feldkamp-Davis-Kress algorithm (FDK) require sufficient projection sampling and a high signal-to-noise ratio to obtain high-quality reconstructions without undersampling artifacts and noise in the reconstructed volumes (Feldkamp et al., 1984).Known as the Crowther criterion, sufficient sampling requires acquiring a number of projections N proj ≥ πN d /2, where N d is the number of projection pixels normal to the sample rotation axis (Crowther et al., 1970;Jacobsen, 2018;Maire & Withers, 2014).This criterion implies that for a 2D detector with 2,000 pixels across, at least 3141 projections must be acquired to satisfy the sampling requirement if the sample fills the entire field of view (FoV).Consequently, with exposure times per projection typically in the range of seconds, regular CT scans in the home-laboratory are slow, lasting from tens of minutes to several hours.
To circumvent the photon flux limitation of home-laboratory CT instruments, considerable research efforts have been devoted to improving and extending both the experimental setups for data acquisition and the reconstruction frameworks.Gantry-based CT setups, denoting instruments where the sample remains fixed and being encircled by source(s) and detector(s), are commonly used in hospitals to allow the patient to rest.For dynamic CT, gantry setups have the advantage that the sample will not be accelerated, in addition to facilitating in situ experiments with complex and bulky sample environments requiring external cables and/or supply lines (Dierick et al., 2014).For materials science with higher demands to spatial resolution, gantry setups are however less common, but custom μCT setups are reported for dynamic μCT (Dierick et al., 2014).Industrial CT scanners usually consist of a fixed source and detector system and a rotating sample, which is necessarily also the overall constraint with a fixed direction of the incoming beam at synchrotron and neutron sources.
On the software side practical iterative reconstruction frameworks that suppress undersampling artifacts and noise are developed.With the advent of compressed sensing (Candès et al., 2006), reconstruction algorithms with few projections have been shown to faithfully recover the three-dimensional sample structure by minimizing the l 1norm of an objective function, provided that the object is sparse (Candès et al., 2006;Chambolle & Pock, 2011;Sidky & Pan, 2008a;Sidky et al., 2006).Regularization methods utilizing total variation (TV) (Rudin et al., 1992) have been demonstrated to provide reliable reconstruction quality under a variety of difficult conditions including few projections, limited-angle views, and noise (Beck & Teboulle, 2009;Chambolle & Pock, 2011;Sidky & Pan, 2008b;Sidky et al., 2006).TV regularization is particularly suitable for samples that are sparse, in practice often implemented as having a sparse image gradient, corresponding to a sample having uniform and homogeneous sub-volumes separated by relatively sharp boundaries.Such uniform and homogeneous regions are often encountered in fluid-flow experiments with immiscible phases and porous media such as idealized glass beadpacks or natural sandstone samples like Fontainebleau and Bentheimer consisting primarily of quartz minerals.
A priori sample information can also be incorporated into iterative reconstruction algorithms.In discrete tomography, the sample is assumed to consist of a few discrete densities (phases) that the algorithm uses as a constraint for the reconstructed density distribution, thereby also effectively doing segmentation as part of the reconstruction framework (Batenburg & Sijbers, 2011;van Aarle et al., 2012).Discrete tomography has been combined with TV regularization to further enhance the image quality (Zhuge et al., 2016).Another approach to utilizing a priori sample information is to incorporate knowledge of the static structure, typically obtained as a high-quality CT scan, as a constraint to the reconstruction algorithm of the sample parts that remain stationary during the dynamics (Myers et al., 2015;Rasmussen et al., 2021;Van Eyndhoven et al., 2015).Prior image constrained compressed sensing (PICCS) is a prime example of this idea (G.H. Chen et al., 2008).The development of iterative reconstruction algorithms has been facilitated by a growing number of open-source toolboxes such as ASTRA (van Aarle et al., 2015), TIGRE (Biguri et al., 2016), TomoPy (Gürsoy et al., 2014), and PyHST2 (Mirone et al., 2014), with GPU-accelerated projection operators readily accessible through high-level programming languages.
The rapidly evolving field of Machine Learning (ML) has been extensively applied to CT during the last decade.Specifically, ML has been used for reconstruction in 2D and 3D, segmentation, classification, and image Water Resources Research 10.1029/2023WR036514 enhancement (J.Chen et al., 2020;Jin et al., 2017;Liu et al., 2020;Park et al., 2018;Seo et al., 2019).CT reconstruction can be performed using a generative adversarial network (GAN) (Yang et al., 2020).This specific GAN is an example of an end-to-end approach, implying that the network was trained to directly retrieve a 3D model from the raw data.As an alternative to end-to-end, the reconstruction can be split into several steps, where for CT reconstruction, the first step could typically utilize a conventional reconstruction algorithm like FDK, followed by a neural network performing image enhancement that essentially reduces noise and undersampling artifacts in the second step (Guo et al., 2022;Liu et al., 2020).In the context of CT, concepts from artificial intelligence have also been used to automatically segment 3D data of porous media (Phan et al., 2021), to assess porous media properties such as porosity, specific surface area and average pore sizes (Alqahtani et al., 2019), and to estimate the permeability of porous media (Tian et al., 2021).
Here we present and discuss our experiences with capturing and reconstructing time-resolved 3D images of quasistatic two-phase flow experiments in porous samples of varying complexity using a conventional commercially available home-laboratory X-ray micro-CT setup.We developed a custom-made experimental rig consisting of syringe pumps, sensors, and a slip ring to enable fast CT acquisitions.We demonstrate our setup for imaging multiphase fluid flow dynamics inside two different porous samples, being a glass bead-pack model system and a Bentheimer sandstone sample.Furthermore, we compare the statistical distribution of pore-filling events in imbibition and drainage experiments in the bead-pack sample to those reported in the literature through simulations and experiments.Our work contributes to the flourishing topic of enabling high-performance 4D-CT with conventional home-laboratory setups, intended to give broader access to scientific experiments formerly primarily conducted at synchrotrons.

Reconstruction Framework for Dynamic CT Scans
As described in the Introduction, "dynamic" or "4D" CT is a balancing act between spatial and temporal resolution (Bultreys et al., 2016).With severe projection undersampling, the much-used non-iterative CT reconstruction algorithms like filtered back-projection (FPB) based on the inverse Radon transform, cannot satisfactorily reconstruct 3D-models.Consequently, additional information or assumptions must be invoked, and to this end, the various iterative reconstruction methods, like simultaneous iterative reconstruction technique (SIRT), and maximum likelihood estimation (MLE), that have been developed over the years are highly useful, as they allow additional constraints to be incorporated into the reconstruction procedure (Kak & Slaney, 2001).
In the current study of reconstructing dynamic data sets, we rely on the iterative reconstruction algorithm known as prior image constrained compressed sensing (PICCS), which is a way to realize sparse-view reconstruction (G.H. Chen et al., 2008;Lauzier et al., 2012).Fundamentally, this algorithm assumes that large regions of the sample are unchanged ("static") during the dynamics, and employs a high-quality 3D model of the static structure as a prior to guide the reconstruction of undersampled dynamic CT data, through minimizing the objective function Here, x prior is the pre-measured high-resolution prior volume, x is the volume to be optimized and α is a weighting parameter between the prior volume and the optimized volume.The first term in Equation 1 is a data fidelity term that ensures consistency with the measured projections p, where A is the forward projection matrix and D is a noise matrix, weighted with a parameter λ to balance the data fidelity and regularization terms.In the present work, D was set equal to the identity matrix.The second term regularizes the gradient-sparsified difference between the a priori information x prior and the target object x, and the third term regularizes the sparsity of the target object x.
The most commonly used sparsifying transform is the total variation (TV), which can be shown to be expressed through the volume gradient as follows (Lauzier et al., 2012), Water Resources Research 10.1029/2023WR036514 TEKSETH AND BREIBY TV is a regularization technique that works most efficiently for images consisting of large homogeneous regions separated by sharp boundaries.The objective function in Equation 1 can be optimized using a non-linear conjugate gradient algorithm.

Samples
Two samples were used in the experiments reported here.The first sample consisted of a glass cylinder of length 9 cm, 4.0 mm outer and 2.5 mm inner diameter filled with a sintered glass bead-pack extending 4.5 cm within the cylinder.The soda-lime glass beads, with a diameter range 250-500 μm, were sintered inside the glass cylinder in an oven at approximately 630°C for 1.5 hr.When soda-lime glass is exposed to X-ray radiation the glass is prone to darkening due to the creation of color-centers (Serrano et al., 2013).While the bead-pack sample was prone to darkening, no structural changes were observed in the experiments reported here.
The second sample consisted of an outcrop of Bentheimer sandstone (Peksa et al., 2015) that was cored to a diameter of 6 mm and a length of 7 mm.The sandstone was rinsed in methanol and left to dry for 3 weeks prior to the experiment.A macro-porosity of 24% ± 2% was estimated for this particular sample based on voxel counting, corresponding to a pore volume of approximately 42 ± 4 μL within the field of view.The sandstone was placed in a Viton sleeve and installed in an RSS Micro core holder (RS Systems, Trondheim, Norway) which is equipped with one inlet at the bottom and one outlet at the top.The core holder was made of polyether ether ketone (PEEK) to ensure high X-ray transparency.To prevent the fluids from permeating between the sandstone sample and the Viton sleeve, a confining pressure of 7 bar with Argon gas was applied to the Viton sleeve, squeezing the sleeve tightly around the sandstone.

Experimental Procedure
The samples were measured using a Nikon XT H 225 ST CT instrument equipped with a Perkin Elmer 1620 AN CS detector with a CsI scintillator.The CT instrument was used as purchased without modifications.To facilitate faster acquisitions, the detector was operated in a 2 × 2 pixel binning mode, yielding 1,000 × 1,000 detector pixel values per projection.The X-ray source was a tungsten anode in reflection geometry with a maximum possible Xray energy of 225 keV.The dynamic CT measurements utilized the in-built continuous scanning mode to maximize the measurement time and to reduce both sample accelerations and the motor positioning overhead time associated with step-and-shoot acquisition schemes.Consequently, each projection thus represents an integration over an angular range.The data sets containing the fluid dynamics were acquired with a batch script triggering a sequence of identical CT measurements.
Thereafter, the syringe pump was triggered to start injecting potassium iodide (KI)-doped water into the sample at a constant flow rate.The capillary and Bond numbers are defined by Ca = μv D /σ and Bo = Δρga 2 /σ, respectively, with ρ being the doped water density (Rumble, 2014), v D the Darcy velocity, σ the surface tension (Ali et al., 2009), g = 9.81 m/s 2 the gravitational acceleration, and a the average pore radius.The low capillary number ensures that we are in a regime dominated by capillary forces.A series of fast CT measurements followed.After each full rotation CT scan, additional overhead time of about 5-6 s had to be added for data transfer.The high-resolution scan and the scans during fluid dynamics were measured at the same acceleration voltage and magnification to ensure similar contrast.The effective voxel size was defined by the geometrical magnification.The experimental parameters for the two experiments are listed in Table 1.

CT Measurements of the Glass Bead-Pack
An a priori "dry" CT measurement of the dry glass bead-pack sample without any liquid was conducted first, to be used as a priori information.Thereafter, the syringe pump was triggered to start injecting potassium iodide (KI)doped water into the bead-pack.A series of fast CT measurements followed (Tekseth & Breiby, 2024b).The highresolution scan and the scans during fluid dynamics were measured at the same acceleration voltage and magnification to ensure similar contrast.

CT Measurements of Bentheimer Sandstone
An a priori scan of the dry Bentheimer sandstone was acquired first, followed by flooding the sample with KI doped water resulting in a partially saturated sample.After flooding, the sample was drained while a series of CT scans were collected (Tekseth & Breiby, 2024b).The low flow rate ensured a capillary-dominated flow regime.

CT Reconstruction and Post Processing
The dry sample scans were reconstructed using the Feldkamp-Davis-Kress (FDK) algorithm combined with a Ram-Lak filter (Feldkamp et al., 1984).The dynamic scans were reconstructed by optimization of the objective function in Equation 1 using a non-linear conjugate gradient algorithm combined with a Polak-Ribière update scheme and a Moré-Thuente line search algorithm (Dunlavy et al., 2010;Tekseth & Breiby, 2024a).Both the forward-and back-projection operators from the TIGRE toolbox (Biguri et al., 2016(Biguri et al., , 2019) ) and the gradient of the TV norm were GPU accelerated for increased reconstruction speed.The reconstructions were performed on a workstation (featuring an AMD Ryzen 9 3900X 12-Core processor and 128 GB of memory) and a GPU card (NVIDIA RTX 3090 with 24 GB of memory).Each function evaluation of Equation 1 had a computation time of approximately 13 s with a full reconstruction lasting approximately 15-20 min.All image processing steps were computed in Matlab TM (The MathWorks Inc., USA), whereas the 3D renderings were visualized with Paraview TM (Lin et al., 2018).
The reconstructed volumes were post-processed using a non-local means filter that efficiently removes noise present in the volumes while preserving edges (Bruns et al., 2017;Buades et al., 2005).The high-resolution sample models were segmented using a manually set threshold value to separate the static pore volumes from the matrix.In the dynamic data sets, the regions corresponding to the static porous matrix were masked out using the segmented high-resolution model, effectively reducing the number of phases present in the dynamic data sets to air and water.These two phases were then segmented using a manually set threshold value.The porosity was computed by counting the number of pore space voxels divided by the total sample volume in the segmented dry volume.Similarly, the saturation S w was defined as the number of voxels identified as "water" in the segmented dynamic scan divided by the total pore volume.The pore space was partitioned into pore bodies and pore throats using pnextract (Dong & Blunt, 2009).A non-linear least-squares method was used to fit a power-law to the pore- TEKSETH AND BREIBY filling and volume-filling events with the uncertainty given as the 95% confidence bounds.A pore body was considered filled if more than 50% of its voxels were classified as "water." We also performed analysis directly on the CT projections by comparisons between subsequent CT scans, somewhat inspired by Armstrong et al. (2014).By subtraction, the difference images of projections obtained 360°a part were obtained.In regions where a saturation change had occurred (water had invaded or retracted), the difference image will be either bright or dark.Counting the number of pixels within such regions provides an estimate of the size of the dynamic event.As noted by Armstrong et al. by segmenting the image before summarizing reduces the noise level within each projection image, making it easier to capture these dynamic events (Armstrong et al., 2014).

A Custom Sample Stage for 4D-CT of Multiphase Flow in Porous Media
A custom sample stage consisting of two syringe pumps, a pressure transducer and a temperature sensor was constructed and used for these experiments, see Figure 1.Both syringe pumps were based on a Newport MFA-PP motorized linear translation stages with a travel range of 25 mm and combined with gas-tight highprecision syringes from Hamilton Company (Nevada, USA).The flow rate ranges associated with the syringe pumps are thus dependent on the cross-sectional area of the syringe used.Here, we used a 250 μL syringe (Hamilton, part #81122).The syringe pumps were controlled with a custom LabView (National Instruments, Texas) program.The liquid supply lines from the syringe pumps were connected to a sample mounting block equipped with a gauge pressure transducer and a temperature sensor.To limit expansion and contraction of the liquid supply lines during pore-scale dynamics, high-modulus PTFE tubes were used (Z.Sun & Santamarina, 2019).The pressure transducer (Honeywell 26PCAFA6G) has a pressure range from 0 to 1 psi (0-69 mbar).The T-type thermocouple (RS Pro #363-0266) has an operational temperature range of 75-250°C . Both sensors can easily be replaced to allow other pressure and temperature ranges with the M5 threads sealed with nitrile rubber O-rings.Measurement series of the two sensors were acquired with a Picotech TM TC-08 thermocouple datalogger operating at approximately 3.3 Hz.The sample stage was outfitted with an 18-circuit through-bore slip ring (SR-EST18, Dynamic Sealing Technologies Inc) to avoid disturbing the liquid supply lines during sample rotation, and to allow continuous unidirectional rotation of the sample stage without rewinding.

Water-Air Dynamics in a Glass Bead-Pack
Here we present the results of a primary (not pre-wetted) imbibition and drainage experiment with 0.5 M KI doped water as the wetting phase and air as the non-wetting phase in a glass bead-pack, with a time-resolution of 26.1 s. Figure 2 shows reconstructed cross sections of the high-resolution (fully sampled) dry bead-pack and of the undersampled wet bead-pack during the imbibition dynamics as obtained with the FDK and PICCS algorithms, respectively.Histograms of the filtered FDK and PICCS cross sections are displayed in Figures 2e and 2f.Note that the X-ray effective attenuation coefficients of doped water and the glass beads are comparable, approx.μ = 14.8 ± 0.8 mm 1 , and we thus relied on difference images between the dry and wet scans to distinguish Clearly, the salt-and-pepper noise seen in the FDK reconstruction caused by the few projections and short exposures significantly reduced the quality, even after filtering.In the reconstruction based on the PICCS algorithm, water and air can easily be discerned.
between the water and the glass.Because the highly undersampled dynamic data sets consisted of only 150 projections with short exposure times, the reconstructed tomographic sections obtained with the conventional FDK algorithm are plagued with reconstruction artifacts appearing as noise and consequently broad peaks in the cross-section histogram.These artifacts are significantly reduced when using the PICCS algorithm with the TV regularization terms in the objective functions, Equation 1, as evidenced by the smooth cross section image and the sharp peaks within the histogram.The analysis clearly demonstrates the superior performance of the PICCS algorithm over FDK, even after the FDK reconstruction had been filtered.A slight tendency of beam hardening can be seen in Figure 2 as the profiles in 2e slightly sag, however, without causing complications for the analysis presented here.
A total of 175 full-revolution consecutive CT data sets were acquired and reconstructed during imbibition.Eight of these reconstructed data sets are visualized in 3D in Figure 3a showing a stable invasion front displacing air from all pores as the water intrudes into the sample.During the entire imbibition stage, only one single pore snapoff event (Blunt, 2017) was observed, albeit indirectly, as an air pocket could be seen surrounded by water after the imbibition front had passed (data not shown).Consequently, the invasion dynamics was dominated by cooperative pore filling events (Blunt, 2017) as the fluid front steadily advanced and got in contact with new dry pore throats.The temporal resolution is sufficient to image the states before and after these cooperative porefilling events, as demonstrated in Figure 4. Events may also be identified through the projection images as shown in Figure 4b, and in particular for Figure 4c, indicating that the event occurred in between the two CT scans shown in Figures 4a2 and 4a3.
As the water-air front necessarily only proceeded from pores and throats already filled with water, and the thoroughly dried sample thus had no pre-wetted layers throughout the porous medium, the process can be modeled as frontal advance (Blunt, 2017).By grouping the number of pore-filling events per timestep, the distribution of pores filled in each event can be deduced, limiting the analysis to those timesteps for which the entire water-air meniscus was contained within the field of view, cf. Figure 3b.According to invasion percolation arguments, the size of the bursts in terms of the pore volume filled is expected to follow a power law distribution, with an exponent of 1.125 in the case of 2D simulations (Martys et al., 1991).To our knowledge, a similar exponent value in 3D has not been reported.By counting the filled pores in each 3D volume, we find an exponent of 1.16 ± 0.19, as shown in Figure 3d.Similarly, for the distribution of burst sizes in terms of volume filled, as presented in Figure 3e, we find an exponent of 1.48 ± 0.39, which compares to the experimental value of 1.27 reported by Singh, Menke, et al. (2017) during imbibition in a Ketton limestone, which has a similar pore structure morphology as bead packs.The differences can be attributed to an insufficient number of measured events and possibly also to the limited lateral size of the sample system.A higher bead-diameter-tube-diameter ratio could reduce wall effects and thus provide a richer playground for effects such as snap-off events.Furthermore, it is well known that there are caveats associated with the numerical fitting of power-law distributions, including bias and inaccuracy, which should be considered in further studies, see for example, (Clauset et al., 2009).Still, being able to capture these dynamic events with an in-house instrument is beneficial as it can be used to verify pore scale models of flow in porous media.
We note with Armstrong et al. (2014) that significant information about the dynamics can also be obtained directly from the projections, giving a temporal resolution which is a factor 10-100 times faster than in the 4D reconstructions.An example is shown in Figures 4b and 4c, where the pore filling events are detected by detector pixels changing intensity.
Following a technical pause of approximately 20 min after the imbibition process, drainage of the injected water from the sample was initiated.In total 140 consecutive CT data sets were acquired, and selected 3D renderings of the drainage front are visualized in Figure 5a.The water phase almost completely saturated the bead-pack at the onset of drainage, owing to the limited presence of snap-off in this system.In contrast to the stable invasion front observed in the imbibition experiment, the intruding non-wetting air phase, having a lower viscosity and density than the liquid water phase, caused an unstable invasion front by capillary fingering (Lenormand et al., 1988).
In this drainage experiment, we were able to observe several significant Haines jumps before the sample reached its residual water saturation S w = 0.4 after approximately 35 min, see Figure 5b.Importantly, the water saturation and the water-air interface could still be observed to fluctuate even after the bulk air phase had passed the imaged field of view.These fluctuations are attributed to the long-ranging non-local effects of Haines jumps occurring below the imaged part of the sample, affecting the water-air meniscus throughout the sample (Armstrong & Water Resources Research 10.1029/2023WR036514 Berg, 2013;Tekseth et al., 2024).While capillary forces are traditionally considered to be a local force, recent studies have shown that the capillary may affect interfaces several pores away from the pore-filling event, depending primarily on the interfacial tension and viscosity.No intermittency was observed during these fluctuations.
As the drainage process proceeded, the displaced volume per timestep was observed to fluctuate with around 0.04 ± 0.02 μL per timestep.However, at certain timesteps the displaced volume per time was observed to exhibit maxima, as noted with arrows in Figure 5c.These peaks of the displaced volume are correlated with the largest Haines jumps observed.The increase in displaced volume for these specific timesteps can be related to the increase in interfacial area as the water is drained.Because more interfacial area becomes available as water is drained from the sample, more energy can be stored in these menisci by increasing their curvature.In particular for the highest peak observed around 1,650 s into the drainage cycle displaces 0.26 μL within two radiographic exposures of 268 ms, consistent with previously reported timescales of Haines jumps (Armstrong & Berg, 2013;Moebius & Or, 2012).3D renderings of the state before and after this jump are visualized in Figure 6.By investigating the projection images each event may be identified as shown in Figure 6b, see Section 3 for procedure.For the particular event in Figure 6a3, the event occurred within projection images #116 and #117 in the corresponding CT scan.
The drainage experiment presented herein can be modeled as an invasion percolation process (Wilkinson & Willemsen, 1983).Here, we find an approximate power-law relationship of pore-filling events with an exponent of 1.54 ± 0.14, as presented in Figure 5d.The value of the exponent matches closely with the values of 1.55 (Blunt, 2017;Strenski et al., 1991) and 1.63 observed by Roux and Guyon (1989), both through invasion percolation simulations.An approximate power-law distribution has also been found experimentally by Bultreys et al, however without presenting an exact value, in a drainage experiment with a Bentheimer sandstone sample using a gantry-based CT setup (Bultreys et al., 2015;Dierick et al., 2014).

Drainage of a Bentheimer Sandstone
Here we present the results obtained from extracting a 1.0 M KI doped water partially saturated Bentheimer sandstone sample.As both the sample diameter was larger and the doped water contained twice as much KI as in the bead-pack experiments, a higher number of projections and a longer exposure time were needed.The resulting temporal resolution in this experiment was approximately 50 s.Figure 7 shows reconstructed cross sections of the Bentheimer sandstone in the fully sampled dry state and the undersampled wet state as reconstructed with the FDK and PICCS algorithms.The wet data set presented was obtained just before air breakthrough.The noise imprinted in both the raw and filtered reconstructed cross sections obtained with the FDK algorithm is clearly visible in Figures 7b and 7c, whereas the presence of noise is reduced in the cross-section obtained with the (c) Zoomed-in plot of the shaded region in (b), also comprising the event captured in (a3) as marked with the red point, showing that this particular pore-filling event occurred within the overhead time in-between two consecutive CT scans.
PICCS algorithm.An estimate of the resolution along an interface profile is shown in Figure 8 given as the distance between the 25%-75% change in the gray value along the interface profile of the white line in Figure 7d (Donnelly et al., 2020).By averaging three interface profiles an estimated resolution of 10 ± 2 μm was obtained for the prior information, 21 ± 2 μm for the dynamic scan reconstructed with FDK and 19 ± 3 μm for the dynamic scan reconstructed with PICCS.Although the resolution improvement is limited, these line profiles clearly show the improved noise-handling of the PICCS algorithm.As time passed, the available water-air surface area increased, causing the pore filling events to increase in size.(d) Number of pore-filling events of a given size.A best-fit power law is displayed with an exponent of 1.54 ± 0.14.(e) Distribution of volume filling.An approximate power law behavior can also be observed for these volume fillings with an exponent of 1.6 ± 0.3.

Water Resources Research
10.1029/2023WR036514 The drainage dynamics in the Bentheimer sample is visualized in Figure 9a.The water is observed to be continuously drained, as shown in the water saturation profile in Figure 9b.However, the pressure signal presented in Figure 9c demonstrates that there were significant pressure fluctuations throughout the dynamics that is associated with Haines jumps (Berg et al., 2013;Haines, 1930).These pressure jumps continued to fluctuate until the air reached breakthrough 1,200 s into the experiment.Figure 10 shows two different views of the evolving airphase as it protruded into the pore-space during the drainage dynamics.We note that the number of pore-filling events was statistically insufficient to provide power-law exponents of these events.In addition, the displacement events observed within the temporal window of 50 s might also constitute one or more Haines jumps, such that multiple subsequent events would be counted as one.While only the drainage dynamics was considered here, this is not a limitation, and a natural extension of the current study could be experiments that include cyclic drainage and imbibition.

Discussion
The adaptations described here for carrying out dynamic CT experiments can be readily implemented in most commercial CT instruments having sufficient space for the complex sample environment.Many further refinements of the presented hardware can be foreseen, including installation of more sensitive pressure sensors combined with a higher frequency data acquisition.Also, a high-quality microphone for acoustical emissions (Moebius et al., 2012) would have provided complementary information to the CT data and could also be used for either triggering fast CT measurements, or as an aid in the reconstruction process.As an example of the latter, if a sufficient number of projections for reconstruction can be chosen before or after (but not across) a large redistribution event like a Haines jump, much smearing of the liquid interfaces in the reconstructed 3D timeframe could conceivably be avoided.Indeed, it is a limitation of our current implementation that the rapid dynamics during Haines jumps results in motion artifacts during reconstruction (Armstrong et al., 2014).In the future, it could also be possible to remedy this problem by actively utilizing sudden pressure changes of readings from the pressure sensor as a gating method to either initiate a scan or as input to the reconstruction algorithm.As indicated in Figures 4 and 6, it could also be possible to specify projection images used in a CT reconstruction after the acquisition by selecting radiographs located before and after the dynamic events as proposed by Armstrong et al. (2014).For the first ∼800 s, the pressure decreased until the local capillary pressure of the water-air meniscus exceeded the capillary entry pressure of the surface pore throats.From this point on, several abrupt changes in the pressure signal, associated with Haines jumps, can be observed, as indicated by arrows.The sudden pressure increase around 1,200 s signifies that the air phase had percolated through the sandstone. Water TEKSETH AND BREIBY The use of a priori information is an obvious and increasingly exploited method for gaining satisfactory reconstructions from undersampled measurements.CT is in this connection one of many instances under the wide umbrella of computational imaging, where the traditional high-performance imaging optics is gradually being redesigned to incorporate active computers as an integral part of the image-forming process.By defining the purpose of the optical setup as providing coded information about the sample, suitably trained software can decode and thereby obtain information that goes well beyond high-resolution imagery (Liebi et al., 2015;Pfeiffer, 2018;Schaff et al., 2015;Wakonig et al., 2019;Zheng et al., 2013).Limiting the scope to CT, we find it natural to expect that future experimental setups with dedicated software will be capable of retrieving information that is currently not accessible without lengthy post-analysis, like detailed local mineralogy, pore surface roughness, and power law exponents related to self-organized critical dynamics (Bak et al., 1987).We note that with diffraction contrast, additional features can be extracted, as recently demonstrated both with X-rays (Mürer et al., 2021) and neutrons (Woracek et al., 2018).
In this article, we have focused on sparsity in the spatial domain using the PICCS algorithm originally developed for medical CT (G.H. Chen et al., 2008).However, as mentioned, there are also efforts reported on exploiting sparsity in the time domain (Goethals et al., 2022;Kaestner et al., 2011).Specifically, during multiphase fluid flow in porous media, at a given pore location, and in the absence of turbulence, the fluctuations between the various phases can be low, thus allowing information from preceding (and/or succeeding) time steps to help constraining the reconstruction.We also mentioned the existence of discrete tomography algorithms like DART (Batenburg & Sijbers, 2011;Zhuge et al., 2016) which with the aim of faster and/or better reconstructions, exploit the presence of only a few phases in certain sample systems.Our experience (not shown) is that these algorithms are sensitive to small density offsets caused by beam hardening and by erroneous classification of intermediate phases across boundaries between high and low attenuation phases.We still consider that it should be possible to combine PICCS and DART.
A limitation of our current implementation is that the iterative reconstructions take a long time to execute due to the many data sets often associated with dynamic CT.Despite our use of GPUs, there is certainly a speed-up to be gained from optimizing the code toward high-performance computing.Moreover, we also suspect that there could be considerable room for reducing the computation load by directing the reconstruction efforts to regions within the sample that de facto exhibit significant changes at the given time step (Heyndrickx et al., 2020).As briefly discussed in the Introduction, machine learning is under rapid development and will prove a highly important tool also for measuring and understanding transport in porous media in the coming years.Once trained, computational neural networks also have the advantage of delivering very fast non-iterative reconstructions.The possibility of extracting wetting angles from pore-scale views of the fluid transport in porous media has gathered much scientific attention because of its importance for understanding the local pore physiochemical conditions (AlRatrout et al., 2017;Blunt et al., 2019;Khanamiri et al., 2020;Mascini et al., 2020;C. Sun et al., 2020).A targeted study toward the wetting angle dynamics would put stringent requirements on the spatial resolution in order to resolve the three-phase contact line accurately.Consequently, a lower temporal resolution would be obtained as either more projections or a higher exposure time would be needed to increase the spatial resolution.
Similarly, from a physics point of view, the observation of power laws has far-reaching implications, and it is satisfying that our setup allows us to estimate the critical exponents for drainage and imbibition in porous media.Still, while it is important to realize that the lateral cross section of the samples studied here is highly limited and perhaps giving artifacts, and the number of events recorded is too low for a careful statistical analysis, the potential for designing such experiments in home-laboratories is certainly demonstrated.

Conclusions
Time-resolved 4D-CT imaging of multiphase fluid flow in porous media has primarily been conducted at synchrotrons, because of the availability of orders of magnitude higher photon flux as compared to laboratory-based CT instruments.As synchrotron access is a scarce resource, methods for conducting two-phase flow experiments with conventional CT equipment are in high demand.In this article, we have outlined our efforts on both hardware and software toward home-laboratory experiments.An experimental stage outfitted with syringe pumps, a pressure transducer and a thermocouple placed onto the rotational stage, combined with a slip ring, enables fast acquisitions without rewinding between subsequent CT scans.Due to the limited photon budget associated with short exposure times, we have implemented a reconstruction scheme relying heavily on the PICCS algorithm to incorporate a priori information, specifically exploiting a high-resolution scan with information about the static regions of the sample.With all these efforts, we were able to achieve a sub-minute temporal resolution and sufficient spatial resolution to resolve the water-air interface.
The developed framework was used to study two-phase flow experiments in a synthetic glass bead-pack sample and in a Bentheimer sandstone.We demonstrated, qualitatively, superior image quality compared to the FDK algorithm.In the experiment with the bead-pack sample we tracked the stable invasion front dynamics during imbibition and the capillary fingering flow regime during drainage with several Haines jumps.Furthermore, during imbibition we found a power-law behavior in the pore-filling events.The burst size distribution, however, was not reminiscent of a power-law behavior.In the drainage experiment we observed a power-law behavior in both the distribution of pore-filling events and the burst sizes, corroborating the exponents reported from both invasion percolation models and through experiments.The drainage experiment with the Bentheimer sandstone showed several pore-filling events, however, the experiment had to be terminated because of an early percolation breakthrough, and a pore-filling event distribution could therefore not be obtained.By sharing these experiences and insights, along with the technical details of the sample environment developed for a conventional CT scanner, we believe that this study can help facilitate entering the increasingly important experimental technique of timeresolved computed tomography, including two-phase flow experiments.The framework can be used to study transient flows, capturing ganglion dynamics, snap-off and Haines jumps.Furthermore, it can be used to study intermittency behavior observed in steady-state liquid-flow experiments.

Figure 1 .
Figure 1.Sketch of the sample stage with a glass bead-pack sample.All equipment necessary to realize the controlled liquid flow was placed onto the sample stage to avoid disturbing the liquid supply lines.A slip ring allowed unidirectional rotation without rewinding.(a) Perspective drawing of the experimental sample environment placed onto the rotational stage.Sample is not to scale.(b) Orthogonal sketch of the experimental setup.

Figure 2 .
Figure 2. Comparison of tomographic cross-sections obtained with the FDK and PICCS algorithms for the glass bead sample.(a) Reconstruction of the high-resolution a priori scan, and a dynamic scan using (b) FDK, (c) FDK with non-local means filtering, and (d) PICCS.Density histograms of (e) non-local means filtered FDK reconstruction and (f) PICCS reconstruction, showing significantly sharper peaks for the latter.(g) Line profiles corresponding to the colored horizontal lines in (a)-(d).Clearly, the salt-and-pepper noise seen in the FDK reconstruction caused by the few projections and short exposures significantly reduced the quality, even after filtering.In the reconstruction based on the PICCS algorithm, water and air can easily be discerned.

Figure 3 .
Figure 3. Dynamics of water imbibing into a glass bead-pack.(a1-a8) 3D perspective renderings of selected timesteps separated by 4 min 22 s as water is imbibed.The bead-pack and the confining cylinder are rendered transparent for clarity.(b) The liquid saturation, increasing linearly with time.The unshaded region highlights the points where the entire water-air interface was within the field of view.(c) The change in measured volume per time step (Δt = 26 s).(d) Distribution of the number of pores filled per event.A non-linear least square fitting procedure provided a best fit power law with an exponent of 1.16 ± 0.19.(e) Distribution of volume fillings, with a power law fit having an exponent of 1.48 ± 0.39.

Figure 4 .
Figure 4. (a1-a4) 3D renderings of four consecutive reconstructed data sets during imbibition of the bead-pack sample.The temporal resolution is sufficient to image the before and after a cooperative pore-filling event, as marked by the arrow.(b) Number of detector pixels filling with water per projection step with the 50 largest events marked with red points.(c)Zoomed-in plot of the shaded region in (b), also comprising the event captured in (a3) as marked with the red point, showing that this particular pore-filling event occurred within the overhead time in-between two consecutive CT scans.

Figure 5 .
Figure 5. Dynamics during water drainage of the bead-pack sample.(a1-a8) 3D visualizations of selected timesteps as water was drained from the bead-pack.Note the increasingly irregular water-air interface with larger surface area and increasingly complex topology with time because of fingering.(b) Water saturation as a function of time, decreasing approximately linearly in the beginning and reaching a residual water saturation value of 0.4 about 2,200 s into the experiment.The unshaded region (400-2,200 s) highlights the drainage part within the field of view.The remaining experiment (2,200 s and onwards) exhibited saturation fluctuations, interpreted as non-local coordinated effects of pore-filling events occurring outside the imaged field of view.(c) Measured change of volume per time step during drainage.As time passed, the available water-air surface area increased, causing the pore filling events to increase in size.(d) Number of pore-filling events of a given size.A best-fit power law is displayed with an exponent of 1.54 ± 0.14.(e) Distribution of volume filling.An approximate power law behavior can also be observed for these volume fillings with an exponent of 1.6 ± 0.3.

Figure 6 .
Figure 6.3D renderings of four consecutive reconstructed data sets separated by ∼26 s during drainage of the bead-pack sample.Here the states before and after a Haines jump event is visualized, as marked by the arrow, occurring during the acquisition of the CT measurement visualized in (a3).The remaining images show that fluid configuration is practically stationary before and after the Haines jump.This particular jump was the biggest observed jump displacing 0.26 μL.(b) Identified porefilling events during drainage obtained from the projection images.The 50 biggest observed events are marked with red points.(c) The gray region of interest marked in (b) corresponds to the projections used in the reconstructions of (a1)-(a4).The sharp point marks the Haines jump event presented in (a3), occurring in the middle of a CT scan.

Figure 7 .
Figure 7.Comparison of tomographic cross-sections obtained with the FDK and PICCS algorithms for the Bentheimer sandstone.(a) Reconstruction of the high-resolution a priori scan, exemplified with a cross section of the dry Bentheimer sample showing the porous structure used as input to the PICCS algorithm.Reconstructed cross-sections of a selected timestep during the dynamics with (b) FDK, (c) FDK and non-local means filtering, and (d) PICCS.The white line in (d) marks the interface profile used in Figure 8. (e) Line profile from the horizontal lines in (a)-(d), showing the superior performance of the PICCS algorithm over FDK.

Figure 8 .
Figure 8. Resolution estimate for the high-quality prior information and the dynamic scan as reconstructed with the FDK followed by non-local means filtering, and PICCS algorithm along an interface.

Figure 9 .
Figure9.Drainage dynamics in a Bentheimer sandstone.(a1-a8) 3D renderings of selected timesteps during drainage with the sandstone rendered semi-transparent.(b) Changes in the water saturation.Note that the invading air phase percolated through the system already after ∼1,200 s, and the experiment was then stopped.(c) Pressure at the inlet as a function of time.For the first ∼800 s, the pressure decreased until the local capillary pressure of the water-air meniscus exceeded the capillary entry pressure of the surface pore throats.From this point on, several abrupt changes in the pressure signal, associated with Haines jumps, can be observed, as indicated by arrows.The sudden pressure increase around 1,200 s signifies that the air phase had percolated through the sandstone.

Figure 10 .
Figure 10.Displacements during drainage in the Bentheimer sandstone with a coloring scheme indicating the timedevelopment.(a and b) Displacements pattern at two different views.The sandstone has been rendered transparent.

Table 1
Parameters Used for the Two Experiments