Hydrodynamic Porosity: A Paradigm Shift in Flow and Contaminant Transport Through Porous Media, Part II

In this work, we build upon our previous finding that hydrodynamic porosity, 𝜃 !"#$%& , is an exponential function of pore-scale 10 flow velocity (or interstitial Reynolds number). We previously discovered this relationship for media with a square cavity geometry – a highly idealized case of the dead-ended pore spaces in a porous medium. Thus, we demonstrate the applicability of this relationship to media with other cavity geometries. We do so by applying our previous analysis to rectangular and non-rectangular cavity geometries (i.e., circular, and triangular). We also study periodic flow geometries to determine the effect of upstream cavities on those downstream. We show that not only does our exponential relationship hold for media with a variety 15 of cavity geometries, but it does so almost perfectly with a coefficient of determination ( R² ) of approximately 1 for each new set of simulation data. Given this high fit quality, it is evident that the exponential relationship we previously discovered is applicable to most, if not all, unwashed media.


Introduction
In this paper, we build on our previous findings, i.e., Young and Kabala (2023), where we discovered and quantified the exponential dependence of hydrodynamic porosity ( !"#$%& ) on pore-scale flow velocity ( '"(& ) or interstitial Reynolds Number ().This relationship is reproduced below:  !"#$%& = ( +  )* , !"#$ )  ⟺  !"#$%& = ( +  )-.& ) ,  =  ℎℎ  (1) Where  !"#$%& is the mobile-zone porosity (what we refer to as the hydrodynamic porosity of the medium),  is the total porosity of the medium (i.e.,  '"(& / #/%0 ), ℎℎ is the height of the through-channel in the dead-end pore model (m),  is the kinematic viscosity (m²/s), and a, b, and c are dimensionless parameters; the exponential parameter, d has units of (s/m) to allow for use of pore-scale velocity in place of Reynolds number.Readers are further reminded that a is the value of the porescale partitioning coefficient () approximated by the boundary-driven, or analogously, Mobile-Immobile Zone model, both of which are discussed at length by Young and Kabala (2023), and originally applied to groundwater flows by Vangenuchten and Wierenga (1976).In the case of the mobile separatrix we study here, it is the value that results from Re → ∞.Finally, the quantity 'a + b,' is the value of  in the creeping flow regime (i.e., the value that results from Re → 0).
Equation ( 1) was discovered and quantified for the dead-end pore geometry reproduced in Figure 1, below.Readers will note that the cavity geometry is square.Thus, we did not previously determine if this relationship holds for more than this single idealized case.This work is motivated by the fact that the dead-ended pore spaces in unwashed porous media (e.g., glacial deposits, fractured rock, and filtration media such as granulated activated carbon, etc.) are poorly approximated by the square cavity geometry we previously studied.The jagged edges of such media result in a random array of pore geometries (as illustrated below).Subsequently, we must be able to answer the question: does the exponential dependence of hydrodynamic porosity,  !"#$%& , on pore-scale flow velocity hold for media with cavity geometries that are not square?If the answer is indeed yes, researchers are justified in utilizing the exponential relationship provided in Eq. ( 1).Further, researchers can expand upon this relationship by determining if  !"#$%& exhibits a dependence on any other aspects of the dead-end pore geometry (e.g., depth, width, hydraulic radius, etc.). 2 Related Work

Flows Past Cavities
In our previous work, we provide a thorough literature review of flows past cavities, citing the major contributions that have led to a comprehensive understanding of cavity vortex structure(s) (Moffatt, 1963;Higdon, 1985;Shen and Floryan, 1985;Fang et al., 1997).Readers are directed to this discussion for further detail.In the following cavity geometries studied by these (and other) authors.Following the precedent set by these authors, we proceed by utilizing the same cavity geometries in our analysis.• Kahler and Kabala (2016) • Kang and Chang (1982) Table 1 illustrates that most research has been carried out on rectangular cavity geometries.This is unsurprising given that sharp corners in flow geometries result in flow separation.Mehta and Lavan (1969), O' Brien (1972), Shen and Floryan (1985), Fang et al. (1997), andHigdon (1985) utilize rectangular geometries and manipulate the cavity width and/or depth to determine the subsequent effect on the cavity vortex structure.Moffatt (1963), in his study on resistive eddies that result from sharp corners, focuses most of his work on triangular structures, while Higdon (1985) and O'neill (1977) study circular cavities.

Our Previous Work
The idealized, square cavity geometry we previously utilized in Young and Kabala (2023) is justified by its foundational use in the study of flows past cavities, as noted by Shankar and Deshpande (2000) and is included in most of the studies referenced in Table 1 above.To exaggerate the dependence of hydrodynamic porosity,  !"#$%& , on pore velocity, we also narrowed the height of the through-channel, relative to the cavity depth and width, which is kept at a one-to-one aspect ratio so that the cavity remains square; readers are referred to Figure 6 of Young and Kabala (2023) for reference.
We then subjected these dead-end pore geometries to inlet flows with an interstitial (channel-based) Reynolds number in the range of 1-100, reasoning that doing so keeps the flow within the laminar regime.To justify this claim, we reference multiple examples that demonstrate a general trend: the onset of turbulence occurs at or above Reynolds numbers of 100 in porous media (Jolls and Hanratty, 1966;Wegner et al., 1971;Latifi et al., 1989;Rode et al., 1994;Bu et al., 2014).For each dead-end pore geometry, we discovered an exponential dependence of hydrodynamic porosity,  !"#$%& , on pore-scale flow velocity, with a coefficient of determination (R²) of approximately 1 for each set of simulation data.As anticipated, this relationship was most exaggerated when the through-channel height of the dead-end pore geometry was reduced to one-quarter of the cavity depth/width.For media with this cavity type,  !"#$%& decreased by 42% over the laminar flow regime.For reference, over the same Reynolds number range,  !"#$%&only decreased by 4% when the through-channel height was equivalent to the cavity depth/width.

Methods
In short, the hydrodynamic porosity,  !"#$%& , of a porous medium with dead-ended pore spaces is ultimately determined by the location of the separatrix (i.e., the streamline that divides the through-channel flow from the recirculatory cavity flow), which we track over a range of flow conditions by varying the magnitude of the interstitial Reynolds number of the inlet flow.
Further discussion on the separatrix is provided by Elderkin (1975), Weiss (1991), Horner et al. (2002), and Kahler and Kabala (2016).We track the location of the separatrix because it defines the magnitude of the mobile zone (i.e., the pore space conducive to through-flow), and therefore the value of the pore-space partitioning coefficient, : As defined above,  is the ratio of the mobile-zone volume ( !"#$%& ) to the total pore space volume ( '"(& ).For an isotropic, or 2D media, like the ones we study in this work,  can be defined in terms of cross-sectional areas, .To relate  to  !"#$%& , we remind readers of the definition of , the total porosity: Similarly,  !"#$%& is defined as: Combining Eq. ( 2) -( 4), we find that the product of the total porosity of the medium and the partitioning coefficient yield the hydrodynamic porosity,  !"#$%& , of the medium: The methods and tools used in this analysis are the same as those used in our previous work.Readers are directed to Young and Kabala (2023) for a description of the numerical flow solver we leverage in this analysis, as well as a detailed outline of the data collection process.Slight modifications to the flow solver program allow for manipulation of the cavity geometry, which we picture in the figures below.With each cavity manipulation, we adjust the refinement region of the solver mesh to fully encapsulate the area in which we assume the separatrix to be.Altogether, we test rectangular cavities with depth-to-width ratios of 2:1, 1.5:1, 1:1, and 0.5:1 (pictured in Figure 2), and 1:2, 1:1.5, 1:1, and 1:0.5 (pictured in Figure 3).The varying channel length in the geometries provided in Figure 3 is an artifact of the decision to extend the length of the through-channel past the cavity by twice the cavity width.We also test two non-rectangular cavities (i.e., a semicircle and an equilateral triangle).Finally, we test a periodic flow geometry to determine the effect upstream cavities have on those downstream.These additional flow geometries are provided below, in Figures 4 and 5.
105    Again, these flow geometries are subjected to an inlet flow with channel-based Reynolds numbers in the laminar flow regime (i.e., Re < 100).For ease of comparison to our previous results, we determine the location of the separatrix at the same Reynolds numbers (i.e., Re = 1.81, 3.61, 5.52, 9.19, 25, 50, 75, 100).Note that we do not test Reynolds numbers in the creeping flow regime, which in this case, is less than 1.This is because we previously found that the separatrix is immobile in this regime (Young and Kabala, 2023).After having subjected each flow geometry to the specified Reynolds numbers, we apply the exponential relationship in Eq. ( 1) that we discovered for the square cavity geometry to each new set of simulation data.We determine the fit quality of the exponential model via the coefficient of determination (R²).Finally, given the sparse sampling we previously used to generate the exponential relationship we define in Eq. ( 1), we also determine the error associated with this low sampling frequency; we do this by testing Reynolds numbers of 1-100 in increments of 5.

Results
Stream plots for each cavity geometry at Reynolds numbers of approximately 1, 10, 50, and 100, are provided in the figures below.In each stream plot, the separatrix is highlighted in red.Again, because the cavity flow is driven by the adjacent through- channel flow, we also provide a stream plot of the entire dead-end pore geometry adjacent to the corresponding cavity flow.
Figure 6 is reproduced from our previous work.

Separatrix Movement: Rectangular Cavity Geometry Manipulation 130
In Figure 7 and 8, we provide the stream plots that result when we adjust the depth and width of the square cavity geometry we previously tested.We also illustrate movement of the separatrix, which is highlighted in red.

Separatrix Movement: Non-Rectangular Cavity Geometries
In Figure 11, we provide the stream plots for non-rectangular cavity geometries (i.e., the semicircle and equilateral triangle) with the separatrix highlighted in red.

Separatrix Movement: Periodic Cavity Geometry
For the reader's edification, a single stream plot example of the periodic cavity flow geometry is provided below.Additional stream plots are not provided because the streamlines in each cavity are effectively identical for each of the cavity geometries pictured in Figure 12.The purpose of testing a periodic flow geometry was to determine the effect of upstream cavities on those downstream.Qualitatively, we can conclude that there is slightly less penetration of the through-flow into the 160 downstream cavity, although this is barely evidenced by the red horizontal line in the figure below.As drawn, this line coincides with the bottom-most point of the separatrix in the upstream cavity.Note that this figure is generated for the periodic flow geometry with cavity spacing one-quarter of the cavity width, at a Reynolds number of 10.

Applying the Exponential Dependence of Hydrodynamic Porosity on Pore-Scale Flow Velocity
The change in hydrodynamic porosity,  !"#$%& , as a function of Reynolds number, is provided in Figure 13 and 14, below.For ease of comparison across cavity types, we normalize  !"#$%& by the static mobile-zone porosity value,  121 , that results from enforcement of the Mobile-Immobile Zone model in the dead-end pore space.We define this quantity as: Where  121 is determined by the relative magnitudes of the through-channel and cavity volumes for each dead-end pore geometry.For example, using Eq. ( 2), we find that for the square cavity,  121 = 4/5.
Applied to each dataset, is the exponential dependence of hydrodynamic porosity,  !"#$%& , on pore flow velocity we previously discovered and described in Eq. (1).We remind readers that although it is the value of the pore-space partitioning coefficient,  =  !"#$%& /, that we numerically quantify in this analysis, we are able to plot the hydrodynamic porosity of the medium,  !"#$%& , because of the direct proportionality between these two quantities.The fitting parameters and coefficient of determination (R²) for each model are summarized in the table below.

Sampling Error
Finally, the error associated with the exponential model derived from the relatively sparse dataset is approximated by fitting the exponential model to a more finely sampled dataset, pictured below.The results provided in Figure 15 are for the square 195 cavity geometry.

Discussion
In this work, we find that the pore-space partitioning coefficient,  =  !"#$%& /, of each of the tested cavity geometries is well-approximated by the exponential relationship in Eq. ( 1).In fact, the coefficient of determination (R²) is approximately 1 for each dataset.Should we test additional cavity geometries in the future, given representative imaging of the pore space in an unwashed media, it would be reasonable to assume that the partitioning coefficient, and therefore the hydrodynamic porosity of that media,  !"#$%& , would be well-approximated by this same exponential relationship.Further, our sampling frequency in the previous study did not yield a significant error in our model fit.In fact, when comparing the exponential models that result from both data sets, the largest difference between them is 0.1%.This error is an order of magnitude smaller than the error associated with the measurement process, which we judge to be small enough to disregard in our analysis.
Of the tested geometries, those that exhibited the largest decrease in the pore-space partitioning coefficient,  =  !"#$%& /, accompanying an increase in Reynolds number from 1 to 100 were rectangular.More specifically, the two largest decreases were for geometries with a cavity width greater than its depth.When the cavity width was less than its depth, however, the decrease in  was the smallest for all the tested geometries.In comparison, adjustment to the cavity depth did not result in such large variations.Although only six rectangular cavities were tested, it appears that  is more sensitive to cavity width than cavity depth.Of course, this conclusion is made for a flow geometry with a through-channel height equivalent to either the cavity width or depth.Reduction of the adjacent through-channel height, as we previously demonstrated, results in more significant decreases to the mobile zone volume.These results are summarized in the table below.Further increasing the applicability of the exponential dependence of hydrodynamic porosity,  !"#$%& , on pore-scale flow velocity, is the fact that the periodic flow geometries do not exhibit a significant difference between the upstream and downstream cavities.Given this finding, we can more confidently apply our hydrodynamic description of porosity onto a periodic medium at the macroscale.As previously mentioned, the next steps for this work are in the experimental quantification of the exponential fit parameters in Eq. (1) (i.e., a, b, and c).The ability to do so for media at sites needing remediation is necessary to demonstrate the ease with which this relationship can be tailored to a specific media.Further, it may be the case that the numerically calculated parameters provided in Table 2 indicate a physical characteristic of the aggregate, effective pore geometry of real media.For example, consider the rectangular cavities with a depth-to-width ratio of 1:2 and 1:0.5.The former has a b-value one order of magnitude larger than the rest, and the latter, a c-value one order of magnitude smaller than the rest.The values in Table 2 also serve as an order of magnitude guide for what the exponential fit parameters should be in column studies.

Conclusions
In this work, we further demonstrate the applicability of the exponential dependence of hydrodynamic porosity,  !"#$%& , on pore-scale flow velocity.We do so by testing a range of cavity geometries that are more representative of what we would expect to see in a real porous medium (i.e., rectangular, non-rectangular, and periodic geometries).Previously, we only tested a single, square cavity geometry, which is a highly idealized case.We show that not only does this exponential relationship hold for media with these new cavity geometries, but it does so almost perfectly with a coefficient of determination (R²) of approximately 1 for each set of simulation data.Given this high fit quality, it is evident that the exponential relationship we previously discovered is applicable to most, if not all, media with dead-ended pore spaces.That said, this work would lend itself well to significant expansion.Although we have highlighted sensitivity to cavity width, and previously, through-channel height to pore depth, there remain other parameters that could affect the exponential fit and perhaps, even inform us about the aggregate, effective nature of the pore spaces of experimentally studied media.

Figure 1 :
Figure 1: Illustration of the pore space and cavities in an unwashed porous media and the idealized dead-end pore geometry used by Young and Kabala (2023) (left and right, respectively).

Figure 5 :
Figure 5: Periodic flow geometries with square cavities.The spacing between the cavities is a fraction of the cavity depth/width (i.e., from left to right, top to bottom: 1, 0.5, 0.33, 0.25).

Figure 12 :
Figure 12: Stream plot of the cavities in the periodic flow geometry with a cavity spacing of one-quarter the cavity depth/width.The horizontal red line marks the bottom-most point of the separatrix of the upstream cavity (left).

Figure 14 :
Figure 14: Normalized decrease in hydrodynamic porosity,   , for media with triangular and semi-circular cavities, compared to media with square cavities.

Figure 15 :
Figure 15: (top) Comparison of the exponential model that results from a sparse dataset and the exponential model that results from 200