Diffusivity¶

The effective diffusivity $D_\text{eff}$ measures how fast a species diffuses through the pore space compared with diffusion in free fluid, $D_0$. Their ratio folds together two effects: how much of the volume is pore (porosity $\varphi$), and how convoluted the paths through it are (the tortuosity factor $\tau$). The same calculation gives the electrical conductivity of an electrolyte-filled pore space with an insulating solid, and heat conduction through a non-conducting matrix.

Solves Laplace's equation for the concentration in the pore space
Backends Tortuosity.jl (default) and TauFactor, on GPU
Reports $D_\text{eff}/D_0$, the tortuosity factor $\tau$ and the formation factor $F$, with the steady concentration field
How to run it Diffusivity in the User guide

The problem being solved¶

Steady Fickian diffusion with a uniform diffusivity $D_0$ in the pore space reduces to Laplace's equation for the concentration $c$:

\[ \nabla^2 c = 0 \qquad \text{in } \Omega_p , \]

with

Boundary	Condition
Inlet face	$c = 1$ (Dirichlet)
Outlet face	$c = 0$ (Dirichlet)
Lateral faces	no flux, $\partial c/\partial n = 0$
Pore–solid surfaces	no flux, $\partial c/\partial n = 0$

Top: the problem both backends solve. The lateral faces are closed, unlike permeability, where they're periodic. Bottom: the backends discretise the inlet and outlet differently, which changes the effective sample length by one voxel.

Closed lateral faces, unlike permeability

Diffusion treats the four faces parallel to the transport direction as impermeable walls, so nothing crosses them. The permeability backends make them periodic instead. A pore that touches a lateral face is a dead end for diffusion, but can be a through-path for flow.

From the steady solution, the total flux $J$ through the sample gives the effective diffusivity, normalised by the whole cross-section $A$ (solid included) and the length $L$:

\[ \frac{D_\text{eff}}{D_0} = \frac{J\,L}{A\,D_0\,\Delta c} \qquad\text{and}\qquad \tau = \frac{\varphi}{D_\text{eff}/D_0} \quad\Longleftrightarrow\quad \frac{D_\text{eff}}{D_0} = \frac{\varphi}{\tau}. \]

$\varphi$ is the total porosity, every pore voxel over every voxel, including pores that carry no flux.
$D_\text{eff}/D_0$ already contains the porosity: it's the bulk effective diffusivity, the quantity to use in a continuum model of the whole material.
The formation factor is its reciprocal, $F = D_0/D_\text{eff} = \tau/\varphi$. Archie introduced it for electrical conduction, as the resistivity of a rock saturated with brine over the resistivity of the brine,¹ which is the same ratio for a conducting fluid in an insulating solid. PoreStore reports it next to $\tau$; here $F$ means the formation factor, not a cumulative fraction.

What "tortuosity" means here¶

PoreStore's tortuosity is the tortuosity factor $\tau$, as defined above and as used by TauFactor and most of the electrochemistry literature.² It's a transport quantity, not a geometric path-length ratio. Some authors define a geometric tortuosity $\tau_g = L_\text{path}/L$ and relate the two by $\tau = \tau_g^2$. Check which convention a paper uses before comparing.³

Because $\varphi$ counts every pore voxel, isolated and dead-end pores raise $\tau$: they add porosity but carry no flux. A sample with no pore path from inlet to outlet has $D_\text{eff} = 0$. Before either backend runs, PoreStore checks for such a path, counting pore voxels as connected only when they share a face (6-connectivity, the links both solvers exchange flux across), and refuses the job when none spans. If the check can't finish within five minutes, the solver runs without it, and a sample with no through-path comes back with $\tau = 0$, a marker for "no through-path", not a physical value.

The backends¶

Tortuosity.jlTauFactor

The default: Tortuosity.jl, a Julia package for GPU-accelerated tortuosity and diffusivity calculations.

Aspect	Setting
Discretisation	Cell-centred finite volume with a 7-point stencil: each pore voxel is linked to its six face neighbours with unit conductance
Dirichlet planes	$c = 1$ and $c = 0$ at the centres of the first and last voxel planes, so $L = N - 1$
Lateral / solid	No flux: links leaving the image or entering solid are omitted
Flux	Measured on the links leaving the inlet plane: $D_\text{eff}/D_0 = J (N-1)/A$
Linear solver	Conjugate gradient, matrix-free, with a two-level coarse-space preconditioner (8³-voxel aggregates) above 10⁵ unknowns
Precision	FP32 on GPU, refined against the FP64 residual
Pre-processing	None. Isolated pores stay in the domain and in $\varphi$.
Settings	Target accuracy $\delta_J$ and Solver tolerance $\delta_c$

Convergence: solve, then solve to the accuracy target. Write the discretised problem as the linear system $\mathbf{K}\mathbf{c} = \mathbf{b}$, with one unknown concentration per pore voxel.

First solve. The linear solve stops when the solver residual, relative to $\lVert \mathbf{b} \rVert$, drops below the Solver tolerance $\delta_c$: $$ r_c = \frac{\lVert \mathbf{b} - \mathbf{K}\mathbf{c} \rVert}{\lVert \mathbf{b} \rVert} < \delta_c . $$ If the solver can't reach $\delta_c$, the job fails.
Flux balance error. On the solved field, the fluxes $J_i$ through the $N-1$ interfaces between neighbouring planes are compared: $$ r_J = \frac{\max_i J_i - \min_i J_i}{\max_i J_i} . $$ At steady state every plane carries the same flux, so the flux balance error measures how far the field is from converged, and it bounds the relative error in $\tau$ and $D_\text{eff}$. The target is $r_J < \delta_J$, where $\delta_J$ is the Target accuracy you ask for: 0.001 accepts a 0.1% error. A small residual doesn't always imply a small flux balance error, which is why the flux balance error is the test that decides.
Accuracy tightening. If $r_J \ge \delta_J$ after the first solve, the solver keeps going from the field it has. Each round computes the residual in FP64, solves for a correction that reduces it further, adds the correction and measures $r_J$ again. It stops when $r_J < \delta_J$, after two rounds, when a round no longer lowers $r_J$, or when its allowance runs out: extra solving time of at most half the time the job had already spent setting up and solving.

A run that stops short of $\delta_J$ still completes, marked not converged, with the flux balance error it reached. Rerun with a smaller $\delta_c$ if you need the target met.

TauFactor 1.2.1 from the TLDR group at Imperial College, on PyTorch,⁴⁵ with one change: before solving, it checks whether any pore path joins the inlet and outlet faces.

Aspect	Setting
Discretisation	Cell-centred finite difference with a 6-neighbour stencil
Dirichlet planes	$c = 1$ and $c = 0$ on the outer faces of the image, half a voxel outside the first and last voxel centres, so $L = N$
Lateral / solid	No flux
Flux	Mean of the fluxes through the $N-1$ interfaces between planes: $D_\text{eff}/D_0 = \bar J\, N / A$
Iterative solver	Red–black successive over-relaxation, $\omega = 2 - \pi/(1.5\,N)$
Precision	FP32 on GPU
Pre-processing	None
Settings	Target accuracy $\delta_J$ and Maximum iterations $n_\text{max}$

Convergence, checked every 100 iterations, so iteration counts are multiples of 100. With $(j)$ marking check $j$:

the flux balance error $r_J$, defined as for Tortuosity.jl, below $\delta_J$, and
a settled tortuosity factor, $\bigl\lvert \tau^{(j)} - \tau^{(j-1)} \bigr\rvert < 2\times10^{-3}$. This threshold is fixed and can't be changed.

The solver stops after $n_\text{max}$ iterations if it hasn't converged, and the job fails. On an image with no pore path from inlet to outlet, it stops at its first convergence checks and reports $D_\text{eff} = 0$ and $\tau = 0$.

Why the two backends differ slightly¶

Both solve the same equation with the same boundary conditions, but they put the fixed concentrations in different places (see the figure). On an open, empty domain both give exactly $\tau = 1$. On a real microstructure the half-voxel shift at each end is a small change in the problem itself, and it shows up as a systematic difference of order $1/N$. In our calibration runs on 250³–650³ samples, the two agreed to within about 0.1% on $\tau$. The largest gaps were in low-porosity samples, where TauFactor ran slightly low. Neither choice is wrong, but pick one backend for a study and use it throughout.

Accuracy and convergence¶

The two backends handle a missed target differently, and both differ from permeability, which keeps a result that didn't converge:

Tortuosity.jl completes with the accuracy it reached when its extra solving runs out (see the backends). It fails only if the linear solver can't reach $\delta_c$.
TauFactor fails if it hasn't converged after $n_\text{max}$ iterations.

The flux balance error in the result is the accuracy you got, whichever backend ran. The saved concentration is stored in single precision, rounded to about 4–5 significant digits.

Reading the result¶

Straight pores give $\tau = 1$. Pores running straight along the transport axis give 1 with either backend. Anything that lengthens or narrows the paths, or adds pore volume that carries no flux, raises $\tau$.
The formation factor is the same number inverted. $F = 1/(D_\text{eff}/D_0)$ carries no information beyond $D_\text{eff}/D_0$, but it's the form to compare with a resistivity measurement on the saturated sample or with Archie's law, $F = \varphi^{-m}$, whose cementation exponent $m$ is often fitted to it. A sample with no through-path has no formation factor, because $F$ would be unbounded, and the result table leaves the row out.
Porosity is in the result. $\varphi = \tau \cdot D_\text{eff}/D_0$ recovers the porosity the job used.
The concentration field shows where transport is hard. In a homogeneous sample the plane-mean concentration falls roughly linearly from inlet to outlet. A steeper drop across a few planes marks a layer that controls $\tau$, which the single number hides. With Tortuosity.jl the first and last planes are exactly 1 and 0; with TauFactor they aren't, because its fixed values sit on the outer faces. Concentration in a pore cluster that touches neither the inlet nor the outlet isn't determined by the problem, so don't interpret it.
Anisotropy. Run all three directions; different $\tau$ along different axes is a property of the material.

TauFactor's saved field for a 64³ synthetic sample, porosity 0.55, with transport along x. The sample was run through PoreStore's own solver code with the app's suggested settings. Left: the concentration on the middle slice, with solid in grey. Right: its mean over the pore voxels of each plane. Because TauFactor fixes the values on the outer faces, the first and last planes aren't exactly 1 and 0.

Choosing a backend¶

Tortuosity.jl for most work: it meets an accuracy target you set rather than an iteration count, completes with the accuracy it reached instead of failing, and takes larger images on the same GPU.
TauFactor to compare with TauFactor-based studies, whose boundary placement it shares.

Conductivity solves the same problem in the solid, with the same backends and settings.

References¶

G. E. Archie, "The electrical resistivity log as an aid in determining some reservoir characteristics," Transactions of the AIME 146 (1942) 54–62. doi:10.2118/942054-G ↩
B. Tjaden, D. J. L. Brett, P. R. Shearing, "Tortuosity in electrochemical devices: a review of calculation approaches," International Materials Reviews 63 (2018) 47–67. doi:10.1080/09506608.2016.1249995 ↩
N. Epstein, "On tortuosity and the tortuosity factor in flow and diffusion through porous media," Chemical Engineering Science 44 (1989) 777–779. doi:10.1016/0009-2509(89)85053-5 ↩
S. J. Cooper, A. Bertei, P. R. Shearing, J. A. Kilner, N. P. Brandon, "TauFactor: An open-source application for calculating tortuosity factors from tomographic data," SoftwareX 5 (2016) 203–210. doi:10.1016/j.softx.2016.09.002 ↩
S. Kench, I. Squires, S. J. Cooper, "TauFactor 2: A GPU accelerated python tool for microstructural analysis," Journal of Open Source Software 8 (2023) 5358. doi:10.21105/joss.05358 ↩

Boundary	Condition
Inlet face	\(c = 1\) (Dirichlet)
Outlet face	\(c = 0\) (Dirichlet)
Lateral faces	no flux, \(\partial c/\partial n = 0\)
Pore–solid surfaces	no flux, \(\partial c/\partial n = 0\)