Permeability¶

Absolute permeability \(k\) measures how easily a single fluid flows through the connected pore space under a pressure gradient. PoreStore computes it by simulating slow, steady flow through the image with a lattice Boltzmann solver, then applying Darcy's law to the result.

Solves the Stokes equations for creeping flow in the pore space, by the lattice Boltzmann method
Backends waLBerla (default) and Kabs Taichi, on GPU
Reports \(k\) in mD and m² (or voxel²), with the steady pressure and velocity fields
How to run it Permeability in the User guide

There are many ways to set up a "permeability simulation", and they don't all give the same number. This page states exactly which problem each backend solves, so you can tell which one matches your experiment, and why two backends can disagree.

The problem being solved¶

In the pore space \(\Omega_p\), the fluid obeys the steady Stokes equations for incompressible creeping flow:

\[ \mu \nabla^2 \mathbf{u} = \nabla p, \qquad \nabla \cdot \mathbf{u} = 0 \qquad \text{in } \Omega_p, \]

with no slip on every pore–solid surface, \(\mathbf{u} = \mathbf{0}\) on \(\partial\Omega_s\). A pressure difference \(\Delta p = p_\text{in} - p_\text{out}\) is applied between two opposite faces of the sample, the inlet and outlet faces normal to the flow direction you choose (X, Y or Z).

Permeability then comes from Darcy's law:

\[ \boxed{\;k = \frac{\mu\, u_D}{|\nabla p|}\;} \qquad\text{with}\qquad u_D = \frac{1}{V}\int_{\Omega_p} u_\parallel \,\mathrm{d}V , \]

where \(u_D\) is the superficial (Darcy) velocity: the flow-direction velocity \(u_\parallel\) summed over the pore voxels and divided by the total number of voxels \(V\), solid included. The mean velocity in the pores is \(u_D/\varphi\), with \(\varphi\) the porosity.

The set-up every backend shares: fixed pressure on the inlet and outlet faces, periodic lateral faces, no-slip on solid. Flow runs along the axis you choose. The image isn't rotated: flow in Y enters at the first plane along the image's second axis.

How the lattice Boltzmann solvers get there¶

Every backend solves the problem with the lattice Boltzmann method (LBM) on the voxel grid itself, one lattice node per voxel, using the D3Q19 velocity set. Each runs until the flow stops changing, and that steady state is the Stokes solution. The common settings are:

Setting	Value	Why
Lattice spacing	1 voxel	Results come out in voxel units and are scaled on read.
Relaxation time	\(\tau = 1\), so the kinematic viscosity is \(\nu = (\tau - \tfrac12)/3 = \tfrac16\)	With single-relaxation-time (BGK) collision, the effective position of a bounce-back wall, and so \(k\), depends slightly on \(\tau\). \(\tau = 1\) is the usual choice that keeps this error small.¹
Viscosity in Darcy's law	\(\mu = \rho_0 \nu = \nu\)	The reference density is \(\rho_0 = 1\), so the backends write Darcy's law with \(\nu\) in place of \(\mu\).
Pressure	\(p = c_s^2 \rho\), with \(c_s^2 = \tfrac13\)	Pressure is imposed through the lattice density.
Driving	\(\Delta\rho = 0.01\) (\(\rho_\text{in} = 1.00\), \(\rho_\text{out} = 0.99\))	Small enough that the flow is creeping: Reynolds and Mach numbers are tiny, so inertia is negligible and the LBM steady state is Stokes flow.

Because Stokes flow is linear, \(k\) doesn't depend on the size of \(\Delta p\) or on \(\mu\). Their only role is to produce a measurable, creeping flow. On this page \(\tau\) is the LBM relaxation time, not the tortuosity factor.

Boundary conditions¶

Inlet and outlet faces¶

Pressure is fixed on the two faces normal to the flow, but each backend fixes it in a slightly different place, and measures the gradient over a different length. On a long, homogeneous sample this doesn't matter. On a short or heterogeneous one it shifts \(k\) by a few percent:

Where pressure is imposed, and how ∇p is measured. waLBerla fits the gradient to the interior pressure field, which excludes entrance and exit effects. Kabs Taichi divides the imposed difference by a length.

Where pressure is imposed, and how ∇p is measured. waLBerla fits the gradient to the interior pressure field, which excludes entrance and exit effects. Kabs Taichi and Palabos divide the imposed difference by a length. Palabos adds open-fluid buffer planes at each end and measures that length between the outermost ones.

Lateral faces: periodic, for every backend¶

The four faces parallel to the flow are periodic

No backend puts walls on the lateral faces. Fluid leaving through the top face re-enters through the bottom face, and the same for front and back. A pore that touches one lateral face is connected to whatever touches the opposite face.

This is the natural choice for a sample cut from a larger, statistically uniform medium: there are no artificial walls to slow the flow near the faces. But your sample isn't really periodic, so the wrap joins pores that aren't neighbours in the real material, and can open or close paths that don't exist:

Left: periodic wrapping turns two half-channels on opposite faces into one full channel. Right: a one-voxel layer of solid on the four lateral faces makes them no-slip walls, which is closer to a sleeved core in a permeameter.

If you want walls, add one layer of solid voxels on the four faces parallel to the flow before uploading, and leave the inlet and outlet faces open. The difference between the periodic and walled answers shrinks as the sample gets wider relative to its pores. A large gap between them means the sample is too narrow to be representative.

Solid surfaces¶

Pore–solid surfaces are no-slip, imposed with bounce-back: a fluid population that would stream into a solid voxel is reflected back the way it came. waLBerla and Kabs Taichi use halfway bounce-back, which puts the wall midway between the last pore node and the first solid node, that is, on the voxel face. The XLB variants use full-way bounce-back. Palabos uses bounce-back nodes on the solid voxels next to pores.

Images that can't flow¶

An image that is all solid or all pore is refused before it runs. Every job, on every backend, also passes PoreStore's spanning-path check, which asks whether some cluster of pore voxels touches both the inlet and the outlet face along the flow axis. Voxels count as connected when they share a face or an edge (18-connectivity, the links a D3Q19 lattice streams along), and a path that only spans by wrapping round a lateral face doesn't count. When no path spans, the launch form refuses the image if it has already been checked, and otherwise the job stops before its solver starts. If the check can't finish within five minutes, the solver runs without it.

The backends¶

waLBerlaKabs Taichi

The default backend: a GPU solver built on the waLBerla framework, with kernels generated by lbmpy.²³

Aspect	Setting
Lattice and collision	D3Q19, single relaxation time (BGK) with lbmpy's incompressible equilibrium, \(\tau = 1\), AA in-place streaming
Inlet / outlet	Fixed density \(\rho = 1.00\) and \(0.99\), imposed on the links that cross the image's inlet and outlet faces (from a ghost layer just outside the image)
Lateral faces	Periodic
Solid	No-slip, halfway bounce-back
Pre-processing	None: no voxels are removed. Before solving, it also checks for itself that some pore path joins inlet to outlet, counting paths that wrap round the lateral faces, and reports \(k = 0\) if none does. PoreStore's spanning-path check normally refuses such an image first.
Pressure gradient	Fitted: \(G\), the least-squares slope of \(\bar p_i\) against \(i\), where \(\bar p_i\) is the pressure \(c_s^2(\rho - 1)\) averaged over the pore voxels of plane \(i\). The \(\max(2, \lfloor N/8 \rfloor)\) planes at each end are left out of the fit.
Permeability	\(k = \nu\, u_D / \lvert G \rvert\)
Convergence	Velocity, permeability and flux-profile criteria, checked every 500 iterations, all passing on 2 consecutive checks. Maximum iterations can be any count up to 10⁶.
Precision	FP32 on GPU, with double-precision reductions
Limits	At most 800 planes along the flow axis, and an inlet and outlet face that each hold some pore
Settings	Maximum iterations \(n_\text{max}\), Velocity tolerance \(\delta_u\), Permeability tolerance \(\delta_k\), Flux-profile tolerance \(\delta_Q\)

Why it fits the gradient: pressure next to the inlet and outlet faces is distorted by entrance and exit effects, where flow converges into or out of the pores on the face. Fitting the slope over the interior measures the gradient the bulk of the sample actually sees. On short samples this is the most faithful measure, and it's why waLBerla can differ by several percent from the backends that divide \(\Delta p\) by the sample length.

A GPU solver from Kabs, an LBM permeability package adapted from Taichi-LBM3D.⁴⁵

Aspect	Setting
Lattice and collision	D3Q19, single relaxation time (BGK), standard second-order equilibrium, \(\tau = 1\)
Inlet / outlet	All 19 populations at pore nodes of the first and last voxel planes are reset each step to equilibrium at \(\rho_\text{in} = 1.00\) / \(\rho_\text{out} = 0.99\), with the velocity copied from the next plane in (a zero-gradient velocity condition)
Lateral faces	Periodic
Solid	No-slip, halfway bounce-back
Pre-processing	None. Isolated pores simply carry no flow.
Pressure gradient	Imposed: \(\lvert\nabla p\rvert = c_s^2 \Delta\rho / N\), where \(N\) is the number of voxel planes along the flow axis
Permeability	\(k = \nu\, u_D / \lvert\nabla p\rvert\)
Convergence	Velocity and permeability criteria, checked every 500 iterations, both passing on 2 consecutive checks
Precision	FP32 on GPU
Settings	Maximum iterations \(n_\text{max}\), Velocity tolerance \(\delta_u\), Permeability tolerance \(\delta_k\). The flux-profile tolerance is ignored.

The density is held at the centres of the first and last planes, which are \(N - 1\) voxel spacings apart, while the gradient divides by \(N\). The imposed gradient also includes the entrance and exit pressure losses. Both effects fade as \(N\) grows.

Palabos¶

A CPU solver built on Palabos,⁶ adapted from LANL's LattEasy. It runs in double precision with MPI. Palabos is available to PoreStore staff only for now, and runs only on a compute tier you name, never on Auto.

Aspect	Setting
Lattice and collision	D3Q19, BGK, \(\tau = 1\)
Inlet / outlet	Two planes of open fluid are added before the inlet and after the outlet. Palabos's local pressure boundary condition is applied on the outermost added planes, with \(\rho_\text{in} = 1\) and \(\rho_\text{out} = 1 - 3\Delta p\), where \(\Delta p = 0.0005\) in lattice units.
Lateral faces	Periodic
Solid	Bounce-back on every solid voxel next to a pore and every solid voxel on the image faces. Solid further inside is inert.
Pre-processing	None. Isolated pores simply carry no flow.
Pressure gradient	Imposed: \(\lvert\nabla p\rvert = \Delta p / (N + 3)\), across the image and the four added planes
Permeability	\(k = \nu\,\bar u_\parallel / \lvert\nabla p\rvert\), with \(\bar u_\parallel\) the mean flow-direction velocity over the whole padded box. At steady state this equals \(u_D\), because the flux through every plane is the same.
Convergence	Permeability only: every 1000 iterations, the relative change in \(\bar u_\parallel\), which equals \(r_k\) because \(k\) is proportional to it, must fall below \(\delta_k\) on two consecutive checks. \(\delta_u\) and \(\delta_Q\) are ignored.
Precision	FP64 on CPU

The fluid buffers give the pressure condition a uniform plane to act on. Their small resistance is counted in the length but not in the geometry, so \(k\) comes out slightly low on short samples.

Earliest stop is iteration 3000

Palabos compares velocity samples 1000 iterations apart, and the first sample is only a baseline, so the earliest possible stop is iteration 3000. A run that ends before two passing comparisons is marked not converged, and so is a run whose mean velocity stays at zero, such as a sample with no flow path.

Kabs XLB-JAX and XLB-Warp¶

Kabs can also run on XLB,⁷ Autodesk's differentiable LBM library, with a JAX or an NVIDIA Warp back end. The physics is Kabs Taichi's, with a different inlet and outlet treatment and a different wall rule:

Aspect	Setting
Lattice and collision	D3Q19, BGK, \(\tau = 1\), FP32 on GPU
Inlet / outlet	Zou–He pressure condition⁸ at \(\rho = 1.00\) and \(0.99\) on the pore nodes of the first and last voxel planes. Pore nodes on the edges of those faces get an equilibrium condition with \(\mathbf{u} = 0\).
Lateral faces	Periodic
Solid	Full-way bounce-back on every solid voxel
Pre-processing	None
Pressure gradient	Imposed: \(\lvert\nabla p\rvert = c_s^2 \Delta\rho / N\), as Kabs Taichi
Permeability	\(k = \nu\, u_D / \lvert\nabla p\rvert\)
Convergence	As Kabs Taichi: velocity and permeability criteria, every 500 iterations, 2 consecutive passes

XLB is in testing and isn't offered on every team yet.

Comparison¶

	waLBerla	Kabs Taichi
Hardware	GPU, FP32	GPU, FP32
Lattice / collision	D3Q19 BGK, incompressible eq.	D3Q19 BGK
Pressure imposed	on the image faces	on the first and last voxel planes
\(\lvert\nabla p\rvert\) from	fit to the interior pressure	\(\Delta p / N\)
Lateral faces	periodic	periodic
Solid	halfway bounce-back	halfway bounce-back
Disconnected pores	kept	kept
Convergence	velocity + \(k\) + flux profile, ×2	velocity + \(k\), ×2

Palabos differs from both: it runs on a CPU in FP64, imposes pressure on two added fluid planes outside each face, divides by \(N + 3\), puts bounce-back nodes on the solid next to pores, and tests only \(k\).

Accuracy and convergence¶

LBM reaches steady state by iterating. PoreStore checks progress every 500 iterations (Palabos: every 1000) and stops when the criteria pass, or after \(n_\text{max}\) iterations. A superscript \((j)\) marks a quantity at check \(j\), and \((j-1)\) at the check before it. The Kabs and waLBerla backends test the relative change in the velocity field and in \(k\):

\[ r_u = \frac{\displaystyle\sum_{\mathbf{x} \in \Omega_p} \sum_{c \in \{x,y,z\}} \bigl\lvert u_c^{(j)}(\mathbf{x}) - u_c^{(j-1)}(\mathbf{x}) \bigr\rvert}{\displaystyle\sum_{\mathbf{x} \in \Omega_p} \sum_{c \in \{x,y,z\}} \bigl\lvert u_c^{(j)}(\mathbf{x}) \bigr\rvert} < \delta_u, \qquad r_k = \frac{\bigl\lvert k^{(j)} - k^{(j-1)} \bigr\rvert}{\bigl\lvert k^{(j)} \bigr\rvert} < \delta_k . \]

waLBerla also tests the flux profile, the flux \(Q_i\) through each flow-normal plane \(i\) of the \(N\) planes:

\[ r_Q = \frac{\displaystyle\sum_{i=1}^{N} \bigl\lvert Q_i^{(j)} - Q_i^{(j-1)} \bigr\rvert}{\displaystyle\sum_{i=1}^{N} \bigl\lvert Q_i^{(j)} \bigr\rvert} < \delta_Q . \]

The result reports these as Velocity change (\(r_u\)), Permeability change (\(r_k\)) and Flux profile spread (\(r_Q\)), each against the tolerance you set: Velocity tolerance \(\delta_u\) and Permeability tolerance \(\delta_k\), both 0.0001–0.001 with a default of 0.001, and Flux-profile tolerance \(\delta_Q\), 0.0001–0.001 with a default of 0.0001. Maximum iterations is the cap \(n_\text{max}\). The check interval is fixed at 500 iterations (Palabos: 1000).

Every criterion must pass on two consecutive checks, so the earliest possible stop is iteration 1500. The iteration cap needn't be a multiple of 500: the iterations after the last check still run, but no check follows them. The flux-profile criterion catches a slow drift where the total flow looks settled but some planes are still filling, which shows up first in samples with dead-end pores.

A run that reaches \(n_\text{max}\) still completes

It's saved and marked not converged, with the last value of each criterion in the result. The value of \(k\) may still be useful, but check the criteria before relying on it. Rerun with more iterations if needed.

Resolution and accuracy¶

LBM accuracy depends on how many voxels span each pore throat. Bounce-back places walls to within about a voxel, so a throat only two or three voxels wide carries a large relative error in its effective width, and its conductance goes roughly as width to the fourth power. In practice:

Results converge as resolution increases. If narrow throats control the flow, check a downsampled copy: if \(k\) moves a lot between resolutions, the finer one isn't fully resolved yet either.
Single-voxel throats pass flow, but their conductance is a lattice artefact, not physics.
Backends differ by a few percent on real samples, mostly because of how they measure \(\nabla p\). On an analytical channel, waLBerla matches the exact solution to within about 1%, and Kabs Taichi agrees with waLBerla to 0.7%. Compare samples within one backend, and say which backend you used when you report a value.
REV. A sample has to be several pore diameters wide in every direction before \(k\) stops depending on where it was cut. Run REV porosity or compare subvolumes to check.

Units¶

Each solver works on a unit lattice, so its permeability is in voxel², \(k_\text{vox}\). Stokes flow is linear and has no intrinsic length scale, so the same geometry at voxel size \(h\) has

\[ k = k_\text{vox}\, h^2 . \]

PoreStore stores \(k_\text{vox}\) and applies \(h\) whenever you read the result, in m² and in millidarcy, with \(1\ \text{mD} = 9.869233\times10^{-16}\ \text{m}^2\). Fixing a wrong voxel size corrects \(k\) immediately, without a rerun. See Units and voxel size.

Reading the result¶

Headline. \(k\) in mD, or in voxel² when the image has no voxel size.
Check convergence first. A run that reached its iteration cap completes with its last criteria in the result; read them before relying on \(k\).
The fields. PoreStore saves the pressure field and the three velocity components on the image grid, in lattice units. waLBerla stores the gauge pressure \(c_s^2(\rho - 1)\) with solid masked. Kabs and Palabos store the lattice density as a pressure proxy. In a converged run the flux through every plane normal to the flow is the same.
Isolated and dead-end pores carry no flow but count in the porosity. The result says how much of the pore space was found isolated or dead-end along the flow before the solver ran.
Anisotropy. Run all three directions; different \(k\) along different axes is a property of the material, or a sign the sample is too small.

Kabs Taichi's saved fields on the middle slice of a 64³ synthetic sample, porosity 0.55, with flow along x from left to right. The sample was run through PoreStore's own solver code with the app's suggested settings. Solid is grey. The pressure proxy is the lattice density, which falls from inlet to outlet. Hover over a panel to read its values.

Choosing a backend¶

waLBerla for most work: GPU-fast, a gradient measure that's robust on short samples, and the strictest convergence test.
Kabs Taichi to compare with Kabs, or with other studies that impose the gradient as \(\Delta p / L\), and for images longer than 800 voxels along the flow.
beta Palabos to compare with Palabos- or LattEasy-based workflows, or when you want double precision.

References¶

C. Pan, L.-S. Luo, C. T. Miller, "An evaluation of lattice Boltzmann schemes for porous medium flow simulation," Computers & Fluids 35 (2006) 898–909. doi:10.1016/j.compfluid.2005.03.008 ↩
M. Bauer et al., "waLBerla: A block-structured high-performance framework for multiphysics simulations," Computers & Mathematics with Applications 81 (2021) 478–501. doi:10.1016/j.camwa.2020.01.007 ↩
M. Bauer, H. Köstler, U. Rüde, "lbmpy: Automatic code generation for efficient parallel lattice Boltzmann methods," Journal of Computational Science 49 (2021) 101269. doi:10.1016/j.jocs.2020.101269 ↩
J. Yang, Y. Xu, L. Yang, "Taichi-LBM3D: A single-phase and multiphase lattice Boltzmann solver on cross-platform multicore CPU/GPUs," Fluids 7 (2022) 270. doi:10.3390/fluids7080270 ↩
Y. Hu et al., "Taichi: a language for high-performance computation on spatially sparse data structures," ACM Transactions on Graphics 38 (2019) 201. doi:10.1145/3355089.3356506 ↩
beta J. Latt et al., "Palabos: Parallel Lattice Boltzmann Solver," Computers & Mathematics with Applications 81 (2021) 334–350. doi:10.1016/j.camwa.2020.03.022 ↩
beta M. Ataei, H. Salehipour, "XLB: A differentiable massively parallel lattice Boltzmann library in Python," Computer Physics Communications 300 (2024) 109187. doi:10.1016/j.cpc.2024.109187 ↩
beta Q. Zou, X. He, "On pressure and velocity boundary conditions for the lattice Boltzmann BGK model," Physics of Fluids 9 (1997) 1591–1598. doi:10.1063/1.869307 ↩