## 2 Materials and Methods
### Lattice-Boltzmann Simulations  
The LBM simulations of the permeability of samples are performed in Palabos [44], an open source computational fluid dynamics software, which is a well-established software package tested in multiple LBM applications [45; 46; 47]. Its modern architecture and storage organization [44] allow efficient parallelization in massive computations. A large library of models and boundary conditions allows for the facile adaptation to a particular phenomenon. The methodology is straightforward: The flow of fluid through the reconstructed digital pore matrix is imposed by a fixed pressure gradient between the inlet and outlet faces of the sample. All simulations were held in LB units, represented as a two-step conversion from physical units to dimensionless, which do not depend on real physical sizes of the system, then from dimensionless to discretized, in order to convert to space and time grid-related units [48]. The gradient is implemented using the common relationship between density and pressure in LBM: \(p_{\mathrm{d}}=\rho c_{\mathrm{s}}^{2}\)[49], where \(p_{\mathrm{d}}\)--is dimensionless pressure, \(c_{\mathrm{s}}\) is the speed of sound in LB units, and \(\rho\)--density in LB units.  
Since the permeability is independent of the pressure, the density at the inlet boundary is constant and equal to \(\rho_{inlet}=1.0\), and the density at the outlet is varied in the interval \(\rho_{outlet}=[0.995;0.95]\) depending on the sample. The resulting gradient values were primarily aimed at the stability and accuracy of the LBM [49]. The methodology implies that the system obeys the Darcy law, which requires \(\mathrm{Re}\leq 1\).  
To avoid the artificial anisotropic velocity slip at the solid boundaries caused by the choice of relaxation time(s), we apply the MRT collision operator [23], which provides a more precise ratio between the kinematic and bulk viscosities of the flow. At the solid wall, the half-bounce-back boundary condition (HBB-BC) [50] is applied. HBB-BC effectively places an impenetrable boundary exactly between the cells. This approach helps improve the momentum exchange algorithm and maintain the Galilean invariance. It also alleviates the effect of the \"staircase\" approximation of complex curved boundaries [49], which is one of the major limitations of LBM in complex geometries. Because all of the samples  
Figure 1: Three-dimensional maps of tight-pore sandstones of the Achimov suite: (**A**) 3_E sample with overall porosity of 7%; (**B**) 3_B sample with overall porosity of 5.3%. Green represents the pore walls; the inner pore space is shown in black. The cube size is 720 μm.  
from the dataset were already voxelized and there was no information about isosurfaces of the wall interfaces, the off-lattice BCs that could potentially improve the accuracy are not applicable. The faces, which are not the inlet or outlet, are padded by symmetry boundary condition.  
Simulations are assumed to converge to a steady state when the mean kinetic energy deviation over 1000 iterative steps decreases below \(\epsilon=10^{-4}\). The relaxation time is \(\tau=1\), which provides optimal stability and accuracy of the simulation according to [PERSON] et al. [51]. In LB simulations of sample 3_C, we experienced problems with convergence because of the substantial number of narrow pores less than 4 voxels in diameter. That is, the CT resolution is insufficient for a permeability evaluation of such tight pores. To avoid the issue, we artificially increased the map resolution by scaling the matrix twofold and interpolating the boundary regions using the nearest-neighbor method (0-order interpolation). Such a technique is quite common [7; 52].