Nonuniqueness of hydrodynamic dispersion revealed using fast 4D synchrotron x-ray imaging

Description

: Scheme of the solute transport experimental setup and X-ray imaging. Note that the dimensions of the experimental apparatus are not scaled with actual size. The injection direction to the sample is from bottom to the top. The pairs of loading and unloading experiments were performed at the same injection rates for 6 difference rates ranging between 0.2 to 6.4 µl/s. The holder was filled with fine sand and the fluids used were 3M potassium iodide (KI)-water solution and de-ionized water.

X-ray image analysis of the concentration field
The sXRCT images were processed with Fiji (38) and Avizo with the assistance of the OpenCV library (39). First, the collected 32-bit images were transformed to 16-bit images in Fiji. Then, the images of water-saturated samples were selected to make a mask of the pore space. The pore space was partitioned out with the use of the watershed algorithm in Avizo. Uaing the OpenCV library, the binary images were masked on grey scale images to remove the solid phase. The masked images were processed in Fiji with the "analyze particles" function to get the mean CT intensity of the pore space.
To establish a relation between CT values and actual concentrations, four calibration experiments were performed. Four solutions with known KI concentrations (0.1, 1, 2, 3 mol/l) were Figure S2: Histogram and cumulative probability distribution of pore radii. Pore size distribution was extracted using image processing. Left y-axis shows the probability of each bin size and the right axis shows the cumulative probability.
injected into the sandpacking and the CT values were obtained. Due to the nearly linear CT response with the actual KI concentration (as shown in Fig. S3), the solute concentrations were then computed based on Eq. 1 (40,41). We did not apply any filters to the grey scale images in order to keep the original CT intensities rather than lose information during the filtration process, where C voxel refers to the solute concentration in the voxel of interest, I voxel , I 3 mol/l and I water indicate the grey scale value of the same voxel, 3 mol/l KI and clean water, respectively.

Estimation of upscaled hydrodynamic dispersion coefficient and pore velocity
To obtain the pore velocity and hydrodynamic dispersion coefficients during loading and unloading, we fitted the average resident concentration curves to the analytical advection-dispersion Initial and boundary conditions: Analytical Solution: where where and in the experiments of this study the Ra number is expected to be smaller, due to the upwards direction of viscous forces, which hamper the downward forces due to the density difference.
Thus, the Rayleigh instability will be even further suppressed in these injections scenarios.
Following (43,44), Ra should be smaller than the critical Rayleigh number (Ra c ) to have the stability in transport: 400 Kg/m 3 ⇥ 9.8 m/s 2 ⇥ 4.8 ⇥ 10 12 m 2 ⇥ 3.47 ⇥ 10 3 m 2.5 ⇥ 10 9 m 2 /s ⇥ 1 ⇥ 10 3 Kg/(m · s) , where ⇢ is the density difference between two fluids, K is the intrinsic permeability of the porous medium, H is the height, g is the gravitational acceleration, is the porosity, D is the diffusion coefficient, and µ is water viscosity. Ra =26.1 is smaller than the critical Rayleigh number (43,44). Thus, we can conclude that that role of density-dependent instability in our observations was negligible.

Effect of the Centrifuge force
To calculate the centrifuge force and its importance in the data analysis, the centrifuge force was calculated at the outmost part of the flow cell.
Since the sample rotated ⇡ radians in 3 seconds, the angular velocity (!) was 1.05 radians per second. Given that the outmost radius of the sample was r = 2.25mm, the centrifugal force acting on a mass in a pore with the diameter of 50µm would be : F r ⇡ ! 2 r⇢D 3 = 4.13⇥10 13 N . For the same control volume using the Hagen-Poisseuille flow, the viscous force for the smallest flow rate would be approximated as F v = 8Qµl/R 2 = 1.7⇥10 11 N . As shown, even for the smallest flow rate the centrifugal force will be negligible. Additionally, the ratio of centrifuge to the gravity acceleration isF = ! 2 r/g = 1.05 2 ⇥ 2.25 ⇥ 10 3 /9.8 = 2.5 ⇥ 10 4 .
This also indicates negligible centrifuge acceleration compared to the gravity.

Diffusion coefficients from molecular dynamics simulation
To assess the change of diffusion coefficients with concentration, we calculated the diffusion coefficients by molecular dynamics simulation in different KI concentration solutions at 0.01 mol/L, 0.1 mol/L, and 1 mol/L. The molecular dynamic simulation was taken in LAMMPS (45). The simulation input were prepared with assistance of Moltemplate (46). The SPC/E water (47) was used and a force field for KI solvent was from Koneshan et al. (48). The nonbonded potential between different molecules was modelled with arithmetic mixing method.
The simulation boxes are shown in Fig. S5. In box A and C, it was initially filled with water molecules while box B was filled with KI solution. The simulation includes two steps. At step one, the three boxes were relaxed in NVE ensemble with a Langevin thermostat at 300K under periodic boundary conditions in x and y directions while a fixed boundary condition was in z condition. A wall was applied in z direction to avoid atoms lost. The three simulations were run for 25 ns to fully relax the system. After the system was relaxed, the three boxes were combined together and ran in a NVT ensemble under periodic boundary conditions in three directions. The dynamic diffusion coefficients were calculated with Einstein equation (Equ. 13) as shown in Fig. S5.
Where r(t) is the position of molecule at time t, D is diffusion coefficient, t is time, hi means the ensemble average.

Impact of concentration-dependent diffusion coefficient on transport time scale
To demonstrate qualitatively the impact of concentration-dependent diffusion coefficient on transport time scale, we performed 1D numerical simulations using the conventional advectiondiffusion equation. @C/@t = v@C/@x + D@ 2 C/@x 2 . For a 1D system with the length of 5⇥10 -3 m, velocity of 10 -5 m/s and D m = 2.5 ⇥ 10 9 m 2 /s, simulations for the loading and unloading scenarios were done (Fig. S6). To demonstrate the difference between loading and unloading due to concentration-dependent diffusion, we defined the deviation factor as C load 1 + C unload shown in Fig. S6b, which show positive values. The average resident concentration profiles for loading and unloading scenarios are shown in Fig.   S6c. Clearly, our results show that even for a 1D case (fixed velocity) the effect of concentrationdependent diffusion coefficient is visible. Note that this mechanisms will be more pronounced in high-concentration gradients such as saline aquifers (27). Additionally, for a 3D system with a large variation of pore velocity, it is expected that delay in unloading scenarios to be more pronounced than the small systems.