Estimation of hydraulic conductivity from 3D µCT images

[hmuhamme]

I'm processing 3D µCT images with dimensions of (435,435,435) pixels. Instead of using generated cubic image shape, I uploaded the 3D µCT image and inserted to the porespy library. After extracting the pore and throat properties, I imported the extracted pores and throats into the OpenPNM library. I should mention that I manually included different methods and functions in the main Python script. However, the Python script (porespy and openPNM) produce hydraulic conductivity, which is much lower than the hydraulic conductivity, measured in the laboratoty. For example, the script results in hydraulic conductivity K=0.0113 cm/second. While, the laboratory-measured hydraulic conductivity using constant head is: hydraulic conductivity K=0.040. The material used for the experiment, is soilless substrate (Peatmoss). I would appreciate it, if any suggestions for improving the script.

I put the script and the output below:

import numpy as np
import scipy.ndimage as ndi
import tifffile as tiff
import porespy as ps
import openpnm as op
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from skimage.filters import threshold_otsu
import scipy.sparse.csgraph as csgraph

image = tiff.imread('peat_flow.tif')
fil_image = ndi.gaussian_filter(image, sigma=1)
print(f"Image intensity range: min={np.min(fil_image):.2f}, max={np.max(fil_image):.2f}")
print(f"Image shape: {image.shape}")

thresh=75 # manual thresholding
im = (fil_image < thresh).astype(np.uint8) # Pores=1, solids=0 for porespy

Apply binary closing with a larger structure to enhance connectivity

im = ndi.binary_closing(im, structure=np.ones((5, 5, 5)))

plt.hist(fil_image.ravel(), bins=256, alpha=0.5, label='Image')
plt.axvline(x=thresh, color='g', linestyle='--', label=f'Threshold={thresh:.2f}')
plt.title('Intensity Histogram')
plt.legend()
plt.show()

Visualize a slice of the binary image

plt.imshow(im[im.shape[0]//2, :, :], cmap='gray')
plt.title(f'Binary Image Slice (Depth={im.shape[0]//2}, thresh={thresh:.2f})')
plt.show()

Estimate porosity

porosity = np.mean(im)
print(f"Porosity: {porosity:.3f}")

pad_z = (0, 1 if im.shape[0] % 2 != 0 else 0) # Z: 435 → 436
pad_y = (0, 1 if im.shape[1] % 2 != 0 else 0) # Y: 435 → 436
pad_x = (0, 1 if im.shape[2] % 2 != 0 else 0) # X: 435 → 436
im_padded = np.pad(im, pad_width=(pad_z, pad_y, pad_x), mode='constant', constant_values=1)
print(f"Padded image shape for snow2: {im_padded.shape}")
snow = ps.networks.snow2(im_padded, boundary_width=0)
net = snow.network
print(f"Shape after snow2 processing: {snow.regions.shape}")

print("Raw snow2 coords range (x, y, z):", net['pore.coords'].max(axis=0), net['pore.coords'].min(axis=0))

pn = op.io.network_from_porespy(net)
project = pn.project

def trim_disconnected_pores(network):
conns = network['throat.conns']
adj_matrix = np.zeros((network.Np, network.Np))
adj_matrix[conns[:, 0], conns[:, 1]] = 1
adj_matrix[conns[:, 1], conns[:, 0]] = 1
n_components, labels = csgraph.connected_components(csgraph.csgraph_from_dense(adj_matrix, null_value=0), directed=False)
if n_components > 1:
print(f"Before trimming: {n_components} connected components")
component_sizes = np.bincount(labels)
largest_component = np.argmax(component_sizes)
print(f"Largest component: {largest_component}, size: {component_sizes[largest_component]}")
pores_to_keep = labels == largest_component
throats_to_keep = np.all(pores_to_keep[conns], axis=1)
network['pore._trimmed'] = ~pores_to_keep
network['throat._trimmed'] = ~throats_to_keep
op.topotools.trim(network=network, pores=np.where(~pores_to_keep)[0], throats=np.where(~throats_to_keep)[0])
print(f"After trimming: Number of pores = {network.Np}, Number of throats = {network.Nt}")
return network

pn = trim_disconnected_pores(pn)
spacing = 46e-6 # Voxel size in meters (this correcponds to the image resolution= 46µm /voxel
dims = snow.regions.shape
coords = pn['pore.coords'].copy()
pn['pore.coords'][:, 0] *= spacing # x
pn['pore.coords'][:, 1] *= spacing # y
pn['pore.coords'][:, 2] *= spacing # z
if 'pore.equivalent_diameter' in pn:
pn['pore.equivalent_diameter'] *= spacing
else:
print("Warning: 'pore.equivalent_diameter' not found. Adding geometry models.")
pn.add_model_collection(op.models.collections.geometry.spheres_and_cylinders)
pn.regenerate_models()
if 'throat.equivalent_diameter' in pn:
pn['throat.equivalent_diameter'] *= spacing
else:
print("Warning: 'throat.equivalent_diameter' not found. Adding geometry models.")
pn.add_model_collection(op.models.collections.geometry.spheres_and_cylinders)
pn.regenerate_models()
pn['throat.equivalent_diameter'] *= spacing
pn.add_model_collection(op.models.collections.geometry.spheres_and_cylinders)
pn.regenerate_models()

print("Scaled network extent (x, y, z):", pn['pore.coords'].max(axis=0), pn['pore.coords'].min(axis=0))

print("Available network properties:", pn.keys())
print(f"Number of pores: {pn.Np}")
print(f"Number of throats: {pn.Nt}")
print("Pore diameter stats (µm): min =", np.min(pn['pore.equivalent_diameter'])*1e6,
", max =", np.max(pn['pore.equivalent_diameter'])*1e6,
", mean =", np.mean(pn['pore.equivalent_diameter'])*1e6)
print("Coordination number stats: min =", np.min(pn['pore.coordination_number']),
", max =", np.max(pn['pore.coordination_number']),
", mean =", np.mean(pn['pore.coordination_number']))

coords = pn['pore.coords']
z_coords = coords[:, 2] # z-axis for flow direction
y_coords = coords[:, 1]
x_coords = coords[:, 0]
tol = 138e-6 # 3 voxels as tolerance

inlet = pn.pores()[z_coords < (coords[:, 2].min() + tol)]
outlet = pn.pores()[z_coords > (coords[:, 2].max() - tol)]
print(f"Number of inlet pores: {len(inlet)}")
print(f"Number of outlet pores: {len(outlet)}")

print("Inlet z-coordinates range:", pn['pore.coords'][inlet, 2].min(), pn['pore.coords'][inlet, 2].max())
print("Outlet z-coordinates range:", pn['pore.coords'][outlet, 2].min(), pn['pore.coords'][outlet, 2].max())

z_max_range = (dims[0] - 1) * spacing
plt.hist(pn['pore.coords'][inlet, 2], bins=np.linspace(0, z_max_range, 20), alpha=0.5, label='Inlet z-coords')
plt.hist(pn['pore.coords'][outlet, 2], bins=np.linspace(0, z_max_range, 20), alpha=0.5, label='Outlet z-coords')
plt.title('Inlet/Outlet Pore Z-Coordinates')
plt.xlabel('Z-Coordinate (m)')
plt.ylabel('Frequency')
plt.legend()
plt.show()

if len(inlet) == 0 or len(outlet) == 0:
raise ValueError("Inlet or outlet pores are empty. Check network extraction or adjust tolerance.")
phase = op.phase.Phase(network=pn, name='water')
phase['pore.viscosity'] = 1.0e-3

phase.add_model_collection(op.models.collections.physics.basic)

phase.add_model(
propname='throat.hydraulic_conductance',
model=op.models.physics.hydraulic_conductance.hagen_poiseuille
)
phase.regenerate_models()

print("Project contents:", [obj.name for obj in project])
print("Phase properties:", phase.keys())
print("Throat diameter stats (µm): min =", np.min(phase['throat.equivalent_diameter'])*1e6,
", max =", np.max(phase['throat.equivalent_diameter'])*1e6,
", mean =", np.mean(phase['throat.equivalent_diameter'])*1e6)
print("Throat length stats (µm): min =", np.min(phase['throat.length'])*1e6,
", max =", np.max(phase['throat.length'])*1e6,
", mean =", np.mean(phase['throat.length'])*1e6)

flow = op.algorithms.StokesFlow(network=pn, phase=phase)
print(f"Phase name: {phase.name}, Flow settings phase: {flow.settings['phase']}")
flow.set_value_BC(pores=inlet, values=100) # Increase inlet pressure to 20 Pa
flow.set_value_BC(pores=outlet, values=0)
flow.run()
phase.update(flow.soln)

Q = flow.rate(pores=inlet, mode='group')[0]

Adjust domain dimensions for z-direction flow

L = dims[0] * spacing # z-axis length (flow direction)
A = (dims[1] * spacing) * (dims[2] * spacing) # y-x plane area (perpendicular to flow)
mu = 1.0e-3
Delta_P = 100.0 # Reflect the increased pressure gradient
K = Q * L * mu / (A * Delta_P)

print("Domain length (z-axis, flow direction):", L)
print("Domain area (y-x plane):", A)

K_mD = K / 0.987e-12 * 1000
print(f'The permeability (z-direction) is: {K_mD:.2f} mD')

Convert to cm/s

rho = 1000
g = 9.81
k_cm_s = K * (rho * g / mu) * 100
print(f'The hydraulic conductivity (z-direction) is: {k_cm_s:.4f} cm/s')

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

op.visualization.plot_connections(pn, ax=ax)
scatter = op.visualization.plot_coordinates(pn, ax=ax, color_by=phase['pore.pressure'])
ax.set_xlabel('X (m)', labelpad=0)
ax.set_ylabel('Y (m)', labelpad=0)
ax.set_zlabel('Z (m)', labelpad=0)
x_max_range = (dims[2] - 1) * spacing
y_max_range = (dims[1] - 1) * spacing
z_max_range = (dims[0] - 1) * spacing
ax.set_xticks(np.linspace(0, x_max_range, 4))
ax.set_yticks(np.linspace(0, y_max_range, 4))
ax.set_zticks(np.linspace(0, z_max_range, 4))
ax.tick_params(axis='x', labelrotation=45, labelsize=8, pad=0)
ax.tick_params(axis='y', labelrotation=45, labelsize=8, pad=0)
ax.tick_params(axis='z', pad=0, labelsize=8)
zticks = ax.get_zticks()
ax.set_zticklabels([f'{z:.4f}' for z in zticks], rotation=0, fontsize=8)
plt.colorbar(scatter, ax=ax, label='Pressure distribution(Pa)')

plt.savefig('3d_pore_network.tif', dpi=400, bbox_inches='tight')
plt.show()

Outputs:
Porosity: 0.425
Padded image shape for snow2: (436, 436, 436)
Shape after snow2 processing: (436, 436, 436)
Raw snow2 coords range (x, y, z): [431.73943662 431.65517241 431.54716981] [2.08695652 2.14285714 2.33333333]
Before trimming: 315 connected components
Largest component: 0, size: 11561
After trimming: Number of pores = 11561, Number of throats = 34062
Scaled network extent (x, y, z): [0.01986001 0.01985614 0.01983758] [0.00010889 0.00011988 0.00011462]
Available network properties: dict_keys(['throat.conns', 'pore.coords', 'pore.all', 'throat.all', 'pore.region_label', 'pore.phase', 'throat.phases', 'pore.region_volume', 'pore.equivalent_diameter', 'pore.local_peak', 'pore.global_peak', 'pore.geometric_centroid', 'throat.global_peak', 'pore.inscribed_diameter', 'pore.extended_diameter', 'throat.inscribed_diameter', 'throat.total_length', 'throat.direct_length', 'throat.perimeter', 'pore.volume', 'pore.surface_area', 'throat.cross_sectional_area', 'throat.equivalent_diameter', 'pore._trimmed', 'throat._trimmed', 'pore.coordination_number', 'pore.max_size', 'throat.spacing', 'pore.seed', 'pore.diameter', 'throat.max_size', 'throat.diameter', 'throat.hydraulic_size_factors', 'throat.diffusive_size_factors', 'throat.lens_volume', 'throat.length', 'throat.total_volume', 'throat.volume'])
Number of pores: 11561
Number of throats: 34062
Pore diameter stats (µm): min = 137.55223017810158 , max = 4731.400432739629 , mean = 663.1749133936644
Coordination number stats: min = 1.0 , max = 36.0 , mean = 5.8925698468990575
Number of inlet pores: 106
Number of outlet pores: 105
Inlet z-coordinates range: 0.00011462295081967213 0.0002511293900184843
Outlet z-coordinates range: 0.01969962699386503 0.01983758041958042
Project contents: ['net', 'water']
Phase properties: dict_keys(['pore.all', 'throat.all', 'pore.temperature', 'pore.pressure', 'pore.viscosity', 'throat.hydraulic_conductance'])
Throat diameter stats (µm): min = 51.90544168639358 , max = 1990.0846888234728 , mean = 370.9344845913096
Throat length stats (µm): min = 92.96858656444171 , max = 10272.198217326604 , mean = 697.972457583606
Phase name: water, Flow settings phase: water
Domain length (z-axis, flow direction): 0.020056
Domain area (y-x plane): 0.00040224313600000004
The permeability (z-direction) is: 11705.16 mD
The hydraulic conductivity (z-direction) is: 0.0113 cm/s

Best,
Hadi