Pore‐Scale Visualization and Quantification of Transient Solute Transport Using Fast Microcomputed Tomography

Solute transport is important in a variety of applications regarding flow in porous media, such as contaminant groundwater remediation. Most recent experimental studies on this process focus on field‐scale or centimeter‐scale data. However, solute spreading and mixing are strongly influenced by pore‐scale heterogeneity. To study this, we developed a novel methodology to quantify transient solute concentration fields at the pore scale using fast laboratory‐based microcomputed tomography. Tracer injection experiments in samples with different degrees of pore‐scale heterogeneity (porous sintered glass and Bentheimer sandstone) were imaged in 3D by continuous scanning at a time resolution of 15 s and a spatial resolution of 13.4 μm. While our calibration experiments indicated a high uncertainty (1σ) on the concentration in single voxels due to imaging noise (± 27% of the total concentration range), we show that coarse gridding these values per individual pore significantly lowers the uncertainty (± 1.2%). The resulting pore‐based tracer concentrations were used to characterize the transport by calculating the solute's arrival time and transient (filling) time in each pore. The average velocities estimated from the arrival times correspond well to the interstitial velocities calculated from the flow rate. This suggests that the temporal resolution of the experiment was sufficient. Finally, the pore‐based transient filling times, the global concentration moment and the global scalar dissipation rate calculated from our experiments, indicated more dispersion in the sandstone sample than in the more homogeneous sintered glass. The developed method can thus provide more insight in the influence of pore‐scale heterogeneity on solute transport.


Introduction
An understanding of solute transport in porous media is key for various applications in engineered and natural porous media, such as building stone performance (Rodriguez-Navarro et al., 2000), remediation of contaminant groundwater (Bekins et al., 2001), and waste management (Ghoraba et al., 2013). Most of the recent experimental studies on solute transport processes focus on field-scale or centimeter-scale data (e.g., Frippiat & Holeyman, 2008;Koestel & Larsbo, 2014). However, solute transport is significantly impacted by pore-scale heterogeneity found in natural porous media . Gathering direct data on pore-scale transient solute concentration fields is complicated by the high temporal and spatial resolutions that are required. Solute transport is primarily controlled by solute spreading and mixing. Mixing refers to the smoothening out of concentration gradients by diffusion, therefore increasing the volume occupied by the solute plume. This phenomenon describes dilution (Bear, 1972;Boon et al., 2017). Mixing processes can be characterized by the dilution index and the scalar dissipation rate (Bolster et al., 2011;Kitanidis, 1994). More specifically, the scalar dissipation rate expresses the decrease of concentration variance (Bolster et al., 2011). Spreading or mechanical dispersion changes the solute plume shape without changing the volume of the tracer plume. It is linked to a nonuniform (pore-scale) velocity field and tortuosity, due to heterogeneous pore space geometry and friction variations (Bear, 1972;Sahimi, 2011). Hydrodynamic dispersion represents the combination of solute mixing and spreading (Bear, 1972;Boon et al., 2016).
To describe flow at the continuum scale, the traditional Advection Dispersion Equation (ADE) or the continuous time random walk formalism are frequently used (Berkowitz et al., 2006). The ADE is a standard approach based on Fick's law and the conservation of mass (Anderson & Cherry, 1979). It uses an average linear velocity for a representative elementary volume, assumes Fickian processes, and uses one hydrodynamic dispersion coefficient in which the effects of solute mixing and spreading are lumped together, assuming perfectly mixed conditions within a representative elementary volume (Anderson & Cherry, 1979). For early times, when the velocity field is not fully sampled, incomplete spreading and mixing at the pore scale exist and the ADE cannot be used (Alhashmi et al., 2015;Bijeljic & Blunt, 2007). Continuous time random walk formulations are believed to be more applicable for such non-Fickian or anomalous processes (Berkowitz et al., 2006). The asymptotic length depends on the transport properties in the rock, which is influenced by pore-scale heterogeneity and flow velocity. This influences the relative effect of advection and diffusion that is characterized by the dimensionless Péclet number Pe (Bijeljic & Blunt, 2007;Péclet, 1827): where v is the interstitial velocity (m/s), L is the characteristic length (m) and D m is the molecular diffusion coefficient (m 2 /s).
For laminar flow, the dependence of transverse hydrodynamic dispersion on the Péclet number is classified in different regimes (Bijeljic & Blunt, 2007;Sahimi, 2011): 1. Restricted diffusion for Pe < 0.3 where mechanical dispersion is so low that molecular diffusion controls the transport processes. 2. Transition for 0.3 < Pe < 5 where mechanical dispersion and diffusion both contribute to hydrodynamic dispersion. 3. Power law for 5 < Pe < 300 where mechanical dispersion dominates in hydrodynamic dispersion, although the effect of diffusion cannot be neglected. 4. Mechanical dispersion for Pe > 300 where the transport processes are defined by mechanical dispersion.
To get a better understanding of solute transport in porous media, numerous experimental studies to measure macroscopic dispersion coefficients have been performed (Frippiat & Holeyman, 2008;Rolle et al., 2009;Ye et al., 2015bYe et al., , 2015a. Koestel and Larsbo (2014) quantified preferential solute transport in soil macropores by microcomputed tomography (micro-CT). Kurotori et al. (2019) and Zahasky et al. (2019) used positron emission tomography (PET) for in situ characterization of fluid/solute transport in geologic media at subcore scales. Boon et al. (2017) use three-dimensional (3D) experimental computed tomography data of solute transport at the centimeter scale to characterize the impact of rock heterogeneity on solute spreading and mixing. However, structural heterogeneity at the pore-scale is assumed to play a major role in solute spreading and mixing . This complicates the definition of the asymptotic regime and upscaling of dispersion coefficients (Bijeljic & Blunt, 2007). Therefore, simulations and experiments investigating the evolution of pore-scale solute concentration fields in consolidated rock materials are very valuable. Optical techniques have been used to evaluate the intermittency of Langrangian velocity and acceleration in transparent porous media in great detail (Holzner et al., 2015). Also mixing and reaction kinetics were quantified on the two-dimensional pore scale by measuring intrapore, conservative concentration fields of a fluorescent tracer with the soft lithography technique (De Anna et al., 2014). Direct micron-scale visualization and quantification of solute concentration fields during spreading and mixing in rock samples have however remained difficult. It is complicated by the high spatial and time resolutions that are required. For example, to capture the concentration change over a characteristic length of 250 μm (corresponding to a coarse grained sandstone) at a Pe of 1, an image has to be taken at least every 35 s. Currently, micro-CT is likely the only technique that can be used for in situ imaging experiments of rock samples at these time resolutions while achieving spatial resolutions in the micrometer range. Bultreys et al. (2016) and Boone et al. (2016) suggested that fast (12-15 s time resolution) laboratory-based micro-CT can be used to image tracer dispersion experiments. Where typical laboratory-based micro-CT scans are taken in about 15 min to 24 h (Wildenschild & Sheppard, 2013), a time resolution up to 15 s was reached in these studies. However, fast lab-based micro-CT typically suffers from low signal-tonoise ratio, due to the short acquisition times (Heyndrickx et al., 2018). This scales with the square root of the number of photons that are detected, thus requiring higher photon fluxes from the X-ray source to maintain the same signal-to-noise ratio. For X-ray tubes used in micro-CT scanners, this would in turn affect the spatial resolution of the images. Synchrotron radiation sources circumvent this problem by providing much higher X-ray fluxes, yet are much less accessible (Berg et al., 2013;Singh et al., 2017).
Due to the low signal-to-noise ratios in laboratory-based micro-CT, this study focuses on developing a method for enhanced image analysis. This is key for the quantitative analysis of the concentration fields in such data and consequently the interpretation of the spreading and mixing behavior. We show the application of this analysis method on two newly acquired experimental data sets (based on the methodology described in Bultreys et al., 2016), taking specific care to quantify the experimental uncertainty and the temporal evolution of micro-CT concentration fields to investigate solute spreading and mixing. First, we performed calibration experiments to investigate the measurement errors on the concentration fields that were extracted from the micro-CT images. Then, fast laboratory-based micro-CT was used to image tracer injection experiments on two different samples under two different flow velocities. The tracer concentration over time in individual pores was measured by coarse gridding the voxel-based concentration images per pore bodies. To analyze the spatial and temporal variation of the tracer concentration, we defined the arrival time and transient (filling duration) time of individual pores. These can be linked to pore-scale velocity distributions and hydrodynamic dispersion. The global concentration moment and scalar dissipation rate were calculated over time for the entire sample. This provided additional information about the global spreading and mixing mechanisms.

Sample Selection and Flow Setup
To facilitate detailed analyses and quantification, we selected two homogeneous porous sintered glass samples for the calibration experiments and the first tracer injection experiment. These two samples had fairly large pores sized between 160 and 250 μm (ROBU P0, Germany), making the micro-CT images more straightforward to analyze. The second injection experiment was performed on Bentheimer sandstone, which has a pore size distribution from 50 to 200 μm (Peksa et al., 2015). It is known for its fairly homogeneous and highly permeable pore structure (Peksa et al., 2015). All samples had a cylindrical core shape with a diameter of 6 mm and a length of 20 mm.
Each sample was placed in a custom-built Hassler-type polymethylmethacrylate flow cell that fits on the micro-CT scanner ( Figure 1). A Viton sleeve was pressed around the sample by a confining pressure of 0.6 MPa. This avoided that the injected aqueous tracer solution bypassed the sample during the tracer injection experiments. CsCl (10 wt%) was chosen as tracer, because of its high solubility in water. The Cs-ions cause the solution to have a high X-ray attenuation coefficient during micro-CT scanning, which is expected to rise approximately linearly with the Cs concentration (Agbogun et al., 2013). A Harvard Apparatus PHD Ultra pump controlled the flow through the sample.

Laboratory-Based Micro-CT
A laboratory-based micro-CT system was used to visualize the solute concentration distribution. Micro-CT is a technique to acquire 3D images of a sample in a nondestructive manner. Therefore, it is a valuable technique to investigate internal structures and dynamic processes in materials (Cnudde & Boone, 2013;Ketcham, 2005;Noiriel, 2015). With micro-CT, a number of radiographic images or projections are taken at various rotation angles of the sample (Ketcham & Carlson, 2001;Wildenschild & Sheppard, 2013). When X-rays pass through the object, X-ray scattering and absorption attenuate the signal (Jackson & Hawkes, 1981). The resulting X-ray attenuation is mainly controlled by the X-ray energy, and the material's density and atomic number (McCullough, 1975;Wildenschild & Sheppard, 2013). For a mono-energetic beam through a homogeneous material, Lambert-Beer's law relates the intensity I of X-ray photons passing through the object with the initial X-ray intensity I 0 , the attenuating material's thickness along the beam d, and the linear attenuation coefficient of the object μ: Equation (2) is only valid for monochromatic X-rays that follow a straight path. In laboratory-based micro-CT, polychromatic X-ray beams are used, causing image artifacts such as beam hardening in the reconstructed images (Cnudde & Boone, 2013;Ketcham & Carlson, 2001).
Micro-CT images consist of voxels, each of which has a grey value corresponding to the local X-ray attenuation coefficient in the sample, therefore allowing to visualize differences in density and composition within its internal structure.
For the experiments presented here, we used a gantry-based micro-CT system specially designed for in situ imaging, with a static sample and a rotating X-ray source-detector system in a horizontal plane Dierick et al., 2014). The system is able to gather fast scans back-to-back with a minimum acquisition time of 12 s per scan . To estimate the uncertainty of concentration fields derived from these images, calibration scans where the pores were saturated with a known CsCl concentration were performed at 12 s per scan. The CsCl tracer injection experiments were imaged at 15 s per scan. Each fast imaging experiment was preceded by a higher quality slower pre-scan to determine the static pore structure. The detailed settings for all scan series are provided in Table 1.
The acquired radiographs were reconstructed to a 3D volume using Acquila (TESCAN-XRE, Belgium). The higher-quality prescans were segmented to extract the pore space by manual thresholding (Avizo 9.5.0, ThermoFisher Scientific). This segmented image was used as a mask to extract the pore space from the fast scans . Because of the high X-ray attenuation coefficient of the dissolved Cs ions, higher tracer concentrations lead to higher grey values in the pore space. Other studies show that the link between grey value and tracer concentration can be assumed to be linear (Agbogun et al., 2013). Therefore, relative concentrations can be calculated based on voxel-based grey values (Zhang et al., 2018): where GV tracer represents the grey value of a pore voxel at a certain time in the tracer injection experiment. GV water and GV 10wt % CsCl are the grey values in the same pore voxel when the sample was, respectively, saturated with water and 10 wt% CsCl.

Calibration Experiments
Calibration experiments were carried out to investigate the quantitative correctness of the concentration fields extracted from the fast scans of the tracer injection experiments. First, a higher-quality micro-CT image of the water-saturated sintered glass sample was acquired. Then, the sample was scanned while it was saturated with different constant CsCl concentrations (0 to 10 wt % CsCl in steps of 2 wt%). While the conditions in the sample remain constant, eight fast scans (eight repeat observations) were taken to assess the repeatability of the grey value measured in each voxel. Although motion artifacts related to the movement of the solute are not taken into account, these scans allow to assess random and systematic errors on the determination of the tracer concentration during the injection experiments, caused by image noise and other artifacts.
The average grey value of the pore volume in the micro-CT image was calculated and plotted against the concentrations. Since the solute concentration was kept constant, the voxel-based grey values should theoretically be constant over the eight repeat scans. However, random noise causes the measured grey value in a voxel to be distributed around this true value. To assess the measurement uncertainty, we estimated the width of this distribution by calculating the standard deviation of each voxel's grey value over the eight repeats (equation (4)): SD is the voxel-based standard deviation over the eight repeat experiments. In a specified voxel i, x it is the grey value in repeat scan t and x is the average grey value in the voxel over the repeat experiments. The average of this value over all pore space voxels is a measure of the uncertainty on the micro-CT based concentration field due to random image noise.
To alleviate random noise, the concentration field can be coarsened by averaging multiple voxels together. To this end, the pore space image was decomposed into individual pore bodies using a pore network extraction algorithm based on a watershed of the pore space's distance map (Raeini et al., 2017). The decomposition of the pore space is visualized in Figure 2, which shows individual pores in different (randomly selected) colors. While this sacrificed sub-pore-scale information by effectively reducing the resolution of the images, the fact that there is a conceptual basis for separating pores at local regions with low hydraulic conductance (Raoof & Hassanizadeh, 2012) means that the method may to an extent aid the interpretation of the experiments. Note that any method that provides a decomposition of the pore space into small parts can be used instead of the employed pore network extraction method, as it is only used here to provide a coarse gridding of the pore space (and thus not for any modeling purposes). The average grey value of the (voxels in a) pore body indicates the pore-averaged tracer concentration in each scan. The time-based standard deviation of this value indicates the pore-averaged concentration uncertainty.

Tracer Injection Experiments 2.4.1. Flow Regimes
For the tracer injection experiments, the sintered glass and Bentheimer samples were first saturated with demineralized water (after CO 2 sparging to remove air). A higher-quality micro-CT image was taken, after which a solution with 10 wt% CsCl was pumped in from below with a constant flow velocity while 40 to 50 fast micro-CT scans were acquired. The injection experiments were carried out with volumetric flow rates of 0.25 and 0.50 μL/s. These flow velocities were estimated to result in different regimes, determined by the Péclet number (equation (1), Boon et al., 2017;Dentz et al., 2011;Péclet, 1827). The characteristic length L was taken to be the pore size (or the average pore diameter based on the network extraction): 1.5 × 10 −4 m for the sintered glass sample and 0.8 × 10 −4 m for the Bentheimer sandstone sample. The molecular diffusion coefficient D m of Cs + at room temperature is 1.77 × 10 −9 m 2 /s (Li & Gregory, 1974). The resulting Péclet numbers (Table 2) indicate that the low and high volumetric flow velocity are, respectively, expected to be in the transition and (close to the) power law regime.
It should be noted that the maximum time resolution that can be achieved imposes an upper limit on the Péclet number for which we can expect the transient concentration profile to be reliably imaged. For example, in a porous medium with a characteristic length of 250 μm, the maximum Péclet number for which the solute movement within 12 s would be smaller than the characteristic length is approximately 3 (though useful data may still be obtained for higher Péclet numbers due to the sample's tortuosity, which increases the solute travel time). There is no fundamental lower limit to the Péclet numbers that can be probed with the method.

Pore-Scale Concentration Fields
The voxel-based and pore-based relative concentrations during tracer injection are calculated similar to the calibration experiments (equation (3)). For the sintered glass sample, network extraction resulted in a decomposition of the pore space into 1,402 pore bodies. In the Bentheimer sample, 2,700 pore bodies were found. The average grey values for each pore body were converted to relative concentrations (equation (3)), allowing to evaluate the concentration change over time. The resulting pore-scale concentration evolution for three randomly selected pore bodies is plotted against time in Figure 2.
To reduce the amount of data into a monotonously increasing relative concentration curve for each pore body, a sigmoidal function was fitted for each pore in Octave (GNU Octave): where A and B are fitting parameters, x is the time in seconds, and the output of the provided function f(x) is the relative concentration (ranging from 0 to 1). The fit seems to be less reliable at the outer parts of the curve ( Figure 3); however, the middle part seems to be fitted accurately. The mean R 2 value was calculated for the fit in each pore. For the sintered glass experiments this was 0.998, for Bentheimer sandstone 0.996, indicating a very good fit for the data sets. We used this fit to characterize the concentration curve for each pore body by defining arbitrary the arrival time of the solute in the pore body as the moment when it reaches a relative concentration of 0.1. The filling duration was defined as the time it takes the relative concentration to rise from 0.1 to 0.9 (Figure 3).

Global Concentration Moment and Scalar Dissipation Rate
To measure global solute mixing in the continuous tracer injection experiments, we determined the decline of concentration variance by using the scalar dissipation rate (Pope, 2000): This can be calculated from the concentration moment (Bolster et al., 2011): Because the pore-based concentration fields determined from the experimental data are more reliable than the voxel-based concentrations, we used discretized equations weighted by the pore volume and summed over all the pore bodies: M j (t) is the concentration moment at time j, c i (t j ) is the concentration in pore i at time j, and v i is the volume of pore i in m 3 . X j (t) is the scalar dissipation rate at time j, calculated as a finite difference approximation (forward difference scheme using the diff function in Matlab).

Calibration Experiments
For the reconstructed micro-CT scans of the eight repeat observations, the voxel size is 13 μm. While there is no single measure that fully characterizes the actual resolution of an experimental image in general, the Fourier ring correlation has been proposed as a fully automatic, quantitative image-based measure without the need for prior information (Nieuwenhuizen et al., 2013). We calculated this measure using a plugin available in Fiji (ImageJ), yielding a Fourier ring correlation resolution of 55 μm. The average grey value of the pore space was plotted against the known relative tracer concentration in the calibration experiment ( Figure 4). While the average grey value shows a clear linear relation with concentration, image noise causes significant random errors on the concentration determined in individual voxels. The minimum noise level is determined by the finite number of detected photons in each voxel (Poisson noise, Cnudde & Boone, 2013). We evaluated the uncertainty on the micro-CT-based concentration measurements with the average standard deviation of a voxel's concentration over eight repeat scans (Figure 4). The average uncertainty range (as indicated by the standard deviation, 1σ) of the concentration is ±27%. The reliability can be increased by averaging over multiple voxels. Assuming similar photon counts in neighboring pore voxels, the signal-to-noise ratio is expected to increase with √N, where N is the number of voxels that are averaged. Given the common approach of describing the pore space as a network of pores separated at local constrictions (regions of low hydraulic conductance) in both theoretical and modeling studies (Mehmani & Tchelepi, 2017), we propose to average voxels per pore body in order to quantify local concentrations during the experiment.
The average volume of the identified pore bodies is 4,214 voxels (1 × 10 7 μm 3 ). We calculated the standard deviation of the average relative concentration for each pore body over the eight repeat experiments. These uncertainty values (1σ) varied from ±0.1% to ±8.5% with an average of ±1.2%. All the pore-based standard deviation values (for all the repeat experiments) are plotted against their pore volume in a 2D histogram ( Figure 5). In general, smaller pore bodies have a larger standard deviation than larger pore bodies. This histogram can also be used to estimate pore-based uncertainties in other samples, provided they have a similar overall X-ray attenuation as the porous glass. For the Bentheimer sample, this indicates an average estimated uncertainty of ±1.6%. The uncertainty on concentrations averaged over pore bodies is thus significantly smaller than those measured on individual voxels (Figure 4). Note that the tracer injection experiments were performed with scanning parameters that are expected to yield a better signal-to-noise ratio than the calibration experiments (15 s versus 12 s acquisition time, respectively).  This uncertainty can be calculated for sub-pore scale down to single voxels; however, then the uncertainty will rise significantly.

Tracer Injection Experiments 3.2.1. Arrival Time and Filling Duration
The arrival time and filling duration in each pore are visualized in Figures 6 and 7. The arrival times show at which time the separate pore bodies start to fill with the tracer solution. The filling duration captures the width of the local concentration curve, which is related to hydrodynamic dispersion (Kanzari et al., 2015;Lapidus & Amundson, 1952).
As expected, the arrival time is linearly related to the distance from the inlet to the pore ( Figure 6). Slightly higher arrival times can be noticed in the middle of the sample compared to the cylindrical borders at the same distance from the inlet. This is likely due to boundary effects. The arrival time reflects pore-scale velocities in the sample. The average linear velocity calculated from the slope of the linear correlation to the data points deviates less than 16.5% from the injected interstitial velocity in all the experiments (Table 2). This indicates that the time resolution of the imaging is sufficient to characterize the transport behavior. The small differences between the different interstitial velocities are probably due to (1) small errors in the measurements of the experimental setup, (2) small errors due to the pump in the setup, and/or (3) errors in the linear fitting. As expected, the average linear velocity from the linear correlation in the sintered glass sample for Pe = 2.68 is approximately half of the linear velocity for Pe = 5.35. For the Bentheimer sandstone the linear velocity is higher than in the sintered glass sample because the porosity is smaller. Here the ratio between the linear velocities for the two Péclet numbers deviates more strongly from the expected value of 2. This  Table 2). may be related to the fact that the variations in pore-scale flow velocities are expected to be higher in more heterogeneous samples with smaller pore sizes.
The filling duration rises approximately linear in function of the distance from the inlet. This reflects longitudinal dispersion. In Figure 8, the filling duration is normalized against a characteristic time (characteristic length/interstitial velocity). This results in a dimensionless filling duration. The histograms show a broader range in filling durations for the Bentheimer sample than for the sintered glass sample. The filling duration in a pore generally depends on the size of the pore and the pore-scale velocity variations. When these variations are small, we expect a small range in dimensionless filling durations. For the Bentheimer sample, this range is broader, which suggests a broader range in pore-scale velocities and consequently more dispersion. The histograms in Figure 8 do not show significant tailing, which indicates the absence of stagnant flow zones. This is in line with the well connected, highly permeable nature of the samples (Bijeljic et al., 2013). Due to boundary effects, pores located in the center of the cylinder have generally smaller filling durations than pores located close to the border of the cylinder (Figure 7). Generally, for lower Péclet numbers (lower flow velocities), a slightly higher spread  in the data of the filling durations was observed (Figure 7). This could be due to the increased importance of diffusion.

Concentration Moment and Scalar Dissipation Rate
The concentration moment was calculated for the whole sample over time, based on the pore-scale concentration data. From this, the scalar dissipation rate, which represents a measure for destruction of the concentration variance, was calculated. In Figure 9, the concentration moment and scalar dissipation rate against time (normalized by the characteristic time) are shown. This shows that the normalized concentration moment and normalized scalar dissipation rate are similar for the different Pe numbers in the same sample. This may be due to the relative small difference in Péclet numbers tested, indicating similar dispersion regimes. In the Bentheimer sample, the scalar dissipation rate has a wider distribution than in the sintered glass, which can be linked to the wider range in filling duration (Figure 8). This means that relative to the characteristic time of the flow, it takes longer before the concentration variance declines by mixing and spreading mechanisms in the continuous injection experiments. This can be linked to the smaller pore sizes and more heterogeneous nature of the sample ( Figure S1 in the supporting information).

Conclusions and Further Research
In this methodological study, fast laboratory-based micro-CT is used to quantify pore-scale tracer concentration fields. Although fast laboratory-based micro-CT images are linked to low signal-to-noise ratios, calibration experiments proved that the concentration in individual pores can be determined much more reliably than in single voxels, at the cost of losing information on the sub-pore scale. This provided a basis for a workflow to quantify local transient concentrations in the tracer injection experiments. From these concentration data, we deduced when the solute arrives in individual pores (arrival time) and how fast the pore-based solute concentration rises (filling duration). An average interstitial velocity was calculated based on the arrival time data for every pore body. The determined average velocities correspond well to the initially calculated velocities, which gives confidence that errors coming from the finite time resolution of the images have a minor influence. The range in dimensionless filling duration is higher for the Bentheimer sandstone, which can be related to increased dispersion due to the more heterogeneous nature of the sample. The principles of the methodology can readily be applied or extended to other experimental setups and to samples with different degrees of heterogeneity, to provide more insight on the influence of pore-scale heterogeneity on solute transport processes. Nevertheless, the samples investigated here had relatively large pore sizes, and future work should thus indicate the reliability of the method in more complex samples, such as carbonate rocks. Furthermore, the maximum Péclet numbers that can reliably be imaged in the experiments are limited by the relation between the achievable micro-CT time resolution and the sample's characteristic length. In such cases, synchrotron imaging may provide more appropriate data. Future work should also indicate if the presented experiments can be adapted to separate spreading and mixing patterns more explicitly (e.g., pulse experiments). There is a significant difference between the sandstone and sintered glass sample, linked to the influence of pore-scale heterogeneity on solute transport.