Received: 2 November 2023

Accepted: 9 May 2024

DOI: 10.1002/vzj2.20354
Vadose Zone Journal

ORIGINAL ARTICLE
S p e c i a l S e c t i o n : Tr i b u t e t o R i e n v a n G e n u c h t e n , R e c i p i e n t o f t h e 2 0 2 3 Wo l f P r i z e f o r A g r i c u l t u r eMultiscale pore-network reconstruction of a fine-textured
heterogeneous soil
Elizabeth M. Pontedeiro1,2Martinus Th. van Genuchten1,3
Mateus G. Ramirez4
Carlos M. P. Vaz5
Silvio Crestana5
Paulo Couto2

William Godoy2
Maira C. O. Lima2

Jian Su3

1 Department of Earth Sciences, Utrecht

University, Utrecht, The Netherlands
2 Department of Civil Engineering, Federal

University of Rio de Janeiro, Rio de Janeiro,
Rio de Janeiro, Brazil
3 Department of Nuclear Engineering,Federal University of Rio de Janeiro, Rio deJaneiro, Rio de Janeiro, Brazil
4 Department of Mechanical Engineering,Federal University of Rio de Janeiro, Rio deJaneiro, Rio de Janeiro, Brazil
5 Embrapa Agricultural Instrumentation, SãoCarlos, São Paulo, Brazil
Correspondence
Elizabeth M. Pontedeiro, Department ofEarth Sciences, Utrecht University,
Princetonlaan 8a, 3584 CB, Utrecht, The
Netherlands. Email:bettinadulley@hotmail.com
Assigned to Associate Editor Jirka Simunek.

Abstract
Digital samples offer many opportunities to study subsurface fluid flow and contaminant transport processes. The pore size distribution of especially fine-textured
porous media often covers many orders of magnitude in the length scale, whichmakes accurate microCT scanning and modeling of the underlying processes difficult. When a single-resolution image is not capable of capturing all relevant details ofa sample, one should scan the sample, or selected parts of it, at different resolutions.
Combining multiple resolutions into one single sample for subsequent pore-scalemodeling is generally not possible due to limitations in computer memory and speed,
thus making it necessary to create a simpler sample containing relevant information from the parent networks. We imaged four samples using different resolutions
to capture the multiscale heterogeneity of a fine-textured soil and combined theminto one overall digital sample based on the original pore networks. The parent networks were characterized using their geometrical properties, correlations between
these properties, and connectivity functions describing the network topologies. Our
approach creates stochastic networks of arbitrary size with the same flow properties
as the parent network. The method, implemented using the PoreStudio pore networkmodel, repeatedly integrates information at two subsequent scales, with the resulting digital sample having the same hydraulic properties as the original samples. The
procedure leads to more useful three-dimensional digital models, facilitating basic
analyses of underlying pore size distributions. Porosity calculations were comparedwith direct measurements, while those for the hydraulic conductivity were comparedwith estimates based on the particle size distribution and nearby field data.

Abbreviations: μCT, computed microtomography; CV, coefficient of variation; DPA, digital pore analysis; PNM, pore network model; PSD, pore size
distribution; REV, representative elementary volume.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original
work is properly cited.
© 2024 The Author(s). Vadose Zone Journal published by Wiley Periodicals LLC on behalf of Soil Science Society of America.
Vadose Zone J. 2024;e20354.
https://doi.org/10.1002/vzj2.20354

wileyonlinelibrary.com/journal/vzj2

1 of 13

2 of 13

1

PONTEDEIRO ET AL.

Vadose Zone Journal

INTRODUCTION

A major motivation of pore-scale imaging using X-ray computed microtomography (μCT) is to estimate properties ofporous media using nondestructive methods. Improvements
in high-resolution digital imaging and computerized image
analysis now enable one to characterize the microstructure ofa range of natural and engineered porous materials, including soils and rocks (e.g., Allaire et al., 2009; Cnudde &Boone, 2013; Godoy et al., 2019; Katuwal et al., 2018; Limaet al., 2020; Salek & Beckingham, 2023; Vaz et al., 2014;
Wildenschild & Sheppard, 2013).
MicroCT has become popular since it provides threedimensional (3D) volumetric imaging data of undisturbed
samples. Pore-scale characterization of porous media is
important for many practical applications involving multiphase fluid flow and the transport of a range of agriculturaland industrial contaminants. Flow and transport processes insuch applications are closely related to pore size, pore shape,
and the connectivity between pores. X-ray μCT provides
an excellent non-invasive means for characterizing dynamic
pore-scale processes (e.g., Van Offenwert et al., 2019). The
digital reconstruction of image stacks permits representationsof the internal pore structure of a sample, the connection ofpores within the pore network, and the underlying flow andtransport processes. Many different approaches can be used tomodel pore-scale processes, such as various lattice Boltzmann
formulations, coupled computational fluid dynamics-discreteelement methods (Elrahmani et al., 2022), and pore networkmodels (PNMs) (Raoof et al., 2013). For our studies, weused a PNM because it is computationally less demanding
than direct methods and permits one to incorporate moreheterogeneity when modeling larger volumes.
Porosity and permeability are key porous media parameters affecting fluid flow and contaminant transport processes.
Porosity can be estimated from 3D imaging data by convertinggrayscale images into binary ones representing solid material and void spaces. Many previous studies indicate that amore complete characterization of the pore size distribution(PSD) of μCT images requires the estimation of pore sizes
taken from multi-resolution images, and then combining theinformation from each resolution to obtain the complete PSD(Cássaro et al., 2017; Jiang et al., 2007; Vaz et al., 2011; Yanget al., 2009). Thus, the first step of our work was to evaluate images of soils generated using μCT and, according to
the property being studied, to define if possible an optimal
imaging resolution for comparison with the measurements.
Soils and rocks generally consist of a range of microand macropores, which together define the overall connected
pore structure of the medium. The differently sized poresare known to contribute to distinct yet interacting water flow
and contaminant transport processes. In saturated media, thelarger pores will favor overall fluid flow processes (includingCore Ideas∙ A stochastic method was used to combine four
computed microtomography samples with different resolutions into one digital sample.
∙ The approach was applied to a fine-textured tropical soil with a very broad pore-size distribution.
∙ The methodology kept the same digital porosity as
the original pore networks and enhanced the flow
properties.

preferential flow) and contribute the most to the permeability of the sample, while the smaller pores are important forcapillarity, fluid retention, diffusion, and chemical reaction
processes. Recent advances in CT facilitate the acquisitionof high-resolution images of 3D porous media, thus permitting more precise studies about the double- or triple-porosity
characteristics of macro-porous soils and fractured rock (Jiang
et al., 2013; Katuwal et al., 2018; Ruspini et al., 2021; Sadeghnejad & Gostick, 2020; Wildenschild & Sheppard, 2013;
among others).
A major limitation of image-based soil and rock porosity measurements derived solely from μCT data is that pores
smaller than the voxel size cannot be detected (Cássaro et al.,2017; Jiang et al., 2012; Pak et al., 2016). Using existingμCT instrumentation to properly characterize heterogeneous
samples having PSDs spanning several orders of magnitude
is especially challenging since the evaluated “digital” porosity/permeability is generally smaller than measured in thelaboratory. The resulting underestimation of porosity canlead to errors in the water retention and hydraulic conductivity/permeability calculations, and hence fluid flow andcontaminant transport processes in porous media. To compensate for this, advanced image processing techniques must
be used for soils and rocks to account for their possible
dual porosity properties. For dual- or triple-porosity media,
which includes many tropical soils (e.g., oxisols and latosols)and carbonate rocks, new imaging and modeling techniques
have been introduced to capture the full range of the pore
structure features. These include approaches using machinelearning (Dragonfly, 2022) or cellular automata and microcomputed tomography (Zubeldia et al., 2016), among others.
The simulation and modeling of physical parameters using
non-destructive techniques such as high-resolution tomography opens new opportunities for dynamic studies at different
scales, like synchrotron-based tomography (Menke et al.,2016), including the dynamics of unsaturated flow in soils.
One challenge is the integration and interpretation of results at
different scales (downscaling and upscaling), such as from theuse of different types of μCTs (Crestana & Cruvinel, 2023).

PONTEDEIRO ET AL.

MicroCT imaging techniques suffer from limitations interms of a trade-off between resolution and sample size.
In practice, the images that are obtained via μCT may beunable to capture all aspects of the pore system of any
realistic medium. The issue is twofold: some porous media
contain noticeable fractions of micro-porosity that cannot
be resolved at micron-scale resolutions, while using a submicron resolution implies that the resulting sample volume is
too small to capture the heterogeneity of the macro-porosity.
One approach is to combine images taken at different resolutions into one global image (Hebert et al., 2014; Jiang et al.,2012; Kastner et al., 2021; Vaz et al., 2014; Wu et al., 2011;
Yao et al., 2023). For example, Jiang et al. (2012) developed
a stochastic approach to combine the observed networks at
different resolutions into a simpler overall network with very
similar overall properties.
While pore network modeling is considered to be an important tool for simulating complex physical processes at the porescale, Pore network model (PNM) predictions often do notagree well with the corresponding laboratory measurements,
especially when complex pore systems are present (Pak et al.,2016). PNMs simplify the geometry of the pore space andgenerally cannot capture all of the characteristics of porous
media due to a lack of resolution. The trade-off between
sample size and μCT image resolution means that a single
μCT volume may not capture all details of a multiscale pore
structure. The realism and precision of the modeling results
depends on the constructed global pore network assumed to
be its complete representation (Xiong et al., 2016), as well as
whether the overall sample can be considered a representative
elementary volume (REV) of its physical properties.
The REV is defined as the minimum volume of a certain property that is large enough to capture a representative
amount of heterogeneity, or the minimum volume in which
the property is insensitive to small changes in volume or sample location. Location here means that the REV is the largest
quantity insensitive to small changes in volume throughout
the entire medium. The REV may vary for different properties(Vik et al., 2013), with the porosity REV usually being smaller
than the permeability REV because of macroscopic variability
in the latter. Selecting a region with proper REV characteristics is as important as the invoked mathematical model toobtain good results with adequate computational cost (Blunt
et al., 2013; Mehmani et al., 2020; Piccoli et al., 2019; Verriet al., 2017). Since fluid flow only occurs in the connectedporosity, 50% of the flow process in many soils and rocks
may well occur in only 10% of the pores (e.g., Mehmani et al.,2020).
Main objective of this paper is to use images at different resolutions to provide global estimates of porosity andpermeability, and the resulting fluid flow process. Digital
soil samples were generated using μCT images and the main
characteristics of the imaged soil (the skeleton components),

Vadose Zone Journal

3 of 13and then compared to those of the original sample. Multiscale imaging and pore-scale modeling, along with statistical
functions based on the original pore networks, were used toobtain the new digital soil sample. The resulting 3D networks
hence were combined and transformed into one multiscalepore network. For this, we used the skeleton of the lower resolution (sample BS4) and added the micro- and meso-porosity
characteristic (pore throats and pore bodies) from the higher
resolution samples such that the multiscale pore system of thedigital sample honored the average connectivity of the sample. The digital soil sample was generated numerically using
the pore-scale network model within the PoreStudio software(PoreStudio, 2023), an update of the PoreFlow software ofRaoof et al. (2013).

2

EXPERIMENTAL

This section describes the measurements, μCT acquisition, thenumerical simulations, and other analyses that were carriedout using the soil samples.

2.1
Sample selection and soil physicalproperties
The μCT and modeling studies were carried out using four
samples taken from an undisturbed block of a Brazilian
tropical oxisol (Rhodic Hapludox). Sample selection, their
imaging with μCT, and the laboratory measurements were carried out by Embrapa/CNPDIA. Complete descriptions of thesoils and experimental procedures are provided by Vaz et al.
(2011, 2014); they are briefly described below.
A cubic sample of undisturbed soil, 20 cm on each side,
was collected at a depth of 20 cm below the surface at the
Embrapa Southeast Cattle Experimental Farm in the city ofSão Carlos, Brazil. The sample was dried at room temperature for 30 days and then cut into two equal parts. One partwas used for measurements of physical properties such as thesoil particle density (using the pycnometer method), the bulk
density, and the particle size distribution (data in Vaz et al.
[2011, 2014] and Cássaro et al. [2017]), and the other partwas used to obtain several samples in the form of cylinders ofdifferent sizes to perform the μCT analyses with a series ofresolutions. Mercury intrusion porosimetry was used to measure the total porosity (0.502 cm3 /cm3 ), while the bulk density
was estimated to be 1.13 g/cm3 . The particle size distribution
by weight was found to be 12% sand, 9% silt, and 79% clay.
Embrapa/CNPDIA used four different sample sizes (cylinders) to obtain the soil μCT images in four distinct resolutions.
Sample sizes in terms of their diameter and height were 0.5 ×
0.5 cm (sample BS1), 1.0 × 2.0 cm (sample BS2), 2.0 × 2.0 cm(sample BS3), and 4.0 × 2.5 cm (sample BS4). The samples

4 of 13

PONTEDEIRO ET AL.

Vadose Zone Journal

2.2

F I G U R E 1 Four cylindrical undisturbed soil samples extracted
from a soil block collected at 0–20 cm of a Brazilian oxisol (Rhodic
Hapludox) with different diameters: 4 cm (BS4), 2 cm (BS3), 1 cm
(BS2), and 0.5 cm (BS1).

F I G U R E 2 Observed (data) and fitted (solid line) water retention
curves of the fine-textured soil.

were wrapped with parafilm to avoid losing soil particles during data acquisition, fixed onto the microtomography sample
stage, and scanned with spatial resolutions (pixel sizes) of
3.6, 6, 12, and 25 μm for samples BS1, BS2, BS3, and BS4,
respectively. Figure 1 shows pictures of the four samples.
Figure 2 shows the measured water retention curve of thesoil, along with the fitted van Genuchten–Mualem function(van Genuchten, 1980). Fitted values of the saturated (θs ) and
residual (θr ) water contents and the semi-empirical shape factors α and n were as follows: θs = 0.497 cm3 /cm3 , θr = 0.284
cm3 /cm3 , α = 0.028 1/cm, and n = 1.750. Unfortunately, no
direct data of the saturated hydraulic conductivity (Ks ) wereavailable of the soil block we used. We used the Rosetta pedotransfer function software of Schaap et al. (2001) to obtain an
approximate estimate of Ks from the particle size distribution
and the bulk density. A 24.7 cm/day value was obtained, consistent with the average of several Ks measurements at field
locations near the sampling site (Vaz et al., 2014).

Multi-scale µCT measurements

Microtomography images were acquired with a cone beam
benchtop X-ray microCT scanner from SkyScan, model 1172
(Bruker). The X-ray tube operated at 100 kV and 100 mA. The
soil samples were rotated stepwise by 0.3˚ over 360˚, while 15
frames were acquired and averaged to represent the shadow
projection of each rotational step. The samples were wrappedwith parafilm to avoid losing soil particles during data acquisition and fixed onto the μCT scanner stage with double-sided
tape and modeling clay to avoid any movement of the samples
during rotation (for more detail, see Cássaro et al. [2017]).
Each soil sample was scanned using four distinct resolutions to permit not only studies of the internal network of
the pores but also to explore the effect of imaging resolution
on the distribution of pore sizes and the equivalent pore
diameters. The scans showed that sample BS1 had a very fine
network, samples BS2 and BS3 had fine to intermediate networks, and sample BS4 had a coarse network. Table 1 shows
the resolution, the porosity obtained from the μCT images, theconnected porosity given by the skeleton, and the hydraulicconductivity obtained for each sample after modeling with
the PoreStudio software using pore network modeling. Ourmodeling results showed that the connected porosity and the
expected hydraulic conductivity obtained using Rosetta could
not be reached using only one resolution of the μCT scans.
Cássaro et al. (2017) noted that the soil showed much variability in terms of its structure because of the presence oflarge, relatively dense aggregates of about 5 mm in diameter, possibly because of the intra-aggregate pores in this heavy
clay soil. The soil total porosity could not be obtained directlyfrom the μCT images since part of the fine-textured (clayey)
medium contained many pores with diameters smaller thanthe best resolution of the scans (3.6 μm for sample BS1). The
loss of part of the porosity (smaller pores or bigger pores) can
impact directly on the modeling of the fluid flow and transport inside the sample, once the pores that are not part of the
microCT pore network can break relevant connections among
percolation clusters at the skeleton of the sample.
The μCT stacks were reconstructed and analyzed with theAvizo 2022.2 software (Thermo Fisher Scientific, 2022). We
developed scripts in Mathematica 12.3 (2021) to track the
diameters and 3D coordinate pore paths, and obtain input datafor pore network modeling. The images of μCT, the porous
media, and the skeleton given by Avizo 2022.2 were obtainedusing 501 slices/stack for the BS1/2/3 samples and 301 slices
for the BS4 stack.

2.3

Pore network modeling

Pore network modeling was carried out using the new
PoreStudio software (Lima et al., 2022; PoreStudio, 2023),

PONTEDEIRO ET AL.

TA B L E 1

SampleVadose Zone Journal

5 of 13

MicroCT porosity and saturated hydraulic conductivity values of the four samples.
µCT resolution(µm)µCT porosity(cm3 /cm3 )Hydraulic conductivity, Ks(cm/day)

Connected porosity(cm3 /cm3 )

BS1

3.6

0.240

0.110

36.1

BS2

6

0.207

0.127

22.41

BS3

12

0.121

0.071

17.34

BS4

25

0.043

0.014

10.98

which comprises a cross-platform pore network code andassociated solver. The software allows simulations of a fullrange of flow and transport processes based on several modules. The different modules permit one to model single andmultiphase flow, tracer transport in single and multiphase systems, and reactive contaminant transport in a single phase.
The software can use the skeleton data extracted from μCT
images or generate a statistically equivalent network using
geometrical parameters (e.g., pore size and pore throat radii)
and morphological parameters (e.g., pore connectivity) basedon analysis of the sample skeleton. Here, we give only a very
short description of the PNM approach used in the PoreStudio software. We refer to Raoof and Hassanizadeh (2010) and
PoreStudio (2023) for more details.
The PNM uses pore bodies (nodes) and pore throats with a
specific topological configuration, with the pore throats prescribing the connections between pore bodies. The larger porebodies are modeled as spheres, and the pore throats are modeled as capillary tubes to serve as connections between thepore bodies. During quasi-static single-phase flow, a pressure gradient is established in the pore network between the
input and output boundaries of the PNM along the flow
direction. Numerical simulations were carried out assuming
connectivity between the top and bottom of the sample (thez-axis), being the direction of flow within the plugs. Volumetric flow through the pore throats is assumed to be laminarand described by the Hagen–Poiseuille equation:
𝑄𝑖𝑗 = 𝑔𝑖𝑗 (𝑝𝑗 − 𝑝𝑖 )(1)

where Qij is the volumetric flow rate through the pore throat
between two adjacent porbodies (i and j), pi and pj are thepressures at the two pore bodies, and gij is the conductanceof the cylindrical shaped pore throat whose value can beobtained using:
𝑔𝑖𝑗 =

𝜋𝑅4𝑖𝑗
8𝜇𝑙𝑖𝑗

(2)

in which Rij is the radius of the pore throat, lij is the lengthof the throat, and μ is the dynamic viscosity of the fluid(water). Considering a well-defined volume sample, the average velocity of fluid in the pore, v¯, is then (Raoof et al.,2013):𝑣̄ =𝑄tot 𝐿
𝑉𝑓(3)

where 𝑄tot is the total flow rate through the pore network,
which can be determined at the inlet or outlet of the pore
network as the sum of all flows, L is the length of the pore
network, and Vf is the total volume of the liquid phase within
the pore network. The absolute permeability Ks of the sample
is then calculated using Darcy’s equation. A small differential of pressure was adopted to guarantee laminar flow (lowReynolds number). The initial moisture condition is assumed
to be full saturation.
In the present study, the software used the skeleton
extracted directly from the μCT images. The program initially gives the number of pore throats and pore bodies. Some
results, such as porosity, permeability, and the average coordination number, are related to the connected pores (clusters).
PNM modeling details are provided in Table 2, while Figure 3
shows the μCT images and the sample porous media at four
different resolutions.
A digital pore analysis (DPA) was further performed to analyze the distribution of several parameters inside the samples.
The PNM approach assumes that the pore space is a network (the skeleton) with a simplified geometry represented byspherical pore bodies and connecting cylindrical pore throats.
Despite being a simplified structure, the skeleton typicallycomprises thousands of points, requiring statistical techniques
for quantitative assessments. A script was developed for the
DPA (Silveira et al., 2022) to make it possible to evaluate
the distribution of total porosity and of the pore-throat andpore-body porosities along the samples. The sample for thispurpose was divided longitudinally into 20 sections. Figure 4
shows the observed skeletons and distributions of the threeporosities along the 20 sections. The data indicate that thepore throat and pore body porosity distributions along samples BS1 and BS3 (3.6 and 12 μm) are very similar, whilesamples BS2 and BS4 show a more heterogeneous porosity
distribution with higher pore throat porosities at sections 9(BS2) and 10 (BS4), respectively.
The frequencies of the coordination number are shown in
Figure 5 for all samples. Sample BS2 (6 μm) showed the best6 of 13

Vadose Zone Journal

PONTEDEIRO ET AL.

F I G U R E 3 Images of the BS soil samples
at different resolutions: (a) a micro slice, (b) the
binarized image, (c) the skeleton image
generated with Avizo, and (d) the pressure
distribution inside the sample as calculated with
the PoreStudio pore network model (PNM).

FIGURE 4

Longitudinal distributions of the different total, pore throat and pore body porosities (Phi) along the four soil samples.

PONTEDEIRO ET AL.

TA B L E 2

Vadose Zone Journal

Pore network model (PNM) parameters of the four subsamples.

SamplePore throats

Pore bodies

Domain dimension(voxels)

Average coordination
numberBS1

118,686

79,5441300 × 1300 × 300

2.98

BS2

341,586

230,934

1600 × 1600 × 300

2.98

BS3

179,735

172,143

1700 × 1700 × 300

2.88

BS4

65,84460,4761550 × 1550 × 150

1.76

FIGURE 5

Coordination number versus frequency of the soil samples.

pore connectivity in the samples, followed by sample BS3(25 μm).
We also calculated the distribution of pore throat andpore body volumes per radius in terms of voxels. Results in
Figure 6 show that sample BS2 has the highest volumes forboth the pore throats and pore bodies, followed by sample
BS3. The DPA data explained that the samples with the best
connectivity and higher volumes of pore throats and pore bodies present higher hydraulic conductivities when modelingusing the PNM.

3

7 of 13

MULTI-RESOLUTION ANALYSISIn order to improve the modeling of a sample having different
orders of magnitude in the size of their pores, we combined
information of the samples with distinct spatial scales andμCT image resolutions into one overall global digital sample.
For this, we used stochastic networks of arbitrary size withthe same fluid low properties as the initial or parent networks.
This was done in four steps. First, we characterized the parent networks in terms of distributions of the geometrical poreproperties and correlations between these properties, as wellas the connectivity function describing the detailed network
topology. Second, to create a stochastic network of arbitrary
size, we generated the required number of nodes and bonds
with the correlated properties of the original network. The
nodes were randomly located in the given network domain
and connected by bonds according to the strongest correlation between node and bond properties, while honoring theconnectivity function. Then, using the generated stochastic
networks, we used PoreStudio to integrate two networks created in the same domain into a single two-scale network bycharacterizing the cross-scale connection structure between
the two networks. Finally, multiscale networks were generated by repeating the above scheme. The overlap/combinationof the networks ensured a better representation of the overall
pore network.
To apply the described method, a REV is required. We
used two samples of our set (BS2 and BS4) to verify if they
represented reasonable REVs. This choice was based on thefact that BS2 gave a hydraulic conductivity closer to the one
estimated, while BS4 had the larger network to which thoseof the finer resolutions could be added. The combination of

8 of 13

PONTEDEIRO ET AL.

Vadose Zone Journal

FIGURE 6

Pore body (top) and pore throat (bottom) volumes versus radius for the four samples.

more than one network is based on the original data fromthe skeletons (pore bodies and throats), as described using
the stochastic approach from PoreStudio. The PNM extracted
from the higher resolution μCT volume (the finer network) ishence combined with a coarser network.
A specific objective of our study was to apply an REV analysis of the hydraulic conductivity of the largest soil sample(BS4) imaged at the lower resolution. We wanted to know
which one of the smaller digital BS1, BS2, or BS3 could bestenhance the pore network of the larger sample. The REV is
used here to define the fundamental size for the hydraulic conductivity measurements, those of the simulations, and theirrespective averages (i.e., xmean ). We further used the coefficient of variation (Cv ) as a dimensionless measure of theheterogeneity of a sample. As reported by Vik et al. (2013),

Cv is defined as the ratio between the standard deviation of aproperty and the property arithmetic mean:√[]
Var (𝑥)
1
1+
𝐶𝑣 =
𝑥mean
4 (𝑛 − 1)(4)

in which the second part of the equation is a correction when
the number of samples (n) is equal to or <10. The hydraulicconductivity class is based on the Cv value, and the sample isconsidered homogeneous when Cv is <0.5. In our study, theCv was evaluated to be 0.28.
The pore network flow simulations were carried out for
all sub-samples. For the REV study, we used modules in
the PoreStudio software to digitally divide the coarser pore
network of sample BS4 into 10 pieces and then estimated

PONTEDEIRO ET AL.

Vadose Zone Journal

9 of 13(a)(b)

F I G U R E 7 Evaluation of hydraulic conductivity by volume: (a)BS2 (6 μm) and (b) BS4 (25 μm).

the porosity and saturated hydraulic conductivity (Ks ) of
an increasingly larger set of combined subsamples until theporosity or conductivity stabilized. For the porosity, the REV
could be easily determined. However, the results in Figure 7
indicate that sample BS4 could not reach a constant plateau,
except for the last two points. This means that the REV for Ks
could not be defined well due to soil heterogeneity, and hence
that the REV size for the proposed methodology should be
about 7000 mm3 or larger (Zhang et al., 2022).
Following Jiang et al. (2012) and others, the geometry andtopology of the individual networks were described statistically in terms of probability distributions of the volume,
radius, and connectivity of the network elements. The next
step was to generate stochastic digital samples based on theoriginal skeletons (Raoof & Hassanizadeh, 2010; Raoof et al.,2013). For this purpose, we used PoreStudio to generate astochastic sample of BS4, as well as of the other samples,
based on their original average connectivity and the original
sets of pore body and pore throat radii. The placement of porebodies and throats in the network was made stochastically by
the program.
The best distributions for the networks were found to belognormal. Figure 8 shows the histogram of the BS4 throat
length and the fitted lognormal distribution.
Using the pore network generation module of PoreStudio(2023), the stochastic BS4 sample was generated assuming
the same average connectivity of the original and stochastic

FIGURE 8

Histogram of the pore throat radius of sample BS4.

skeletons. This is done by statistically placing the pore bodies and pore throats in the new skeleton. We verified results
obtained with the stochastic sample of BS4 with the original
digital sample. Figure 9a compares pore throat volumes (byradius) obtained using PoreStudio with those obtained withthe software for the stochastic sample (along with the sample correlations), while results for the pore body volumes areshown in Figure 9b. The pore connectivity had a variation of
6% comparing the stochastic with the original digital sample. The data show good agreement, including the generated
curves. Table 3 provides an overview of all BS4 pore networkmodeling results.
To capture the multiscale PSD of the rock, we imaged four
different soil samples at different resolutions and integrated
the data to produce a pore-network model that combines information at several length scales that cannot be recovered froma single tomographic image. We used the network integration
approach of Jiang et al. (2007, 2013) to generate PNMs ofarbitrarily large volumes that can incorporate multiscale poresystems extracted from images obtained at different resolutions. The larger volume sample in our study was BS4, imaged
with a resolution of 25 μm, in which we then incorporated
features extracted from the other systems (3.6, 6, and 12 μm
scans). Sample BS4 showed less pores and throats, comparedwith the other samples, and also had a lower hydraulic conductivity. The conductivity, connectivity, and total porosity
all showed increasing trends with the imaging resolution.
Sample BS1 had the finest network and displayed a calculated hydraulic conductivity value higher than the values
estimated from the particle size distribution and approximate
field estimates.
We explored three different combinations of samples togenerate the scaled stochastic samples, that is, by combining BS4 and BS1, BS4, and BS2, and BS4 + BS3. In allcases, we used BS4 since this sample had the largest (andcoarsest) network with the largest pores and combined its
network with those of the smaller sample with higher resolution to better capture the smaller pore bodies and pore

10 of 13

PONTEDEIRO ET AL.

Vadose Zone Journal

FIGURE 9
TA B L E 3

Correlation of the digital and stochastic pore throat (a) and pore body (b) volumes of sample BS4.
Pore network model (PNM) parameters and results.

SamplesThroats

Pore bodies

Volume (mm3 )Porosity

Conductivity(cm/day)

BS4 original

65,84460,4767711.3

0.135

10.96

BS4 stochastic

66,73661,4477732.0

0.132

10.80

throats. Combining the pore systems of two networks facilitated having enough overlap to include the representation ofintermediate-size pores.
PNM results obtained with the combinations of 25 +3.6 μm, 25 + 6 μm, and 25 + 12 μm are reported in Table 4.
The best results were obtained by combining the μCT images
with resolutions of 6 μm (BS2) and 25 μm (BS4). The PNM
modeled hydraulic conductivity of sample BS2 imaged at
a resolution of 6 μm was closest to the estimated value of
24.7 cm/day. The pore system of sample BS4 received thefiner part of the BS2 network and hence also had a higher
Ks value. The combination of BS4 and BS2 pore throats andpore bodies showed an increase relative to the BS4 skeleton
but remained less than Ks of the BS2 skeleton. This meansthat the two samples had a good area of overlap in their poresystems.
The combination of BS4 and BS1 skeletons gave thesecond-best estimate of the hydraulic conductivity, but with
the number of pore bodies and pore throats being very close
to the sum of the amounts of both skeletons. This means

PONTEDEIRO ET AL.

Vadose Zone Journal

11 of 13F I G U R E 1 0 PoreStudio generated water retention data and fitted curve of the fine-textured soil (left) and comparison of the optimized
PoreStudio generated retention curve (black line) with the PoreStudio curve adjusted with the residual water content (blue line) and thelaboratory-measured retention curve (green line) (right).

T A B L E 4 Pore networking modeling results using the combined
skeletons of two samples.

Pore bodies

Conductivity(cm/day)

SamplePore throats

BS4 + BS1

158,021

128,549

20.1

BS4 + BS2

277,953

204,449

22.3

BS4 + BS3

120,953

100,729

18.0

that the pore systems of the two samples had a very narrow overlap range. Finally, the combination of BS4 and BS3
gave the smallest improvement in the BS4 conductivity; theskeletons of the generated pore network (pore bodies andpore throats) showed a much broader overlap in their poresystems.
Much of the above analysis focused on the saturated
hydraulic conductivity, Ks , which is determined mostly by
the larger connected pores. We also still ran PoreStudio forthe various networks to obtain predictions of the soil water
retention curve. Figure 10 shows the results for the pore network of sample BS1. The PoreStudio calculations were basedon a pore network made up of mostly connected pore bodies
and pore throats and do immediately account for the presenceof a residual water content (θr ), as is assumed in most standard formulations of the water retention curve. The residual
water content of especially fine-textured soils can be quite
large, as is the case for the selected clay soil in our study.
One could correct for this reason the PoreStudio results withthe residual water content of 0.284 cm3 /cm3 (Figure 1). This
would produce the middle curve in Figure 10, which is still
about 0.08 units below the measured curve near saturationand has a higher value of the van Genuchten G n parameter(4.66 vs. 1.75 for the laboratory-measured curve). These differences point to a need to measure water contents in the very
dry range, including portions of the curve where corner, film,
and vapor flow become important (e.g., Peters et al., 2015).

This issue appears is important for fine-textured (clay) soils,
such as those used in this study, since their nano-size pores are
below the resolution of most or all μCT methods (e.g., Godoyet al., 2024).

4

CONCLUSIONSIn this study, we showed one way of combining the μCT
results obtained for samples having different sizes and imaged
with different resolutions, with as ultimate objective to obtain
more representative conductivity/permeability estimate of thelarger sample imaged with a lower resolution. The idea is toobtain simpler pore networks, based on μCT stacks of the soil
or rock samples, that carry information from distinct samples imaged at higher resolutions and is implemented in pore
networks obtained from lower resolution imaging. The technique provides samples with a less-refined pore system that if
scanned at very high resolution, thus permitting one to perform other modeling studies such as contaminant transportand multiphase flow in less time and with less computational
demands. The procedure was illustrated here by combining
the images of two skeletons obtained with different resolutions. The approach can be extended to more than two samples
by using the same steps as described in this paper.
AU T H O R C O N T R I B U T I O N S
Elizabeth M. Pontedeiro: Conceptualization; methodology; writing—original draft. Martinus Th. van Genuchten:
Supervision; writing—review and editing. William Godoy:
Formal analysis; software. Mateus G. Ramirez: Conceptualization; formal analysis. Carlos M. P. Vaz: Data curation;
writing—review and editing. Silvio Crestana: Supervision;
writing—review and editing. Maira C. O. Lima: Writing—review and editing. Paulo Couto: Supervision. Jian Su:
Supervision; writing—review and editing.

12 of 13

Vadose Zone Journal

C O N F L I C T O F I N T E R E ST STAT E M E N T
The authors declare no conflicts of interest.
ORCID
Elizabeth M. Pontedeiro https://orcid.org/0000-00031418-3314
Martinus Th. van Genuchten https://orcid.org/0000-00031654-8858
William Godoy https://orcid.org/0000-0003-1879-9386
Mateus G. Ramirez https://orcid.org/0000-0002-22683274
Carlos M. P. Vaz https://orcid.org/0000-0003-1336-8958
Paulo Couto https://orcid.org/0000-0002-3847-8638
REFERENCES
Allaire, S. E., Roulier, S., & Cessna, A. J. (2009). Quantifying preferential flow in soils: A re-view of different techniques. Journal of
Hydrology, 378, 179–204. https://doi.org/10.1016/j.jhydrol.2009.08.
013
Blunt, M. J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi,
P., Paluszny, A., & Pentland, C. (2013). Pore-scale imaging and modelling. Advances in Water Resources, 51, 197–216. https://doi.org/10.
1016/j.advwatres.2012.03.003
Cássaro, F. A. M., Durand, A. N. P., Gimenez, D., & Vaz, C. M. P.
(2017). Pore size distributions of soils derived using a geometrical approach and multiple resolution microCT images. Soil Science
Society of America Journal, 81, 468–476. https://doi.org/10.2136/
sssaj2016.09.0291
Cnudde, V., & Boone, M. N. (2013). High-resolution X-ray computed
tomography in geosciences: A review of the current technology andapplications. Earth-Sciences Review, 123, 1–17. http://doi.org/10.
1016/j.earscirev.2013.04.003
Crestana, S., & Cruvinel, P. E. (2023). Developing spectroscopic
and imaging techniques for advanced studies in soil physics basedon results obtained at Embrapa instrumentation. Brazilian Journalof Physics, 52(6), Article 200. https://doi.org/10.1007/s13538-02201202-8
Dragonfly. (2022). Dragonfly 2022.2 [Computer software]. Comet
Technologies Canada Inc. https://www.theobjects.com/dragonfly
Elrahmani, A., Al-Raoush, R. I., Abugazia, H., & Seers, T. (2022). Porescale simulation of fine particles migration in porous media using
coupled CFD-DEM. Powder Technology, 398, 117130. https://doi.
org/10.1016/j.powtec.2022.117130
Godoy, W., Pontedeiro, E. M., Barros, R. A. B. R., de Vries, E. T.,
Raoof, A., van Genuchten, M. Th., & Couto, P. (2024). Modeling subresolution porosity of a heterogeneous carbonate rock sample. Vadose
Zone Journal. https://doi.org/10.1002/vzj2.20348
Godoy, W., Pontedeiro, E. M., Hoerlle, F., Raoof, A., van Genuchten, M.
Th., Santiago, J., & Couto, P. (2019). Computational and experimental
pore-scale studies of a carbonate rock sample. Journal of Hydrology
and Hydromechanics, 67(4), 372–383. https://doi.org/10.2478/johh2019-0009
Hebert, V., Garing, C., Luquot, L., Pezard, P. A., & Gouze, P. (2014).
Multi-scale X-ray tomography analysis of carbonate porosity. In S.
M. Agar & S. Geiger (Eds.), Fundamental controls on fluid flow incarbonates: Current workflows to emerging technologies (Vol. 406,
61–79). Geological Society Special Publication.

PONTEDEIRO ET AL.

Jiang, Z., van Dijke, M., Sorbie, K., & Couples, G. (2013). Representation of multiscale heterogeneity via multiscale pore networks.
Water Resources Research, 49, 5437–5449. https://doi.org/10.1002/
wrcr.20304
Jiang, Z., van Dijke, M., Wu, K., Couples, G., Sorbie, K., & Ma, J. (2012).
Stochastic pore network generation from 3D rock images. Transport
in Porous Media, 94, 571–593. https://doi.org/10.1007/s11242-0119792-z
Jiang, Z., Wu, K., Couples, G., van Dijke, M., Sorbie, K., & Ma,
J. (2007). Efficient extraction of networks from three-dimensional
porous media. Water Resources Research, 43, W12S03. https://doi.
org/10.1029/2006WR005780
Kastner, J., Zaunschirm, S., Rendl, S., & da Silva, E. P. (2021). Influence of rheocasting and rare earth metals on the 3D-microstructureof ZK60 Mg-alloy measured by multiscale X-ray computed tomography. Materials Characterization, 174, 110996. https://doi.org/10.
1016/j.matchar.2021.110996
Katuwal, S., Hermansen, C., Knadel, M., Moldrup, P., Greve, M. H.,& de Jonge, L. W. (2018). Combining X-ray computed tomography
and visible near-infrared spectroscopy for prediction of soil structural properties. Vadose Zone Journal, 17, 160054. https://doi.org/10.
2136/vzj2016.06.0054
Lima, M. C. O., Pontedeiro, E. M., Ramirez, M. G., Favoreto, J., Santos,
H. N., van Genuchten, M. Th., Borghi, L., Couto, P., & Raoof, A.
(2022). Impacts of mineralogy on petrophysical properties. Transport
in Porous Media, 145, 103–125. https://doi.org/10.1007/s11242-02201829-w
Lima, M. C. O., Pontedeiro, E. M. B. D., Ramirez, M. G., Boyd, A., vanGenuchten, M. Th., Borghi, L., Couto, P., & Raoof, A. (2020). Petrophysical correlations for permeability of coquinas (carbonate rocks).
Transport in Porous Media, 135, 287–308. https://doi.org/10.1007/
s11242-020-01474-1
Mehmani, A., Kelly, S., & Torres-Verdín, C. (2020). Leveraging digitalrock physics workflows in unconventional petrophysics: A review ofopportunities, challenges, and benchmarking. Journal of Petroleum
Science and Engineering, 190, 107083. https://doi.org/10.1016/j.
petrol.2020.107083
Menke, H. P., Andrew, M. G., Blunt, M. J., & Bijeljic, B. (2016).
Reservoir condition imaging of reactive transport in heterogeneous
carbonates using fast synchrotron tomography; Effect of initial pore
structure and flow conditions. Chemical Geology, 428, 15–26. https://doi.org/10.1016/j.chemgeo.2016.02.030
Pak, T., Butler, I. B., Geiger, S., van Dijke, M. I. J., Jiang, Z., & Surmas,
R. (2016). Multiscale pore network representation of heterogeneous
carbonate rocks. Water Resources Research, 52, 5433–5441. https://doi.org/10.1002/2016WR018719
Peters, A., Iden, S. C., & Durner, W. (2015). Revisiting the simplified
evaporation method: Identification of hydraulic functions considering
vapor, film and corner flow. Journal of Hydrology, 527, 531–542.
https://doi.org/10.1016/j.jhydrol.2015.05.020
Piccoli, I., Schjønning, P., Lamandé, M., Zaninid, F., & Moraria, F.
(2019). Coupling gas transport measurements and X-ray tomographyscans for multiscale analysis in silty soils. Geoderma, 338, 576–584.
https://doi.org/10.1016/j.geoderma.2018.09.029
PoreStudio. (2023). PoreStudio. https://porestudio.com
Raoof, A., & Hassanizadeh, S. M. (2010). A new method forgenerating pore-network models of porous media. Transport in
Porous Media, 81(3), 391–407. https://doi.org/10.1007/s11242-0099412-3

PONTEDEIRO ET AL.

Raoof, A., Nick, H. M., Hassanizadeh, S. M., & Spiers, C. J. (2013).
PoreFlow: A complex pore-network model for simulation of reactive transport in variably saturated porous media. Computers &Geosciences, 61, 160–174. https://doi.org/10.1016/j.cageo.2013.08.
005
Ruspini, L. C., Øren, P. E., Berg, S., Masalmeh, S., Bultreys, T.,
Taberner, C., Sorop, T., Marcelis, F., Appel, M., Freeman, J., &Wilson, O. B. (2021). Multiscale digital rock analysis for complex
rocks. Transport in Porous Media, 139, 301–325. https://doi.org/10.
1007/s11242-021-01667-2
Sadeghnejad, S., & Gostick, J. (2020). Multiscale reconstructionof vuggy carbonates by pore-network modeling and image-based
technique. SPE Journal, 25(01), 253–267. https://doi.org/10.2118/
198902-PA
Salek, M. F., & Beckingham, L. E. (2023). Effect of X-ray computed
tomography imaging parameters on quantification of petrophysicalproperties. Earth and Space Science, 10, e2023EA003095. https://doi.org/10.1029/2023EA003095
Schaap, M. G., Leij, F. J., & van Genuchten, M. Th. (2001). Rosetta:
A computer program for estimating soil hydraulic parameters withhierarchical pedotransfer functions. Journal of Hydrology, 251(3-4),
163–176. https://doi.org/10.1016/S0022-1694(01)00466-8
Silveira, T. M., Hoerlle, F., Rocha, A. S., Lima, M. C. O., Ramirez, M. G.,
Pontedeiro, E. M., Genuchten, M., Cruz, D. O. A., & Couto, P. (2022).
Effects of carbonated water injection on the pore system of a carbonate
rock (coquina). Journal of Hydrology and Hydromechanics, 70(2), 1–
12. https://doi.org/10.2478/johh-2022-0001
Thermo Fisher Scientific. (2022). Avizo (Version 2022.2) [Computer
software]. https://www.fei.com/software/amira-avizo/
van Genuchten, M. Th. (1980). A closed form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Science Societyof America Journal, 44, 892–898. https://doi.org/10.2136/sssaj1980.
03615995004400050002x
Van Offenwert, S., Cnudde, V., & Bultreys, T. (2019). Pore-scalevisualization and quantification of transient solute transport using
fast microcomputed tomography. Water Resources Research, 55(11),9279–9291. https://doi.org/10.1029/2019WR025880
Vaz, C. M. P., de Maria, I. C., Lasso, P. O., & Tuller, M. (2011). Evaluation of an advanced benchtop micro-computed tomography systemfor quantifying porosities and pore-size distributions of two Brazilian oxisols. Soil Science Society of America Journal, 75(3), 832–841.
https://doi.org/10.2136/sssaj2010.0245
Vaz, C. M. P., Tuller, M., Lasso, P. R. O., & Crestana, S. (2014).
New perspectives for the application of high-resolution benchtop
X-ray MicroCT for quantifying void, solid and liquid phases insoils. In W. G. Teixeira, M. B. Ceddia, M. V. Ottoni, & G. K.
Donnagema (Eds.), Application of soil physics in environmental analyses: Measuring, modelling and data integration, progress in soil
science (pp. 261–281). Springer. https://doi.org/10.1007/978-3-31906013-2_12

Vadose Zone Journal

13 of 13

Verri, I., Torre, A. D., Montenegro, G., Onorati, A., Duca, S., Mora,C. A., Radaelli, F., & Trombin, G. (2017). Development of a digital rock physics workflow for the analysis of sandstones and tightrocks. Journal of Petroleum Science and Engineering, 156, 790–800.
https://doi.org/10.1016/j.petrol.2017.06.053
Vik, B., Bastesen, E., & Skauge, A. (2013). Evaluation of representative elementary volume for a vuggy carbonate rock, Part 1:
Porosity, permeability and dispersivity. Journal of Petroleum Scienceand Engineering, 112, 36–47. https://doi.org/10.1016/j.petrol.2013.
03.029Wildenschild, D., & Sheppard, A. P. (2013). X-ray imaging and analysis
techniques for quantifying pore-scale structure and processes in subsurface porous medium systems. Advances in Water Resources, 51,
217–246. https://doi.org/10.1016/j.advwatres.2012.07.018
Wu, K., Jiang, Z., Ma, J., Couples, G. D., van Dijke, M. I. J., & Sorbie,
K. S. (2011). Multiscale pore system reconstruction and integration.
International Symposium Society of Core Analysts.
Xiong, Q., Baychev, T. J., & Jivkov, A. P. (2016). Review of pore
network modelling of porous media: Experimental characterisations,network constructions and applications to reactive transport. Journalof Contaminant Hydrology, 192, 101–117. http://dx.doi.org/10.1016/
j.jconhyd.2016.07.002
Yang, Z., Peng, X.-F., Lee, D.-J., & Chen, M.-Y. (2009). An image-based
method for obtaining pore-size distribution of porous media. Environmental Science & Technology, 43(9), 3248–3253. https://doi.org/10.
1021/es900097e
Yao, J., Liu, L., Yang, Y., Sun, H., & Zhang, L. (2023). Characterizing
multi-scale shale pore structure based on multi-experimental imaging
and machine learning. Natural Gas Industry, 10, 361–371. https://doi.
org/10.1016/j.ngib.2023.07.005
Zhang, H., Abderrahmane, H. A., Arif, M., Kobaisi, M., & Sassi, M.
(2022). Influence of heterogeneity on carbonate permeability upscaling: A renormalization approach coupled with the pore networkmodel. Energy Fuels, 36, 3003–3015. https://doi.org/10.1021/acs.
energyfuels.1c04010
Zubeldia, E. H., Ozelim, L. C. de S. M., Cavalcante, A. L. B., & Crestana,S. (2016). Cellular automata and x-Ray microcomputed tomography
images for generating artificial porous media. International Journalof Geomechanics, 16(2), 04015057. https://doi.org/10.1061/(ASCE)
GM.1943-5622.0000527

How to cite this article: Pontedeiro, E. M., vanGenuchten, M. Th., Godoy, W., Ramirez, M. G., Vaz,
C. M. P., Crestana, S., Lima, M. C. O., Couto, P., &Su, J. (2024). Multiscale pore-network reconstructionof a fine-textured heterogeneous soil. Vadose Zone
Journal, e20354. https://doi.org/10.1002/vzj2.20354