# Pore network models to determine flow statistics and structural controls in
variably saturated porous media

Ilan Ben-Noah  Juan J. Hidalgo  Marco Dentz

###### Abstract

We study the abilities of pore network models of different complexity to
determine the flow statistics and their structural controls in partially
saturated porous media. The medium permeability and hydraulic tortuosity are
the basic parameters for upscaling flow problems from the pore to the Darcy
scale. They represent average flow properties. However, upscaling and
predicting dispersion and anomalous solute transport from the pore to the
continuum scale requires knowledge of the velocity distribution, not only its
mean values. Considering four different network models of increasing
complexity, we analyze the statistical and structural properties of the fluid-
filled pore space that determines the flow statistics. We consider statistical
network models based on regular lattices with the same statistical properties
as the porous medium regarding coordination number and pore-size distribution.
We consider regular lattices which are characterized by uniform coordination,
and diluted lattices, and random lattices, which are characterized by a
distribution of coordination numbers. Furthermore, we consider a detailed
network model, which accounts for the spatial location of pores, their
coordination numbers, and the sizes of pore bodies and throats. The flow
behaviors estimated from these network models are compared to direct numerical
single-phase flow simulations in the digitized images of a fully and partially
saturated two-dimensional porous medium and different saturation degrees. We
find that the statistical network models can capture the saturation dependence
of permeability and tortuosity but are not able to reproduce velocity
statistics of even the velocity range observed in the direct flow simulations.
The detailed network models, in contrast, provide excellent estimates for all
flow statistics. This indicates that the configuration and correlation of the
fluid phase are crucial structural controls of the observed distribution of
flow velocities.

## Plain Language Summary

Conceptualizing a porous media as a network of conductors sets a compromise
between the oversimplifying conceptualization of the media as a bundle of
capillary tubes and the computationally expensive and unobtainable detailed
description of the media’s geometry needed for direct numerical simulations.
These models are abundantly being used to evaluate single and multiphase flow
characteristics. The different flow characteristics are valuable in evaluating
phenomena that may or may not be relevant for different applications. Here, we
evaluate how different information about the pore space affects the ability of
the network model to evaluate different flow characteristics.

We found that the resistance of a media to the fluid flow can be estimated by
the general stochastic features of the media (its size and connectivity).
However, to account for more complex phenomena, such as solute transport and
dispersion through the media, a piece of detailed information about the
spatial location of the fluids is needed.

††journal: Advances in Water Resources

[inst1]Institute of Environmental Assessment and Water Research (IDAEA),
Spanish National Research Council (CSIC), Barcelona, Spain. [inst2]Department
of Environmental Physics and Irrigation, Institute of Soil, Water and
Environmental Sciences, The Volcani Institute, Agricultural Research
Organization, Rishon LeZion, Israel

Flow through porous media can be conceptualized using Network models

In partially saturated media, the pore scale velocities distribution and
tortuosity are governed by the phase configuration and not by pore geometry.

Network properties to match DNS flow statistics are identified

##  1 Introduction

Network models are commonly used to evaluate and analyze quasi-static [1] and
dynamic [2], single and multiphase flow through porous media. These models are
used to predict the flow and transport through the media or to upscale pore-
scale features to the Darcy scale hydraulic properties, such as retention and
permeability. Conceptualization of the media as a pore network (rather than a
bundle of capillaries) goes back to the work of Fatt [3], where at its
simplest form, it conceptualizes the media as a lattice of conductors,
disregarding the pore shape. Another conceptualization involves an
interconnected pore network with assumed circular pores [4]. A more complex
(and realistic) modeled pore network is comprised of both bonds (relating to
the pore throats) that represent most of the resistance and sites (relating to
the pore bodies) that represent most of the capacitance [5]. In this context,
many pore network models use a dual probability density function for the pore
bodies’ and throats’ size distributions [6]. More detailed pore network
models, accounting better for pores shapes and connections, such as
equilateral triangles, squares, or stars [7, 8], do not capture the full
complexity of the actual pores but represent a compromise between the over-
simplicity of a cylindrical capillary tube network and the directly-imaged
pore shape models [9].

Reduction in the water saturation degree decreases the volume fraction in
which the flow can take place and the mean and size distribution of water-
filled pores. The pore system connectivity is also reduced, which increases
the structural tortuosity of the medium. In turn, the reduction of the
saturation degree reduces the media’s permeability and affects the flow
velocity distribution. The velocity distribution broadens because of the
formation of dead-end regions with slow velocities. These broad velocity
distributions give rise to dispersive and non-Fickian solute transport [10,
11, 12].

The pore structure also affects the velocity distribution [13, 4]. The impact
stems from the effect of the pore system geometry with wide pore bodies (low
velocity) and narrow throats (high velocity) and the Newtonian fluid’s
viscosity with lower velocities near the grains surfaces and higher velocities
in the pore throats centers. In this context, the capability of pore-networks
conceptualization to evaluate the complex velocity fields needs to be
investigated. Moreover, the importance of the different phenomena for
different saturation degrees is unknown. Here, we use different network models
based on different levels of information about the pore space and assumptions
about the pore geometry to evaluate the importance of the different physical
phenomena.

We first describe the experimental results of steady flow through variably-
saturated milifluidic devices, the image analysis method used for pore
segmentation and pore space characterization, and the numerical simulation of
the flow on the exact pore space geometry. Second, we describe the different
network models evaluated in this work (i.e., regular lattice, diluted lattice,
random lattice, and detailed network representation of the pore system) and
the different pore geometries under consideration (i.e., strip and parabolic).
Lastly, we present the ability of the different models to evaluate the media’s
permeability, the spatial distribution of pressures and flow rates, the
hydraulic tortuosity of the media, and the complex velocity distribution. The
flow characteristics are compared to the values obtained by the numerical
simulations.

##  2 Porous media and pore-scale flow

In this section, we describe the porous media and pore-scale flow
characteristics used in this work and to test the performance of the different
network models given in Section 3.

###  2.1 Pore space characterization

Porous media geometry is obtained from the binary images, i.e., the
representation of the pores and the matrix as zeros and ones, of a milifluidic
device at different saturation degrees (Figure 1) from the experimental work
of Jiménez-Martínez et al. [12]. The experiments established a steady flow of
a water and glycerol mixture in variably saturated (i.e., containing also a
stagnant air phase) hydrophilic millifluidic device. Four effective saturation
degrees of $S_{e}=1.00$, $0.80$, $0.73$, and $0.65$ (Figure 1a,b,c, and d in
accordance) are considered, where $S_{e}$ is defined as the fraction of the
porosity ($\phi$) filled with the continuous and percolated liquid phase. Flow
was controlled so that the phase distributions were not altered during the
flow experiments, that is, steady-state imbibing conditions were attained. The
flow rates, corresponding to a capillary number of about $3\cdot 10^{-5}$, are
considered capillary-dominated.

Figure 1: Binary images of the percolated flowing phase (in white) for a)
saturated media ($S_{e}=1$) b) $S_{e}=0.80$, c) $S_{e}=0.73$, and d)
$S_{e}=0.65$. Images were taken from the experimental work of Jiménez-Martínez
et al. [12]. The solid medium comprises cylinders (pillars) with a mean
diameter of $0.83$ mm (and standard deviation of $0.22$ mm), a height of
$b=0.5$ mm, and a porosity of $0.72$. The pixel size of the 2D images is
$\Delta x=0.032$ mm. The images do not include the entire original flow device
and depict a flow domain of length $L=105$ mm ($3292$ pixels) along the main
flowing axis ($x$-axis) and width $W=70$ mm ($2182$ pixels)

####  2.1.1 Image analysis and pore space partitioning

The image analysis used in this work follows the methodology presented in
[14]. In a nutshell, a “$\lambda$-based” medial axis method is used to segment
the pore space into distinct pores. This method uses the binary images of the
pore space (Figure 1) to form the distance map, i.e., the Euclidean distance
of each point from the closest phases’ interface (Figure 2a). In the distance
map, positive and negative values are assigned for the pore space and stagnant
phase in accordance. Then, the distance map curvatures and the eigenvalues of
the Hessian matrix of the distance map are calculated at each point (Fig. 2b).
A new binary image is then derived from the small eigenvalue ($\lambda_{1}$)
map, in which negative values are assigned with 1 (i.e., the blue colored
areas in Fig. 2b) and the rest with 0. The union of the lines formed by the
skeletonization of the new binary image and the complementary to the original
binary image gives a binary image of the throats (Fig. 2c). Thus, the union of
the complementary to the throats image and the original binary image gives a
new binary image of the segmented pores (Fig. 2d). The pore bodies’ sizes, in
terms of the maximal disk that can fit within the pore, are given by the
maximal value of the original distance map in each pore (blue circles in Fig.
2d), i.e., not accounting for the throats for the confinement of the sphere.
Similarly, the throats’ sizes are the maximal values of the distance map on
each distinct throat line, i.e., disregarding the tortuosity of the throat
line. The eigenvalue ($\lambda_{1}$) map is processed using a Gaussian filter
as discussed in the appendix of Ben-Noah et al. [14]. For the saturated case,
a single-layer filter with a standard deviation ($\sigma=5$) is used.
Conversely, a Gaussian pyramid filter is used for the partially saturated
cases with a base level standard deviation ($\sigma_{0}=0.25$), standard
deviation increments of ($\Delta\sigma=0.005$), and a scale space factor of
($\gamma=0.75$).

####  2.1.2 Saturation degree and pore space topology

The effect of the saturation degree on the pore size distributions and
connectivity was discussed in Ben-Noah et al. [14], and the effect of $S_{e}$
on the velocity distribution in Ben-Noah et al. [11], for reader’s
convenience, we summarize it here. The reduction in $S_{e}$ reduces the mean
water-filled pore size distribution because capillary forces preferentially
drive air to invade large pores. Moreover, air, being the non-wetting phase,
invades water-filled pores through their center. Thus breaking up large pores
into small water films surrounding the grains. The latter mechanism forms a
dual pore-size distribution (Figure 3a,b). The throats’ size distributions
(Fig. 3b) of this granular media are of the same order of magnitude as the
pore bodies, slightly decreasing with the reduction in $S_{e}$. The
coordination number (Z, Fig. 3c) is the number of interconnected pores. The
reduction in the saturation degree increases the number of dead-end pores
(Z=1) and of pores connected in series ($Z=2$). In turn, the mean coordination
number ($\langle Z\rangle$) reduces with the reduction in $S_{e}$, with a
maximal value of $3.14$ and minimal $2.97$ for the saturated ($S_{e}=1$) and
$S_{e}=0.65$ cases in accordance.

Figure 2: a) distance map b) small eigenvalues ($\lambda_{1}$) of the Hessian
matrix map, c) binary image of the pore throats, and d) binary image of the
segmented pore space in a zoomed-in section of the milifluidic device in the
$S_{e}=0.65$ case. Blue circles in (d) present the pore body in terms of a
maximal disc. Figure 3: Probability density functions of the a) pore bodies
radius ($r_{t}$), b) pore throats radius ($r_{t}$), and c) coordination number
distributions of the different effective saturation degrees. The mean values
are given in the figure labels. The pores bodies’ and throats’ sizes and
connectivity are determined using the image analysis specified in section
2.1.1

###  2.2 Direct numerical simulations of pore-scale flow

The flow field is characterized in terms of its velocity distribution,
stagnant region volume, permeability, and hydraulic tortuosity. The two-
dimensional water flow field for each analyzed image is determined by solving
the steady-state Stokes equation assuming an incompressible laminar flow

$\mu\nabla^{2}\mathbf{u}-12\mu\mathbf{u}/b^{2}=\nabla p$ (1)
$\nabla\cdot\mathbf{u}=0,$

where $\mathbf{u}$ [L/T] is the flow velocity, $\mu$ [ML-1T-1] the fluid’s
dynamic viscosity and $b$ [L] the device height. The term
$-12\mu\mathbf{u}/b^{2}$ in (1) accounts for the viscous forces induced by the
millifluidic device’s top and bottom plates [15].

####  2.2.1 Numerical setup

We use the COMSOL [16] multiphysics platform for the numerical simulation. A
unit pressure gradient was imposed between the left (inlet) and right (outlet)
boundaries, i.e., $p(0,y,t)=\rho g_{c}L$ and $p(L,y,t)=0$, with $\rho$ [ML-3]
the fluid density, $g_{c}$ [LT-2] the gravitational constant, and L [L] the
device length. The linearity of the flow equation (1) implies that both
pressure and velocity (normalized by their mean values) are independent of the
fluid properties (i.e., density and viscosity). A no-slip boundary condition
was set at the phases’ interfaces. Using the no-slip assumption on the liquid-
gas interface during the liquid (wetting), phase flow was evaluated and
justified in Ben-Noah et al. [11].

####  2.2.2 Flow statistics

##### Permeability

The permeability $k$ [L2] of a medium is defined as the flux of a viscosity
unit that a unit pressure (P, [ML-1T-2]) gradient will invoke [17]

$k=\frac{q_{m}\mu}{|\nabla P|}=\frac{q_{m}\mu}{\rho g_{c}}.$ (2)

The flux $q_{m}$ is given in terms of the volumetric flow rate, or discharge
rate $Q$ as $q_{m}=\frac{Q}{Wb}$ [L T-1], which can be measured at any medium
cross-section.

##### Distribution of flow speeds

We determine the probability density function (PDF) of Eulerian flow speeds
$u=|\mathbf{u}|$, denoted by $f(u)$, through areal sampling from the direct
numerical flow simulations as

$\displaystyle f_{u}(u_{k})=\frac{\sum_{i}\mathbb{I}(u_{i},u_{k})A_{i}}{\Delta
u_{k}\sum_{i}A_{i}},$ (3)

where $A_{i}$ [L2] is the area of the $i$-th mesh element, and $u_{i}$ is the
flow speed at the Gaussian point of the $i$-th mesh element. The indicator
function $\mathbb{I}$ is one if the velocity $u_{i}$ is within the interval
$(u_{k},u_{k}+\Delta u_{k}]$ and zero otherwise.

##### Hydraulic tortuosity

The hydraulic tortuosity $\chi$, which measures the ratio of averaged
trajectory length to the linear distance, accounts for the divergence of the
streamline from the macroscopic flow direction. Thus, it can be calculated by
the ratio of the mean flow speed $\overline{u}_{e}$ [LT${}^{-}1$] and mean
flow velocity (in the flow direction) $\overline{u}=q_{m}/S_{e}\phi$ [18]

$\chi=\frac{\overline{u}_{e}}{\overline{u}}.$ (4)

##### Stagnant regions

The media permeability, hydraulic tortuosity, and velocity distribution of the
different saturation degrees of the numerical simulation are considered the
“ground truth” and the basis for estimating the applicability of the different
network models.

##  3 Pore network models

We present here pore network models of different complexity to represent the
porous media discussed in the previous section. We first summarize the basic
methodology for solving flow on networks. Then, we discuss two representations
of the pore geometry and the corresponding pore-scale flow distribution and
target flow statistics. With these preparations, we define the network models
used to represent the flow properties in the model porous media presented in
the previous section.

Pore network models represent a porous medium as a network of sites and bonds,
which connect them. Sites are typically located at the centers of pore bodies,
the bonds align with the pore throats. Both are identified during the pore-
space partitioning process. These basic network features can be supplemented
with information on the pore shape, diameter, and other relevant geometric and
hydraulic characteristics. In the following, we use the terminology ”sites”
and ”bonds” when we refer to the network model and ”body” and ”throat” when we
refer to the actual porous medium.

###  3.1 Flow on networks

The flow rate $Q_{ij}$ along the bond connecting the sites $i$ and $j$ is
given by

$\displaystyle Q_{ij}=-g_{ij}(P_{j}-P_{i}),$ (5)

where $P_{i}$ and $P_{j}$ are the pressures at sites $i$ and $j$. The bond
conductance is $g_{ij}=\kappa_{ij}/\ell_{ij}$ [M-1L4M-1] with $\ell_{ij}$ [L]
the bond length and $\kappa_{ij}$ [M-1L3M-1] the bond conductivity. The bond
conductivity depends on the pore geometry as discussed in detail in Subsection
3.2. Note that $Q_{ij}=-Q_{ji}$. If $Q_{ij}>0$, flow is directed from node $i$
to $j$ and vice versa. Volume conservation requires that the sum of the flow
rates in the bonds meeting at a site $j$ is zero,

$\sum_{[ji]}Q_{ij}=0,$ (6)

where the notation $\sum_{[ji]}$ denotes summation over the sites $j$ that are
connected to site $i$. The coordination number of site $i$, that is, the
number of sites $j$ connected to site $i$ is $Z_{i}=\sum_{[ji]}$. The flow
rate $\mathcal{Q}_{i}$ at an internal site $i$ of the network is evaluated
from the sum of the flow rates $Q_{ij}$ through the bonds that are connected
to the site $i$ as

$\mathcal{Q}_{i}=\frac{1}{2}\sum_{[ji]}|Q_{ij}|.$ (7)

In words, $\mathcal{Q}_{i}$ equals half the sum of the absolute values in the
rows of the flow rate matrix $\mathbf{Q}$. Combining equations (5) and (6)
gives the following system of equations at all internal sites of the network,

$\mathbf{L}\cdot\mathbf{P}=\mathbf{0},$ (8)

where $\mathbf{P}$ is the vector of pressures at all sites of the network,
$\mathbf{L}=\mathbf{D}-\mathbf{G}$ is the weighted Laplacian matrix. It is
given by the difference of the degree matrix $\mathbf{D}$ and the conductance
matrix $\mathbf{G}$. The degree matrix is a diagonal matrix that contains the
degrees of each node, that is $D_{ik}=\delta_{ik}\sum_{[ji]}g_{ij}$. The
weight matrix contains the conductances between connected sites
$G_{ij}=g_{ij}$ and $0$ else. The pressure boundary conditions used here are
implemented by defining $\mathbf{L}^{\prime}$ such that
$L^{\prime}_{ij}=L_{ij}$ for all sites $i$ that are not, and
$L^{\prime}_{ij}=\delta_{ij}$ for all sites $i$ that are in the inlet or
outlet boundaries such that

$\mathbf{L}^{\prime}\cdot\mathbf{P}=\mathbf{b},$ (9)

where the vector $\mathbf{b}$ on the right side of Eq. (9) contains the
boundary conditions. We assign the pressure $b_{j}=p_{0}=\rho g_{c}L$ at the
sites $j$ at the inlet boundary and $b_{j}=0$ else.

###  3.2 Pore geometries and conductivities

Pore network models can also include representations of the geometry of the
pore space [4, 19, 20], in which the bodies represent most of the phase
capacitance and throats most of the resistance to the flow [5].

####  3.2.1 Rectangular pore shape

The rectangular (”strip-work”) pore shape is characterized by the pore channel
height $b$ and width $w$. Recall that from the pore-space partitioning, we
know the sizes of the pore bodies and pore throats and their distributions
(see Fig. 3). Since most of the resistance to the flow is in the throats, we
set $w=2r_{t}$ (black line in Figure 4b). For a rectangular channel, the
conductivity is given by [21]

$\kappa=\frac{b^{3}w}{12\mu}-\frac{16b^{4}}{\mu}\sum_{n=1}^{\infty}\frac{\cosh(\frac{2\beta_{n}w}{b})-1}{\beta_{n}^{5}\sinh(\frac{2\beta_{n}w}{b})},\qquad\beta_{n}=(2n-1)\pi,$
(10)

The distribution of bond conductivities is then given by

$\displaystyle
p_{\kappa}(k)=\int\limits_{0}^{\infty}drf_{t}(r)\delta[k-\kappa(r)].$ (11)

Note that the pore throat, in principle, has zero length as it marks the
boundary between two pores. Thus, assigning a finite length to the bond
underestimates the conductance of the pore. The velocity profile in a
rectangular pore is approximated by the Hagen-Poiseuille equation

$\displaystyle u(y)=\frac{3Q}{2bw}\left(1-\frac{y^{2}}{w^{2}}\right),$ (12)

where the pore flow rate is $Q=\kappa|\Delta P|/\ell$, where $\Delta P$ is the
pressure drop along the length $\ell$ of the pore.

####  3.2.2 Parabolic pore shape

The second geometry under consideration assumes a parabolic shape of the pore
(red line in Figure 4b). This parametrization of the aperture is agile enough
to give a good approximation to the actual aperture, for example, given by the
minimal distance from the pore-media interface of the medial axis (blue line
in Figure 4a and b). where the aperture ($w$) of a pore changes from
approximately the pore throat diameter ($2r_{t}$) at one end ($\mathbb{x}=0$)
to the pore body diameter ($2r_{b}$) at its other end ($\mathbb{x}=l$), where
$\mathbb{x}$ is the axis of the main flow direction taken to be perpendicular
to the pore throat facing away from the pore center. Moreover, the parabola
has a minimal value at the pore throat, which means that the derivative of
$w(x)$ at the throat location $\mathbb{x}=l$ is zero. Then, the pore aperture
is given by

$w(\mathbb{x})=2r_{t}+2(r_{b}-r_{t})\frac{\mathbb{x}^{2}}{l^{2}}$ (13)

Figure 4: a) A zoomed-in section in the $S_{e}=0.65$ case, the blue line
depicts the bond medial-axis ($\mathbb{x}$) from the throat to the pore
centers, b) the bond aperture evaluated from twice the distance map (blue)
value (i.e., twice the distance from the nearest solid-pore interface), the
strip geometry assumption (black), and the parabolic geometry (red). And bonds
conductance probability density function ($f(g)$) for c) strip (Eq. 11) and d)
parabolic pore geometries (Eq. 15).

For laminar flow, the effective pore conductance $\kappa_{e}$ is evaluated by
the harmonic sum of the local conductivities along the pore channel with
variable width $w(x)$ by using Eq. 10 for $\kappa(\mathbb{x})$

$\frac{1}{\kappa_{e}}=\frac{1}{\ell}\int\limits_{0}^{\ell}\frac{dx}{\kappa(x)}.$
(14)

Note that $\kappa_{e}=\kappa_{e}(r_{b},r_{t})$ is a function of the body and
throat radii because $\kappa(x)$ depends on $w(x)$. Assuming that $r_{t}$ and
$r_{b}$ are independent variables, that is, disregarding the correlation
between the bodies and the throats and also disregarding the correlation
between $r_{b}$ and the coordination number ($Z$), the distribution of
conductivity can be determined as

$\displaystyle
p_{\kappa}(k)=\int\limits_{0}^{\infty}dr\int\limits_{0}^{r}dr^{\prime}f_{b}(r)f_{t}(r^{\prime})\delta[k-\kappa_{e}(r,r^{\prime})].$
(15)

Note that we use the condition that $r_{t}<r_{b}$. This condition causes a
bias towards smaller throats and larger bodies. The magnitude of this bias
depends on the overlap of the pore body and throat size distributions (Figure
3). For a small overlap, the bias is small. On the other hand, a significant
overlap means that the pore bodies and throats are of the same size. In this
case, the pore geometry is essentially rectangular. In the model media, the
rectangular pore geometry (Fig. 4c) has many more small conductors for the
unsaturated cases than the parabolic geometry (Fig. 4d). The saturated cases
(blue lines) are very similar in both magnitude and distribution because the
differences between $f(r_{b})$ and $f(r_{t})$ are not very large (Fig. 3).

For a more detailed representation, the pore space is reconstructed by joining
semi-parabolas at the pore throat such that they have the same $r_{t}$, but in
general, different $r_{b}$, see Figure 5. The pore geometry also affects the
intra-pore velocity distribution with higher velocity in the narrow section of
the pore throat due to volume conservation. That is, the mean cross-sectional
velocity is inversely proportional to the $w(x)$. The no-slip condition at the
pore-grain boundary gives a parabolic velocity profile for laminar flow
conditions. Here, for simplicity, we approximate the velocity profile across
the pore by the Hagen-Poiseuille equation

$u(\mathbb{x},\mathbb{y})=\frac{3Q}{2bw(\mathbb{x})}\left(1-\frac{\mathbb{y}^{2}}{w(\mathbb{x})^{2}}\right),$
(16)

where $Q=\kappa_{e}|\Delta p|/\ell$ is the pore flow rate. Note that we set
the average cross-sectional velocity to
$\overline{u_{\mathbb{x}}}={Q_{b}}/{bw(\mathbb{x})}$. The maximum flow
velocity is $u_{m}=3\overline{u_{\mathbb{x}}}/2$ at the center of the conduit.
The resulting velocity distribution is illustrated in Figure 5.

Figure 5: The velocity distribution in a single bond in the detailed network.
The axis $\mathbb{y}=0$ is the center of the bond, which is assumed
symmetrical around this axis.

###  3.3 Flow statistics

##### Permeability

The permeability $k$ of the network is determined from the flow rates at the
inlet sites, which are given by $\mathcal{Q}_{j}=\sum_{[ji]}Q_{ji}$, where $j$
is a site at the inlet boundary. The flux $q_{m}$ is then given by
$q_{m}=\sum_{j}Q_{j}/(Wb)$. From $q_{m}$ we obtain the network permeability
according to the expression in Eq. 2. For regular and diluted networks, $k$
can be evaluated using the Effective Medium Approximation (EMA), see A.

##### Distribution of flow speeds

The distribution of flow speeds in a single bond is denoted by $p_{ij}(u)$.
For the parabolic profile (12) it is given by

$\displaystyle p_{ij}(u)=\frac{bw_{ij}}{3Q_{ij}\sqrt{1-2ubw_{ij}/Q_{ij}}}$
(17)

for $0<u<3Q_{ij}/2bw$. For the velocity profile (16) it is accordingly given
by

$\displaystyle
p_{ij}(u)=\frac{1}{\ell_{ij}}\int\limits_{0}^{\ell_{ij}}dx\frac{bw_{ij}(x)}{3Q_{ij}\sqrt{1-2ubw_{ij}(x)/3Q_{ij}}}.$
(18)

The speed distribution across the network is then obtained by the weighted sum
of the $p_{ij}$, that is,

$\displaystyle f_{u}(u)=\frac{\sum_{i>j}p_{ij}(u)A_{ij}}{\sum_{i>j}A_{ij}}.$
(19)

##### Distribution of flow rates

The distribution of flow rates $Q_{ij}$ is determined in analogy to the
distribution of flow speeds as

$\displaystyle f_{Q}(Q_{k})=\frac{\sum_{i>j}\mathbb{I}(Q_{ij},Q_{k})}{\Delta
Q_{k}}.$ (20)

##### Hydraulic tortuosity

We suggest to estimate the hydraulic tortuosity of the network as the ratio
between the average flow rate magnitude and the average flow rate in the
network, that is,

$\tau=\frac{\sum_{i>j}|Q_{ij}|}{\sum_{i>j}Q_{ij}\cos(\alpha_{ij})},$ (21)

where $\alpha_{ij}$ is the bond angle with the $x$-axis. We can define an
average bond velocity as

$\displaystyle v_{ij}=\frac{Q_{ij}}{A_{ij}},$ (22)

where $A_{ij}$ is the average surface area of the bond cross-section. If we
assume that the bond lengths are not very variable and can be replaced by the
network average, $\ell_{ij}=\langle\ell\rangle$, definition (21) becomes

$\displaystyle\tau=\frac{\sum_{i>j}|v_{ij}|V_{ij}}{\sum_{i>j}v_{ij}\cos(\alpha_{ij})V_{ij}}=\frac{\langle
v\rangle}{\langle v_{x}\rangle}=\chi,$ (23)

where we set the bond volume $V_{ij}=A_{ij}\langle\ell\rangle$.

###  3.4 Network models

####  3.4.1 Regular networks

Figure 6: Regular (hexagonal) network.

A regular network (Figure 6) refers to a lattice with a constant coordination
number ($Z$). Regular networks are the most basic networks under consideration
in this study. They require only a pore size distribution as input parameters
and an estimate of the mean coordination number to decide the lattice
geometry. The coordination number is considered the single most important
topological feature of the medium [22]. It has been shown that regular
networks and topologically disordered networks variable $Z$ have similar
transport properties if they share the same mean coordination number $\langle
Z\rangle$ [23]. The mean coordination number of the study media presented in
Section 2.1 is close to $\langle Z\rangle=3$ for all cases, see Fig. 3c.
Therefore, a regular hexagonal (honeycomb) lattice with $Z=3$ and a constant
bond length $\ell$ is used. Moreover, for each saturation degree, the number
of sites of the corresponding honeycomb network model is equal to the number
of pores obtained from the pore space segmentation in Subsection 2.1.1, while
the physical dimensions of the network are kept constant. This implies that
the bond lengths of the different network models are slightly different. The
longest diagonal of the network is aligned with the main flow direction. The
bond conductances are assigned randomly from the distribution $f(g)$ (Fig. 4).
As discussed in Subsection 3.2, we consider both rectangular and parabolic
pore geometries. This process disregards any spatial correlation in the pore
size distribution, which is in line with the pore structure of the porous
media considered in this study. For correlated networks, we refer readers to
Sahimi [24].

####  3.4.2 Diluted regular networks

Figure 7: Diluted regular network.

The diluted regular networks comprise a more detailed variation of the regular
networks (Figure 7). In these networks, a regular lattice with
$Z_{lattice}>\langle Z\rangle$ is imposed. Then bonds are eliminated until the
mean coordination number of the diluted network is equal to the mean
coordination number of the porous medium, $\langle Z_{diluted}\rangle=\langle
Z\rangle$. That is, the coordination number is variable but does not
necessarily follow the distribution of coordination numbers in the porous
medium. As for the regular networks, the bond conductances are assigned
randomly from the distribution $f(g)$ (Fig. 4) for both rectangular and
parabolic pore geometries. For each medium, the mean coordination number can
be evaluated from the number of bonds ($N_{b}$) and sites ($N_{s}$) through
the relation $\langle Z\rangle=2{N_{b}}/{N_{s}}$. It should be noted that for
irregular $3D$ media, this ratio does not necessarily hold because a throat
might connect more than two pores [14]. The regular and diluted models (square
and hexagonal) has the same number of bonds ($N_{b}$) equal to the number of
pore throats taken from the image analysis. Since we keep the original device
dimension, this means that the bond lengths are slightly different
($l_{sq}>l_{hex}$).

For the diluted regular networks, we use a regular square lattice
($Z_{grid}=4$) for all cases ($S_{e}$) and a hexagonal lattice ($Z_{grid}=3$)
only for the lowest saturation degree ($S_{e}=0.65$ for which the mean
coordination number is $\langle Z_{media}\rangle=2.97$). Comparing this case
with different lattices allows us to separately evaluate the role of the
lattice dimensionality and dilution. We start with a fully saturated grid,
i.e., with a number of bonds equal to the number of throats obtained from the
pore segmentation using the image analysis. Then, the conductors are randomly
assigned from the distribution $f(g)$. As in the regular network model, we
compare the strip and parabolic pore geometries (described in Subsection 3.2)
by accounting for the different conductances support ($g$) and distribution
($f(g)$). Then, the grid is diluted to obtain the $\langle Z\rangle$ of each
case (of different $S_{e}$), and in consequence, $N_{s}$ equals the number of
pores in the device. The dilution method applied here follows a simplification
to the protocol described by [25]. Each bond is randomly assigned with a
dilution (or elimination) number. Then, an elimination threshold value is used
so that all bonds with an elimination number larger than the threshold are
assigned zero conductivity. The elimination threshold is related to the mean
coordination number through

$f_{et}=\langle Z_{\text{media}}\rangle/Z_{\text{grid}},$ (24)

where $Z_{\text{grid}}$ is the coordination number of the regular lattice
before dilution. The conductance distribution is adjusted to account for the
dilution, i.e., adding $g=0$ with a dilution probability of
$p(g=0)=1-N_{b}/N_{s}$.

The random elimination of bonds by assigning zero conductivity $k=0$ can
create singular points or regions, that is, subregions of the network that are
not connected to the inlet and outlet. Such regions are artifacts of the
network generation procedure and do not reflect the real medium structure.
From a numerical point of view, disconnected regions may cause numerical
problems since they lead to a singular or ill-conditioned Laplacian matrix.
The challenge of addressing this singularity issue was raised in several
diluted network models [25]. Usually, this issue is solved by preconditioning
the matrix or by using a search algorithm to find and remove all ill-connected
sites and clusters. In this contribution, we use the Singular Value
Decomposition (SVD) method as an efficient way to eliminate the singular
points of the matrix [26]. Regarding the flow behavior, diluted regular
networks may display more stagnant regions than regular networks, that is,
sites of very low or zero velocity due to the occurrence, for example, of
dead-end pores.

####  3.4.3 Random network

Figure 8: Random network.

In this section, we describe the generation of a random network (Figure 8)
that has the same coordination number distribution as the porous medium (not
just the same mean). We start from a regular square lattice of sites that are
also connected along the diagonals, i.e., with a constant coordination number
of $Z_{lattice}=8$. This value is equal to the largest observed coordination
number ($Z=8$, Fig. 3c). Unlike for the previous networks, in this case, the
length is not constant because the diagonals are longer than the horizontal
and vertical bonds by a factor of $\sqrt{2}$. As for the previous networks
(sections 3.4.1 and 3.4.2), conductance values are assigned randomly to all
the bonds from the distribution $f(g)$. Adjusting the distribution of
coordination numbers of the lattice to a prescribed distribution is done
iteratively through diluting and concentrating schemes as described in the
following.

The algorithm is initialized by assigning a target coordination number
$\mathbb{Z}_{j}\leq Z_{lattice}$ to each site $j$ drawn randomly from the
distribution $f(Z)$ (Fig. 3c), and by assigning a random elimination number to
each bond of a site from a uniform distribution between $0$ and $1$. All site
coordination numbers at initialization are $Z_{j}^{(0)}\geq Z_{lattice}$. At
the next step (dilution step), if a site coordination number is larger than
the target value ($Z_{j}^{(0)}>\mathbb{Z}_{j}$), the bond with the smallest
elimination number is removed by assigning zero conductance. After the
dilution step, all coordination numbers are updated. The dilution step is
followed by a concentration step. If the grid coordination number at the site
is smaller than the target value ($Z_{j}^{(1)}<\mathbb{Z}_{j}$), then a random
bond is added (to one of the available locations, keeping the initial
conductance value) and a new random elimination number is assigned to the
bond. After the concentration step, all coordination numbers are updated, and
a dilution step follows. These steps are repeated until the distribution of
coordination numbers in the random network has converged to the target
distribution. As a convergence criterion, we use a threshold value for the
difference between the target and lattice mean coordination number. That is,
the algorithm has converged after the $k$th iteration if
$|\langle\mathbb{Z}\rangle-\langle Z^{(k)}\rangle|<10^{-4}$. For the media
under consideration, the coordination number distribution converges to the
prescribed distribution after less than 30 iteration steps.

####  3.4.4 Detailed network representation of the pore system

Figure 9: Detailed network representation of the pore system. Taken in a
Zoomed-in section of the $S_{e}=0.73$ case.

The detailed representation of the pore system accounts for the spatial
location of the pores, their coordination number, and the sizes of the pore
bodies and throats (Figure 9). In order to generate the detailed network, the
segmented pore label image (attained from the image analysis), in which each
pore is assigned a different index (label), is used to form a graph. The graph
holds information about the connectivity of each pore and its neighbors (on an
8-connectivity directions basis). The site location is set equal to the
location of the maximal distance map value in the pore, i.e., to the center of
the maximal disk (Fig. 2d). The bond length between each connected pores
couple is set as the distance between pore centers. The pore geometry in the
detailed network representation is parabolic, see Subsection 3.2.2.

##  4 Results and discussion

In this section, we study the estimates of the different network models for
permeability, hydraulic tortuosity, and the distributions of pressure and flow
rate and compare them to the results of the direct number simulations, which
represent the ground truth. These are the key quantities for the understanding
of pore-scale flow and solute dispersion in saturated and partially saturated
porous media [27, 28, 11].

###  4.1 Permeability

The ability to estimate the permeability ($k$) of a representative elementary
volume (REV) is a basic prerequisite for using network models to upscale pore-
scale flow phenomena to the Darcy and field scales. Figure 10 shows the $k$
estimated from the different network models and $k$ based on the effective
medium approximations (EMA) versus the permeabilities obtained from the direct
numerical simulations (DNS). The estimates from the pore network models are
based on averages of 25 network realizations for each pore geometry. The
regular and diluted networks give a relatively narrow span of permeabilities
for different realizations, suggesting that the domain is large enough to be
considered a REV. As expected, for all network models, the permeability
monotonically increases with the saturation degree $S_{e}$. Even if the
permeabilities from the network models and direct simulation do not coincide,
the network models may give a good estimate for the relative permeability
$k_{r}=k/k_{s}$, where $k_{s}$ is the saturated (intrinsic) permeability. The
effect of $S_{e}$ on $k_{r}$ for the different networks is presented in figure
SF1 in the supplementary material.

Figure 10: The permeability ($k$, Eq. 2) of the network compared to the
numerically evaluated permeability for the a) regular, b) diluted, c) random,
and d) detailed networks. The different values of k are for the different
saturation degrees under consideration (Fig. 1). The boxplots of the regular
(a), diluted (b), and random (c) models are evaluated from 25 realizations for
each pore geometry (strip or parabolic). The extreme values exceeding the
$95\%$ confidence level (box) are depicted with triangles.

##### Regular networks

The regular networks (and also the diluted and random) estimate better the
relative permeability ($k_{r}$, Fig. SF1) than the absolute one ($k$, Fig.
10). This is in accordance with the literature [22]. Here, the regular
networks slightly underestimate the effect (slope) of $S_{e}$ on the relative
permeability $k_{r}=k/k_{s}$ (Fig. SF1a). The increase in the standard
deviations of the permeabilities (seen in the box plots in Fig. 10a,b) with
the reduction of the saturation $S_{e}=1$ in the hexagonal grid is due to the
changes in the pore size distribution with a broader distribution and many
more small pores accounting for the narrow water films.

The $k$ evaluated by the EMA for the rectangular pore geometry (red
rectangular markers in Fig. 10a) and parabolic (blue circle markers) pore
geometries are very close. This is because of their relatively similar
conductance probability density functions (Fig. 4). As expected, the $k$
evaluated by the EMA is very close to the mean $k$ evaluated by the regular
network models (box plots, Fig. 10a). Interestingly, the strip pore geometry
(red) gives better agreement between the EMA and the network. Maybe because of
the smaller variability in the bond conductances.

##### Diluted networks

The diluted networks (Fig. 10b) provide a better estimate for $k_{r}$ as a
function of $S_{e}$ than the regular networks (Fig. SF1b), but give worse
estimates for $k$, mostly underestimating the data from the direct numerical
simulations. The strip pore geometry in the diluted network (red boxes in Fig.
10b) presents higher variability between realizations (wider box plots)
compared to the parabolic geometry. This is because of the combined effect of
the larger fraction of low conductances (Fig. 4a) and the network’s dilution,
which may form singular and dead-end sites. Interestingly, the $k$ values
obtained from the diluted EMA are lower than the ones evaluated by the network
models.

##### Random networks

The estimates by the random networks (Fig. 10c) for $k$ and kr are similar to
those of the regular and diluted networks. This suggests that the added
information held in the distribution of $Z$ as opposed to the mean
coordination number does not improve the estimation of $k$.

##### Detailed networks

The permeability estimates of the detailed network model (Fig. 10d) are in
excellent agreement with the direct numerical simulations. This confirms the
ability of the detailed network model to represent the mean flow behavior
through porous media. However, this model relies on detailed pore space
information that may often not be available.

###  4.2 Tortuosity

Figure 11 shows the hydraulic tortuosity estimates from the network models
versus the tortuosity obtained from the direct numerical flow simulations
using the digitized images of the porous media. Recall that the definition of
tortuosity used for the network models in Eq. (21) is different from the
definition in Eq. (4) for the direct numerical simulations. The two
definitions coincide only if the bond lengths are constant, which is the case
for regular networks. We observe for the networks and direct numerical
simulations that tortuosity increases monotonically as the saturation degree
decreases, which is due to the increasing complexity of the fluid phase
configuration. Generally, the network models overestimate the values obtained
from the direct flow simulations. A possible cause is the geometry of the
lattice underlying the network models and its angle with respect to the
pressure gradient. We define the relative tortuosity
$\tau_{r}={\tau}/{\tau_{s}}$, where $\tau_{s}$ is the tortuosity of the fully
saturated case $S_{e}=1$. The effect of $S_{e}$ on the tortuosity can be
evaluated from the slope of the markers or centers of the boxes compared to
the dashed line in Figure 11. It can be seen that all the models, except the
regular networks, give a reasonable estimation of $\tau_{r}$. This suggests
that the changes in the conductance distribution alone can not explain the
increase of tortuosity with the reduction of $S_{e}$.

Figure 11: The hydraulic tortuosity on the network ($\tau$, Eq. 21) compared
to the numerically evaluated hydraulic tortuosity ($\chi$, Eq. 4, for the a)
regular, b) diluted, c) random, and d) detailed networks.

##### Regular networks

In the regular lattice models, tortuosity ($\tau$) originates from the
geometry of the grid and from the distribution of the conductances. The
conductances in this model are randomly assigned and not correlated. Thus,
these models may underestimate preferential pathways and, therefore, slightly
overestimate the actual hydraulic tortuosity $\chi$, see Figure 11a.

##### Diluted networks

As outlined above, diluted regular networks may be characterized by singular
regions, that is, disconnected subnetworks, and may also display stagnation
points, such as dead-end pores. These features increase the hydraulic
tortuosity so that the $\tau$ for the diluted networks (Figure 11b) are
generally higher than for the regular networks (Figure 11a).

##### Random networks

The random networks (Figure 11c) have the largest tortuosity, which suggests
that the number of diluted bonds has a more significant effect on the
tortuosity than the formation of stagnation and singular points.

##### Detailed networks

The tortuosity estimated from the detailed network model gives the best
approximation for $\chi$. This can be traced back to the effect of the dead-
end regions on the preferential flow.

###  4.3 Velocity and flow rate distribution

The velocity and flow rate distributions play a key role in the understanding
and prediction of solute transport from the pore and continuum scales.
Variations in the micro-scale flow velocities give rise to the phenomenon of
hydrodynamic dispersion [29]. The distribution of pore-scale velocities
determines the evolution of solute dispersion and the tailing behavior of
solute breakthrough curves [27, 30, 11]. The velocity distribution is a key
property for the prediction of large-scale solute transport. It is not clear a
priori that network models capture the velocity distribution observed in
porous media because 1) networks do not account for the intra-pore velocity
distribution, 2) networks do not reproduce the exact topological features of
the pore space, and 3) the networks may display more stagnation sites, that is
dead-end pores than one finds in the real medium.

The effect of the saturation degree on the spatial pressure and flow
distribution is depicted in Figure 12 for the saturated case and the partially
saturated medium with $S_{e}=0.65$. The saturated case has a relatively
uniform flow velocity at the device scale and local low flow velocity regions
between grains that are aligned in the main flow direction. For the partially
saturated case in contrast, the local pressure gradients diverge from the main
flow direction, causing a more tortuous flow. It is also visible that the
unsaturated case has significant regions of low flow velocities and small
pressure gradients. Figure 13a shows the probability density functions of the
flow velocity from the direct numerical flow simulations with flow velocities
that range over almost ten orders of magnitude. The distributions are
characterized by significant tails toward low velocities. As the complexity of
the fluid-filled medium portion increases, the probability of low velocities
increases, which is due to the low-velocity regions depicted in Figure 12. The
low-velocity end of the velocity distribution determines the long-time tails
of solute breakthrough and the evolution of dispersion.

Figure 12: Spatial distribution of the water pressure (a and c) and flow
velocity (b and d) for the fully saturated (a and b) and the $S_{e}=0.65$ (c
and d) cases, evaluated by the numerical simulations.  Figure 13: a) The
velocity (normalized by the mean flow rate) probability density function (PDF)
of the numerical simulation for the different saturation degrees ($S_{e}$), b)
the flow rate at the bonds (Qb) and the velocity distribution (u, Eq. 16) PDFs
evaluated by the detailed representation of the pore space as a network model,
c) the velocity PDF of the parabolic and strip pore geometries evaluated by
the random network, and d) the bonds flow rate PDFs of the strip pore geometry
evaluated by the regular and diluted network models. In these plots, $\Delta
u_{k}$ intervals are chosen by dividing the logarithmic span of the velocities
arbitrarily into 40 (equal on a logarithmic scale) intervals.

####  4.3.1 Regular, diluted and random networks

Figure 14 shows the spatial pressure and flow rate distributions for the
medium with saturation $S_{e}=0.65$ obtained from the (top panels) regular and
(center panels) diluted and (bottom panels) random networks. The three network
models have the same conductance distribution and mean coordination numbers of
$\langle Z\rangle\approx 3$. The random network matches the distribution of
coordination numbers of the porous medium.

Stagnation points, that is, sites of zero flow rate, are relatively frequent
in the square lattice (diluted network, Figure 14d) but not very common in the
hexagonal lattice (regular network, Figure 14b). This suggests that the
orientation of the pore throats and not just the connectivity of the pore
space is important for pressure and flow rate distributions. For the saturated
case (not shown) both network models give spatial pressure and flow rate
distributions that look qualitatively similar to those from the direct
numerical simulations. Also, the velocity and flow rate PDFs estimated from
the networks compare well to those obtained from the direct flow simulations
(see Figures 13c and d). For decreasing saturation degrees, however, the three
network models are not able to describe either the pressure and flow rate maps
or their PDFs. While none of the network models is able to reproduce the
increase of low velocities with decreasing saturation, the random network
picks up the wide range of flow velocities observed in the direct flow
simulations. These observations indicate that the preferential flow patterns
observed in the unsaturated medium in Figure 12 cannot be explained by the
variability in conductance or coordination number alone.

For the diluted networks, the grid dilution to achieve the desired mean
coordination number of $\langle Z\rangle=3$ increases the number of stagnation
and dead-end points compared to the regular square grid. These points give
rise to the white regions in the $Q_{s}$ map (Figs. 14d), and blue regions in
the pressure map (Fig. 14c). However, dilution also gives rise to a more
pronounced preferential flow as compared to the regular network. The pressure
distributions of the square diluted lattice (figure 14c) are more disordered
than the regular hexagonal pressure distribution (figure 14a), for the same
mean coordination number. This is because of the larger amount of dead-ends
(white regions in figure 14b and 14d) caused by the higher absolute dilution
(starting from a higher saturated $Z$). This illustrates the limitation of a
grid’s mean coordination number to generally represent the medium’s
connectivity.

The random network (Figures 14e and f) shows a smaller number of singular
points and dead-end pores than the diluted network (Fig. 14c and d), despite
the larger degree of dilution (from $Z_{grid}=8$ for the base lattice to
$\langle Z_{model}\rangle\approx 3$). This suggests that a wide distribution
of the coordination number increases the pore-space connectivity. For
illustration, consider the situation of ten pores, nine of which have a
coordination number of $1$ and one of them a coordination number of $9$. Then,
they are all connected.

In general, the size and spatial configuration of the stagnation areas are not
well described by the lattice networks. This is to be expected because these
network models (regular, diluted, and random) do not account for the spatial
correlation of the porous space, only mean properties such as the mean
coordination number and the one-point distributions of conductance and
coordination number (only random network). In these network models, the
changes in the saturation degree of the real media manifest in the total size
assigned to the sites and bonds, or in other words, in the porosity
represented by the network, that is, the ratio of the area of sites and bonds
to the area of the medium. Note that the singular regions in the diluted
networks are merely artifacts of the network generation process and do not
contribute to the estimation of saturation of porosity. Moreover, the air
phase is typically spatially correlated, as shown in Figure 12, which is not
accounted for in these network models. Nevertheless, the random network models
can reproduce some of the features obtained from the direct numerical
simulations, such as the preferential flow paths (Fig. 14f) and the broad
spectrum of flow in flow rates (Figure 13).

Figure 14: Spatial distribution of the pressure (P) at the sites (a, c, and e)
and flow rate at the sites ($Q_{s}$, Eq. 7) (b, d and f) for the regular
(hexagonal) network (a and b), the diluted (square) network (c and d), and the
random network (e and f) of the $S_{e}=0.65$ case with a parabolic pore
geometry. The width of the black line represents the conductance (g) of each
bond in the depicted realization.

####  4.3.2 Detailed networks

Figure 15 shows the pressure and flow rate maps for $S_{e}=1$ and $0.65$
obtained from the detailed network representations of the media. The detailed
networks show similar pressure and flow rate distributions as the numerical
simulations (Figure 12) for the two saturation degrees. Also, the PDFs of flow
rates closely match the ones obtained from the direct flow simulations, as
shown in Figure 13b. The detailed network model accounts only for the
percolated connected fluid phase. Thus, it does not generate singular points
or regions. In this model, the number of stagnation sites with zero velocity
is small and comparable to the number of sites with $Z=1$. This is in
accordance with the numerical simulations in which zero velocity prevails only
at the solid-fluid interface and does not exist at any finite volume (mesh
element). As seen in the map of flow rates both from the direct flow
simulations (Fig. 12), the air phase encloses extensive connected dead-end
regions (of many pores) in which the flow is very slow. These regions do not
participate in the conductance of fluids, which reduces permeability and
increases preferential flow. Moreover, these regions give rise to non-Fickian
solute transport due to diffusion-based solute entrapment and release [31, 32,
11]. All of these features are present in the detailed network models shown in
Figure 15. This observation shows the importance of accounting for the
morphology of the connected fluid phase and pore space and demonstrates the
viability of using network models to quantify not only the mean flow but also
the flow distribution in porous media.

The detailed model is able to reproduce all features of the PDF of flow rates
for the different media (Figure 13b). It accounts correctly for the increase
of low flow rates with decreasing saturation and captures the full range of
flow rates. The flow rate PDFs are very similar in magnitude and shape to the
volumetric flow rate PDFs, indicating that there is no large variability in
the pore cross-sections. For the saturated case, however, the range of flow
rates is much larger than the range of volumetric flow rates, indicating that
the distribution of intrapore velocities is important. However, with
increasing heterogeneity of the pore space, i.e., with decreasing saturation,
this difference vanishes and both flow rates and volumetric flow rates show a
similar range and shape. Thus, the intrapore velocity distribution
reproduction is less important for these cases.

Figure 15: Spatial distribution of the pressure (P) at the sites (a and c) and
flow rate at the sites ($Q_{s}$) (b and d) for the detailed network (c and d)
for the $S_{e}=1$ (a and b) and $S_{e}=0.65$ cases with a parabolic pore
geometry. The width of the black line represents the conductance (g) of each
bond.

##  5 Conclusion

The conceptualization of a porous medium as a network of conductors is
valuable and effective for quantifying flow and transport phenomena and
evaluating the media’s hydraulic properties. The amount of detailed
information required for pore-space characterization depends on the complexity
of the flow phenomenon and the flow and transport processes of interest. In
this study, we demonstrated how different levels of information about the
phase configuration in the pore space can be used to evaluate different
properties of single-phase flow in variably saturated media. The media
permeability, even for the partially saturated cases, is well quantified even
by relatively simple regular network models. For the saturated medium under
consideration, the regular network model gives a good approximation of the
hydraulic tortuosity and velocity distribution, albeit not for the velocity
range. However, for the partially saturated cases, the more complex features,
such as the media tortuosity and velocity distribution, are governed by the
spatial configuration of the fluid phase. Thus, the accurate evaluation of
these factors relies on a relatively detailed description of the pore space,
its spatial distribution and connectivity. This information is available for
the pore space in saturated porous media but not necessarily for the more
complex fluid distribution and fluid-filled pore space in a partially
saturated medium. Addressing this issue requires the quantification of the
spatial configuration of the fluid phase and identifying its relevant
quantifiers.

## Acknowledgments

I.B.N, M.D. acknowledge funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement No. HORIZON-MSCA-2021-PF-01 (USFT). I.B.N. , J.J.H. and M.D.
acknowledge the support of the Spanish Research Agency
(10.13039/501100011033), Spanish Ministry of Science and Innovation through
grants CEX2018-000794-S and HydroPore PID2019-106887GB-C31. J.J.H.
acknowledges the support of the Spanish Research Agency
(10.13039/501100011033), the Spanish Ministry of Science and Innovation and
the European Social Fund “Investing in your future” through the “Ramón y
Cajal” fellowship (RYC-2017-22300).

##  6 Open Research

The data used in this work is provided as supplementary material for the
review process and later will be made available in a data repository with a
permanent identifier.

##  Appendix A Effective medium approximation (EMA)

The effective medium approximation (EMA) states that the network of variable
conductors can be represented by a regular uniform network with an effective
conductance $g_{e}$ . The EMA automatically obeys the smooth field
approximation [33], i.e., the macroscopic pressure field is a smooth function
of position in the porous medium [23], or in other words, the pressure
gradient at each pore is the projection of the macroscopic pressure gradient
in the direction of the pore axis [34]. The effective conductivity $g_{e}$ is
calculated from the conductance distribution $f_{g}(g)$ by requiring that the
average of the flow perturbation caused by the EMA is zeroed [35]

$\int\frac{g-g_{e}}{g+(\frac{Z}{2}-1)g_{e}}f(g)dg=0$ (25)

Then the EMA permeability of the media ($k_{EMA}$ [m2]) can be evaluated from
$g_{e}$

$k_{EMA}=\frac{\mu N_{W}}{bW}g_{e}$ (26)

where $N_{W}$ is the number of junctions on the lattice cross-section
perpendicular to the main flow direction. For a square grid
$N_{W}=\frac{W}{l_{sq}}$ while for the hexagonal case, where the hexagons are
oriented with the long diagonal in the direction of the flow,
$N_{W}=\frac{W}{\sqrt{3}l_{hex}}$.

For the diluted network, we describe the conductance distribution in a binary
form as a summation of two terms [34], a fraction $p$ of diluted bonds with
zero conductance and the fraction of the surviving bonds ($1-p$) so that
$f_{1}(g)=p\delta(g)+(1-p)f(g)$, where $\delta(g)$ is the Dirac delta function
and $f(g)$ is the conductance probability density function (4).

## References

  * Blunt [2001] M. J. Blunt,  Flow in porous media—pore-network models and multiphase flow,  Current opinion in colloid & interface science 6 (2001) 197–207.
  * Joekar-Niasar and Hassanizadeh [2012] V. Joekar-Niasar, S. Hassanizadeh,  Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: A review,  Critical reviews in environmental science and technology 42 (2012) 1895–1976.
  * Fatt [1956] I. Fatt,  The network model of porous media,  Transactions of the AIME 207 (1956) 144–177.
  * Koplik [1982] J. Koplik,  Creeping flow in two-dimensional networks,  Journal of Fluid Mechanics 119 (1982) 219–247.
  * Blunt [2017] M. J. Blunt, Multiphase flow in permeable media: A pore-scale perspective, Cambridge university press, 2017.
  * Dullien [1975] F. A. Dullien,  New network permeability model of porous media,  AIChE Journal 21 (1975) 299–307.
  * Hoogland et al. [2016] F. Hoogland, P. Lehmann, R. Mokso, D. Or,  Drainage mechanisms in porous media: From piston-like invasion to formation of corner flow networks,  Water Resources Research 52 (2016) 8413–8436.
  * Mason and Morrow [1991] G. Mason, N. R. Morrow,  Capillary behavior of a perfectly wetting liquid in irregular triangular tubes,  Journal of Colloid and Interface Science 141 (1991) 262–274. doi:https://doi.org/10.1016/0021-9797(91)90321-X.
  * Ben-Noah et al. [2023] I. Ben-Noah, S. P. Friedman, B. Berkowitz,  Dynamics of air flow in partially water-saturated porous media,  Reviews of Geophysics (2023) e2022RG000798.
  * Van Genuchten and Cleary [1979] M. T. Van Genuchten, R. Cleary,  Movement of solutes in soil: Computer-simulated and laboratory results,  in: Developments in soil science, volume 5, Elsevier, 1979, pp. 349–386.
  * Ben-Noah et al. [2023] I. Ben-Noah, J. J. Hidalgo, J. Jimenez-Martinez, M. Dentz,  Solute trapping and the mechanisms of non-fickian transport in partially saturated porous media,  Water Resources Research 59 (2023) e2022WR033613.
  * Jiménez-Martínez et al. [2017] J. Jiménez-Martínez, T. Le Borgne, H. Tabuteau, Y. Méheust,  Impact of saturation on dispersion and mixing in porous media: Photobleaching pulse injection experiments and shear-enhanced mixing model,  Water Resources Research 53 (2017) 1457–1472.
  * de Anna et al. [2017] P. de Anna, B. Quaife, G. Biros, R. Juanes,  Prediction of the low-velocity distribution from the pore structure in simple porous media,  Phys. Rev. Fluids 2 (2017) 124103. doi:10.1103/PhysRevFluids.2.124103.
  * Ben-Noah et al. [2024] I. Ben-Noah, J. J. Hidalgo, M. Dentz, Efficient pore space characterization based on the curvature of the distance map, 2024\. arXiv:2403.12591.
  * Homsy [1987] G. M. Homsy,  Viscous fingering in porous media,  Annual Review of Fluid Mechanics 19 (1987) 271–311. doi:https://doi.org/10.1146/annurev.fl.19.010187.001415.
  * COMSOL® [2018] COMSOL®,  Introduction to comsol multiphysics®,  v. 5.4. www.comsol.com., COMSOL AB, Stockholm, Sweden (2018).
  * Muskat and Botset [1931] M. Muskat, H. Botset,  Flow of gas through porous materials,  Physics 1 (1931) 27–47.
  * Koponen et al. [1996] A. Koponen, M. Kataja, J. Timonen,  Tortuous flow in porous media,  Physical Review E 54 (1996) 406–410. doi:https://doi.org/10.1103/PhysRevE.54.406.
  * Gostick et al. [2016] J. Gostick, M. Aghighi, J. Hinebaugh, T. Tranter, M. A. Hoeh, H. Day, B. Spellacy, M. H. Sharqawy, A. Bazylak, A. Burns, et al.,  Openpnm: a pore network modeling package,  Computing in Science & Engineering 18 (2016) 60–74.
  * Liu et al. [2024] Y. Liu, W. Gong, H. Xiao, M. Wang,  A pore-scale numerical framework for solute transport and dispersion in porous media,  Advances in Water Resources 183 (2024) 104602. URL: http://dx.doi.org/10.1016/j.advwatres.2023.104602. doi:10.1016/j.advwatres.2023.104602.
  * Boussinesq [1868] J. Boussinesq,  Mémoire sur l’influence des frottements dans les mouvements réguliers des fluids,  Journal de mathématiques pures et appliquées 13 (1868) 377–424.
  * Hunt and Sahimi [2017] A. G. Hunt, M. Sahimi,  Flow, transport, and reaction in porous media: Percolation scaling, critical-path analysis, and effective medium approximation,  Reviews of Geophysics 55 (2017) 993–1078.
  * Sahimi et al. [1990] M. Sahimi, G. R. Gavalas, T. T. Tsotsis,  Statistical and continuum models of fluid-solid reactions in porous media,  Chemical Engineering Science 45 (1990) 1443–1502.
  * Sahimi [1998] M. Sahimi,  Non-linear and non-local transport processes in heterogeneous media: from long-range correlated percolation to fracture and materials breakdown,  Physics Reports 306 (1998) 213–395.
  * Raoof and Hassanizadeh [2010] A. Raoof, S. M. Hassanizadeh,  A new method for generating pore-network models of porous media,  Transport in porous media 81 (2010) 391–407.
  * Ben-Noah et al. [2024] I. Ben-Noah, J. J. Hidalgo, M. Dentz, Note on using singular value decomposition to solve diluted pore networks, 2024. arXiv:2403.13462.
  * Dentz et al. [2018] M. Dentz, M. Icardi, J. J. Hidalgo,  Mechanisms of dispersion in a porous medium,  Journal of Fluid Mechanics 841 (2018) 851–882. doi:https://doi.org/10.1017/jfm.2018.120.
  * Puyguiraud et al. [2021] A. Puyguiraud, P. Gouze, M. Dentz,  Pore-scale mixing and the evolution of hydrodynamic dispersion in porous media,  Physical Review Letters 126 (2021) 164501.
  * Bear [1988] J. Bear, Dynamics of fluids in porous media, Courier Corporation, 1988.
  * Gouze et al. [2021] P. Gouze, A. Puyguiraud, T. Porcher, M. Dentz,  Modeling longitudinal dispersion in variable porosity porous media: control of velocity distribution and microstructures,  Frontiers in Water 3 (2021) 766338.
  * Liu and Kitanidis [2012] Y. Liu, P. K. Kitanidis,  Applicability of the dual-domain model to nonaggregated porous media,  Groundwater 50 (2012) 927–934. doi:https://doi.org/10.1111/j.1745-6584.2011.00909.x.
  * Gjetvaj et al. [2015] F. Gjetvaj, A. Russian, P. Gouze, M. Dentz,  Dual control of flow field heterogeneity and immobile porosity on non-fickian transport in berea sandstone,  Water Resources Research 51 (2015) 8273–8293. doi:https://doi.org/10.1002/2015WR017645.
  * Burganos and Sotirchos [1987] V. Burganos, S. V. Sotirchos,  Diffusion in pore networks: effective medium theory and smooth field approximation,  AIChE journal 33 (1987) 1678–1689.
  * Friedman et al. [1995] S. P. Friedman, L. Zhang, N. A. Seaton,  Gas and solute diffusion coefficients in pore networks and its description by a simple capillary model,  Transport in porous media 19 (1995) 281–301.
  * Kirkpatrick [1973] S. Kirkpatrick,  Percolation and conduction,  Reviews of modern physics 45 (1973) 574.


# Covariant Fermionic Fields of Space-like Particles

Radhakrishnan Balu  Army Research Laboratory Adelphi, MD, 21005-5069, USA
<EMAIL_ADDRESS>
(Received: date / Accepted: )

###### Abstract

We develop covariant fermionic fields of space-like particles. As an
application of the formalism we discuss the example of superluminous tachyons
with imaginary rest mass and spin 1/2 forming fermionic ensembles that is
relativistically covariant. We use the tachyonic single particle Hilbert space
to build the quantum stochastic fields as this notion is more consensual,
instability in vacuum quantum field than faster than light particles, in field
theory with applications to superconductivity via tachyonic condensation. We
also establish the localizability and as a consequence the causality of
tachyons (in terms of anticommutators) using the systems of imprimitivity
framework and derive the Dirac equation for the space-like particle that has
an additional helicity term.

##  1 Introduction

In our earlier work we treated covariant Quantum Fields via Lorentz Group
Representation of Weyl Operators [1], [2] constructed Fock spaces for massless
Weyl fermions and covariant white noises. Here, we develop methods for second
quantizing imaginary mass fermions conforming to special theory of relativity
that are space-like particles, that is based on systems of imprimitivity, to
field theoretic context. Systems of imprimitivity [5], [3], [4] is a powerful
characterization of dynamical systems, when the configuration space of a
quantum system is described by a group, from which infinitesimal forms in
terms of differential equations ($Shr\ddot{o}dinger$, Heisenberg, and Dirac
etc), localizability, and the canonical commutation relations can be derived
[6] and [7]. In Varadarajan’s work [4] space-like particles are not considered
as imaginary mass particles with speed greater than light and the non
existence of finite unitary representation for the little group created a lack
of interest in general. However, there are some recent work in studying
tachyons [10] and overcoming the group representation problem that inspired
our work. Let us look at the 2 + 1 spacetime case considered in this work as
the little group in this case is O(2, 1) with the Lie algebra generators
satisfying the following conditions:

$[J_{z},K_{x}]=-ik_{y};[J_{z},K_{y}]=-ik_{x};[K_{y},K_{x}]=-iJ_{z};$ (1)

Where the generators are given by the Pauli matrices as:

$J_{z}=1/2\sigma_{3};K_{x}=i/2\sigma_{1};K_{y}=i/2\sigma_{2};$ (2)

An operator of the form $U=e^{it(\alpha J_{z}+\beta K_{x}+\gamma K_{y})}$ is
in general not unitary. However, as Schwartz [10] observed the relation

$\displaystyle U^{\dagger}$ $\displaystyle=U\sigma_{3}U^{-1}\text{ in
conjunction with the metric }$ (3)
$\displaystyle\langle\psi,\sigma_{3}\phi\rangle\rightarrow\langle
U\psi,\sigma_{3}U\phi\rangle$
$\displaystyle=\langle\psi,U^{\dagger}\sigma_{3}U\phi\rangle=\langle\psi,\sigma_{3}\phi\rangle.$
(4)

Another way to look at the above expression is to use $\sqrt{\sigma_{3}}$
which is also unitary and can be absorbed in to U as

$\displaystyle\langle U\psi,\sigma_{3}U\phi\rangle$
$\displaystyle=\langle\sqrt{\sigma_{3}}^{\dagger}U\psi,\sqrt{\sigma_{3}}U\phi\rangle.$
(5) $\displaystyle\sqrt{\sigma_{3}}$ $\displaystyle=\begin{bmatrix}1&0\\\
0&i\end{bmatrix}.$ (6)

Borrowing ideas from quantum information processing it is like rotating all
the components of the fibers, each fiber would be represented by two qubits,
by the same angle that leaves the systems invariant.

As we will see later a family of metrics including the above one arise
naturally in describing the covariant spinor fields in our framework.
Formally, a spinor field is a map from a G-space where Poincar$\grave{e}$
group acts to the spin group representation. Schwartz considers the operator
$\sigma_{3}$ as helicity and as a way to distinguish particles against anti-
particles as the mass shell hyperbola in the case of space-like particles is a
single sheet one. We will take this consideration in building our fiber bundle
to describe the Hibert space of the space-like particle.Our frame is based on
orthomoduar lattices, quantum logic, where the above inner product can arise
as the duality $\eta$ which is a $\theta-$bilinear form satisfying:

$\displaystyle\langle x_{1}+x_{2},y\rangle$ $\displaystyle=\langle
x_{1},y\rangle+\langle x_{2},y\rangle.$ (7) $\displaystyle\langle
x_{1},y_{1}+y_{2}\rangle$ $\displaystyle=\langle x,y_{1}\rangle+\langle
x,y_{2}\rangle.$ (8) $\displaystyle\langle cx,dy\rangle$
$\displaystyle=c\langle c,y\rangle d^{\theta}.$ (9)

To elaborate further on this, by choosing different values of $\theta$
different inner product can be defined. A polarity $\eta$ is duality which is
also involutive, $(\eta(\eta(x))=x)$, and orthocomplemented lattices are
endowed with a polarity. Finally, we note that covariant metrics are required
to describe spinor fields, as they are not topological, of both superluminous
(faster than light) and subluminous (slower than light) particles.

Our quantization is based on SI and an important theorem by Mackey that
characterizes such systems in terms of induced representations, and so let us
recollect key notions in Clifford algebras, spinor fields, and Schwartz spaces
[4] before discussing our main results on covariant tachyonic fields and
localizability of the particles by allowing imaginary values for momentum [9].

##  2 Little groups (stabilizer subgroups)

Some of the different systems of imprimitivity that live on the orbits of the
stabilizer subgroups are described below. It is good to keep in mind the
picture that SI is an irreducible unitary representation of Poincar$\grave{e}$
group $\mathscr{P}^{+}$ induced from the representation of a subgroup such as
$SO_{3}$, which is a subgroup of homogeneous Lorentz, as
$(U_{m}(g)\psi)(k)=e^{i\\{k,g\\}}\psi(R^{-1}_{m}k)$ where g belongs to the
$\mathscr{R}^{4}$ portion of the Poincar$\grave{e}$ group, m is a member of
the rotation group, and the duality between the configuration space
$\mathbb{R}^{4}$ and the momentum space $\mathbb{P}^{4}$ is expressed using
the character the irreducible representation of the group $\mathbb{R}^{4}$ as:

$\displaystyle\\{k,g\\}$
$\displaystyle=k_{0}g_{0}-k_{1}g_{1}-k_{2}g_{2}-k_{3}g_{3},p\in\mathbb{P}^{4}.$
(10) $\displaystyle\hat{p}:x$ $\displaystyle\rightarrow e^{i\\{k,g\\}}.$ (11)
$\displaystyle\\{Lx,Lp\\}$ $\displaystyle=\\{x,p\\}.$ (12) $\displaystyle\
\hat{p}(L^{-1}x)$ $\displaystyle=\hat{Lp}(x).$ (13)

In the above L is a matrix representation of Lorentz group acting on
$\mathbb{R}^{4}$ as well as $\mathbb{P}^{4}$ and it is easy to see that
$p\rightarrow Lp$ is the adjoint of L action on $\mathbb{P}^{4}$. The
$\mathbb{R}^{4}$ space is called the configuration space and the dual
$\mathbb{P}^{4}$ is the momentum space of a relativistic quantum particle.

The stabilizer subgroup of the Poincar$\grave{e}$ group $\mathscr{P}^{+}$
(Space-like particle) can be described as follows: [12]: The Lorentz frame in
which the walker is at rest has momentum proportional to (0,1,0,0) and the
little group is again SO(3, 1) and this time the rotations will change the
helicity.

##  3 Smooth Orbits in Momentum Space

We need to described few ingredients to construct the Hilbert space of space-
like fermions namely, the fiber bundle, the fiber vector space, an inner
product for the fibers, and an invariant measure. The 3+1 spacetime Lorentz
group $\hat{O}(3,1)$-orbits of the momentum space $\mathscr{R}^{4}$, where the
systems of imprimitivity established will live, described by the symmetry
$\hat{O}(3,1)\rtimes\mathbb{R}^{4}$. The orbits have an invariant measure
$\alpha^{+}_{m}$ whose existence is guaranteed as the groups and the
stabilizer groups concerned are unimodular and in fact it is the Lorentz
invariant measure $\frac{dp}{p_{0}}$ for the case of forward mass hyperboloid.
The orbits are defined as:

$\displaystyle Y^{+,1/2}_{im}$
$\displaystyle=\\{p:p^{2}_{0}-p^{2}_{1}-p^{2}_{2}-p^{2}_{3}=-m^{2},p_{0}>0\\},\text{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}
single sheet hyperboloid}.$ (14)

Each of these orbits are invariant with respect to $\hat{O}(3,1)$ and let us
consider the stabilizer subgroup of the first orbit at p=(0,1,0,0). Now,
assuming that the spin of the particle is 1/2 and mass $m>0$ let us define the
corresponding fiber bundles (vector) for the positive mass hyperboloid that
corresponds to the positive-energy states by building the total space as a
product of the orbits and the group $SL(2,C)$.

$\displaystyle\hat{B}^{+,1/2}_{im}$ $\displaystyle=\\{(p,v)\text{
}p\in{\hat{Y}^{+1/2}_{im},}\text{
}v\in\mathscr{C}^{4},\sum_{r=0}^{3}p_{r}\gamma_{r}v=mv\\}.$ (15)
$\displaystyle\hat{\pi}$ $\displaystyle:(p,v)\rightarrow{p}.\text{ Projection
from the total space }\hat{B}^{+,2}_{im}\text{ to the base
}\hat{Y}^{+,1/2}_{im}.$ (16)

It is easy to see that if $(p,v)\in B_{im}^{+,1/2}$ then so is also
$(\delta(h)p,S(h^{*-1})v)$. Thus, we have the following Poincar$\grave{e}$
group symmetric action on the bundle that encodes spinors into the fibers:

$h,(p,v)\rightarrow(p,v)^{h}=(\delta(h)p,S(h^{*-1})v).$ (17)

###### Proposition 1

There is a fiber representation of spin 1/2 space-like particles with spinor
fields supporting a family of covariant inner products.

For $m>0$ the fiber of $B_{im}^{+,1/2}$ at $p^{(}m)=(0,m,0,0)$ is spanned by
the vectors $\\{e_{1}+ie_{4},e_{2}+ie_{3}\\}$ and so $B_{im}^{+,1/2}$ is two
dimensional with the covariant form $-im^{-1}\langle\gamma_{3}v,v\rangle$.
This positive definite due to skew symmetric $\gamma_{3}$ and
$-im^{-1}\langle\gamma_{3}v,v\rangle$ is purely imaginary. This 1-form is
covariant with respect to the spin representation
$S(h)^{*}\gamma_{3}S(h)=\gamma_{3}$. As we mentioned earlier, both
superluminal and subluminal particles support such a family of covariant
1-forms and thus metrics in describing their spinor fields. The same set of
arguments can be applied to $B_{-im}^{+,1/2}$ implying that the bundles
$B_{-im}^{+,1/2}$ also converge to $B_{0}^{+}$. The endomorphism (chirality or
helicity operator) $\Gamma=i\gamma_{0}\gamma_{1}\gamma_{2}\gamma_{3}$
transforms $B_{im}^{+,1/2}(p)\rightarrow B_{-im}^{+,1/2}(p),\forall p\in
X_{im}^{+},m>0$ as it anticommutes with all the gammas. In the limit $\Gamma$
leaves the fibers of $B_{0}^{+}$ invariant leading to higher degeneracies with
$\Gamma=\begin{bmatrix}1&0\\\ 0&-1\end{bmatrix}$. This means $\Gamma$ commutes
with all of $S(h)$ implying that $(p,v)\in B_{0}^{+}\Rightarrow(p,\Gamma v)\in
B_{0}^{+}$.

$\Gamma$ has eigen values $\pm 1$ at fiber (0, 1, 0, 0) and hence is true of
all fibers. The stability group $E^{*}$ at (0, 1, 0, 0) given by the
generators mentioned during the introduction. is still invariant and we have
$E_{2}$ is the stabiliser group at (0,1,0,0).

Now, we can define the states of the space-like particles on the Hilbert space
$\hat{\mathscr{H}}^{+,1/2}_{im}$, square integrable functions on Borel
sections of the bundle $\hat{B}^{+,1/2}_{im}=\\{(p,v):(p,v)\in
B_{im}^{+,1/2},\Gamma v=\mp v\\}$.

with respect to the invariant measure $\beta^{+,1}_{im}$, whose norm induced
by the inner product is given below:

$\left\|\phi\right\|^{2}=\int_{X^{+}_{m}}p_{3}^{-1}\langle\phi{p},\phi{p}\rangle.{d\beta}^{+,1/2}_{im}(p).$
(18)

The invariant measure and the induced representation of the Poincar$\grave{e}$
group from that of the Weyl fermion are given below:

$\displaystyle{d\beta}^{+,1/1}_{im}(p)$
$\displaystyle=\frac{dp_{1}dp_{2}dp_{3}}{2(-m^{2}+p_{1}^{2}+p_{2}^{2}+p_{3}^{2})},m>0,\text{imaginary
values for p are allowed}.$ (19) $\displaystyle(U_{h,x}\phi)(p)$
$\displaystyle=\text{exp i}\\{x,p\\}\phi(\delta(h)^{-1}p)^{h}.$ (20)

We need Schwartz spaces, and their duals the tempered distributions to
guarantee Fourier transforms, to move with ease between configuration and
momentum space descriptions for Dirac equation as well as localizability and
so let us define them.

###### Definition 2

Let V be an n-dimensional Euclidean space with the positive definite product
$\langle.,.\rangle$ and Lebesgue measure dv. for any translation invariant
differential operator D and any complex polynomial q on V, the function
$\phi\rightarrow sup_{v\in V}\left|q(v)D(\phi)(v)\right|$ is a seminorm on
$C_{c}^{\infty}(V)$ and the collection of these seminorms induces a locally
convex topology for $C_{c}^{\infty}(V)$. Its completion may be identified with
$\mathscr{S}(V)$, the space of $C^{\infty}$ functions on V for which
$sup_{v\in V}\left|q(v)D(\phi)(v)\right|<\infty$ for all q and D. Intuitively,
it is the function space of all infinitely differentiable functions that are
rapidly decreasing at infinity along with all partial derivatives

Let $\\{v_{1},\dots,v_{n}\\}$ be an orthonormal basis for V and
$D_{v_{1},\dots,v_{n}}=(\partial/\partial
x_{1})^{\nu_{1}}\dots(\partial/\partial x_{n})^{\nu_{n}}$ linearly span the
algebra of translation invariant differential operators, and
$x_{1},\dots,x_{n}$ be the linear coordinate functions on V associated with
the chosen basis. A topology on $\mathscr{S}(V)$ can be induced by the
collection of seminorms
$sup_{x_{1}\dots,x_{n}}abs{(1+x_{1}^{2}+\dots+x_{n}^{2})^{k}(D_{v_{1},\dots,v_{n}}\phi)(x_{1},\dots,x_{n})}$
for various k=0,1,2,…, and
$\nu_{1},\dots,\nu_{n}\geq 0$. A tempered distribution E is a complex valued
linear functional on $C_{c}^{\infty}(V)$ which is continuous in the topology
defined. By extending these functionals to $\mathscr{S}(V)$ we may regard a
tempered distribution as a continuous linear functional on the Schwartz space
$\mathscr{S}(V)$.

Fourier transform is an an automorphism on Schwartz space as

$\hat{\phi(x)}=(2\pi)^{-n/2}\int_{V}exp[-i\langle
x,v\rangle]\phi(v)dv,\forall\phi\in\mathscr{S}(V),x\in V.$ (21)

The Fourier transform of a tempered distribution $E(\phi\rightarrow
E(\phi),\phi\in\mathscr{S}(V))$ is given by $\hat{E}(\phi)=E(\hat{\phi})$. We
denote the dual space of $\mathscr{S}(V)$ Schwartz functions that consists of
the tempered distributions as $\mathscr{S}^{\prime}(V)$ and they form the
Gelfand nuclear triple
$\mathscr{S}(\mathscr{H}_{0}^{+2})\subset\mathscr{H}_{0}^{+2}\subset\mathscr{S}^{\prime}(\mathscr{H}_{0}^{+2})$.

Let us now state and establish the main result for the case of massive Weyl
fermion with space-like momentum from the representation of a subgroup (space-
like little group) of Poincarè.

###### Theorem 3

Space-like Weyl representation of Poincar$\grave{e}$ group is a transitive
system of imprimitivity.

Proof: Let $g:\mathscr{G}\rightarrow U_{g}(\hat{\mathscr{H}}^{+,1/2}_{im})$ be
a homomorphism from the two dimensional Euclidean group $\mathscr{G}=E_{2}$ to
the unitary representation of the group in $\hat{\mathscr{H}}^{+,1/2}_{im}$.
We note that it is a stabilizer subgroup which is also closed at the momentum
(0,0,m,0) and so $H/\mathscr{G}$ is a transitive space.

Consider a map $v(g):\mathscr{G}\rightarrow\hat{\mathscr{H}}^{+,1/2}_{im}$
satisfying the first order cocycle relation
$v(gh)=v(g)+U_{g}v(h),g,h\in\mathscr{G}$. An example of such a map is the
following: [11]

$\displaystyle\mathscr{H}$
$\displaystyle=\oplus_{j=0}^{\infty}\mathscr{H}_{j}.$ $\displaystyle H$
$\displaystyle=1\oplus\oplus_{j=1}^{\infty}H_{j}.$ $\displaystyle U_{t}$
$\displaystyle=e^{-itH},t\in\mathscr{G}.$ $\displaystyle v(t)$
$\displaystyle=tu_{0}\oplus\oplus_{j=1}^{\infty}(e^{-itH_{j}}u_{j}-u_{j}).$

Now, we can define the Weyl operator $V_{g}=W_{g}(v(g),U_{g})$ where
$g\in\mathscr{G}$ in the Fock space
$\Gamma_{s}(\hat{\mathscr{H}}^{+,1/2}_{im})$.

This is a projective unitary representation satisfying the commutator relation
$V_{g}V_{h}=e^{iIm\langle v_{g},U_{g}v_{h}\rangle}V_{h}V_{g}$ and let us
denote the homomorphism from $\mathscr{G}$ to $V_{g}$ by m. This guarantees a
map b that satisfies $b(gh)=b(g)m(h),g\in\mathscr{P},h\in\mathscr{G}$ and such
map can be constructed by considering the map $c(x\rightarrow c(x)$ as Borel
section of $\mathscr{P}/\mathscr{G}$ (the choice of this section not a
canonical one but immaterial to our purpose here) with $c(x_{0})=e$. The map
$\beta$ maps $g\in\mathscr{P}\rightarrow g\mathscr{G}$

$\displaystyle a(g)$ $\displaystyle=c(\beta(g))^{-1}.$ (22) $\displaystyle
b(g)$ $\displaystyle=m(a(g)).$ (23)

Then the strict cocycle $\phi$ satisfies
$\phi(g_{1},g_{2})=b(g_{1}g_{2})b(g_{2})^{-1}$.

We can now set the SI relation using the above cocycle as follows:

$\displaystyle U_{h,x}\phi(p)$
$\displaystyle=expi\\{x,p\\}\phi(g,g^{-1}x)f(g^{-1}x),f\in\mathscr{H}$
$\displaystyle\text{ character representation is defined in equation
}\eqref{eq: poissonBr}.$ $\displaystyle P_{E}F$ $\displaystyle=\chi_{E}f\text{
Position operator}.$

We can construct the conjugate pair of field operators for the Fock space
$\Gamma_{s}(\hat{\mathscr{H}}^{+,1/2}_{0}))$ as follows:

Let $p_{g}$ be the stone generator for the family of operators
$P_{gt,p},g\in\mathscr{P},t\in\mathbb{R}$ and q(g) = p(ig) and we get the
creation and annihilation operators as $a(g)^{\dagger}=\frac{1}{2}(q(g)-ip(g)$
and $a(g)=\frac{1}{2}(q(g)+ip(g)$.

Strictly speaking we need to build fermionic Fock space and set up the
corresponding annihilation and creation operators. Instead, we have
constructed the bosonic Fock space, as there is no counterpart for Weyl second
quantized operators in fermionic Fock space, and using the unification of the
[13] two spaces, $dB=(-1)^{\Lambda}dA$, we can construct the fermionic
processes. This involves fashioning reflection processes and we leave out the
details as our main focus here is the white noise process. $\blacksquare$

##  4 Localizability

We establish the localizability of tachyons of spin 1/2 by constructing a
projection valued measure (PVM) as part of the SI relation. As our main
theorem establishes an SI that comes with a PVM on spacetime
$\mathscr{R}^{4}$. We would like to have a PVM on $\mathscr{R}^{3}$ to signify
localization on space only degrees of freedom. We modify the proof from
Varadarajan’s book for the particle with real mass to tachyons of theorem
9.16. This theorem is basis for establishing the fact that photons are not
localizable. The main line of arguments run like this: we consider a Lorentz
transformatiom from an inertial frame that leads to a configuration separated
space-only coordinates. We can denote that unitary transformation by
$U_{h,a},a=a_{0},\hat{a}$ and call this subgroup of Poincar$\grave{e}$ as M.
As series of lemmas on equivalent induced representations on M and necessary
and sufficient conditions for the existence of PVM on $\mathscr{R}^{3}$
culminates in the following theorem.

###### Theorem 4

[4] Let the representation U associated with a system $\mathscr{S}$ be defined
by a pair $(X^{+}_{m},\pi)$, where $m\geq 0$, and $\pi$ is a representation of
the stability subgroup of $L^{*}$ at the point p = (m, 0, 0 ,0) when m ¿ 0 and
(1, 0, 0, 1). The system is localizable when m ¿ 0. If m = 0, $\mathscr{S}$ is
localizable if and only if there exists a representation $\pi^{\prime}$ of
$K^{*}$ such that $\pi$ and $\pi^{\prime}$ have equivalent restrictions to the
group $T^{\prime}$ of all diagonal matrices of $K^{*}$.

As we mentioned earlier we want to keep in mind there is a finite unitary
representation for the space-like little group using the generalized inner
product with Hermitian matrix $\sigma_{3}$. This theorem can be modified for
the case of tachyons with non zero mass by considering the stability subgroup
at p = (0, 0, 0, m) and selecting a strict $(H^{*},X_{im})$-cocycle C with
values in the Hilbert space $\mathscr{K}$ in which $\pi$ acts such that C
defines the representation $\pi$ at p. The unitary representation U acts in
$\mathscr{H}=\mathscr{L}^{2}(X_{im},\mathscr{K},\beta^{+}_{im})$, and for
$(h,a)\in M$ we have

$(U_{h,a}f)(p)=exp[-i\langle
a,p\rangle]C(h,\delta(h)^{-1}p)f(\delta(h)^{-1}p).$ (24)

The proof for the above theorem depends upon a lemma 9.12 [4] that provides
the necessary and sufficient conditions for existence of projection valued
measures.

###### Lemma 5

Let V be a representation of M in some Hilbert space $\mathscr{H}$. Let m be a
real number greater or equal to zero. Then, the existence of a projection
valued measure P based on $R^{3}$ such that (V, P) is a system of
imprimitivity, it is necessary and sufficient that there exists a
representation s, $s\rightarrow s(h)$ of $K^{*}$ in Hilbert space
$\mathscr{K}$ such that V is equivalent to the representation $V^{*}$ of M in
$\mathscr{L}^{2}(R^{3},\mathscr{K},\alpha_{m})$ defined by

$\displaystyle(V^{s}_{h,a}f)(y)$ $\displaystyle=exp[-i\langle
a,y\rangle]s(h)f(\delta(h)^{-1}y).$ (25) $\displaystyle d\alpha_{m}$
$\displaystyle=\frac{dy}{2(m^{2}+\langle y,y\rangle)^{1/2}}.$ (26)

##  5 Dirac equation

Let us rewrite the equation 15 in the single particle Hilbert space to match
the final form available in the literature [9] and then apply the Fourier
transform as we are working with Schwartz spaces the transform is well
defined.

$\displaystyle\sum_{r=0}^{3}p_{r}\gamma_{r}v$ $\displaystyle=mv\\},\Gamma
v=\mp v.$ (27) $\displaystyle(\sum_{r=0}^{3}ip_{r}\gamma_{r}-\Gamma m)v$
$\displaystyle=0.\Rightarrow$ (28)
$\displaystyle(\sum_{r=0}^{3}i\frac{\partial}{\partial x_{r}}\gamma_{r}-\Gamma
m)\Phi$ $\displaystyle=0.\text{ Schwartz spaces allow Fourier transform}$ (29)

##  6 Causality for Dirac Tachyons

Let us now look at how causality arise from SI relations on space-like
particles lifted to second quantization for a simple situation. As we have
constructed Weyl operators that form a system of imprimitivity (V, P) on the
space $L^{2}((E^{*}),\mu)\cong\Gamma_{s}(\hat{\mathscr{H}}^{+,1/2}_{im})$ by
induced an representation of the
Poincar$\grave{e}\mathscr{P}^{*}=H^{*}\rtimes\mathbb{R}^{4}$ group we have

$\displaystyle\alpha\rightarrow U_{\alpha}$ $\displaystyle(U_{\alpha}f)(x)$
$\displaystyle=x(a)f(x),(x\in\mathbb{R}^{4}).$ $\displaystyle(W_{ah}f)(x)$
$\displaystyle=x(a)(V_{h}f)(x),(x\in\mathbb{R}^{4}).$ $\displaystyle
U_{\alpha}\phi(x)U^{-1}_{\alpha}$ $\displaystyle=\phi(\alpha
x),\phi(x)\in\Gamma_{s}(\hat{\mathscr{H}}^{+,1/2}_{im}).$

Let us start with the anti-commutator of Schwartz distributions as
$[\phi(x),\psi(x)]_{+}=\delta(x,y),\phi,\psi\in\Gamma_{s}(\hat{\mathscr{H}}^{+,1/2}_{im})$
and apply the SI relation to get the anticommutation relation [10]

$\displaystyle[\phi(x),\psi(x)]_{+}$ $\displaystyle=\delta^{3}(x,y).$
$\displaystyle\Rightarrow[U_{\alpha}\phi(x)U^{-1}_{\alpha},U_{\alpha}\psi(x)U^{-1}_{\alpha}]_{+}$
$\displaystyle=\delta^{3}(\alpha x,\alpha
y)=[\phi(x),\psi(x)]_{+}=\delta^{3}(x,y).$
$\displaystyle\Rightarrow[\phi(x),\psi(x)]_{+}$
$\displaystyle=\delta^{3}(x-y).$

##  7 Summary and conclusions

We derived the covariant field operators for space-like particles using
induced representations of groups and expressed them in terms of systems of
imprimitivity. We established the results for the massless Weyl fermion case
by inducing a representation of Poincar$\grave{e}$ group from the subgroup
that is a stabilizer at the momentum (0,0,1,0). Finally, we derived the Dirac
equation and established localizability for tachyons.

## References and Notes

  * [1] Radhakrishnan Balu: Covariant Quantum Fields via Lorentz Group Representation of Weyl Operators, Proceedings of the International Conference on Quantum Probability and Related Topics QP38 (2020).
  * [2] Radhakrishnan Balu: Covariant Quantum White Noise from Light-like Quantum Fields, Special issue on Quantum Probability and Related Topics, Journal of Stochastic Analysis, JSOA (Accepted (2020)).
  * [3] Infinite dimensional group representations, Bulletin of the American Mathematical Society, 69, 628 (1963).
  * [4] V.S. Varadarajan: Geometry of quantum theory, Springer (1985).
  * [5] Wigner: Unitary representations of the inhomogeneous Lorentz group, Annals of mathematics, 40, 149 (1939).
  * [6] Newton T D and Wigner E P 1949 Localized states for elementary systems Rev. Mod. Phys. 21 400
  * [7] Radhakrishnan Balu: Kinematics and Dynamics of Quantum Walk in terms of Systems of Imprimitivity, J. Phys. A: Mathematical and Theoretical (2019).
  * [8] Hida, T.; Streit, L. On quantum theory in terms of white noise. Nagoya Math. J. 68 (1977), 21–34
  * [9] U. D. Jentschura and B. J. Wundt, Eur. Phys. J. C 72, 1894 (2012).
  * [10] Charles Schwartz, ”Toward a quantum theory of tachyon fields,” Int. J. Mod. Phys. A 31, 1650041 (2016).
  * [11] K. R. Parthasarathy: An Introduction to Quantum Stochastic Calculus, Birkhauser, Basel (1992)
  * [12] Y.S. Kim and E.S. Noz: Phase Space Picture of Quantum Mechanics: Group Theoretical Approach, Lecture Notes in Physics Series, World Scientific Pub Co Inc (1991).
  * [13] R. L. Hudson and K. R. Parthasarathy, Unification of boson and fermion quantum stochastic calculus, Commun. Math. Phys. 104 (1986) 457-470.