## 2 Materials and Methods
The general LBE is  
\[f_{i}(x+c_{i}\Delta t,t+\Delta t)=f_{i}(x,t)+\mathbf{\Omega}_{i}(f) \tag{6}\]  
In the Bhatnagar-Gross-Krook (BGK) model, the evaluation of the particle distribution function at each time step is separated into two primary processes, collision and streaming. The collision mechanism is mathematically represented as follows  
\[\Omega_{i}(f)=-\frac{f_{i}-f_{i}^{eq}}{\tau}\Delta t \tag{7}\]  
where \(f_{i}^{eq}\) is the equilibrium distribution function and \(s\) is the relaxation time. Macroscopic quantities such as density and momentum were computed from the distribution functions as follows  
Figure 1:0 \begin{table}
\begin{tabular}{c c} \hline \hline
**Velocity Direction** & **Weight** \\ \hline i = 0 & \(\frac{1}{3}\) \\ i = 1, \(\ldots\), 6 & \(\frac{1}{18}\) \\ i = 7, \(\ldots\), 18 & \(\frac{3}{36}\) \\ \hline \hline \end{tabular}
\end{table} 
Table 2: The weighting factor for the D3Q19 model. Figure 3: Discrete velocity in the D3Q19 model scheme. The simulation consisted of streaming and collision steps, with collisions handled via the Single Relaxation Time (SRT) approach based on the BGK model. This method, known for its numerical stability and reduced sensitivity to viscosity, helps suppress spurious velocity fluctuations [38]. Simulations ran for 30,000 time steps ensured steady-state flow with no-slip conditions at the top and bottom walls enforced using the bounce-back scheme. A pressure gradient along the \(x\)-direction drove the fluid through the porous domain. The halfway bounce-back condition was applied at solid surfaces as validated in previous studies [39; 40]. To minimize numerical slip inherent to BGK models. After reaching convergence, permeability was computed from the velocity and pressure distributions using a Darcy-based formulation adapted for LBM  
\[q=-KA\frac{\Delta P}{\mu L} \tag{9}\]  
Effective porosity was estimated from LBM velocity fields using a velocity-thresholding image processing method. Voxels with non-zero velocities were classified as effective pores, while regions with velocities below a threshold were considered ineffective. The threshold was determined as the 10 th percentile of the lowest velocity distribution, isolating flow-contributing pores from stagnant zones. A binary mask was generated from this segmentation, and effective porosity was computed using  
\[\phi_{eff}=\frac{V_{eff}}{V_{total}} \tag{10}\]  
where \(V_{eff}\) is the volume of effective pore space and \(V_{total}\) is the total pore volume. This effective porosity accounts solely for interconnected pores that are actively traversed by the fluid. Velocity and pressure gradients were obtained from the LBM simulation and visualized using ParaView 5.9.1 with the stream tracers feature for each rock type. Figure 4 illustrates the fluid flow within each rock sample measuring 256 \(\times\) 256 \(\times\) 256 voxels and displays the 3D microstructure of Sandstone and Brazilian pre-salt carbonate rocks, highlighting significant variations in pore distribution and connectivity. Figure 4: Velocity field distributions for each sample of rock types. (**a**) SW01, (**b**) SW02, (**c**) SW03, (**d**) SW04, (**e**) SW05, (**f**) SW06, (**g**) SW07, (**h**) SW08, (**i**) SW09, (**j**) SW10, (**k**) Bandera Brown, (**l**) Bandera Gray, (**m**) Buff Berea, (**n**) Bentheimer, (**o**) Berea, (**p**) BSG, (**q**) BUG, (**r**) Castlegate, (**s**) Kirby, and (**t**) Leopard. These microstructural differences form the basis for the quantitative analysis of the relationship between rock physical properties using 3D data. Geometric tortuosity was computed using [PERSON]'s algorithm, a deterministic weighted shortest-path method applied to the segmented 3D pore structure obtained from micro-CT scans. This algorithm accurately identifies the minimal flow path between inlet and outlet nodes without heuristic approximations, making it well-suited for characterizing tortuous pathways within complex pore networks.