Article

Hydrodynamic Porosity: A New Perspective on Flow through

Porous Media, Part I
August H. Young 1,2,* and Zbigniew J. Kabala 3,*

1 Duke Center for WaSH‐AID, Duke University, Durham, NC 27701, USA

2 Mechanical Engineering and Materials Science, Duke University, Durham, NC 27710, USA
3 Civil and Environmental Engineering, Duke University, Durham, NC 27710, USA

* Correspondence: ahf12@duke.edu (A.H.Y.); zbigniew.kabala@duke.edu (Z.J.K.)

Abstract: Pore‐scale flow velocity is an essential parameter in determining transport through porous
media, but it is often miscalculated. Researchers use a static porosity value to relate volumetric or
superficial velocities to pore‐scale flow velocities. We know this modeling assumption to be an over‐
simplification. The variable fraction of porosity conducive to flow, what we define as hydrodynamic
porosity, 𝜃 , exhibits a quantifiable dependence on the Reynolds number (i.e., pore‐scale flow
velocity) in the Laminar flow regime. This fact remains largely unacknowledged in the literature. In
this work, we quantify the dependence of 𝜃 on the Reynolds number via numerical flow sim‐
ulation at the pore scale for rectangular pores of various aspect ratios, i.e., for highly idealized dead‐
end pore spaces. We demonstrate that, for the chosen cavity geometries, 𝜃 decreases by as
much as 42% over the Laminar flow regime. Moreover, 𝜃 exhibits an exponential dependence
on the Reynolds number, Re = 𝑅. The fit quality is effectively perfect, with a coefficient of determi‐
nation (R2) of approximately 1 for each set of simulation data. Finally, we show that this exponential
dependence can be easily fitted for pore‐scale flow velocity through use of only a few Picard itera‐
tions, even with an initial guess that is 10 orders of magnitude off. Not only is this relationship a
more accurate definition of pore‐scale flow velocity, but it is also a necessary modeling improve‐
ment that can be easily implemented. In the companion paper (Part 2), we build upon the findings
reported here and demonstrate their applicability to media with other pore geometries: rectangular
and non‐rectangular cavities (circular and triangular).
Citation: Young, A.H.; Kabala, Z.J.
Hydrodynamic Porosity: A New Keywords: hydrodynamic porosity; cavity; dead‐end pore; pore velocity; volumetric velocity;
Perspective on Flow through Reynolds number; groundwater remediation
Porous Media, Part I. Water 2024, 16,
2158. https://doi.org/10.3390/
w16152158

Academic Editor: Constantinos V.

1. Introduction
Chrysikopoulos Fluid flow and transport through porous media are ubiquitous to natural and engi‐
neered systems including groundwater remediation, surface water treatment, and various
Received: 7 May 2024
Revised: 27 July 2024
other industrial processes. Recent studies underscore the importance of considering pore‐
Accepted: 27 July 2024
scale flow dynamics and structures, namely dead‐end pores, in elucidating fluid transport
Published: 30 July 2024 mechanisms within porous media and describing macroscale trends.
In the study of porous media, the pore space is typically broken down into two re‐
gions: mobile and “immobile” zones, as described by van Genuchten and Wierenga [1] in
application to groundwater flows. In the mobile zone, solute transport occurs via advec‐
Copyright: © 2024 by the authors.
Submitted for possible open access
tion and dispersion. The immobile zone is defined by isolated volumes of cavities or dead‐
publication under the terms and
ended pore space adjacently located to well‐connected, mobile regions. In these zones,
conditions of the Creative Commons fluid recirculates in eddies, and solute transport is limited to the mechanism of vortex‐
Attribution (CC BY) license enhanced diffusion. By this definition, the “immobile” label is a misnomer—fluid in the
(https://creativecommons.org/license cavity space is technically mobile; it does not, however, move through the pore space.
s/by/4.0/). Thus, the fluid in this zone remains immobile relative to the flow in the mobile zone. We
illustrate these zones for an arbitrary matrix subject to an imposed flow in Figure 1, below.

Water 2024, 16, 2158. https://doi.org/10.3390/w16152158 www.mdpi.com/journal/water

Water 2024, 16, 2158 2 of 27

Magnification “A” provides an example of a mobile zone, and magnification “B” contains
an example of an immobile zone in the form of a dead‐end pore. The model we use to
study this dead‐end pore volume is pictured in Figure 2.

Figure 1. Mobile zone composed of well‐connected pore space (A), and an immobile zone in the
form of a poorly connected/dead‐end pore (B).

Figure 2. Mobile–immobile zone model of the dead‐end pore space (left), boundary‐driven condi‐
tion, and vortex‐enhanced diffusion (middle), shear‐driven boundary condition yielding a deform‐
able, mobile separatrix (right).

The relevance of the dead‐end pores pictured in Figure 1, though media‐specific, can
be profound, especially in unwashed media such as glacial till or fractured rock. Indeed,
the presence of dead‐end pores has been shown to be abundant in the subsurface. For
example, Lee et al. [2] found the fraction of immobile water content in an undisturbed soil
core to range from 0.42 to 0.82, while Casey et al. [3] measured the average fraction to be
0.62 in the field. Testing in the 1960s revealed the significance of dead‐end pores in reser‐
voir rock: Fatt et al. [4] estimated a total volume of 20% in limestone and shellstone core
samples, whereas Coats and Smith [5] estimated a volume of roughly 10% in sandstone
core samples.
Typically, a substantial immobile zone volume results in solutes favoring preferential
flow paths (i.e., well‐connected pore volumes), bypassing immobile zones and dead‐end
pores. Models neglecting immobile zones and preferential flow paths may underestimate
solute movement by half, as indicated by Jaynes et al. [6], highlighting the significance of
considering dead‐end pores in fluid transport dynamics. Dead‐end pores have also been
implicated in other macroscopic phenomena. Leismann et al. [7] discuss tailing in large‐
scale propagation processes, attributing it to pollutant persistence in immobile zones or
dead ends of the pore space. Lake [8] identifies viscous fingering instabilities and stagnant
areas, or dead‐end pores, as critical limiting factors in miscible displacement processes in
porous media, impacting remediation, CO2 sequestration, and energy extraction.
Water 2024, 16, 2158 3 of 27

Recent studies have delved further into the impact of dead‐end pores on fluid flow
and macroscopic transport phenomena (e.g., see Bordoloi et al. [9]). Gao et al. [10] find
that immobile water in dead‐end pores not only affects solute transport processes but also
crucially regulates breakthrough and tailing in soil columns. Yuan et al. [11] observe that
dead‐end pores impede the efficient removal of non‐aqueous phase liquids from the pore
space due to slow mass transfer rates between the mobile and immobile zones. Nguyen
et al. [12] discuss the effect of dead‐end pores on porous battery electrode performance,
underscoring the importance of understanding their role in optimizing electrode design.
Following an abundance of such evidence, Khuzhayorov et al. [13] conclude that zones of
immobile liquid significantly influence transport in porous media.
Building upon the insights gained from studies on dead‐end pores and pore‐scale
flow dynamics, the concept of effective porosity emerges as a central focus in understand‐
ing fluid transport in porous media. However, the interpretation of effective porosity has
been mired in ambiguity across different disciplines, leading to varying definitions and
applications. Take, for example, the definition of effective porosity assigned by the textile
industry. In the context of hernia meshes, effective porosity is meant to define changes to
the pore morphology after implantation of the mesh in situ (Jacombs et al. [14]). This is
quite different from the definition used by Sevee [15] in a study on the effective porosity
of marine clay. In this study, effective porosity describes the void space in the clay that
participates in advective transport. Still, other definitions describe it as the difference be‐
tween the total porosity minus the soil water content at 0.33 bar (Helalia [16]; Timlin et al.
[17]). Readers are directed to Hapgood et al. [18], Flint and Selker [19], Cartwright et al.
[20], and Ma et al. [21] for additional, alternative definitions. This ambiguity necessitates
a more nuanced approach to characterize the portion of porosity influencing fluid
transport.
In this work, we introduce the concept of hydrodynamic porosity, denoted as 𝜃 ,
to specifically capture flow‐driven fluctuations in the pore volume conducive to through‐
flow. The objective of our analysis is to provide a more precise understanding of the dy‐
namics underlying fluid and contaminant transport within porous media. The concept of
hydrodynamic porosity not only differs from the many definitions of effective porosity,
but also from that of dynamic porosity, which describes the shrinking or swelling of po‐
rous media (Sheng et al. [22]; McDonald et al. [23]), and dynamic effective porosity, which
describes variably saturated porous media around the capillary fringe (Luo et al. [24]). By
delving into the physics of pore‐scale flow, we lay the theoretical groundwork for our
study, aiming to provide a comprehensive analysis of hydrodynamic porosity and its im‐
plications for modeling fluid flow and transport in porous media, with particular empha‐
sis on the role of dead‐end pores.

1.1. Effective Porosity

Among the authors who define effective porosity as the portion of constant porosity
used to transmit fluid through porous media are Li et al. [25], Kabala and Kim [26], Kim
[27], Lindsay [28], and Werth [29]. Here, we note that while these authors qualitatively
consider its hydrodynamic nature, they do not quantify it. For example, Li et al. [25] qual‐
itatively note the dependence of effective porosity on flow velocity in a study of sedimen‐
tary rock flows. The authors attribute the observed decrease to the presence of selective
pathways (i.e., fissures and cracks) in the studied rock formation. They do not discuss an
explicit relationship between what they term effective porosity and flow velocity. Further,
the morphology of the rock formation studied in this work is starkly different than that of
granular media.
In his doctoral dissertation at Duke and a publication with his advisor, Kim [27] and
Kabala and Kim [26] provide the most thorough discussion on the hydrodynamic quality
of effective porosity; the authors state a dependence of effective porosity on both pore
geometry and Reynolds number, Re = 𝑅. They demonstrate this dependence in the same
idealized pore space we study in this work using the FIDAP software. The authors show
Water 2024, 16, 2158 4 of 27

that effective porosity varies by at least an order of magnitude for creeping flows and
postulate that it may “vary by as much” for non‐creeping flows but do not develop an
explicit relationship between effective porosity and Reynolds number (Kabala and Kim
[26]).
Other authors have qualitatively considered the impact of dead‐end pores on effec‐
tive porosity. Lindsay [28] notes that in water‐saturated paper, flow can be restricted by
mechanical obstructions in the form of isolated and dead‐end pores. The author also dis‐
cusses the development of stagnant zones in high Reynolds number flows, and the ab‐
sence of these zones in creeping flows due to the ability of the flow to closely follow abrupt
changes in medium geometry. In effect, the author contemplates the dependence of effec‐
tive porosity on the Reynolds number, Re = 𝑅, but does not account for it in his analysis.
The concept of flow‐dependent porosity also appears in a study on baleen, though in a
much different context. Werth [29] describes a linear relationship between fringe porosity
and incident flow velocity until a certain limiting velocity is reached, at which point the
effective porosity of the baleen decreases. Given that the application of this work is to
learn more about the mechanical properties of baleen to better understand the feeding
behaviors of various whale species, the spatial scale and flow path morphology are fun‐
damentally different than those studied in this work.

1.2. Cavity Flows

As we discuss in the next section, flow past cavities (also known as dead‐end or blind
pores), which serve as the theoretical basis for this research, are a popular research topic
and comparatively well explored. There are ample publications concerning contaminant
transport and flow manipulation and instabilities in dead‐end pores. The bulk of this work
pertains to the geometric manipulation of pore geometry to determine the effect on the
induced flow (Moffatt [30]; Coats and Smith [5]; Higdon [31]; Shen and Floryan [32]; Fang
et al. [33]).
Three review articles cover the field well (Meier et al. [34]; Worthington [35]; Yan et
al. [36]), but there are many other notable works, e.g., Li et al. [37], Yuan and Rezaee [38],
Foroughi et al. [39], Verbovšek [40], Fenni et al. [41], Kango et al. [42], and Yao et al. [43].
Studies on flow modulation have also been published, though they are relatively less
explored (Jana and Ottino [44]; Howes and Shardlow [45]; Horner et al. [46]; Kahler and
Kabala [47]). Most of this work is geared at industry for the guise of expediting rate‐lim‐
ited manufacturing processes (e.g., etching, finishing, cleaning, etc.) that are applied to
surfaces with cavities (Chilukuri and Middleman [48]; Alkire et al. [49]; Fang et al. [50]).
In this work, we leverage the same physical theory that describes flows past cavities
to define a hydrodynamic porosity function 𝜃 𝑣 𝜃 𝑅 . Previously, Kahler
and Kabala [47,51,52] used the same approach to describe contaminant transport in po‐
rous media—likening duct or surface flow to flow through a series of well‐connected
pores, and flow past grooves and cavities on surfaces to flow past dead‐end pores in gran‐
ular media. Such results are crucial to understanding how phenomena like contaminant
rebound in groundwater reservoirs post‐remediation can be mitigated. Although Kahler
and Kabala implied the hydrodynamic nature of porosity, they did not quantify it. In this
paper, we demonstrate and quantify, explicitly, this relationship for the first time. We also
illustrate the ease with which this relationship can be incorporated into flow and contam‐
inant transport models. To understand why contaminant rebound after traditional pump‐
and‐treat groundwater remediation takes place, and how it could be mitigated, research‐
ers need to account for hydrodynamic porosity, 𝜃 .

2. Pore Scale Flow Modeling

Here, we are solving the continuity equation (mass conservation) and the Navier–
Stokes equation (the evolution of momentum) for incompressible Newtonian fluid (water)
subject to no‐slip and no‐penetration boundary conditions on the walls of the channel and
Water 2024, 16, 2158 5 of 27

cavity. The pressure is specified at the outlet. The imposed inlet velocity distribution is
parabolic, as in a fully developed channel flow (analogous to Hagen–Poiseuille flow).
𝛻⋅𝑢 0

𝜕𝑢
𝜌 𝑢 ⋅ 𝛻𝑢 𝛻𝑝 𝜇𝛻 𝑢
𝜕𝑡

𝑢 𝑥, 𝑦 𝑤𝑎𝑙𝑙 0

𝑝 𝑥 𝑐ℎ𝑎𝑛𝑛𝑒𝑙 𝑙𝑒𝑛𝑔𝑡ℎ, 𝑦 0
where u is a 2D velocity vector, ρ is density, t is time, p is pressure, and μ is dynamic
viscosity. The Navier–Stokes equation is solved for steady state, i.e., the time derivative
of the velocity is zero. Since we do not refer to these equations later, we list them unnum‐
bered.
We assume a (quasi) steady‐state, laminar flow in our model. It is typical for ambient
groundwater conditions, which may vary on a time scale of weeks, as well as for pump‐
and‐treat remediation conducted predominantly with a constant/fixed pumping rate. The
ambient flow is usually characterized by Re = 𝑅 < 1, whereas forced‐gradient flows typi‐
cally do not exceed Re = 𝑅 < 100 even near a pumping well and thus are also laminar. For
our numerical experiments, we chose Reynolds numbers that span this range: 1, 10, 50,
and 100; they help us demonstrate the variability of hydrodynamic porosity. We further
discuss our selection of Reynolds numbers in Section 3 Methods.

2.1. Total Porosity

The total porosity of a medium, 𝜃, is defined by the cumulative volume of the mobile
and immobile zones; and more specifically, as the ratio of the total pore space volume to
the total volume of the media, 𝑉. In an isotropic or 2D medium, like the ones we study in
this work, this definition can be written in terms of cross‐sectional areas:
𝑉 𝐴
𝜃
𝑉 𝐴
We note that 𝑉 𝑉 𝑉 and 𝐴 𝐴 . In the analysis that follows, we rep‐
resent the pore volume (or area) by the dead‐end pore model, which we illustrate in Figure
2. As a result, the pore volume can be defined in terms of mobile and immobile zones.
Again, we provide the cross‐sectional area expression for an isotropic medium or 2D me‐
dia.
𝑉 𝑉 𝑉 →𝐴 𝐴 𝐴 (2)
Thus, we can define the total porosity of the medium, 𝜃, as a function of mobile and
immobile zone porosities, 𝜃 and 𝜃 , respectively, by dividing Equation (2) by
the total volume of the medium, 𝑉 (or cross‐sectional area, 𝐴):
𝑉 𝑉
𝜃 𝜃 𝜃 (3)
𝑉

2.2. Pore‐Scale Flow Velocity

An expedient (and reductive) way to define the pore‐scale flow velocity is to use the
volumetric velocity, also known as the superficial or Darcy/Forchheimer velocity. This
quantity is equivalent to the flow velocity of a fictitious fluid flowing through an entire
cross section of the medium rather than just through the void space conducive to flow:
𝑄
𝑣 𝑣 𝑞 (4)
𝐴
Water 2024, 16, 2158 6 of 27

where 𝑄 is the volumetric flow rate, 𝐴, is the cross‐sectional area of the medium, and 𝑞,
is the “flux.” Neglecting inertial effects, Darcy’s law relates the volumetric velocity to the
pressure gradient applied to the medium (Brutsaert [53]; Muljadi et al. [54]; Bear [55]):
𝜇 𝜇 1
𝛻𝑝 𝑣 ⇔ 𝛻ℎ 𝑣 ⇔ 𝛻ℎ 𝑣 (5)
𝑘 𝑘𝛾 𝐾
where 𝑝 is the pressure, ℎ is the pressure head, 𝑘 is the permeability of the medium, 𝐾
is the hydraulic conductivity, 𝛾 is the specific weight of the fluid, and 𝜇 is its viscosity.
The equivalences we show here are to illustrate the preferred forms in oil and gas reser‐
voir modeling (left) and groundwater hydrology (right). When inertial effects cannot be
neglected, as is the case for high Reynolds number flows, we must utilize the quadratic
correction term, introduced by Èstudes [56] and Forchheimer [57]. Given this adjustment,
Darcy’s law becomes the Forchheimer–Dupuit law:
𝜇 1 1
𝛻𝑝 𝑣 𝐵𝜌𝑣 𝑛 ⇔ 𝛻ℎ 𝑣 𝐵 𝑣 𝑛 (6)
𝑘 𝐾 𝑔
where 𝑛 is a unit vector in the direction of the volumetric velocity, 𝜌 the flow density,
and 𝐵, a coefficient that can be found experimentally (Chen et al. [58]). Depending on the
flow conditions, the volumetric velocity can be estimated experimentally from Equation
(5) or Equation (6) (i.e., from Darcy’s law or the Forchheimer–Dupuit law, respectively).
To attain the true pore‐scale flow velocity, which would be needed to determine quantities
such as contaminant transport time, the volumetric velocity must be modified by the me‐
dium’s porosity (Bear [55]):
𝑣
𝑣 (7)
𝜃
Back‐of‐the‐envelope calculations may simply use the total porosity of the medium.
If the mobile zone porosity is known, then 𝜃 would instead be used, but as a static
quantity. This would also be suitable for washed media without any cavities or other ef‐
fectively immobile zones. Equation (7) can be derived from a simple conservation of mass
analysis:
𝐴 1
𝑣 𝐴 𝑣𝐴 ⇔ 𝑣 𝑣 ⇔ 𝑣 𝑣
𝐴 𝜃
Given the fact that immobile zones do not contribute to through‐flow, we know that
use of the medium’s total porosity is an oversimplification. Pore‐scale flow velocity
should instead be defined by the total volume that is conducive to flow—a quantity that
is itself dependent on pore‐scale flow velocity. Equation (7) should instead read
𝑣 𝑣
𝑣 (8)
𝜃 𝑣 𝜃 𝑅

where, rather than being modeled as static quantity, as is typically the case, 𝜃 is
itself a function of pore‐scale flow velocity, and thus, we proceed with the following no‐
tation: 𝜃 𝑅 . The implicit nature of Equation (8), while seemingly more difficult to
solve than Equation (7), is quickly resolved by a few Picard iterations. Not only is Equa‐
tion (8) a more accurate description of pore‐scale flow velocity, but it is also a necessary
improvement in the modeling of induced subsurface flows that can be easily imple‐
mented.

2.3. Defining the Pore‐Space Partitioning Coefficient

Given that the pore space can be broken down into mobile and immobile zones, we
can define a pore‐space partitioning coefficient, 𝜉, to describe the ratio of pore space condu‐
cive to through‐flow:
Water 2024, 16, 2158 7 of 27

𝑉 𝐴
𝜉 (9)
𝑉 𝐴
We note that the pore‐space partitioning coefficient is related to the hydrodynamic
porosity of the medium, 𝜃 𝑅 , by
𝑉 𝑉 𝜃 𝑅
𝜉 𝜃 𝑅 (10)
𝑉 𝑉 𝜃
In the analysis that follows, it is the behavior of the pore‐space partitioning coefficient
that we numerically quantify as a function of the Reynolds number. We are able to use
our results to describe the hydrodynamic porosity, 𝜃 𝑅 , of the medium because of
the direct proportionality between these two quantities.

2.4. The Dead‐End Pore Model

To geometrically simplify the pore space of a porous medium, we study a single cav‐
ity or dead‐end pore, as pictured in Figure 1. The idea of a poorly connected, or dead‐end
pore was first explored by Turner [59], who studied channel flow with distributed pockets
of stagnant fluid. Although Turner admitted such pore spaces would play a role in diffu‐
sion throughout the pore space, other researchers such as Fatt [60] and Goodknight et al.
[61] initially regarded dead‐end pore spaces as regions through which diffusion could not
occur. Deans [62] noted that the division of pore space into flowing and stagnant regions,
separated by a “resistance to mass transfer” is an “extreme limit” that can only be justified
on the grounds of simplicity. Following this conclusion, Coats and Smith [5] relaxed the
definition of the dead‐end pore to account for diffusion, but still referred to the dead‐end
pore volume as stagnant. Physically, we know this enforcement to be an oversimplifica‐
tion of the recirculatory flow within the dead‐end pore space. Chilukuri and Middleman
[48] corrected for this oversimplification by describing mass transport from dead‐end
pores as a result of vortex‐enhanced diffusion—a conclusion that coincides with a series
of publications that detail the vortex structures within dead‐ended pores (Moffatt [30];
Mehta and Lavan [63]; O’ Brien [64]; Shen and Floryan [32]; Kang and Chang [65]; Fang et
al. [33]). The evolution of the dead‐end pore model is illustrated in Figure 2, below.
Separation of the mobile and immobile zone volumes is described by the boundary‐
or shear‐driven flow models; we refer to the boundary between these zones as the sepa‐
ratrix. In application to flow, the idea of a separatrix was first postured by Elderkin [66],
who describes the boundary as a trajectory that is topologically abnormal in comparison
to nearby trajectories. Weiss [67] later used this concept to describe the defining limit be‐
tween free and trapped fluid regions. Other publications refer to this boundary as a di‐
viding (Moffatt [30]; O’ Brien [64]; Higdon [31]) or separating streamline (Shen and
Floryan [32]; Alkire et al. [49]). The modeling and experimental work on which we build,
i.e., Kahler and Kabala [47,51,52], and even earlier publications such as Horner et al. [46],
use this same terminology (the separatrix) to describe the fluidic boundary between the
mobile and immobile zones in the idealized dead‐end pore space.
In the case of the boundary‐driven model (which is essentially the commonly studied
driven‐lid problem), the geometric boundary between the mobile and immobile zones
also serves as the fixed location of the separatrix. With an increase in Reynolds number of
the adjacent through‐channel flow, the vortex structures within a cavity translate and
smear in the direction of the imposed boundary condition movement. Such results have
been illustrated by many and summarized by Shankar and Deshpande [68]. In shear‐
driven flows, the separatrix is free to move about the cavity space. As discussed by Fang
et al. [50] and Kahler and Kabala [47], the exact location of the separatrix depends on the
Reynolds number of the adjacent through‐channel flow in the mobile zone. This means
that, unlike in the case of the boundary‐driven flow condition, the mobile and immobile
zones cannot be defined based on the geometry of the pore space alone. Instead, the vol‐
umes of these zones must be defined as flow dependent. In application to square cavity
Water 2024, 16, 2158 8 of 27

flow, a comparison between the enforcement of the boundary‐driven and shear‐driven

flow conditions is provided in Figure 3 below.

Figure 3. Vortex location for the mobile and immobile separatrix (left and right, respectively). The
location of the mobile separatrix is determined from a shear‐driven flow condition. The location of
the immobile separatrix is determined from a boundary‐driven flow condition. These results are
generated for the idealized pore space provided in Figure 4 for Re = 𝑅 = 10, see Section 3.2.2.

In the study of flow past cavities, it is standard practice to enforce the shear‐driven
boundary condition. Researchers who initially studied these flows, e.g., Moffatt [30], im‐
mediately identified through‐flow penetration into the cavity space upon investigation.
O’Neill [69] and Wakiya [70] found that attachment of the separatrix to the cavity wall
occurs at some depth into the cavity and as Higdon [31] states, not at the sharp, leading
edge of the cavity. In fact, for rectangular cavities exceeding a given depth ratio, research‐
ers found that the downstream attachment of the separatrix, or dividing streamline, oc‐
curs at the bottom of the cavity wall, mimicking the behavior of a sudden‐expansion flow
(Shen and Floryan [32]; Alkire et al. [49]). Enforcement of the boundary‐driven flow model
would, in this application, yield significant error.
With these findings considered, we again refer to the discussion presented by Li et
al. [25] wherein the effective porosity of the studied medium is found to be flow depend‐
ent. The authors of this study explain that this dependence is a result of physical macro‐
structures that act as preferential flow paths during high flow volumes. What the authors
do not discuss, is the hydrodynamic behavior of the immobile zones within the medium—
a behavior that is driven by separatrix movement in to and out of each effectively dead‐
end pore. As flow volumes increase, the separatrix moves toward its neighboring through
channel, and the immobile zone it defines grows. The result of an increase in flow volume
is a decrease in hydrodynamic porosity, 𝜃 𝑅 . It is only when the shear‐driven flow
condition is applied to the dead‐end pore that this behavior is observed. If the mobile and
immobile zones are improperly defined by the simplified boundary‐driven flow condi‐
tion, this behavior is missed.

3. Methods
To determine the hydrodynamic porosity, 𝜃 𝑅 , of a porous medium, we study
the medium at the pore scale and assume an idealized dead‐end pore geometry. The
movement of the separatrix is tracked over a range of interstitial Reynolds numbers to
determine the relative magnitudes of the mobile and immobile zones, which we then use
to calculate the value of the pore‐space partitioning coefficient, 𝜉 𝜃 𝑅 ⁄𝜃.
Water 2024, 16, 2158 9 of 27

3.1. Numerical Flow Solver

To observe movement of the separatrix in the idealized pore space, we use Mathemat‐
ica’s numerical differential equation solver, NDSolve, to solve the continuity equation
(mass conservation), Navier–Stokes equations (momentum evolution), and associated
boundary conditions. The solver domain is a replica of the idealized pore space utilized
by Kahler and Kabala [47] and is similar in geometry to the domain commonly used in
the study of flow past cavities (Chilukuri and Middleman [48]; Higdon [31]; Fang et al.
[33]). For the sake of simplicity, the flow is modeled as being two‐dimensional. The height
of the through channel is the same as the depth and width of the dead‐end pore (i.e., the
dead‐end pore has a depth ratio of 1:1). In the study of flow past cavities, this geometry is
by far the most prevalent, as noted by Shankar and Deshpande [68]. The through channel
is extended past the dead‐ended pore by twice the channel height to eliminate any end
effects associated with the outflow boundary condition. Finally, the dead‐ended pore is
located one‐fourth of the way into the through‐channel given the need to input a fully
developed flow profile at the through‐channel inlet. To exaggerate the movement of the
separatrix as a function of the Reynolds number, the solver domain is manipulated such
that the through channel becomes much narrower than the depth and width of the cavity
space.
The idealized flow geometry is discretized through use of the ElementMesh function,
which, by default, generates a second‐order, triangular element mesh. The interior and
boundary mesh elements are further refined by specifying upper limits on the Max‐
CellMeasure and MaxBoundaryCellMeasure. A brief convergence analysis of the interior and
boundary mesh cell sizes is provided in the Supplemental Materials. Further, a refinement
region is specified at the geometric boundary of the channel–cavity interface to ensure
proper resolution of the separatrix.
The solver itself is defined by the system of equations that describe steady‐state flow
through the idealized pore space (i.e., the incompressible form of the continuity and Na‐
vier–Stokes equations), as well as the boundary conditions that are assigned to the solver
domain. These equations are normalized by the channel height, ℎ, and the average inlet
flow velocity, 𝑈. The flow is assumed to be a steady state and restricted to the laminar
flow regime. A set of Dirichlet conditions are applied to the boundaries of the solver do‐
main (i.e., the no‐slip condition at the domain walls and a uniform pressure condition,
wherein the pressure is arbitrarily set to zero, at the domain outlet). The inlet velocity
profile is defined by the Hagen–Poiseuille model for fully developed channel flow. Fi‐
nally, flow is assigned to the entire idealized pore space given that the application of this
work is to fully saturated porous media. The properties of water at standard conditions
are assigned to the fluid.
We note that the non‐dimensional form of the Navier–Stokes equations is used in this
analysis. The scaling on the pressure term is appropriate for flows that are dominated by
convective action (i.e., flows in which viscous effects are relatively negligible). This choice
was made to replicate the scaling utilized by Kahler and Kabala [47], who studied flows
with channel‐based Reynolds numbers of Re = 𝑅 = 0.01–10. Fang et al. [33] used the same
scaling, though for admittedly higher channel‐based Reynolds numbers in the range of
50–1600. For our case, and the work on which we build, this scaling is justified because
these studies aim to capture flow phenomena driven by convective action (i.e., changes in
momentum to the bulk flow). When the Reynolds number approaches 0, and the flow can
be approximated as creeping, the scaling on the pressure term can be achieved through
use of the flow viscosity.
For this system of equations, NDSolve utilizes the Finite Element Method (Wolfram
Finite Element Method [71]) to arrive at a solution. In general, the solver method is auto‐
matically determined by the results of symbolic analysis. The use of the Finite Element
Method is triggered by specific user inputs. For example, specification of the Navier–
Stokes equation using the “Inactive” operator or boundary conditions defined by the Di‐
richletCondition function prompt the use of this solver method. Implementation of this
Water 2024, 16, 2158 10 of 27

solver method can be verified by validating that the solution contains an ElementMesh
(Wolfram Symbolic and Numerical Computation [72]). The outputs of the solver are three
interpolating functions that describe the pressure and velocity fields within the flow do‐
main. Streamlines are visualized through the use of the built‐in StreamPlot command. The
error associated with each solver output is determined by the solver mesh (i.e., domain
discretization), assigned working precision, and solver method.
The above‐mentioned functions are very well documented in the Wolfram Language
(Mathematica) Documentation [72,73] and in Wolfram Numerical Solutions of PDEs [74].
Furthermore, in Supplemental Materials, we provide the code used in this study. The
main errors are numerical, but our convergence analysis, also included in Supplemental
Materials, demonstrates that they are negligible. The adopted 2D flow simplifying as‐
sumption is consistent with flow through fractured‐rock medium with wide fractures, and
thus is not severely limiting—moreover, it is invoked by numerous studies in the litera‐
ture (see further discussion at the end of Section 6).

3.2. Data Collection

In this work, we vary the Reynolds number of the through‐channel flow (i.e., the flow
in the mobile zone), and the depth ratio of the idealized dead‐end channel–cavity geome‐
try. As noted by Fang [75] and others, the location of the separatrix is a function of both
Reynolds number and geometry (Mehta and Lavan [63]; O’ Brien [64]; Higdon [31]; Kim
[27]).

3.2.1. Reynolds Numbers

Below a Reynolds number of 1, the separatrix remains stationary (Kahler and Kabala
[47]). Flows of this nature are classified by the creeping flow regime, where viscous effects
dominate. It is not until we study Reynolds numbers within the inertial flow regime that
we observe a mobile separatrix. This is because the location of the separatrix is dictated
by the inertia of the adjacent through‐flow. For this reason, we impose the following Reyn‐
olds number, Re, ranges to the through‐flow:
 Re = 𝑅 = 0.01–1
(to verify the stationary nature of the separatrix in the creeping flow regime)
 Re = 𝑅 = 1–100
(to illustrate the mobility of the separatrix in the laminar flow regime)
To remain within the laminar flow regime, we limit our Reynolds number to a max‐
imum of 100. This choice is admittedly arbitrary, given that transition to turbulence in
pipe flow typically occurs over a diameter‐based Reynolds number of 2000 and at least
one order of magnitude above the particle‐based Reynolds number at which the deviation
from Darcy’s Law occurs in porous media (Bear [55]); deviation from Darcy’s Law gener‐
ally occurs between a particle‐based Reynolds number of 1 and 10.
In terms of particle‐based Reynolds numbers, there is ample evidence that the onset
of transitionary behavior occurs around 100. For columns of packed spheres, Jolls and
Hanratty [76] report the onset of transitionary behavior within the range of 110–150.
Wegner et al. [77] found a slightly lower range of 90–120 for beds of packed spheres. Latifi
et al. [78] encountered transitionary behavior at 110, also for a bed of packed spheres;
however, the authors did note unsteady laminar flow behavior until 370. A similar study
conducted by Rode et al. [79] reports transitionary behavior in the range of 110–150. Fi‐
nally, Bu et al. [80] define a critical particle‐based Reynolds number of 100 as the cutoff
for laminar flow, with the onset of turbulence occurring between 230–400.
If we instead consider the interstitial Reynolds number, which is based on average
pore size and average pore‐scale flow velocity, we encounter the commonly cited Reyn‐
olds number ranges provided by Dybbs and Edwards [81], summarized in Table 1 below.
Water 2024, 16, 2158 11 of 27

Table 1. Reynolds number ranges corresponding to pre‐turbulent flow regimes, as provided by

Dybbs and Edwards [81].

Flow Regime Reynolds Number Range

Creeping/Darcy 0–1
Inertial 1–10
Laminar, non‐linear 10–150
Laminar, unsteady 150–300

We can easily imagine representing our system in terms of the interstitial Reynolds
number. Given our assumption that the medium is homogeneous, we know the average
pore size. Pore‐scale flow velocity is typically calculated by dividing the flux through the
medium by the porosity of the medium, but in our case, we will impose it directly by
assigning a mean velocity to the through‐channel flow. Capping our Reynolds number at
a value of 100 keeps our analysis within the steady laminar flow regime.

3.2.2. Flow Geometries

The effect of pore geometry manipulation has been extensively studied in the litera‐
ture. In these studies, the authors vary the type of cavity (i.e., rectangular, circular, etc.),
the depth ratio of the cavity, and the size of the cavity relative to the size of the adjacent
through channel. To replicate the results obtained by Kahler and Kabala [47], we use an
idealized geometry wherein the height of the through channel is equivalent to the depth
and width of the cavity geometry, as provided in the top‐left quadrant of Figure 4. To
exaggerate the mobility of the separatrix in the laminar flow regime, the depth and width
of the cavity, relative to the through channel, are equally increased in magnitude. These
geometries are provided in Figure 4, below, and referred to by channel–cavity depth ratios
(i.e., 1:1, 3:4, 1:2, and 1:4). For example, the depth ratio 1:2 corresponds to the geometry in
which the cavity depth and width are twice that of the through‐channel height.

Figure 4. Study geometries 1:1 3:4, 1:2, and 1:4; specified by their depth ratio (i.e., channel height to
cavity depth/width) in non‐dimensional space.

We note that in this paper (Part 1), we consider rectangular pores of various aspect
ratios, i.e., a highly idealized case of the dead‐ended pore spaces. In the companion paper
(Part 2), we build upon our findings from this paper and demonstrate their applicability
to media with other geometries: rectangular and non‐rectangular cavities (i.e., circular,
and triangular).
Water 2024, 16, 2158 12 of 27

3.3. Measurement Method

To determine the value of the pore‐space partitioning coefficient, 𝜉 𝜃 𝑅 ⁄𝜃 ,
and therefore the hydrodynamic porosity, 𝜃 𝑅 , for each inlet flow condition, we
use the following procedure:
1. Generate a monochromatic stream plot of the flow.
2. Draw the separatrix in the area between the bulk flow in the through channel and the
recirculatory flow in the dead‐end pore space in a contrasting color, using the stream‐
lines in the stream plot for guidance.
3. Use an interpolating function to mathematically describe the location of the separa‐
trix.
4. Define the area below the separatrix as the immobile zone, and the area above as the
mobile zone and use numerical integration to quantify the magnitude of these re‐
gions.
To draw the separatrix, we use a DynamicModule in Mathematica to generate an inter‐
polating function that includes five points (or more) of our choosing between the mobile
and immobile zones. This process is pictured below in Figure 5. To determine the sizes of
the mobile and immobile zones, we apply numerical integration to the resulting interpo‐
lating function.

Figure 5. A DynamicModule automatically places five points on a monochromatic stream plot in the
vicinity of the separatrix(left). Given user input (i.e., movement of these five points to the approxi‐
mate location of the separatrix), the DynamicModule produces an interpolating function that can be
used to describe the location of the separatrix (right).

Given that the shape and location of this interpolating function are a direct result of
user input, there is an inherent error built into the measurement process that we are una‐
ble to precisely quantify. Additional errors in the measurement process result from the
chosen resolution of the stream plots which is in turn limited by the quality of the solver
mesh and working precision assigned to the numerical solver method.

4. Results
4.1. Separatrix Movement
Because the cavity flow is driven by the adjacent through‐channel flow, we start by
providing a stream plot of the entire dead‐end pore geometry adjacent to the correspond‐
ing cavity flow in Figure 6. We then provide stream plots for each cavity geometry at
Reynolds numbers of 1, 10, 50, and 100 in Figure 7.
Water 2024, 16, 2158 13 of 27

Figure 6. Stream plot magnification example.

Figure 7. Dead‐end pore stream plots for each depth ratio (1:1, 3:4, 1:2, and 1:4 from top to bottom)
for Reynolds numbers (1, 10, 50, and 100 from left to right).

With the flow in the flow‐through channel (shown in Figure 6) from the left to the
right, note in Figure 7 clear inertial effects: with the increasing Reynolds number the center
of the cavity vortex moves to the right and rises up.
The corresponding movement of the separatrix as a function of Reynolds number
was first explored by Kahler and Kabala [47]. In their work, the authors track the bottom‐
Water 2024, 16, 2158 14 of 27

most point of the separatrix to determine its maximum penetration depth into the dead‐
end pore. To confirm the numerical accuracy of the results produced in this work, we
replicate this plot, which is provided in the Supplemental Material. We expand upon this
plot by tracking the movement of the separatrix for three additional flow geometries, as
pictured below.
The results obtained by Kahler and Kabala [47] reveal an immobile separatrix in the
creeping flow regime (i.e., Re < 1), and a mobile separatrix in the inertial flow regime. In
the former, the bottom‐most point of the separatrix does not exceed 25% of the depth of
the idealized pore space. Our results replicate the separatrix behavior observed by Kahler
and Kabala [47] and are within, at most, a 2% difference. Over the range of Reynolds num‐
bers plotted in Figure 8, the maximum penetration depth of the separatrix diminishes by
20%.

Figure 8. Maximum relative penetration depth of the separatrix into the dead‐end pore space as a
function of Reynolds number, Re.

Building upon the test conditions utilized by Kahler and Kabala [47], we observe
movement of the separatrix toward the geometric boundary of the cavity space for Reyn‐
olds numbers approaching 100. Movement of the separatrix within the cavity space is fur‐
ther exaggerated by manipulation of the flow geometry, as illustrated in Figure 9 below.

Figure 9. Immobile zone vortex streamlines, corresponding to Re = 1, bounded by the separatrix

(highlighted in yellow) for each channel–cavity depth ratio, 1:1, 3:4, 1:2, and 1:4 (from left to right).
Additional separatrix locations for Reynolds numbers corresponding to 10, 50, and 100 are plotted
in orange, purple, and blue, respectively.
Water 2024, 16, 2158 15 of 27

Movement of the separatrix toward the adjacent through channel (and out of the cav‐
ity space) results in a decrease in the pore‐space partitioning coefficient, 𝜉
𝜃 𝑅 ⁄𝜃, and therefore, the hydrodynamic porosity, 𝜃 ⬚ ,of the medium. For
example, media with cavities described by the dead‐end pore model pictured in the top‐
left quadrant of Figure 4 experience roughly a 4% reduction in 𝜃 over the tested
Reynolds number range. By comparison, when the channel–cavity depth ratio is 1:4 (i.e.,
the cavity depth is four times that of the through‐channel height), 𝜃 𝑅 decreases
by approximately 42%. See Table 2 for a summary of the changes associated with each
channel–cavity depth ratio; note that given Equation (10), the provided percent‐decrease
values are the same for the pore‐space partitioning coefficient, 𝜉, and the hydrodynamic
porosity, 𝜃 𝑅 , of the medium.

Table 2. Maximum and minimum pore‐space partitioning coefficient, 𝜉 𝜃 𝑅 ⁄𝜃, values cor‐
responding to an increase in Reynolds number in each tested flow geometry.

Pore‐Space Partitioning Coefficient, 𝝃 𝜽𝒎𝒐𝒃𝒊𝒍𝒆 /𝜽

Depth Ratio Max. Min. % Decrease
1:1 0.84 0.8 4.34
3:4 0.81 0.75 7.47
1:2 0.79 0.67 15.81
1:4 0.87 0.51 41.50

4.2. Exponential Dependence of Hydrodynamic Porosity on Pore‐Scale Flow Velocity

When plotted, the pore‐space partitioning coefficient, 𝜉 𝜃 𝑅 ⁄𝜃, and there‐
fore the hydrodynamic porosity, 𝜃 𝑅 , of the medium, approaches the value sug‐
gested by the boundary‐driven model we previously discussed for Reynolds numbers ap‐
proaching the upper limit of the laminar regime (Re = 𝑅 = 100). See Figure 10 below for
evidence of this behavior. In this figure, we also demonstrate that these quantities exhibit
a linear dependence on Reynolds number, in the inertial flow regime (Re = 𝑅 = 1–10).
Extrapolation of this relationship past Re = 𝑅 = 10 is provided to illustrate the error that
would result from not utilizing the exponential relationship provided in Equation (12).
Again, we remind readers of the direct proportionality between the partitioning coeffi‐
cient, 𝜉 𝜃 𝑅 ⁄𝜃, and the hydrodynamic porosity of the medium, 𝜃 .

Linear and Exponential Model Fits

Pore- Space Partitioning Coefficient

0.84
ξ = θ mobile/θ

Re < 10
0.82
Re > 10

0.8

0 20 40 60 80 100
Reynolds Number (Re)
Water 2024, 16, 2158 16 of 27

Figure 10. The pore‐space partitioning coefficient, 𝜉 𝜃 𝑅 ⁄𝜃 , exhibits a linear dependence

on Reynolds number, Re = 𝑅, for Re = 1–10 and an exponential dependence for Re = 𝑅 = 1–100. At
Reynolds numbers approaching 100, the partitioning coefficient of the flow geometry approaches
the value predicted by the boundary‐driven flow model. The results pictured here are for the 1:1
channel–cavity depth ratio.

The fit parameters for each channel–cavity depth ratio are provided in Table 3. The
quality of each fit is measured by the coefficient of determination (R2), which is approxi‐
mately 1 for each tested depth ratio. In the inertial flow regime, the pore‐space partitioning
coefficient, 𝜉 𝜃 𝑅 ⁄𝜃, exhibits a linear dependence on Reynolds number. How‐
ever, as the Reynolds number increases, the partitioning coefficient deviates from this lin‐
ear model. Instead, a nonlinear dependence on Reynolds number explains the calculated
partitioning coefficient values at higher Reynolds numbers, and more significantly, over
the entire range of Reynolds numbers, Re, in the laminar regime (Re = 1–100). This de‐
pendence is well fitted by an exponential function:
𝑐ℎ𝑒𝑖𝑔ℎ𝑡
𝜉 𝑎 𝑏𝑒 ⟺𝜉 𝑎 𝑏𝑒 , 𝑤ℎ𝑒𝑟𝑒𝑑
𝜈
This expression can be easily re‐written using the help of Equation (10) to define the
hydrodynamic porosity of the medium, 𝜃 :
𝜃 𝑣 𝑎 𝑏𝑒 𝜃⟺𝜃 𝑅 𝑎 𝑏𝑒 𝜃 (12)
where a is the pore‐space partitioning coefficient, 𝜉 𝜃 𝑅 ⁄𝜃, value approximated
by the boundary‐driven model (i.e., Re → ∞), and the quantity ‘a + b’ is the value in the
creeping flow regime (i.e., Re → 0). We note that these values, as written, are the same for
Equations (11) and (12).

Table 3. Hydrodynamic porosity, 𝜃 𝑅 , parameters for Reynolds number dependence (d) and
pore‐flow velocity (c) defined in Equations (11) and (12).

Equation (12) Exponential Fit Parameters

Depth Ratio 𝑹𝟐 a b c (s/m) d
1:1 1.0000 0.80 3.80 x 10 25.91 3.19 10
3:4 1.0000 0.75 6.20 x 10 33.62 4.14 10
1:2 1.0000 0.67 1.26 x 10 49.79 6.13 10
1:4 0.9993 0.50 3.66 x 10 77.32 9.52 10

Mathematically, the fit parameter, b, drives the exponential behavior of our fit. When
on the order of magnitude of 10 the exponential behavior of our fit becomes most ex‐
aggerated. This is exemplified in Figure 11, below, for depth ratios 1:4 and 1:2.
Water 2024, 16, 2158 17 of 27

Normalized Hydrodynamic Porosity (θmobile /θMIM ):

Function of Reynolds Number

1.6 Depth
Ratio
θ mobile/θ MIM

1:1
1.4
3:4
1:2
1.2
1:4

1
0 20 40 60 80 100
Reynolds Number (Re)

Figure 11. The exponential decay of hydrodynamic porosity, 𝜃 𝑅 , as a function of Reynolds

number, Re, for media with cavities of varying channel–cavity depth ratios. For ease of comparison,
𝜃 𝑅 is normalized by the value that corresponds to the mobile–immobile zone model, 𝜃 .

In Figure 11, we normalize the hydrodynamic porosity, 𝜃 𝑅 , by the static mo‐

bile‐zone porosity value, 𝜃 , that results from enforcement of the mobile–immobile
zone model (i.e., the boundary‐driven flow condition at the cavity–channel interface) in
the dead‐end pore space, defined as
𝜃 𝜉 𝜃 (13)
The geometric pore‐space partitioning coefficient, 𝜉 , is determined by the relative
magnitudes of the through‐channel and cavity volumes of each dead‐end pore geometry.
For the 1:1 depth ratio pictured in Figure 4, 𝜉 = 4/5. For the 1:2 depth ratio pictured in
Figure 4, however, 𝜉 = 8/12. We define these quantities to compare the relative change
in hydrodynamic porosity across the four geometries provided in Figure 4.

5. Calculating Pore‐Scale Flow Velocity and Hydrodynamic Porosity Parameters

5.1. Pore‐Scale Flow Velocity
We can use this newfound exponential relationship outlined in Equation (12) to fill
in the details of Equation (8):
𝑣
𝑣 (14)
𝑎 𝑏𝑒 𝜃
Here, we remind readers that 𝑣 is the volumetric velocity used in Darcy’s or Forch‐
heimer–Dupuit’s law (Equations (5) and (6), respectively). To illustrate the ease at which
the exponential relationship in Equation (12) can be incorporated into current models, we
provide a brief example utilizing the idealized pore geometry provided in the top left
quadrant of Figure 4, the associated exponential fit coefficients (i.e., 𝑎 0.80 , 𝑏
3.80 10 , and 𝑐 25.91 𝑠⁄𝑚), and the properties of water assumed by Kahler and Ka‐
bala [47]. In this example, we assume the volumetric velocity has a magnitude 10 (m/s)
and a total porosity (𝜃) of 0.4—an arbitrary total porosity value in the range typically
measured for unconsolidated, unwashed media (i.e., 20–45%) (Woessner and Poeter [82]).
Water 2024, 16, 2158 18 of 27

Using only 3 Picard iterations with an initial guess of 2 10 (m/s), we converge

to 3 decimal places and a pore‐scale flow velocity of 2.99 10 (m/s). Below, we illus‐
trate the power of Mathematica’s built‐in Nest function, which will, for the same parame‐
ter set, converge for an initial guess spanning twenty orders of magnitude. In the code
provided below, we use the sister, NestList function to quickly determine the number of
iterations needed for convergence.
Incidentally, why only 3 iterations? This is not surprising once one realizes that our
exponential function f(x) fulfills the Lifschitz condition: |xn+1 − xn| < |f’(xn)| |xn − xn−1| or
|f’(xn)| < 1 that is necessary for the mapping, f(x), to be a contraction and thus have a fixed
point—see the Banach Fixed‐Point Theorem [83].

In Table 4 below, we repeat this process for an initial guess 10 orders of magnitude
larger and smaller than our original guess to illustrate the robust nature of this method.

Table 4. The pore‐scale flow velocity (m/s) can be determined in only a few Picard iterations.

Pore‐Scale Velocity (m/s)

Iteration Number
Initial Guess 1 2 3 4
2 10 0.00299016 0.00299344 0.00299345 0.00299345
2 10 0.00298329 0.00299342 0.00299345 0.00299345
2 10 0.003125 0.00299388 0.00299346 0.00299345

From this demonstration, we conclude that implementation of Equation (12) is no

more mathematically burdensome to implement than the constant porosity model, which
results in unnecessary error.

5.2. Hydrodynamic Porosity Fit Parameters

In this work, we illustrate how to solve for the exponential fit parameters in Equa‐
tions (11) and (12) through numerical simulation. But this is an idealized case—porous
media have an array of randomly distributed and sized pore spaces, not to mention drastic
changes in pore geometry. If we wanted to use the relationship provided in Equation (12),
we would need to quantify the volumetric and pore‐scale flow velocities to then be able
to determine the corresponding exponential fit parameters. This could be achieved
through column experiments with a known pressure gradient (to solve for the volumetric
flow velocity), and a measurable tracer pulse (to solve for the pore‐scale flow velocity).
Given that Equation (12) contains three unknown parameters, we would need to conduct
this experiment at least three times, each at a different, judiciously selected flow condition.
We note in passing that measuring the pore‐scale flow velocity is not trivial and will re‐
quire care as demonstrated or implied by Wood et al. [84], Berkowitz and Scher [85],
Haggerty and Gorelick [86], Medina and Carrera [87], Carrera et al. [88], Bolster et al. [89],
Water 2024, 16, 2158 19 of 27

and others. Since this issue is not central to the focus of this work, we will address it in the
follow‐up paper.
Here we offer only a numerical example and work with an artificial dataset, gener‐
ated from use of Equation (14) and the exponential fit parameters we used in our last
example (i.e., 𝑎 0.80, 𝑏 3.80 10 , and 𝑐 25.91 𝑠⁄𝑚). We assume we have meas‐
ured volumetric velocities that correspond to the laminar flow regime (Re = 1–100); for the
characteristic height, we assign to the dead‐end pore geometry (i.e., 0.001 m), the corre‐
sponding pore‐scale velocity range is approximately 0.001–0.1 (m/s). We assume we meas‐
ure the volumetric flow velocities 0.01, 0.05, and 0.09 (m/s). Given this assumption, we can
calculate the corresponding pore‐scale flow velocities:
Moving forward, we assume that we have measured these values experimentally and
that we actually do not know the values of our exponential fit parameters, a, b, and c. We
can use the synthetic dataset provided in Table 5 to calculate them. This is easily achieved
through the use of Mathematica’s FindFit function. Specifying a Newton solution method,
we arrive at the anticipated parameter values, exactly in 100 iterations (the default value
used by the FindFit function).

Table 5. Artificial volumetric and pore‐scale flow velocity (m/s) data set, generated from Equation
(13).
𝟐 𝟐
Volumetric Velocity (m/s) 𝟏𝟎 Pore‐Scale Flow Velocity (m/s) 𝟏𝟎
1 3.05939
5 15.6122
9 28.1241

Here we note that we have utilized the arbitrary precision assigned by Mathematica to
generate a set of pore‐scale velocity values. Realistically, the precision of our measure‐
ments would be restricted by our measurement device. Thus, we re‐run the above calcu‐
lation, but this time with only three significant digits rather than the six that we previously
used. This time, we arrive at a, b, and c values of 0.801, 0.0366, and 26.4, respectively. The
error in our estimated values is 0.09%, 3.65%, and 1.73%. If we increase the precision of
our measurements to four significant digits, the error in our estimated values reduces to
0.14%, 0.10%, and 0.25%, respectively. Clearly, accurate estimation of a, b, and c will re‐
quire numerical fine‐tuning and surplus velocity data. This process will need to be con‐
ducted via numerical and experimental column tests, which we are currently pursuing.
Water 2024, 16, 2158 20 of 27

6. Discussion
At scale, the implications of these results are easily observed. For example, consider
a periodic medium well approximated by the dead‐end pore model with a measurable
total porosity of 0.4. If we assume the non‐dimensional idealized pore space pictured in
the top left quadrant of Figure 4, we know that the entire area of the through channel
contributes to the mobile porosity of the medium. The immobile zone is occupied by the
vortex, which resides within the dead‐end pore. For a Reynolds number of 1, the pore‐
space partitioning coefficient (𝜉) of the dead‐end pore geometry is approximately 0.84.
Given this value, we can calculate the hydrodynamic porosity of the medium using Equa‐
tion (10):
𝜃 𝜉𝜃
As a result, the medium has a hydrodynamic porosity of approximately 0.34. If we
had used the boundary‐driven flow condition to determine 𝜃 𝑅 (provided in Fig‐
ures 2 and 3), we would have under‐approximated the value at 0.32 by roughly 5%; simi‐
larly, we would have over‐approximated the immobile porosity by roughly 20%. For the
most exaggerated channel–cavity depth ratio we test (1:4), use of the boundary‐driven
model would have resulted in a 42% error in 𝜃 . This is due to the deep impingement
of the through‐flow into the dead‐end pore for this configuration. We summarize these
calculations in Table 6 below:

Table 6. Error in the mobile–immobile zone model for a periodic medium with a total porosity of
0.4 and cavity geometry pictured in the top left quadrant of Figure 4.

Pore‐Space Partitioning Mobile Zone Po‐ Immobile Zone Po‐

Coefficient, rosity, rosity,
𝝃 𝜽𝒎𝒐𝒃𝒊𝒍𝒆 ⁄𝜽 𝜽𝒎𝒐𝒃𝒊𝒍𝒆 𝜽𝒊𝒎𝒎𝒐𝒃𝒊𝒍𝒆
Hydrodynamic po‐
0.84 0.336 0.064
rosity model
Mobile–immobile
0.8 0.32 0.08
zone model
Mobile–immobile zone model error (%) 5% 20%

Comparing the shear‐driven and boundary‐driven models, we see that the error in
the latter (in terms of the pore‐space partitioning coefficient) is largest for flows approach‐
ing the creeping flow regime and does not become less than 1 until a Reynolds number of
roughly 49 for the channel–cavity depth ratio of 1:1. For the exaggerated flow geometries
(i.e., 3:4, 1:2, and 1:4 depth ratios), the Reynolds number must exceed 51, 48, and 45, re‐
spectively. Referring to Figure 3, we see a 17% error associated with the use of the bound‐
ary‐driven model.
We note that the exponential relation in (12) along with its excellent fit quality is not
surprising, as its analogs show up in several areas of mathematical physics. Exponential
forms play an important role in the solution of differential equations and are common in
groundwater flow modeling—see the low‐velocity non‐Darcy flow model in shale and
tight reservoirs (e.g., Wang and Sheng [90]) and the various renditions of the well function
(e.g., Theis [91], Hantush [92]); also, see their applications in the percolation theory (e.g.,
Hunt et al. [93]); Gardner’s equation (Gardner [94]) relating hydraulic conductivity and
matric potential in flows through unsaturated porous media (vadose zone); and even the
relationship between soil water content and electrical resistivity (Pozdnyakov et al. [95]).
Furthermore, in nature, an exponential decline in hydraulic conductivity with depth is
considered a hallmark of catchment hydrology (Ameli et al. [96]).
As already mentioned in Section 3.2.2 Flow Geometries, in this paper (Part 1) we con‐
sider rectangular pores of various aspect ratios, i.e., highly idealized dead‐end pore
spaces. Our main finding is that hydrodynamic porosity can be described by an
Water 2024, 16, 2158 21 of 27

exponential function of pore‐scale flow velocity (or interstitial Reynolds number). In the
companion paper (Part 2), we build upon this finding and demonstrate its applicability to
media with other geometries: rectangular and non‐rectangular cavities (i.e., circular, and
triangular). Furthermore, we show that not only does the exponential relationship hold
for media with a variety of cavity geometries, but it does so almost perfectly with a coef‐
ficient of determination (R2) of approximately 1 for each new set of simulation data. In this
way, the research results have, hopefully, practical value in applications where unwashed
porous media play an important role.
The fit parameters to our exponential relationship are likely to be affected slightly by
numerical errors in the least‐square fitting algorithm, which can be minimized by appro‐
priately choosing tolerance parameters. However, the main source of errors in their esti‐
mation is likely to be anomalous (non‐Fickian) tracer transport in column experiments.
How often such transport will be encountered remains to be seen—we will investigate it
in column lab experiments in a follow‐up paper.
Naturally, the next steps for this work are an application to groundwater flow models
at the macro‐scale and laboratory column experiments. The ability to quantify the expo‐
nential fit parameters for sites needing remediation, as previously discussed, is also nec‐
essary to demonstrate the ease with which this relationship can be tailored to any given
media. We are currently pursuing these research avenues.
As mentioned earlier, we have verified our numerical results by conducting conver‐
gence analysis, i.e., refining the numerical grids used by the Mathematica superfunction
NDSolve—see the Supplemental Materials. Our numerical errors related to finite‐element
discretization do not exceed 2% and are likely significantly smaller. Given that our analy‐
sis exactly aligns with the analytical solution for channel flow (i.e., the Hagen–Poiseuille
flow profile), we are confident that our numerical results are excellent approximations of
the true (but unknown) solutions of the Navier–Stokes and continuity equations involved
in our boundary‐value problems.
We note that several articles and a dissertation on flow over cavities present results
analogous to ours. Although none of them introduce the concept of hydrodynamic porosity,
their results confirm ours qualitatively or even quantitatively (for specific Reynolds num‐
bers and specific cavity geometries). These papers are mostly limited to a short range of
Reynolds numbers. They include studies of flow in tubes with circumferential cavities and
in 2D channels with cavities. The numerical studies confirming our results include Take‐
matsu [97]; Friedman [98]; Stevenson [99]; O’Brien, [64]; Higdon [31]; Driesen et al. [100];
and Young et al. [101]. The experimental studies confirming our results include Taneda
[102]; Shankar and Deshpande [68]; and particularly the dissertation Laskowska [103].
There are also mixed studies (numerical and experimental) confirming our results: Shen
and Floryan [32]; and Pan and Acrivos [104].

7. Study Limitations
The location of the separatrix, the streamline separating the flow in the flow‐through
channel from the vortex in the dead‐end cavity, is determined by an approximate, yet
relatively accurate, hybrid graphical‐computational method illustrated in Figure 5. From
the stream plot, the separatrix location, and ultimately the hydrodynamic porosity,
𝜃 𝑅 , is determined by means of an interpolation function between 5 (or more)
points placed (graphically) on the smooth separatrix. The area below the separatrix (and
above) is then calculated by analytic integration. The repeatability of this measurement
process is reasonably good—for the idealized pore space pictured in Figure 2 (i.e., a chan‐
nel–cavity depth ratio of 1:1), and a Reynolds number of 0.01, the maximum penetration
depth of the separatrix was found to vary by 0.005 across five measurements. Relative to
an average penetration depth of 0.25, this variation accounts for less than 2% of the mag‐
nitude of the mean.
Our algorithm for determining the separatrix location is robust and is presented in
our highly legible code available in the Supplemental Materials. Its errors are smaller than
Water 2024, 16, 2158 22 of 27

the numerical errors of solving the PDEs. An alternative method could be based on image
processing, but with our robust algorithm, there is no need for it.
Additionally, we consider only the immobile zones generated by dead‐end pores and
cavity‐like structures that generate flow separation. We do not consider the immobile
zones that result from molecular forces between the media and the flowing solution, nor
do we consider surface tension effects. However, we assume that these immobile zones,
just as the ones studied in this work, would behave similarly so long as the Reynolds
number of the interstitial flow is appropriate, i.e., up to about 100.
Finally, this work focuses on groundwater flow, i.e., flow through fully saturated
porous media with dead‐end pores (unwashed sand and gravel, fractured rock, etc.). Con‐
sequently, our results do not apply to multiphase flow such as water and air in unsatu‐
rated (vadose) zone. Also, our results are derived for Newtonian, incompressible fluids
and thus do not apply to non‐Newtonian, rheological flows.

8. Conclusions
In this work, we explore a physical phenomenon that has been largely neglected in
the literature. Generally, the porosity of a medium that is conducive to through‐flow, what
we define as hydrodynamic porosity, 𝜃 , is still thought to be a static parameter. Re‐
searchers define this quantity using the cumulative volume of the flow through channels
(i.e., mobile zones); cavities and other effectively immobile zones are not considered to
significantly contribute to the through‐flow volume. We find this approximation to be a
notable oversimplification in the creeping and inertial flow regimes, resulting in the mis‐
representation of porous media flows and associated transport processes.
Although a few researchers have previously acknowledged the dependence of what
they analogously define as effective porosity on the flow velocity in rock formations, these
researchers did not quantify this relationship. In this work, we demonstrate and quantify,
explicitly, the dependence of hydrodynamic porosity, 𝜃 , on pore‐scale flow velocity,
for the first time. We begin with a direct replication of the results provided by Kahler and
Kabala [47], wherein the boundaries between the mobile and immobile zones of porous
media are shown to be hydrodynamic in nature. We then study the movement of this
boundary over a range of interstitial, or channel‐based, Reynolds numbers in the laminar
flow regime in porous media. The movement of this boundary defines the pore‐space par‐
titioning coefficient, 𝜉 𝜃 𝑅 ⁄𝜃 (i.e., the fraction of the pore space conducive to
through‐flow), and therefore, the hydrodynamic porosity, 𝜃 𝑅 , of the medium.
Given the direct proportionality between these quantities, we find both to have an expo‐
nential dependence on Reynolds number or pore‐scale flow velocity, which we provide
in Equations (11) and (12). Finally, we show that this dependence can be easily incorpo‐
rated into porous media flow modeling using only a few Picard iterations, even with an
initial guess that is over 10 orders of magnitude off.
The flow‐dependent nature of hydrodynamic porosity, 𝜃 𝑅 , plays an unmis‐
takable role in the transport of contaminants in porous media. Moving forward, our re‐
search opens avenues for further exploration into the impact of dead‐end pores on fluid
flow dynamics and macroscopic transport phenomena. By integrating insights from re‐
cent studies and leveraging advances in modeling techniques, we can advance our under‐
standing of porous media behavior and contribute to the development of more accurate
predictive models for practical applications in groundwater remediation, surface water
treatment, and various industrial processes.
The exponential relationship and our results reported here can find broad applica‐
tions: from groundwater flow to flow‐through filters (for example, our work: Young et al.
[101]) and flows for cleaning grooved surfaces, such as in the production of computer
motherboards (e.g., Fang et al. [33,50]; Fang [75]).
Finally, our segmenting of the material into two papers has a guiding logic: Part 1
(this paper) presents the hydrodynamic porosity as an exponential function of pore‐scale
flow velocity (or interstitial Reynolds number) for a pore with a square cavity—a highly
Water 2024, 16, 2158 23 of 27

idealized case of the dead‐ended pore spaces in a porous medium. In contrast, Part 2 (com‐
panion paper: Young and Kabala [105]) demonstrates the applicability of the exponential
relationship to media with other cavity geometries: rectangular and non‐rectangular (e.g.,
circular and triangular). It appears this relationship is widely applicable to porous/frac‐
tured media.

Supplementary Materials: The following supporting information can be downloaded at:

https://www.mdpi.com/article/10.3390/w16152158/s1, The supporting materials include Codes:
MDPI_Water_Paper1_DataFile.nb and MDPI_Water_Paper1_SolverModules.wl and convergence
analysis: MDPI_Water_MeshConvergenceAnalysis.nb and MDPI_Water_MeshConvergenceAnaly‐
sis.pdf.
Author Contributions: Conceptualization (Z.J.K.), data curation (A.H.Y.), formal analysis (A.H.Y.),
investigation (A.H.Y. and. Z.J.K.), methodology (A.H.Y. and. Z.J.K.), project administration (Z.J.K.),
supervision (Z.J.K.), validation (A.H.Y.), visualization (A.H.Y.), writing—original draft (A.H.Y. and.
Z.J.K.), writing—review and editing (A.H.Y. and. Z.J.K.). All authors have read and agreed to the
published version of the manuscript.
Funding: This work was, in part, supported by the Bill & Melinda Gates Foundation through Duke
University’s Center for WaSH‐AID by a grant, OPP1173370 (Sanitation Technology Cluster). And
the APC was funded by the same source. All opinions, findings, and conclusions or recommenda‐
tions expressed in these works are those of the author(s) and do not necessarily reflect the views of
the Foundation, Duke, or the Center.
Data Availability Statement: The simulation data that support the findings of this study, the corre‐
sponding Mathematica/Wolfram language code files (notebooks), and convergence analysis note‐
book are available in Supplemental Materials.
Acknowledgments: We acknowledge the academic editor, assistant editor, and three anonymous
reviewers for their careful reading of our initial manuscript and for offering constructive and in‐
sightful critique. Their comments and suggestions allowed us to significantly improve the manu‐
script.
Conflicts of Interest: The authors declare no competing conflicts of interest.

References
1. van Genuchten, M.T.; Wierenga, P.J. Mass‐Transfer Studies in Sorbing Porous‐Media. I. Analytical Solutions. Soil Sci. Soc. Am.
J. 1976, 40, 473–480. https://doi.org/10.2136/sssaj1976.03615995004000040011x.
2. Lee, J.; Horton, R.; Noborio, K.; Jaynes, D.B. Characterization of Preferential Flow in Undisturbed, Structured Soil Columns
Using a Vertical TDR Probe. J. Contam. Hydrol. 2001, 51, 131–144. https://doi.org/10.1016/s0169‐7722(01)00131‐0.
3. Casey, F.X.M.; Horton, R.; Logsdon, S.D.; Jaynes, D.B. Immobile Water Content and Mass Exchange Coefficient of a Field Soil.
Soil Sci. Soc. Am. J. 1997, 61, 1030–1036. https://doi.org/10.2136/sssaj1997.03615995006100040006x.
4. Fatt, I.; Maleki, M.; Upadhyay, R.N. Detection and Estimation of Dead‐End Pore Volume in Reservoir Rock by Conventional
Laboratory Tests. Soc. Petrol. Eng. J. 1966, 6, 206–212. https://doi.org/10.2118/1441‐pa.
5. Coats, K.H.; Smith, B.D. Dead‐End Pore Volume and Dispersion in Porous Media. Soc. Petrol. Eng. J. 1964, 4, 73–84.
https://doi.org/10.2118/647‐pa.
6. Jaynes, D.B.; Logsdon, S.D.; Horton, R. Field Method for Measuring Mobile/Immobile Water Content and Solute Transfer Rate
Coefficient. Soil Sci. Soc. Am. J. 1995, 59, 352–356. https://doi.org/10.2136/sssaj1995.03615995005900020012x.
7. Leismann, H.M.; Herding, B.; Krenn, V. A Quick Algorithm for the Dead‐End Pore Concept for Modeling Large‐Scale Propa‐
gation Processes in Groundwater. In Developments in Water Science; Celia, M.A., Ferrand, L.A., Brebbia, C.A., Gray, W.G., Pinder,
G.F., Eds.; Elsevier: Amsterdam, The Netherlands, 1988; pp. 275–280. https://doi.org/10.1016/S0167‐5648(08)70101‐1.
8. Lake, L.W. Enhanced Oil Recovery; Prentice Hall: Englewood Cliffs, NJ, USA, 1989.
9. Bordoloi, A.D.; Scheidweiler, D.; Dentz, M.; Bouabdellaoui, M.; Abbarchi, M.; de Anna, P. Structure Induced Laminar Vortices
Control Anomalous Dispersion in Porous Media. Nat. Commun. 2022, 13, 3820. https://doi.org/10.1038/s41467‐022‐31552‐5.
10. Gao, G.; Feng, S.; Zhan, H.; Huang, G.; Mao, X. Evaluation of Anomalous Solute Transport in a Large Heterogeneous Soil Col‐
umn with Mobile‐Immobile Model. J. Hydrol. Eng. 2009, 14, 966–974. https://doi.org/10.1061/(ASCE)HE.1943‐5584.0000071.
11. Yuan, Q.; Ma, Z.; Wang, J.; Zhou, X. Influences of Dead‐End Pores in Porous Media on Viscous Fingering Instabilities and
Cleanup of NAPLs in Miscible Displacements. Water Resour. Res. 2021, 57, e2021WR030594.
https://doi.org/10.1029/2021WR03059.
Water 2024, 16, 2158 24 of 27

12. Nguyen, T.‐T.; Demortière, A.; Fleutot, B.; Delobel, B.; Delacourt, C.; Cooper, S.J. The Electrode Tortuosity Factor: Why the
Conventional Tortuosity Factor Is Not Well Suited for Quantifying Transport in Porous Li‐Ion Battery Electrodes and What to
Use Instead. npj Comput. Mater. 2020, 6, 123. https://doi.org/10.1038/s41524‐020‐00386‐4.
13. Khuzhayorov, B.K.; Makhmudov, Z.M.; Zikiryaev, S.K. Substance Transfer in a Porous Medium Saturated with Mobile and
Immobile Liquids. J. Eng. Phys. Thermophys. 2010, 83, 263–270. https://doi.org/10.1007/s10891‐010‐0341‐3.
14. Jacombs, A.S.W.; Karatassas, A.; Klosterhalfen, B.; Richter, K.; Patiniott, P.; Hensman, C. Biofilms and Effective Porosity of
Hernia Mesh: Are They Silent Assassins? Hernia 2020, 24, 197–204. https://doi.org/10.1007/s10029‐019‐02063‐y.
15. Sevee, J.E. Effective Porosity Measurement of a Marine Clay. J. Environ. Eng. 2010, 136, 674–681.
https://doi.org/10.1061/(ASCE)EE.1943‐7870.0000205.
16. Helalia, A.M. The Relation between Soil Infiltration and Effective Porosity in Different Soils. Agric. Water Manag. 1993, 24, 39–
47. https://doi.org/10.1016/0378‐3774(93)90060‐N.
17. Timlin, D.J.; Ahuja, L.R.; Pachepsky, Y.; Williams, R.D.; Gimenez, D.; Rawls, W. Use of Brooks‐Corey Parameters to Improve
Estimates of Saturated Conductivity from Effective Porosity. Soil Sci. Soc. Am. J. 1999, 63, 1086–1092.
https://doi.org/10.2136/sssaj1999.6351086x.
18. Hapgood, K.P.; Litster, J.D.; Biggs, S.R.; Howes, T. Drop Penetration into Porous Powder Beds. J. Colloid Interface Sci. 2002, 253,
353–366. https://doi.org/10.1006/jcis.2002.8527.
19. Flint, L.E.; Selker, J.S. Use of Porosity to Estimate Hydraulic Properties of Volcanic Tuffs. Adv. Water Resour. 2003, 26, 561–571.
https://doi.org/10.1016/S0309‐1708(02)00182‐3.
20. Cartwright, N.; Nielsen, P.; Perrochet, P. Behavior of a Shallow Water Table under Periodic Flow Conditions. Water Resour. Res.
2009, 45, W05408. https://doi.org/10.1029/2008WR007306.
21. Ma, C.; Shi, W.; Zhan, H. On the Vertical Circulation Wells in a Leaky‐Confined Aquifer. J. Hydrol. 2022, 608, 127676.
https://doi.org/10.1016/j.jhydrol.2022.127676.
22. Sheng, G.; Javadpour, F.; Su, Y. Dynamic Porosity and Apparent Permeability in Porous Organic Matter of Shale Gas Reservoirs.
Fuel 2019, 251, 341–351. https://doi.org/10.1016/j.fuel.2019.04.044.
23. McDonald, P.J.; Istok, O.; Janota, M.; Gajewicz‐Jaromin, A.M.; Faux, D.A. Sorption, Anomalous Water Transport and Dynamic
Porosity in Cement Paste: A Spatially Localised 1H NMR Relaxation Study and a Proposed Mechanism. Cem. Concr. Res. 2020,
133, 106045. https://doi.org/10.1016/j.cemconres.2020.106045.
24. Luo, Z.; Kong, J.; Yao, L.; Lu, C.; Li, L.; Barry, D.A. Dynamic Effective Porosity Explains Laboratory Experiments on Water Table
Fluctuations in Coastal Unconfined Aquifers. Adv. Water Resour. 2023, 171, 104354. https://doi.org/10.1016/j.advwa‐
tres.2022.104354.
25. Li, H.; Ohtsuka, Y.; Mori, N.; Inagaki, T.; Misawa, S. Effective Porosity and Specific Yield of a Sedimentary Rock Determined
by a Field Tracing Test Using Tritium as a Tracer. Environ. Geol. 1996, 27, 170–177.
26. Kabala, Z.J.; Kim, Y.‐W. Dynamic Effective Porosity: Numerical Simulations. J. Res. Inst. Eng. Technol. 2011, 30, 91–95.
27. Kim, Y.‐W. Pore‐Scale Flow and Contaminant Transport in Porous Media; Ph.D. Dissertation, Duke University, Durham, NC, USA,
2006.
28. Lindsay, J.D. Relative Flow Porosity in Fibrous Media‐Measurements and Analysis, Including Dispersion Effects. Tappi J. 1994,
77, 225–239.
29. Werth, A.J. Flow‐Dependent Porosity and Other Biomechanical Properties of Mysticete Baleen. J. Exp. Biol. 2013, 216, 1152–1159.
https://doi.org/10.1242/jeb.078931.
30. Moffatt, H.K. Viscous and Resistive Eddies near a Sharp Corner. J. Fluid Mech. 1963, 18, 1–18.
https://doi.org/10.1017/s0022112064000015.
31. Higdon, J.J.L. Stokes Flow in Arbitrary Two‐Dimensional Domains: Shear Flow over Ridges and Cavities. J. Fluid Mech. 1985,
159, 195–226. https://doi.org/10.1017/s0022112085003172.
32. Shen, C.; Floryan, J.M. Low Reynolds Number Flow over Cavities. Phys. Fluids 1985, 28, 3191–3202.
https://doi.org/10.1063/1.865366.
33. Fang, L.C.; Cleaver, J.W.; Nicolaou, D. Hydrodynamic Cleansing of Cavities. In Proceedings of the 8th International Conference
on Computational Methods and Experimental Measurements (CMEM 97), Rhodes, Greece, 21–23 May 1997; pp. 391–401.
34. Meier, H.; Zimmerhackl, E.; Zeitler, G. Modeling of colloid‐associated radionuclide transport in porous groundwater aquifers
at the Gorleben site, Germany. Geochem. J. 2003, 37, 325–350. https://doi.org/10.2343/geochemj.37.325.
35. Worthington, S.R. Estimating effective porosity in bedrock aquifers. Groundwater 2022, 60, 169–179.
https://doi.org/10.1111/gwat.13171.
36. Yan, S.; Yang, M.; Sun, C.; Xu, S. Liquid Water Characteristics in the Compressed Gradient Porosity Gas Diffusion Layer of
Proton Exchange Membrane Fuel Cells Using the Lattice Boltzmann Method. Energies 2023, 16, 6010.
https://doi.org/10.3390/en16166010.
37. Li, B.; Li, X.‐S.; Li, G.; Jia, J.‐L.; Feng, J.‐C. Measurements of Water Permeability in Unconsolidated Porous Media with Methane
Hydrate Formation. Energies 2013, 6, 3622–3636. https://doi.org/10.3390/en6073622.
38. Yuan, Y.; Rezaee, R. Comparative Porosity and Pore Structure Assessment in Shales: Measurement Techniques, Influencing
Factors and Implications for Reservoir Characterization. Energies 2019, 12, 2094. https://doi.org/10.3390/en12112094.
Water 2024, 16, 2158 25 of 27

39. Foroughi, S.; Bijeljic, B.; Gao, Y.; Blunt, M.J. Incorporation of Sub‐Resolution Porosity into Two‐Phase Flow Models with a Mul‐
tiscale Pore Network for Complex Microporous Rocks. Water Resour. Res. 2024, 60, e2023WR036393.
https://doi.org/10.1029/2023WR036393.
40. Verbovšek, T. Variability of Double‐Porosity Flow, Interporosity Flow Coefficient λ and Storage Ratio ω in Dolomites. Water
2024, 16, 1072. https://doi.org/10.3390/w16081072.
41. Fenni, M.; Guellal, M.; Hamimid, S. Influence of Porosity Properties on Natural Convection Heat Transfer in Porous Square
Cavity. Phys. Fluids 2024, 36, 056108. https://doi.org/10.1063/5.0206797.
42. Kango, R.; Chandel, A.; Shankar, V. A Statistical Model for Estimating Porosity Based on Various Parameters of Flow Through
Porous Media. Water Pract. Technol. 2024, 19, 1936–1947. https://doi.org/10.2166/wpt.2024.114.
43. Yao, P.; Zhang, J.; Qin, Z.; Fan, A.; Feng, G.; Vandeginste, V.; Zhang, P.; Zhang, X. Effect of Pore Structure Heterogeneity of
Sandstone Reservoirs on Porosity‐Permeability Variation by Using Single‐Multi‐Fractal Models. ACS Omega 2024, 9, 23339–
23354. https://doi.org/10.1021/acsomega.3c09957.
44. Jana, S.C.; Ottino, J.M. Chaos‐Enhanced Transport in Cellular Flows. Philos. Trans. R. Soc. A 1992, 338, 519–532.
https://doi.org/10.1098/rsta.1992.0018.
45. Howes, T.; Shardlow, P.J. Simulation of Mixing in Unsteady Flow through a Periodically Obstructed Channel. Chem. Eng. Sci.
1997, 52, 1215–1225. https://doi.org/10.1016/s0009‐2509(96)00361‐2.
46. Horner, M.; Metcalfe, G.U.; Wiggins, S.; Ottino, J.M. Transport Enhancement Mechanisms in Open Cavities. J. Fluid Mech. 2002,
452, 199–229. https://doi.org/10.1017/s0022112001006917.
47. Kahler, D.M.; Kabala, Z.J. Acceleration of Groundwater Remediation by Deep Sweeps and Vortex Ejections Induced by Rapidly
Pulsed Pumping. Water Resour. Res. 2016, 52, 3930–3940. https://doi.org/10.1002/2015wr017157.
48. Chilukuri, R.; Middleman, S. Cleaning of a Rough Rigid Surface‐Removal of a Dissolved Contaminant by Convection‐Enhanced
Diffusion and Chemical‐Reaction. J. Electrochem. Soc. 1984, 131, 1169–1173. https://doi.org/10.1149/1.2115772.
49. Alkire, R.C.; Deligianni, H.; Ju, J.B. Effect of Fluid‐Flow on Convective‐Transport in Small Cavities. J. Electrochem. Soc. 1990, 137,
818–824. https://doi.org/10.1149/1.2086562.
50. Fang, L.C.; Nicolaou, D.; Cleaver, J.W. Transient Removal of a Contaminated Fluid from a Cavity. Int. J. Heat Fluid Flow 1999,
20, 605–613. https://doi.org/10.1016/s0142‐727x(99)00050‐8.
51. Kahler, D.M.; Kabala, Z.J. Rapidly Pulsed Pumping Accelerates Remediation in a Vertical Circulation Well Model. Water 2018,
10, 1423. https://doi.org/10.3390/w10101423.
52. Kahler, D.M.; Kabala, Z.J. Acceleration of Groundwater Remediation by Rapidly Pulsed Pumping: Laboratory Column Tests. J.
Environ. Eng. 2019, 145, 06018009. https://doi.org/10.1061/(asce)ee.1943‐7870.0001479.
53. Brutsaert, W. Hydrology: An Introduction; Cambridge University Press: Cambridge, UK, 2005.
54. Muljadi, B.P.; Blunt, M.J.; Raeini, A.Q.; Bijeljic, B. The Impact of Porous Media Heterogeneity on Non‐Darcy Flow Behaviour
from Pore‐Scale Simulation. Adv. Water Resour. 2016, 95, 329–340. https://doi.org/10.1016/j.advwatres.2015.05.019.
55. Bear, J. Dynamics of Fluids in Porous Media. Soil Sci. 1975, 120, 162–163.
56. Èstudes, D.J. Thèoriques Et Pratiques Sur Le Mouvement Des Eaux; Dunod: Paris, France, 1863.
57. Forchheimer, P. Wasserbewegung Durch Boden. Zeit. Ver. Deut. Ing. 1901, 45, 1781–1788.
58. Chen, Y.‐F.; Zhou, J.‐Q.; Hu, S.‐H.; Hu, R.; Zhou, C.‐B. Evaluation of Forchheimer Equation Coefficients for Non‐Darcy Flow in
Deformable Rough‐Walled Fractures. J. Hydrol. 2015, 529, 993–1006. https://doi.org/10.1016/j.jhydrol.2015.09.021.
59. Turner, G.A. The Flow‐Structure in Packed Beds—A Theoretical Investigation Utilizing Frequency Response. Chem. Eng. Sci.
1958, 7, 156–165. https://doi.org/10.1016/0009‐2509(58)80022‐6.
60. Fatt, I. Pore Structure of Sintered Glass from Diffusion and Resistance Measurements. J. Phys. Chem. 1959, 63, 751–752.
https://doi.org/10.1021/j150575a031.
61. Goodknight, R.C.; Klikoff, W.A.; Fatt, I. Non‐Steady‐State Fluid Flow and Diffusion in Porous Media Containing Dead‐End
Pore Volume. J. Phys. Chem. 1960, 64, 1162–1168. https://doi.org/10.1021/j100838a014.
62. Deans, H.A. A Mathematical Model for Dispersion in the Direction of Flow in Porous Media. Soc. Petrol. Eng. J. 1963, 3, 49–52.
https://doi.org/10.2118/493‐pa.
63. Mehta, U.B.; Lavan, Z. Flow in a Two‐Dimensional Channel with a Rectangular Cavity. J. Appl. Mech. 1969, 36, 897.
https://doi.org/10.1115/1.3564799.
64. O’Brien, V. Closed Streamlines Associated with Channel Flow over a Cavity. Phys. Fluids 1972, 15, 2089–2097.
https://doi.org/10.1063/1.1693840.
65. Kang, I.S.; Chang, H.N. The Effect of Turbulence Promoters on Mass‐Transfer‐Numerical‐Analysis and Flow Visualization. Int.
J. Heat Mass Transf. 1982, 25, 1167–1181. https://doi.org/10.1016/0017‐9310(82)90211‐3.
66. Elderkin, R.H. Separatrix Structure for Elliptic Flows. Am. J. Math. 1975, 97, 221–247. https://doi.org/10.2307/2373669.
67. Weiss, J.B. Transport and Mixing in Traveling Waves. Phys. Fluids A 1991, 3, 1379–1384. https://doi.org/10.1063/1.858068.
68. Shankar, P.N.; Deshpande, M.D. Fluid Mechanics in the Driven Cavity. Annu. Rev. Fluid Mech. 2000, 32, 93–136.
https://doi.org/10.1146/annurev.fluid.32.1.93.
69. O’Neill, M.E. Separation of a Slow Linear Shear‐Flow from a Cylindrical Ridge or Trough in a Plane. Z. Angew. Math. Phys. 1977,
28, 439–448. https://doi.org/10.1007/bf01601625.
70. Wakiya, S. Application of Bipolar Coordinates to 2‐Dimensional Creeping Motion of a Liquid. III. Separation in Stokes Flows.
J. Phys. Soc. Jpn. 1975, 45, 1756–1763. https://doi.org/10.1143/jpsj.45.1756.
Water 2024, 16, 2158 26 of 27

71. Wolfram Finite Element User Guide. Available online: https://reference.wolfram.com/language/FEMDocumentation/tuto‐

rial/FiniteElementOverview.html (accessed on 1 June 2023).
72. Wolfram Symbolic and Numerical Computation. Available online: https://reference.wolfram.com/language/#Symboli‐
cAndNumericComputation (accessed on 1 December 2022).
73. Wolfram Language (Mathematica) Documentation. Available online: https://reference.wolfram.com/language/ (accessed on 15
June 2024).
74. Wolfram Numerical Solutions of Partial Differential Equations. Available online: https://reference.wolfram.com/language/tu‐
torial/NDSolvePDE.html (accessed on 24 June 2024).
75. Fang, L.‐C. Effect of Mixed Convection on Transient Hydrodynamic Removal of a Contaminant from a Cavity. Int. J. Heat Mass
Transf. 2003, 46, 2039–2049. https://doi.org/10.1016/s0017‐9310(02)00507‐0.
76. Jolls, K.R.; Hanratty, T.J. Transition to Turbulence for Flow through a Dumped Bed of Spheres. Chem. Eng. Sci. 1966, 21, 1185–
1191. https://doi.org/10.1016/0009‐2509(66)85038‐8.
77. Wegner, T.H.; Karabelas, A.J.; Hanratty, T.J. Visual Studies of Flow in a Regular Array of Spheres. Chem. Eng. Sci. 1971, 26, 59–
63. https://doi.org/10.1016/0009‐2509(71)86081‐5.
78. Latifi, M.A.; Midoux, N.; Storck, A.; Gence, J.N. The Use of Micro‐Electrodes in the Study of the Flow Regimes in a Packed‐Bed
Reactor with Single‐Phase Liquid Flow. Chem. Eng. Sci. 1989, 44, 2501–2508. https://doi.org/10.1016/0009‐2509(89)85194‐2.
79. Rode, S.; Midoux, N.; Latifi, M.A.; Storck, A.; Saatdjian, E. Hydrodynamics of Liquid Flow in Packed Beds: An Experimental
Study Using Electrochemical Shear Rate Sensors. Chem. Eng. Sci. 1994, 49, 889–900. https://doi.org/10.1016/0009‐2509(94)80025‐
1.
80. Bu, S.; Yang, Y.; Zhang, X.; Tao, W. Experimental Study of Transition Flow in Packed Beds of Spheres with Different Particle
Sizes Based on Electrochemical Microelectrodes Measurement. Appl. Therm. Eng. 2014, 73, 1525–1532.
https://doi.org/10.1016/j.applthermaleng.2014.03.063.
81. Dybbs, A.; Edwards, R.V. A New Look at Porous Media Fluid Mechanics—Darcy to Turbulent. In Fundamentals of Transport
Phenomena in Porous Media; Bear, J., Corapcioglu, M.Y., Eds.; Springer: Dordrecht, The Netherlands, 1984; pp. 199–256.
https://doi.org/10.1007/978‐94‐009‐6175‐3_4.
82. Woessner, W.W.; Poeter, E.P. Total Porosity. In Hydrogeologic Properties of Earth Materials and Principles of Groundwater Flow;
Groundwater Project: Guelph, ON, Canada, 2020.
83. Banach Fixed‐Point Theorem, Wikipedia. Available online: https://en.wikipedia.org/wiki/Banach_fixed‐point_theorem (ac‐
cessed on 18 July 2024).
84. Wood, W.W.; Kraemer, T.F.; Hearn, P.P., Jr. Intragranular Diffusion: An Important Mechanism Influencing Solute Transport in
Clastic Aquifers? Science 1990, 247, 1569–1572. https://doi.org/10.1126/science.247.4950.1569.
85. Berkowitz, B.; Scher, H. On characterization of anomalous dispersion in porous and fractured media. Water Resour. Res. 1995,
31, 1461–1466. https://doi.org/10.1029/95WR00483.
86. Haggerty, R.; Gorelick, S.M. Multiple‐rate mass transfer for modeling diffusion and surface reactions in media with pore‐scale
heterogeneity. Water Resour. Res. 1995, 31, 2383–2400. https://doi.org/10.1029/95WR10583.
87. Medina, A.; Carrera, J. Coupled Estimation of Flow and Solute Transport Parameters. Water Resour. Res. 1996, 32, 3063–3076.
https://doi.org/10.1029/96WR00754.
88. Carrera, J.; Sánchez‐Vila, X.; Benet, I.; Medina, A.; Galarza, G.; Guimerà, J. On matrix diffusion: Formulations, solution methods
and qualitative effects. Hydrogeol. J. 1998, 6, 178–190. https://doi.org/10.1007/s100400050143.
89. Bolster, D.; Dentz, M.; Le Borgne, T. Solute dispersion in channels with periodically varying apertures. Phys. Fluids 2009, 21,
056601. https://doi.org/10.1063/1.3131982.
90. Wang, X.; Sheng, J.J. Effect of low‐velocity non‐Darcy flow on well production performance in shale and tight oil reservoirs.
Fuel 2017, 190, 41–46. https://doi.org/10.1016/j.fuel.2016.11.040.
91. Theis, C.V. The Relation between the Lowering of the Piezometric Surface and the Rate and Duration of Discharge of a Well
Using Ground‐Water Storage. Eos Trans. Am. Geophys. Union 1935, 16, 519–524. https://doi.org/10.1029/TR016i002p00519.
92. Hantush, M.S. Modification of the Theory of Leaky Aquifers. J. Geophys. Res. 1960, 65, 3713–3725.
https://doi.org/10.1029/JZ065i011p03713.
93. Hunt, A.; Ewing, R.; Ghanbarian, B. Percolation Theory for Flow in Porous Media; Springer: Berlin/Heidelberg, Germany, 2014;
Volume 880.
94. Gardner, W.R. Some Steady‐State Solutions of the Unsaturated Moisture Flow Equation with Application to Evaporation from
a Water Table. Soil Sci. 1958, 85, 228–232.
95. Pozdnyakov, A.I.; Pozdnyakova, L.A.; Karpachevskii, L.O. Relationship between Water Tension and Electrical Resistivity in
Soils. Eurasian Soil Sci. 2006, 39, S78–S83. https://doi.org/10.1134/S1064229306130138.
96. Ameli, A.A.; McDonnell, J.J.; Bishop, K. The Exponential Decline in Saturated Hydraulic Conductivity with Depth: A Novel
Method for Exploring Its Effect on Water Flow Paths and Transit Time Distribution. Hydrol. Process. 2016, 30, 2438–2450.
https://doi.org/10.1002/hyp.10777.
97. Takematsu, M. Slow viscous flow past a cavity. J. Phys. Soc. Jpn. 1966, 21, 1816–1821. https://doi.org/10.1143/JPSJ.21.1816.
98. Friedman, M. Flow in a circular pipe with recessed walls. J. Fluid Mech. 1970, 37, 5–8. https://doi.org/10.1115/1.3408490.
99. Stevenson, J.F. Flow in a tube with a circumferential wall cavity. J. Appl. Mech. Trans. ASME 1973, 40, 355–361.
https://doi.org/10.1115/1.3422987.
Water 2024, 16, 2158 27 of 27

100. Driesen, C.H.; Kuerten, J.G.M.; Streng, M. Low‐Reynolds‐number flow over partially covered cavities. J. Eng. Math. 1998, 34, 3–
21. https://doi.org/10.1023/A:1004235021527.
101. Young, A.H.; Hotz, N.; Hawkins, B.T.; Kabala, Z.J. Inducing Deep Sweeps and Vortex Ejections on Patterned Membrane Sur‐
faces to Mitigate Surface Fouling. Membranes 2024, 14, 21. https://doi.org/10.3390/membranes14010021.
102. Taneda, S. Visualization of separating Stokes flows. J. Phys. Soc. Jpn. 1979, 46, 1935–1942. https://doi.org/10.1143/JPSJ.46.1935.
103. Laskowska, A. Experimental Studies of Flows in Porous Media and Selected Models of the Pore Space. Ph.D. Dissertation, Strata
Mechanics Research Institute Polish Academy, Krakow, Poland, 1996; pp. 179.
104. Pan, F.; Acrivos, A. Steady flows in rectangular cavities. J. Fluid Mech. 1967, 28, 643–655.
https://doi.org/10.1017/S002211206700237X.
105. Young, A.H.; Kabala, Z.J. Hydrodynamic Porosity: A New Perspective on Flow Through Porous Media, Part II. Water 2024, 190,
41–46.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual au‐
thor(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.