Iterative Algorithm for 4D Tomography Reconstruction Using a Single Projection per Time Step

In this paper, a new method for the mathematical tomographic reconstruction of time-varying objects (4D tomography) is presented. The method is applicable under the conditions that the amount of substance at each spatial point does not decrease, and only one projection image is available per time point during the process of the object changing. The temporal resolution of the proposed method is obtained by using an iterative approach to reconstruction. Prior knowledge about object morphology is exploited to improve the convergence of the algorithm by reducing the number of processed spatial points. Verification of the algorithm operation was carried out using model examples. The presented method was experimentally tested on the tomographic reconstruction of the dynamics of the capillary rise of a liquid. Our experiments demonstrated that monitoring fast processes, when the speed of the process does not allow to measure more than one tomographic projection, is potentially possible.


I. INTRODUCTION
Today, the computed tomography (CT) method is an all-powerful, non-destructive tool for studying the internal structure of objects. CT uses a set of two-dimensional images that are registered while the studied object is illuminated by X-rays at several projection angles [1], [2]. The object may either be stationary or non-stationary. The setups for X-ray tomographic imaging are constantly being improved [3]- [5] due to the new statements of the tasks to be solved. Traditionally, in tomographic studies, it is assumed that the object remains unchanged, and the data obtained during the measurements are used to reconstruct the three-dimensional structure of that object.
The associate editor coordinating the review of this manuscript and approving it for publication was Yongqiang Zhao .
The assumption about the stationarity of the object ceases to be valid when studying dynamic processes. Now, the object is described not by one three-dimensional digital image, but by a series of such images. Each of them corresponds to the state of the object at some point in time. Thus, the method of 3D tomography turns into a 4D tomography method (3D + time). One of the first areas in which the task of monitoring non-stationary objects, including tracking of anatomical structures, arose was medicine [6]. Observations were carried out over the dynamics of the development of organ anomalies. The observation time for the object, which can be continued for several hours, days, or weeks, was divided into bins. The time of observation of the object (total duration of time bins) and the time of the tomography scan (timeframe) were two different times. The study was a retrospective study. If the object is non-stationary, the tomography problem is formulated as follows: To reconstruct a digital image of the object at different timeframes. Within the timeframe, the object is considered to be in a quasi-stationary state. To reconstruct a quasi-stationary object from a set of tomographic projections, algorithms developed for traditional tomography are used.
In the early 2000s, real-time high-speed conical computed tomography [7] was implemented in hardware. The time of observation of the object and the tomography scan time become one period of time. Experiments with a bouncing golf ball and magic hand phantom have demonstrated the potential for capturing projections from a non-stationary object in realtime. Previously, a quasi-static 3D image was reconstructed by the Feldkamp method [8]. The demonstrated possibility of a hardware solution for collecting projections of a highspeed tomographic method has opened up new prospects for its application in studies of dynamic processes. The concept of a model of the process itself has been introduced. Thus, the model of human organ movement, based on the property of preserving the local tissue volume, was taken into account when reconstructing quasi-stationary states of the organ in [9].
Another type of dynamic process is periodically repeating processes [10]- [12]. For such processes, projections are measured over several periods. The sequentially captured projections are rearranged by phase of the process to belong to the same quasi-stationary state, after which one of the reconstruction techniques developed for the stationary case is applied [13].
The study of dynamic processes at synchrotron stations can also be considered a separate class of problems. With a synchrotron source, it is possible to greatly reduce the projection acquisition time [14]- [16]. Tomographic experiments with insects [17]- [19] have demonstrated the ability to monitor the movement of their limbs without involving special reconstruction algorithms. Laboratory experiments to study the work of muscle groups responsible for various functional actions have also been carried out with small animals [20].
At the same time, the improvement of the algorithmic part of tomographic methods, applied to the study of nonstationary objects, was begun. An algorithm for calculating the field of direct and reverse deformation to study the periodic processes occurring in an object has been proposed [21]. It allows the possibility to clarify the result of the reconstruction of quasi-static objects in the series. A spatial-temporal regularization approach to studying non-periodic changes in the state of an object has been presented [22]. The projections measured during all time frames are used together (not one by one) in reconstruction algorithms. The result of the reconstruction is a function of space and time. The usage of a non-uniform projection acquisition scheme [23] in measurements has made it possible to increase time-space resolution. The terms of the space and time penalty are combined into one regularizing term. The further addition of a piecewise constant function for the separate analysis of stationary and non-stationary areas of an object in the reconstruction algorithm has made it possible to increase the accuracy of the reconstruction, which has been demonstrated using the data of neutron tomography [24].
To complete this review, we draw a few conclusions. There are two approaches to the implementation of 4D tomography. The hardware approach consists of ultra-fast measurement of several complete sets of tomographic projections, during which the object is considered to be quasi-stationary. Not one, but several rotations of the gantry around the object or rotations of the object are performed if the emitter-detector system is stationary. 3D reconstruction algorithms are applied independently for each set of projections. The measurements require high flux probing and very fast and expensive detectors. In the second approach, most of the work is assigned to the reconstruction algorithms used. We propose a fast, novel algorithm to solve the optimization task. As in previous works [24], we used in the algorithm the information about stationary areas but instead of the regularization term usage, we update the parts of the sinogram calculated with the current solution in agreement with the measured values. The algorithm works fast enough since a new solution describing the next quasi-stationary state of the object is built using only one projection. The algorithm is fast but imposes more stringent restrictions on the relationship between the rate of the process under study and the measurement time of one projection. The time resolution is determined by the exposure time. We assumed that changes in the object under study occur in such a way that the amount of matter (and the absorption of X-ray radiation as well) at each point of the object does not decrease with time. These processes include fluid flow in some porous media [25], [26], crystal growth, and 3D printing. Measurements may be carried out on a laboratory X-ray source with a parallel tomographic scheme and a low-dose load. The algorithm uses a digital image of the object obtained before the start of the dynamic process. The whole method is performed in several stages: tomography of the initial state of the object, reconstruction of its digital image, binarization of the reconstructed image, tomography of a dynamic process, and slice-by-slice reconstruction of a dynamic process using an initial state of the object as a binary mask.
The remainder of this paper is organized as follows: Section 3 presents a detailed description of the slice-byslice reconstruction algorithm. Section 4 describes the verification of the algorithm on model 2D data (IVA), and experimental testing of the algorithm with the tomography of the capillary rise of a liquid (IVB). It is followed by the discussion (Section V). Section VI summarizes the results. A description of the liquid rise model is presented in Appendix A. Appendix B presents the results of the liquid rise simulation. In conclusion, the results obtained and the prospects for the development of the proposed approach are discussed.

II. STATEMENT OF THE PROBLEM
In this paper, we consider the filling of a porous object with a liquid as an example of a dynamic process. However, all the reasoning is valid for other processes as well, provided that the conditions listed below are met. Our proposed algorithm for the reconstruction of dynamic processes images is based on the following assumptions: 1. We assumed that in the process of changing the object under study, the amount of matter at each spatial point does not decrease with time. For example, let us consider the filling of a sample with a liquid. This means that in each elementary volume (voxel) of the object, the amount of liquid cannot decrease. The voxel value can either remain unchanged or increase in the corresponding digital image (image of the distribution of the absorption coefficient, which is in our case an equivalent to the amount of liquid in voxel). 2. We can obtain only one projection image of the object per time point while registering the dynamic process (filling the object with a liquid) 3. The unchangeable structure of the object is known in advance, which means we have a complete tomography of the initial (empty) state in the case of filling the object with a liquid. We supposed the object substance will be unchangeable, while the pores will be a changeable area. A changeable area is an area in which the dynamic process will take place. This condition is not necessary, but data on the unchangeable structure of the object improves the convergence of the algorithm by reducing the number of processed spatial points. 4. In the experiment and simulation, the plane wavefront approximation is used, i.e., the source is far enough away to consider the incident radiation beams parallel to one another. In this case, we can consider a 3D object as a set of horizontal 2D slices, and in the process of the algorithm operation, each such spatial slice is processed independently of the others. Figure 1 shows a scheme for obtaining tomographic projections in a parallel beam.
Based on these data, it is necessary to restore the dynamics of the process (filling the object with liquid) -the image of a 3D object for each time point.

III. DESCRIPTION OF THE ALGORITHM
Before the algorithm begins, the 3D object under study is separated into an array of horizontal 2D slices. For each such slice, the algorithm reconstructs the dynamics of changes in the slice in time. After the algorithm is finished, all the reconstructed slices are collected into 3D objects for each time point. A set of such 3D objects in time provides the desired 4D image of the dynamic process.
To describe the steps of the algorithm, we introduce the following notations (see Figure 2): N is the number of registered projection images of the object. By projection, we mean the linearized [27], [28] value of the detector reading. In the proposed approach, N also equals both the number of rotation angles of the object and the number of time points in the experiment. I n is the projection of the object in a dynamic experiment, obtained at an angle of rotation a n , at time t n , where n = 1 . . . N .

FIGURE 2.
Scheme for obtaining and designating experimental data: t 1 < t j < t i < t N are time points; a 1 , a j , a i , a N and I 1 ,I j , I i , I N are corresponding rotation angles and projections, respectively; S e is the experimental sinogram. The top row shows the states of one slice of the object during the experiment at the mentioned time points. The experimental sinogram S e is shown on the right. Numbers 1-4 on the first (left) image of the object slice and on the experimental sinogram denote pores and their projection traces, respectively. Below, under the states of the object slice, the corresponding projection images of the slice state are shown. There is a gradual, uneven filling of the pores with liquid from time t 1 to time t N . The lighter the pore, the more liquid it contains (the absorption coefficient is higher).
The set of projections for each spatial slice of the object forms the experimental sinogram S e . The slice has a thickness of one voxel; voxels in the slice are parametrized by a set of 2D vectorsr. This set is divided into two parts: unchangeable, corresponding to the solid structure, and changeable, corresponding to the pores to be filled. The set of reconstructionsR consists of N elements R n that are the reconstructions of the state of the object at all time points t n . R 0 is the reconstruction of the initial state of the object. We denote the reconstructed values of the absorption coefficient (amount of liquid in the voxel) in the given voxelr at time t n by R n (r) and by R 0 (r) for the initial state.
Based on our assumption about the process (the amount of material on any point does not decrease), the condition R n r ≥ R n−1r must be met for any n and anyr. We define the VOLUME 10, 2022 mean L2-norm to the single reconstruction ||R n || 2 and to the set of reconstructions ||R|| 2 : where ch.a. stands for the changeable area.
The sequence of steps of the algorithm ( Figure 3, Algorithm 1) is as follows: Step 0. Data initialization.
Create N identical sinograms S n ← S e (n = 1 . . . N ) corresponding to the initial state of the object (with no liquid). During further steps of the algorithm, the sinograms S n , (n = 1 . . . N ) will change.
Step 1. Use of experimental data.
From the experiment, we know that n-th row of S e corresponds to the projection angle a n , and n-th row in each sinogram S n should be equal to the experimental data corresponding to the projection angle a n . For every n = 1 . . . N , we replace the n-th row in the sinogram S n by the n-th row in the experimental sinogram S e . Other N − 1 remaining rows stay unchanged.
Calculate the set of reconstructionsR by use of the algebraic SIRT method [29] for every updated sinogram S n .
Step 2b. Use of initial conditions. For every r of the unchangeable area, replace the value of R n (r) with R 0 (r), if the reconstruction of the initial state of the object is known.
Step 3. Use of prior information about the process. For every n and every pointr of the changeable area of the object, check the condition of the non-decreasing amount of material. If it is not satisfied, correct the corresponding values for each R n , and save them into theR iter , which consists of R iter n and will be used at the next step of the algorithm: R iter n (r) = min(R n (r), R n+1 (r)) Step 4. Loopback or exit. Calculate the sinograms S n from eachR iter and go back to Step 1. If the result of Steps 1 -3 leads to an insignificant change of the set of reconstructions (the L2 norm change by a value less than 10 −5 ) we exit the algorithm. A minor change of the set of reconstructions means that the replacements of rows in sinograms (Step 1) and the corrections due to the information about the process (Step 2) do not noticeably affect the result. Thus, the calculated sinograms S n are close to the experimental sinogram S e and describe the observed process.
Remark 1: Step 2a can also employ other methods of tomographic reconstruction, e.g., FBP [29]. The choice of the specific method depends on the parameters of the experiment

Algorithm 1 Tomography Reconstruction of Time-Dependent Data
Input: experimental 2D sinogram S e , each n-th (n = 1 . . . N ) row S n e corresponds to the experimental data at the angle a n (at the time moment t n ), R 0 (r) is 2D reconstruction of the empty porous structure,r unchangeble -a region on reconstruction where the reconstruction values shouldn't change due to prior information about the sample. ε is the precision. Prior information: Concentration of liquid can't decrease R n (r) ≤ R n+1 (r) , for n = 1 . . . N − 1 Output: R n (r)-2D reconstructions of the object at each moment of the time t n for n = 1 . . . N Initialization: Create N sinograms for the iterative procedure S n ← S e for n = 1 . . . N Create N reconstructions to calculate the first iteration of the reconstruction set R n = SIRT (S n ) Step 1: Use of experimental data do Update n-th row in each S n sinograms to be consistent with experimental data for n in 1 . . . N : Step 2a: Reconstruction.
Step 2b: Use of initial conditions for every r of the unchangeable arear unchangeble .
for n in 1 . . . N : R n (r) ← R 0 (r) , ∀r ∈r unchangeble Step 3: Use prior information about the process.
for n in 1 . . . N − 1 : (signal/noise ratio, the total number of projections, etc.). The algebraic method works better at a low signal/noise ratio while FBR is fast for good enough data.
Remark 2: If the initial state of the object is not known the Step 2b is skipped.
Remark 3: Our calculations showed that the algorithm either converges within a reasonable number of iterations or does not converge at all. Thus, the exit criterion can be replaced by reaching a large enough number of iterations (e.g. 10 000) and control of the L2 norm change at the last iteration. In Section IVA the equivalence of these conditions is shown.

IV. EXPERIMENTAL VERIFICATION
The most evident field of application of the proposed algorithm is the study of the dynamics of fluid movement in porous media. To test the algorithm, we perform two experiments on the filling of a porous structure when the ground truth is known. First, we used a model object with artificially generated data. In this case, we set the ground truth ourselves. The second experiment deals with a physical object. We observe a rise of oil in a vertical capillary. The simple geometry allows comparing the tomographic results with a theoretical model from hydrodynamics as well as visual observations. The following subsections give details of these experiments.

A. VIRTUAL EXPERIMENT 1) SETUP
For the virtual experiment, we generate the object: a twodimensional model of a porous structure (one horizontal slice of a three-dimensional object; see Figure 4). We will refer to this object as the model object. It consists of 1024 (32 × 32 × 1) voxels. Some of them form the unchangeable area representing the matrix, the others are pores. The pores are located in the center of the object to avoid tomography reconstruction artifacts and form the changeable area. They are empty at t = 0. The matrix and pores are shown as white and black in Figure 4, respectively. Further, we simulate the filling of the space by liquid. We generate the state of the model object in 100 time points with an equal interval between them: at any time t n for each voxel r in the changeable area, we set a number c n (r) between 0 and 1 representing the absorption coefficient (which corresponds to the amount of liquid). For fixed r, the sequence c n (r) monotonically grows with n. Figure 5 shows the model object for several time points; the color of pores transforms from black to white while the liquid fills the pore, and the absorption coefficient grows. We calculate the projections of the model object and use them as input information for our algorithm.
We generated sets of N projections I n for N = 100, 50, 25 considering every first, second, and fourth-time point. The projection angles are uniformly distributed across the circle for each case. A comparison of the results for different values of N gives estimations for the effect of the number of projections.

2) DATA ANALYSIS
First, we analyzed data for N = 100 and after that apply the same procedure for lower time resolution. In this subsection, the results correspond to the large value of N , unless other stated explicitly.
The algebraic SIRT method from the ASTRA toolbox software package was used as a reconstruction method from Step 2 of the described algorithm. We perform 10,000 iterations of our algorithm. Figure 6 shows intermediate reconstructed images for the middle of the whole time interval (50 th point in time). By the 1000 th iteration of the algorithm, the resulting image is in good agreement with the ground truth (RMS deviations of the reconstructed values from those specified in the model for iterations 1, 200, and 1000 were 0.0322, 0.0095, and 0.0049, respectively). Figure 7 shows the dependence of the reconstruction norm on the iteration number (Figure 7a) and the change in the differences between the norms of successive iterations (Figure 7b). In this case, the total number of iterations of 1000 seems to be a reasonable compromise between the reconstruction accuracy and the algorithm running time.

B. EXPERIMENT WITH A REAL OBJECT
We tested the proposed algorithm for the tomographic reconstruction of the rise of oil in a capillary due to surface tension.
In this case, we were able to check the obtained data of 4D tomography visually using shadow X-ray projections, as well as compare them to simple hydrodynamic calculations.

1) SETUP AND DATA ACQUISITION
The experiment was carried out on a laboratory tomograph ''TOMAS'' at the FSRC ''Crystallography and Photonics'' RAS [30], [31]. A conventional X-ray tube with a molybdenum anode was used as a source. We used a pyrolytic graphite monochromator crystal to isolate the characteristic line. The energy of the probe radiation was 17.5 keV. The distance from the source to the sample was 1250 mm, and from the sample to the detector was 15 mm. The ratio of these distances allowed using the plane wavefront approximation and working in the approximation of a parallel tomographic scheme (see Item 4 from Section II and Figure 1). Radiation was recorded with a XIMEA xiRAY11 detector with a pixel size of 9 × 9 µm and a field of view of 36 × 24 mm. The characteristic time for a complete tomographic experiment is 0.5-2 hours, depending on the absorption level of the object under study. The characteristic time of the experiment and the size of the detector's field of view define the capillary size ranges and the working liquid's viscosity. We took a glass capillary with a height of 33.4 mm with an internal channel of 200 µm in diameter and the silicone oil (polymethilsiloxane) PMS-1000 with a viscosity of 1000 cSt, surface tension of 0.019 N/m, and a contact angle with a glass of 0 degrees. Appendix A provides details on the calculation of these parameters.
The capillary was installed to a special holder with a container for working liquid (see Figure 9). The scanning protocol of the experiment consisted of two parts. First, we performed a tomography study of an empty capillary (without liquid). After that, oil was poured into the capillary holder, and in the second part of the experiment, we repeated the tomography study during the rise of the liquid in the capillary. In each tomographic scan, we recorded 444 projections with a sample rotation step of 0.45 degrees. Each projection exposure time was 1.5 seconds, and the pause between projections due to the sample rotation and detector initialization was 3.1 seconds. The total duration of one scan was 34 minutes. Thus, for each angle of rotation of the sample, we obtained projections of both an empty ( Figure 10a) and a partially filled capillary (Figure 10b). For each pair of images, we calculated the logarithm of the ratios of the normalized signal at each point (the difference in the images), as displayed in Figure 10c. The obtained differences were used as the initial data for the algorithm's operation (see Section III). The calculations were performed with the Python-based implementation of the algorithm on Intel Core i7-5930K, 64 Gb RAM, Nvidia GTX 980Ti. The algebraic SIRT method was used as a reconstruction method. We made 100 iterations of the proposed algorithm for each tomography slice since further iterations do not significantly improve an image's quality due to the noise in experimental data. The calculation of each slice took 20 minutes.  After processing the data with the proposed algorithm, we obtained reconstructed images of the horizontal slices of the object at all time points (Figure 11). The subsequent connection of the slices into a single object allowed obtaining a 4D image of the capillary rise.

2) COMPARISON WITH THEORY AND MODELING
To assess the accuracy of the dynamic process reconstruction by the presented algorithm, it is necessary to compare the VOLUME 10, 2022 reconstructed 4D image with the experimental data. In the case of a simple object (e.g., a capillary), the independent results can be obtained from the shadow projection or the calculations.
Differences between the shadow projections of the empty and partially filled capillary provide distinct liquid boundaries and give the accuracy of the liquid column height determination of the order of a few pixels (∼10 µm). Theoretical dependence of the liquid column height h over the liquid level far from the capillary on time t follow from the balance between the lifting force of surface tension, pulling down gravity, and decelerating viscous friction (see Appendix A). It is expressed by the formula [32] where H = 2σ cos θ ρgR , T = 16µσ cos θ ρ 2 g 2 R 3 are reference height and time. They are determined by the surface tension coefficient and the contact angle σ, θ, the radius of the capillary R, the density and viscosity of the liquid ρ and µ, and the immersion depth of the lower edge of the capillary h 0 ; t 0 is a free parameter that allows considering the initial interval when inertia is relevant (see Appendix B for details).
In the tomographic reconstruction, we consider horizontal layers may be reconstructed independently of each other. We consider a layer at the height h filled the mean adsorption coefficient of the voxel inside the capillary reaches a specified value. The time when the filling criterion is satisfied is less sensitive to the percentage of the lower layer since the liquid fills the lower layers quite quickly. For the upper layers, we have a small number of the projections with the filled layer, and signal to noise ratio is less than for the lower layers.  Figure 12 shows a comparison of the experimental data of visual observations, a calculation using the non-inertial capillary rise model (see Appendix B for details), and the reconstruction of process dynamics by the proposed algorithm. For the latter, three filling criteria are used.
If the object to be filled has a more complex shape (e.g., a porous ceramic filter or rock core), it can be difficult to visually determine the quantitative characteristics of filling the porous medium with a liquid from the projection data. Approximate analytical methods for calculating the flow are also not applicable. For comparison, one can use, for example, the dependence of the mass of a liquid inside of the sample on time and estimate the position of the interface calculated from the tomographic reconstruction. The former dependence can be measured in the experiment independently, the position of the interface can be found numerically, e.g., by use of the Lattice-Boltzmann equation (LBE) [33], which results in the same data as the algorithm -the distribution of the liquid concentration over voxels depending on time. An example of the application of the LBE to the problem of filling a capillary is given in Appendix B.

V. DISCUSSION
As can be seen from the results obtained, the proposed algorithm for tomographic reconstruction makes it possible to reconstruct the dynamics of non-stationary monotonic processes. In this case, the temporal resolution is potentially limited only by the exposure time of one projection. It leads to a decrease in the total experiment time, and, consequently, to a lower accumulated radiation dose for the object under study. It also makes it possible to apply this method for studying dynamic processes on laboratory X-ray radiation sources, where experiments on 4D tomography with the reconstruction of quasi-static states by a large set of projections are severely limited in their capabilities due to the relatively long exposure time. At the same time, the application of the proposed algorithm to synchrotron sources will make it possible to study extremely fast and non-periodic processes.
In our virtual experiment, we observe the monotonic dependence of the reconstruction quality on the number of projection images. In real experiments, the accuracy is limited by noise and at some point, the increase in the number of projection images does not lead to better reconstruction. The total duration of the process is a given value. The scanning protocol provides a tradeoff between a number of projection images that defines the accuracy of reconstruction and exposure time per projection that corresponds to signal to noise ratio.

VI. CONCLUSION
Herein, it has been demonstrated that the presented algorithm for the mathematical reconstruction of 4D tomography makes it possible to reconstruct dynamic processes with only one projection image available per time point while the process of object change is non-decreasing in terms of the intensity of the voxels of reconstructed images. We proved the feasibility of the proposed algorithm using the data from two experiments: a virtual one with a complex structure of the porous object when the input data were generated artificially and the ground truth is known a priori, and a real one with the data obtained by the X-ray tomography of a simple object. The design of the latter experiment allowed independent monitoring of the filling process by shadow pictures and simple modeling. Experimental verification by studying the process of filling a capillary with a liquid demonstrated the possibility of the practical application of the algorithm. An example of the application of the proposed method is the study of fluid flows in porous objects. Tomography of stationary porous objects today are performed on different scales [34], [35]. This is often connected with the tasks of mineralogy [36], but not always alone. The study of hierarchically structured porous morphology is requested in the manufacturing of batteries [37], for example. Since it is no longer enough to visually analyze the reconstructed images of porous objects, methods for the automatic binarization of tomographic images of porous structures, used as a mask in our algorithm, have also been developed [38]. We understand that the trade-off between the rotation speed and the signal-to-noise ratio in the monitoring of dynamic processes is very important [39]. Our experiments demonstrated that monitoring fast processes, when the speed of the process does not allow to measure more than one tomographic projection, is potentially possible. The second advantage of such an approach is its possible application in medical studies due to extremely low doses during CT. An extra small exposition time produces an extremely bad signal-to-noise ratio; however, the development of new approaches for working with tomographic data collected under ultra-low-dose conditions [40], [41] is in progress.
At the moment, we are working to improve the proposed 4D CT algorithm. We plan to adapt it to such dynamic processes where changes in the objects under study can be non-monotonic.

APPENDIX A MODELING OF CAPILLARY RISE
The reference values for the capillary rise are the height in equilibrium H = 2σ cos θ ρgR and the time is given by Washburn's law [41] T = 16µσ cos θ ρ 2 g 2 R 3 . We used silicone oil (polymethylsiloxane) as the working fluid since it wets the glass (the contact angle is close to zero) and allows selecting a viscosity from a wide range. The surface tension coefficient of silicone oil weakly depends on its type (viscosity) and equals 0.019 N/m. The size of the detector's field of view determines the value of H as approximately 4 cm and fixes a capillary radius of approximately 0.1 mm. For a capillary radius of 0.1 mm and oils with a kinematic viscosity of ν = µ ρ = 10 −3 m 2 /s, the reference time is T = 3260 s, which is comparable to the typical tomography scan time.
For non-inertial capillary rise when the surface tension force, gravity, and viscous friction are in balance, the following differential equation describes the dependence of the liquid column height above the oil level h(t) on time [41]: R 2ḣ + ρgh where σ, θ are the surface tension and the contact angle, R is the radius of the capillary, ρ and µ are the density and viscosity of the liquid, the h 0 is the immersion depth of the lower edge of the capillary (see Figure 13) The numerical coefficient in the viscous friction term reflects the assumption that a parabolic velocity profile takes place in all sections.
A solution to this equation is: The non-inertial approximation is valid if the dimensionless number = µ 128σ cos θ ρ 3 g 2 R 5 > 2 proportional to the viscosity is large enough (exceeds 2 [42], [43]).
The parameters of the present experiment correspond to = 5.1 · 10 4 .

APPENDIX B SIMULATION OF THE FLOW IN AN ARBITRARY CHANNEL
An analytical solution exists only for channels of a simple shape -for example, a gap between planes, a cylindrical channel with a circular cross-section. The formulation of the problem requires explicit expressions for the meniscus curvature and the ratio between the mean velocity and viscous friction. In other cases, particularly for the flow in a porous structure, it is necessary to numerically simulate the flow within the framework of the Navier-Stokes equation, the Stokes approximation, or other approaches for viscous fluid dynamics modeling. The use of the Lattice-Boltzmann equation (LBE) [44] seems an effective tool for the numerical simulation of the flow in a complex domain. Solving Navier-Stokes equations for two-phase flow by volume-offluid method with OpenFOAM software [45] turned out to be enormously time-consuming [46]. The main idea of the LBE method is the following. For a uniform grid in the physical space and in the space of velocities, one introduces quasiparticles and solves the problem with respect to the distribution function for the density of these quasiparticles. At each time step, a particle can move VOLUME 10, 2022 to a finite number of states in the neighborhood of the initial one, and the probability of transition to one or another state depends on the external forces and the forces of interparticle interaction. The parameters of the equations define the magnitude of the surface tension, contact angle, and viscosity. The LBE method allows one to calculate the flow of multiphase media in a region with arbitrary geometry [31]. We used the Palabos code [47] for calculations with the LBE method. The applied procedure neglects the effect of air and does not perform simulation in the points occupied by the air or the capillary.
We simulated a flow in a vertical channel of D = 2 mm ( Figure 13). There was a capillary in the center of the channel, and the capillary's inner diameter was 2R = 0.2 mm, while the wall thickness was d = 0.3 mm. The configuration used in the experiment and shown in Figure 12 (left) was set with the following initial conditions: The liquid formed a layer 19 mm deep, and the empty capillary was immersed to a depth of 13.3 mm. We performed the calculation of unsteady threedimensional flow, despite the fact that the initial and boundary conditions had axial symmetry. The grid step in space was 1/50 mm and the capillary diameter equaled 10 grid steps; this spatial resolution is sufficient to reproduce the maximum capillary rise height value with an accuracy of 2.7% [44]. The time step was equal to 10 −5 s, which is close to the maximum value that ensures the stability of the computational scheme. The shape of the liquid surface is shown qualitatively in Figure 13 (right), while Figure 14 shows the results of calculations using the LBE method. The calculations showed that the time of rising to a certain height linearly depends on the kinematic viscosity for ν from 10 −6 to 10 −3 m 2 /s, which is consistent with the non-inertial approximation. To reduce the computation time, we set and then rescaled the obtained time multiplied by 10 −3 .
The calculation of 1 s of physical time took 54 s on CPU Intel (R) Core (TM) i7-5930K, 64 Gb RAM A comparison of the experimental results with the calculation and analytical solution (A1) is shown in Figure 15, where the value h = 0 corresponds to the initial oil level outside of the capillary. Calculation using the LBE method improved the accuracy in the initial stage when the liquid entered the capillary, although the assumption of a parabolic velocity profile was not fulfilled. The same conclusion was obtained in [43]. The discrepancy between the results of the calculations by the LBE method at large values of time can be explained by the influence of the outer boundaries of the computational domain. In other words, the capillary pressure in the gap between the channel and the outer boundaries reduces the pressure drop that drives the fluid.
The time for the capillary to rise to a certain height based on the non-inertial modeling and the experiments differed by a value close to constant (42 s for the parameters used). This difference corresponds to the time required for the initial acceleration of the liquid. The thin lines in Figures 11 and 14 show the curve obtained with a shift of 42 s as a result of modeling