Measuring porous media velocity fields and grain bed architecture with a quantitative PLIF-based technique

Porous media flows are common in both natural and anthropogenic systems. Mapping these flows in a laboratory setting is challenging and often requires non-intrusive measurement techniques, such as particle image velocimetry (PIV) coupled with refractive index matching (RIM). RIM-coupled PIV allows the mapping of velocity fields around transparent solids by analyzing the movement of neutrally buoyant micron-sized seeding particles. The use of this technique in a porous medium can be problematic because seeding particles adhere to grains, which causes the grain bed to lose transparency and can obstruct pore flows. Another non-intrusive optical technique, planar laser-induced fluorescence (PLIF), can be paired with RIM and does not have this limitation because fluorescent dye is used instead of particles, but it has been chiefly used for qualitative flow visualization. Here, we propose a quantitative PLIF-based methodology to map both porous media flow fields and porous media architecture. Velocity fields are obtained by tracking the advection-dominated movement of the fluorescent dye plume front within a porous medium. We also propose an automatic tracking algorithm that quantifies 2D velocity components as the plume moves through space in both an Eulerian and a Lagrangian framework. We apply this algorithm to three data sets: a synthetic data set and two laboratory experiments. Performance of this algorithm is reported by the mean (bias error, B) and standard deviation (random error, SD) of the residuals between its results and the reference data. For the synthetic data, the algorithm produces maximum errors of B & SD = 32% & 23% in the Eulerian framework, respectively, and B & SD = −0.04% & 3.9% in the Lagrangian framework. The small-scale laboratory experimental data requires the Eulerian framework and produce errors of B & SD = −0.5% & 33%. The Lagrangian framework is used on the large-scale laboratory experimental data and produces errors of B & SD = 5% & 44%. Mapping the porous media architecture shows negligible error for reconstructing calibration grains of known dimensions.


Introduction
Below the surface water of rivers and streams, there exists a subsurface zone where porous media flows impact ecological and biological processes by exchanging water-laden solutes with the surface water (Buffington andTonina 2009, Tonina andBuffington 2009a).This exchanging phenomenon is referred to as hyporheic exchange, and it not only supports benthic life (Freixa et al 2016), it also harbors biogeochemical reactions that affect global warming rates (Quick et al 2016).
One example of organisms relying on hyporheic exchange is the survival of salmonids (Stuart 1953).Female salmonids create a riverbed structure-called a redd-that is shaped like a dune and lay eggs within the structure (Deverall et al 1993).
The redd is shaped in such a way that downwelling hyporheic fluxes deliver dissolved oxygen to the roe, as well as regulate the surrounding temperature (Tonina and Buffington 2009b).Upwelling fluxes evacuate any embryo respiratory byproducts, causing a continuous cycle of fresh oxygen supply to the eggs (Cooper et al 1965, Cardenas et al 2016).
Experimentally investigating these porous media flows can lead to a better understanding of their processes and aid conservation and restoration efforts.
In a laboratory setting, flows can be visualized with nonintrusive optical techniques, such as particle image velocimetry (PIV) (Adrian 2005, Daaboul et al 2009, Li et al 2017), particle tracking velocimetry (PTV) (Malik et al 1993, Ohmi andLi 2000), or planar laser-induced fluorescence (PLIF) (Crimaldi 2008).PIV and PTV methods rely on neutrally buoyant micron-sized seeding particles to map flow fields, whereas PLIF methods utilize soluble dyes that are fluoresced with light.These techniques can be paired with the refractive index-matching (RIM) technique (Budwig 1994), where the refractive index (RI) of the liquid is matched to that of the solid to study flow around and within complex objects, like grains (Rubol et al 2018) or vertical cylinders (Weitzman et al 2014).While RIM-PIV and RIM-PTV have been used to quantitatively map porous media flow fields (Cenedese and Viotti 1996, Voermans et al 2017, Carrel et al 2018, Hilliard et al 2021), seeding particles adhere to the solid medium over time, causing a loss of transparency and obstruction of the porescale flows.PLIF techniques avoid these drawbacks by using a soluble dye instead of seeding particles.A few studies have employed PLIF techniques to quantitatively study freestream flows, such as Gendrich and Koochesfahani (1996, molecular tagging velocimetry, or MTV) and Dahm et al (1992, scalar image velocimetry, or SIV).However, the literature is lacking on employing PLIF techniques to map porous media flow fields quantitatively.
MTV uses phosphorescent molecules that turn into longlifetime tracers when excited by UV-laser light.These molecules are mixed into the working fluid, and regions of interest are successively 'tagged' with a pulsed UV-laser grid resulting in a Lagrangian displacement vector.Studies have shown this method's exceptional performance in mapping freestream flows (Hu andKoochesfahani 2006, Bohl andKoochesfahani 2009).However, it has not been used in porous media flows because the laser grid will not appear within the solid phase, where the phosphorescent molecules do not exist.
SIV applies the scalar transport equation to spatiotemporal scalar images to extract the underlying velocity field.High-speed planar imaging captures laser-induced fluorescence whose concentration is a conserved scalar variable.Care is taken to ensure the spatial and temporal resolutions are less than the molecular diffusion and advection time scales, respectively.A few studies have exemplified the application of this method in turbulent freestream flows (Dahm et al 1991, 1992, Su and Dahm 1996).Similarly to MTV, SIV does not work in porous media flows when using planar measurements because in porous media flows, the scalar field is intrinsically three-dimensional and thus the scalar transport equation requires a volumetric measurement of the scalar quantity.
When PLIF techniques are employed to study flows, intensity gradients exist between empty and saturated regions of the fluorescent dye.Instead of using MTV or SIV methods, one could simply track the movement of the intensity gradient to obtain a velocity field.Injecting a fluorescent dye plume into a porous media flow creates a sharp gradient which can be thought of as a leading edge, or front, of the injected fluorescent dye plume that can be tracked to obtain a velocity field (figure 1).The sharp gradient also appears between the solid and liquid phases when the porous medium is completely saturated with a fluorescent dye, allowing investigations into the architecture of the porous medium (Hilliard et al 2021).Here, we propose a methodology to obtain two-dimensional velocity field measurements within porous media by tracking the advective movement of a dye front, and revisit a methodology to reconstruct the porous media architecture, all with PLIF techniques.

Methods
We use synthetic and experimental data sets to quantify the performance of the proposed methodology to automatically map a porous media flow field from image analysis at the pore and Darcian (several pores) scales.The methodology is also applied manually to obtain a reference velocity field for the automatic tracking algorithm when the true velocity field is unknown.Performance of the methodology is quantified by the mean (bias error) and standard deviation (random error) of the residuals between the measured and reference velocity fields.We present these values as relative errors by normalizing the residuals with the reference velocity, providing a percent error rather than presenting nominal residual velocities because those are very small within pores.

Synthetic data set.
We numerically generate a synthetic data set containing the explicit movement of a circular dye plume (figure 2(a)).This synthetic data set allows for the performance evaluation of manually tracking an advective dye front as well as the automated tracking algorithm.The synthetic binary images are created in MATLAB ® by explicitly defining the center location and diameter (40 pixels, or 1.6 mm) of a non-deforming circle using the 'poly2mask' function over a black background.The images are then successively added so that the movement of the circular dye plume appears as a continuous dye injection.Though this data set is not experimentally obtained, experiment parameters such as the recording rate and image resolution need to be specified to extract velocities in mm s −1 .Therefore, the synthetic recording rate and image resolution is 4 Hz and 25.6 px mm −1 , respectively, which are values expected for the small-scale laboratory experiments.The movement of the circular dye plume is specified in three regionsdownwelling (frames 1-166), crosswelling (frames 167-336), and upwelling (frames 337-481).In the downwelling and upwelling regions, the movement of the dye circle center is set to five pixels per frame (0.78 mm s −1 ) in the positive and negative y-directions, respectively.In the crosswelling region, a parabolic equation specifies the trajectory of the circle center, being defined by three points: the last coordinate of the downwelling zone, the first coordinate of the upwelling zone, and a third point between and below the two aforementioned points.The x-coordinate shift of the circle center in this crosswelling region is set to five pixels per frame, resulting in an average velocity magnitude of 0.89 mm s −1 .

Small-scale laboratory experimental setup.
A thin flow cell is constructed to simulate two-dimensional porous media flow (figure 3(a)).The front, left, and right walls are of 3 mm thick pane glass to allow optical access, and a 1 cm thick piece of acrylic plastic painted matte black creates the back wall.A 3 mm thick HDPE divider (also painted matte black) is installed so that the flow cell and part of the porous medium is split into two sides, creating a porous media flow with strong overall movements in both the x-and ydirections.A small peristaltic pump (figure 3(b)) keeps the fluid on the left side approximately 1 mm deeper than the fluid on the right side.This maintains a constant head differential   flow cell are hereafter referred to as the small-scale laboratory experiments.
The porous medium is made of non-spherical 3 mm mean diameter transparent polymer grains (Tetrafluoroethylene hexafluoropropylene vinylidene fluoride, or THV 221-GZ, produced by 3M) with a specific gravity of 1.97 and a RI of 1.365.We use an aqueous 24.2% by-weight glycerin solution as the non-toxic fluid, which matches the RI of the grains and allows optical access through the solid (Hilliard et al 2021).
Rhodamine B fluorescent dye (554 nm peak excitation wavelength, 576 nm peak emission wavelength, both in water; Kristoffersen et al 2014) is mixed with the RIM fluid and injected into the porous medium.532 nm 50 mW continuous waveform diode lasers shining through sheet optics (creating a ∼0.5 mm-thick laser light sheet) are placed on both sides of the flow cell (figure 3(b)) to fluoresce the injected dye and illuminate the flow in a planar two-dimensional fashion (figure 2(b)).An orange Schott glass filter (OG550) is placed between the camera lens and the flow cell so that only the fluoresced light from the dye is observable.Once the dye is injected near the surface of the porous bed on the left side of the flow cell, live recording begins using a 5-MPx LaVision Imager MX camera, and data is collected until the dye is no longer visible in the porous medium.Only the downwelling portion of the flow (left side of the flow cell) is analyzed from the small-scale laboratory experiment data because the dye plume becomes much more diffuse and untraceable as it moves further through the bed (see supplementary material 1).
Raw images are preprocessed before extracting velocities to improve the contrast between the empty and saturated regions of dye.For the small-scale laboratory experiments, dark-frame subtraction is performed to improve the contrast in the images.Afterwards, an Otsu threshold filter (Otsu 1979) is applied to create binary images which are finally cleaned up using the 'Despeckle' function in the ImageJ software.
2.1.3.Large-scale laboratory experimental setup.We also perform a case study involving a salmon redd model that is constructed in the Aquatic Imaging Flume (AIF) at the University of Idaho Center for Ecohydraulics Research (figure 4(a)).The AIF is a tilting and recirculating flume that is 7 m long, 0.5 m wide, and 0.7 m deep (Moreto et al 2022).The redd model has a length of 1 m, from the start of the pit to the end of the tailspill (DeVries 2012), and is made up of THV, identical to what is used in the small-scale laboratory experiments (figure 5).The other 6 m of the flume is filled with crushed glass of similar size to the THV (figure 4(b)).The flow through the model redd is studied in a similar planar two-dimensional fashion as the small-scale laboratory experiments.Rhodamine B dye is mixed with a non-toxic RIM fluid of 15% by-weight MgSO 4 solution, which is used as the working fluid in the AIF.The dye is injected at specific points along the windward surface of the model redd, and green laser light (532 nm) fluoresces the dye as it moves through the porous structure.Images have a resolution of 16.5 px mm −1 and are captured at 2 Hz on the same camera as the small-scale laboratory experiments, also with an orange Schott glass filter in place (figure 2(c)).For the data presented in this study, the average streamwise velocity is 0.19 m s −1 with a depth of 0.103 m for the surface flow.This case study will hereafter be referred to as the large-scale laboratory experiment.Like the small-scale laboratory experiments, the raw images from the large-scale laboratory experiments are preprocessed to improve the dye contrast, but with a few added steps to improve signal quality.Data recording began before dye is injected, so the first few images are treated as a background-averaged over time and then subtracted from the data set, increasing the visibility of the dye movements.This is done in lieu of the dark-frame subtraction process.Next, a Gaussian blur filter is applied in the ImageJ software to reduce noise for the binarization of the data, which is done by an Otsu threshold filter (Otsu 1979).Finally, any small objects present in the planar binary images that appear to be discontiguous from the main dye injection plume are eliminated based on an area threshold (<100 000 px 2 ).This is done because portions of the injected dye plume can sometimes move out of and back into the illumination plane, causing the appearance of disconnected dye objects that affect the velocity readings.While this three-dimensional motion is real, and the disconnected dye objects are likely connected to the dye plume in 3D space, the planar setup of these experiments is incapable of properly tracking these movements.

Velocity extraction method
The proposed methodology assumes that an injected dye plume moves chiefly advectively within a porous medium such that diffusion has negligible impact on its displacement.At the spatial scales we are investigating, advective dispersion is indistinguishable from molecular diffusion.However, at the velocity magnitudes that we are investigating, the dye front remains distinct, unlike other studies where the advective velocity could be comparable to molecular diffusion rates (Wirner et al 2014).
This assumption suggests that tracking the evolution of the plume front over time provides the displacements ∆x i in the ith direction due to the flow field.Thus, the local velocity components, u i , are quantified by dividing the displacement of the front along the i-th direction by the change in time, T Here, we focus on two-dimensional measurements such that the i-th direction is either horizontal (x, i = 1) or vertical (y, i = 2).Examples of injected dye plumes and their fronts can be seen in figures 2(b) and (c).
The application of this method is at both the pore scale and larger Darcian scale, depending on the dye injection conditions and image resolution.Tracking the dye front at these scales can be done either manually or automatically.Manual tracking involves a simple point-and-click process to track the movement of the dye front from frame to frame in the ImageJ software package (Ershov et al 2022).However, to remove the tedious process of manually tracking the dye front, we develop an automatic algorithm to track the dye front and map the flow field.Automatically tracking the dye front requires interrogation cells superimposed on the data, but the size of these cells depends on the flow scale.At the pore scale, a regular grid of square interrogation cells (e.g.32 × 32 pixels, 64 × 64 pixels, or larger) is superimposed on the preprocessed, binary images to create partitions for individual analysis.Sizing the interrogation cells depends on the scale of the analyzed flow and should be chosen such that the cell size is 1-2 pore diameters.This framework is advantageous when a uniform dye front cannot be identified (figure 2(b)) because it splits the dye movement into a set of small individual dye fronts whose movements can be tracked within each interrogation cell, resulting in a velocity distribution measured in the Eulerian framework.At this scale, the method maps the velocity field from the fingering of the dye around grains to measure local pore velocity.Conversely, large dye injections and/or lower image resolution will cause the dye plume to span several pores (figure 2(c)).This results in a dye plume with a regular front, which can be better tracked with a single interrogation cell that is of the size of the images and results in a velocity distribution measured in the Lagrangian framework.This framework is advantageous when a uniform dye front can be identified because no partitioning is needed to track the dye front and therefore requires much less computational effort.

Automatic velocity field reconstruction
The automatic tracking algorithm works by leveraging the properties of the 'regionprops' MATLAB ® function (R2021a; MATLAB, Image Processing Toolbox ™ , 1993; pp 1-2790; last accessed on 23rd Jan, 2023) (see supplementary material 2).The function input is a binary image and the output is a set of properties of the input binary image, including a property entitled 'BoundingBox'.This property contains the x-coordinate, y-coordinate, width, and height of a rectangle that envelops all binary objects present in the input image (figure 6).First, we quantify the velocity of each box side, u i,n,j , where i indicates the direction (1, horizontal, or 2, vertical), n indicates the temporal location (i.e.frame number), and j indicates the first or second side for the respective direction (1 for the left or top side and 2 for the right or bottom side).This results in four bounding box side velocities for a single temporal location (2) However, we are interested in the velocity of the dye front as a rigid body and not in its deformation.Thus, the final flow velocity components, u i,n , which move the front as a rigid body at the n-th location, are quantified based on the following assumptions: when the velocities of opposing sides (left-right or top-bottom) have opposite directions, we take the vectoral sum between u i,n,1 and u i,n,2 .When they have the same direction, we select the velocity of the front side, determined by the signs of u i,n,1 and u i,n,2 (see coordinate axes in figure 6): with the velocity resultant, U n , being: This algorithm can be applied in the Eulerian or Lagrangian framework.In the Eulerian framework, each interrogation cell provides one velocity reading.Once a velocity reading is acquired, it stops analyzing the interrogation cell and moves on to the next cell providing a local velocity point.In the Lagrangian framework, a velocity reading is acquired for each successive frame pair (n, n + 1, … m where m is the last frame).
Spatially locating the velocity vectors differs slightly between the two frameworks.Another property of the 'regionprops' MATLAB ® function is the 'Extrema' property, which contains x-and y-coordinates of the object's extrema.In the Eulerian framework, the chosen extrema for the reading location in each interrogation cell is dictated by the final velocity components.In the Lagrangian framework, the final velocity components are only used on the first reading for an initialization, with locations for each subsequent reading determined in a particle path-line fashion, where u 1,n and u 2,n dictate new x n+1 and y n+1 coordinates, respectively.All vectors are spatially located by their tail end.
Similar to other image-based analyses like PIV, a minimum velocity threshold equal to the smallest detectable velocity is specified as one pixel per frame.This is implemented in the algorithm because any single pixel movement seen in the data may not be true dye movement-it could be illumination fluctuations or background camera noise.If the value of u i,n,j is equal to this threshold, that respective bounding box side velocity is set to zero, but the resultant velocity is still calculated if other bounding box side velocities are greater than this threshold.
Multiple processing steps are applied in the Eulerian framework before applying the bounding box.One is to filter out any interrogation cells that do not see a significant portion of the dye plume.We refer to this filtration as edge thresholding (ET).Any cell that does not have a specified percentage of its area filled with dye at some point in time is not used to obtain a velocity.Results of tracking the dye front in the smallscale laboratory experiment data indicate an optimal value of ET to be around 90%.Another processing step is applying an object size threshold (OST) within each analyzed interrogation cell.In some cases, discontiguous dye objects are present at a single temporal location within an interrogation cell.This complicates the tracking of the dye front because the bounding box envelops all objects within the interrogation cell, even if they are not contiguous.The OST counteracts this downfall by eliminating any small dye objects that are smaller than a specified percentage of the entire interrogation cell area.The optimal value of OST can vary between 1% and 5%, depending on the nature of the data and selected interrogation cell size.Other processing steps include eliminating dye objects that are not in contact with an interrogation cell edge (the dye object seems to appear out of nowhere, but in reality, this is out-of-plane motion and untraceable for this planar two-dimensional analysis) and eliminating any dye objects that span the entire width and/or height of the interrogation cell as this would limit the velocity readings to only one component.
A volume of dye moving through a porous medium will leave a trail of dye that can affect the bounding box velocity readings because the tail may have a different velocity than the front.To eliminate the dye tail, we develop an adaptive masking process for the Lagrangian framework in which a masking line moves and rotates with the dye front so that any dye behind this line is removed.This is the only processing step applied in the Lagrangian framework, aside from the bounding box analysis.
Rotating the masking line depends on the ratio of the velocity components of a previous reading: θ t,n is a temporary angle of the masking line for frame n and u i,n−1 is the i-th velocity component at frame n − 1.The actual angle, θ n , applied to the masking line at the n-th frame is calculated from the temporary angle, its value for the previous frame (n − 1), and a relaxation coefficient, R c : Equation ( 6) is used to smoothly rotate the masking line between frames by damping sharp changes in direction of the dye front.Optimal R c values typically range between 0.01 and 0.05 as found from automatically tracking the large-scale laboratory experiment data.Translating the masking line with the dye front is based on u i,n−1 and the time between frames, T: This masking operation requires initialization, so n begins with the second frame of the data set and the first masking line parameters, θ 1 and X i,1 , are user-specified.
Finally, a postprocessing filter can be applied to remove potential outliers for both frameworks.The filter removes any vector with a magnitude less than or greater than the mean velocity ±F standard deviations, with F being a suggested value of 2.

Automatic tracking algorithm performance evaluation
The performance evaluation of the automatic tracking algorithm compares its quantified velocity field with the known velocity field of the synthetic data and those quantified from manually tracking the dye front from the two laboratory experiments.To verify that manually tracking a dye front provides an adequate reference velocity field, the synthetic data is also manually tracked and results are compared to the explicitly specified velocity field.
Performance is quantified by analyzing the relative residuals (Res n , expressed in %) between the automated (au) and reference (ref) velocity (being either the true or manually tracked velocity field), the bias or mean error (B, expressed in %), and the standard deviation (SD, expressed in %) of the relative residuals: Res n (9) where m is the total number of velocity locations.Positive bias (B) indicates higher values for the automated method compared to the reference velocity field.
For the small-scale laboratory experiment, validating the automatic tracking algorithm with the manual velocities involves placing 32 × 32-pixel interrogation cells at each manually-tracked location.The automatic tracking algorithm then extracts velocities from the same successive frames as when the manual tracking occurred.This ensures proper validation for the automatic tracking algorithm when compared to the manually-tracked velocity field because the readings occur at spatially and temporally identical locations.We plot the automated versus reference U n and quantify the linear trend of the point cloud, along with the residual analysis.A perfect match results in all data falling on a 1:1 line and an R 2 of one.Deviation from a regression line slope of one suggests less accuracy with over-or under-estimations indicated by values greater or lower than one, respectively.Values of R 2 less than one suggest lower precision.Finally, we plot the cumulative distribution function (CDF) of U n for the automated and reference flow fields.Comparing the distributions helps identify which velocities have different statistical representations.

Porous media architecture reconstruction
We develop a solid model of the porous bed structure in a similar fashion done by Hilliard et al (2021) to show how PLIF techniques can be used to analyze a porous medium architecture and develop an intricate model that could be used in computational fluid dynamics (CFD) studies.Knowing the characteristics of the bed architecture is also necessary to better understand the experimentally measured flow fields because the reconstructed model can offer insight to certain pore-scale flow phenomenon.
After the dye injections are performed in the small-scale laboratory experiment, the porous bed is saturated with dye.The fluorescent dye offers contrasting boundaries between the solid (dye absence) and liquid (dye presence) phases of the porous medium when illuminated with the laser light sheet.These boundaries are utilized to recreate a model of the porous structure, done in the VGStudio MAX software.Images are captured every 0.1 mm in the camera's viewing axis direction by traversing the flow cell on a moving platform with 0.01 mm resolution.The entire bed reconstruction is about 24 × 24 × 14 mm.
Two perfectly spherical artificial grains are sketched in the images for model reconstruction calibration.The spheres have a radius of 0.59 mm (15 pixels) and are treated as if they exist in the real grain bed during image processing.Care is taken to place these spheres in empty pore spaces so that they would not appear to merge with actual grains.

Method performance with synthetic data
Results of the synthetic data analysis reveal the merits and limitations of the Eulerian and Lagrangian frameworks for the automatic tracking algorithm (figure 7 and supplementary material 3).These results are compared to the true, explicitly specified velocities that are used to numerically generate the synthetic data set (section 2.1.1).
A regular grid of 64 × 64-pixel interrogation cells is used to evaluate the Eulerian framework.This interrogation cell size is selected because the circle's diameter (40 pixels) fits well within an interrogation cell for the downwelling and upwelling regions.The circle's movement in these regions are tracked well, with errors of B & SD = 0.3% & 0.7% and 0.4% & 0.8%, respectively, indicating a slight overestimation of the true, explicitly specified velocity.The performance of the algorithm in the crosswelling portion is lower (B & SD = 32% & 23%) because of how the superimposed grid splits the circular dye architecture into multiple fronts.This causes some interrogation cells to only see lateral portions of the front, and while they do read the correct x-direction movement, they also detect a spurious y-direction movement due to the fringe of the dye front entering the cell.This fringe effect also occurs twice in both the downwelling and upwelling regions, but the dye front fringe grows in the xdirection rather than the y-direction.While the errors for the crosswelling region are significantly higher when compared to the other welling regions, they are still acceptable when investigating natural porous media flows because these errors indicate the tracking results are well within 100% of the true velocities.
Because of the nature of the synthetic data, the Lagrangian framework produces much better results-it does not split the circular dye architecture into multiple fronts and therefore does not encounter the movement of the dye front fringe.The downwelling and upwelling regions have negligible error and the crosswelling region gives errors of B & SD = −0.04%& 3.9%, indicating a slight underestimation of the true velocity.
The plume is also manually tracked to ensure that a manually-tracked velocity field can be used as reference data for the laboratory experiments.Results give errors of B & SD = 1.1% & 13% in the downwelling and upwelling regions, and 2.8% & 12% in the crosswelling region (see supplementary material 4).These results verify that manually tracking an advective dye front can give an accurate reference velocity field.

Method performance with small-scale laboratory experiments
Manually tracking the downwelling dye front imaged in the small-scale laboratory experiments results in 1157 spatial points where velocity is extracted.Figure 8 shows the spatial locations of 48 points, or 4% of the manually tracked points.This amount was chosen to offer a clearer view of the comparison.
Evaluating the automatic tracking algorithm with the manually tracked data in the small-scale laboratory experiments involves placing 32 × 32-pixel interrogation cells at the spatial locations where a manual reading occurs and reading the dye movement at the same temporal location (section 2.4).ET is not used because the interrogation cells are spatially located at exact locations as opposed to a regular grid.The final postprocessing filter (±F) is also not used because comparisons are desired for all velocity readings.All other processing filters are used (section 2.3), which eliminated dye objects that met the filtration criteria and prevented the algorithm from obtaining a velocity reading at some manual locations.
Values for both the manually-and automatically-tracked velocities are similar, ranging from 0.1-3.2mm s −1 to 0.3-2.9mm s −1 , respectively (figure 8).Velocity directions are also similar, with an average angle between the manual and automatic vectors of 18 • .This direction divergence stems from how the two methods (manual and automatic) differ in tracking the dye front.If the bounding box side velocities indicate movement in exactly opposite directions for the automatic tracking algorithm, the vectoral sum is taken between the two readings and determines the orientation of the final vector.In the manual approach, either the equivalent resultant velocity (U n ) is tracked, or the velocity components (i.e.u 1,n and u 2,n ) happen to be recorded as two separate movements, and the spatially closer of the two to the automatic reading is used for comparison.
Out of the 1157 manually-tracked velocities, 687 locations obtained results from the automatic tracking algorithm due to the processing steps taken in the Eulerian framework discussed in section 2.3.Comparison between the automatically and manually tracked velocities has a regression line slope of around 0.86 and R 2 of 0.75 (figure 9(a)).This slope indicates that the automatic tracking algorithm underpredicts the true velocity, and although the y-intercept is greater than zero, it is less than half the minimum detectable velocity component of the automatic tracking algorithm for this experiment (0.16 mm s −1 ).The error analysis gives B & SD = −0.5% & 33%, respectively.
The minimum U n that the algorithm can produce is 0.31 mm s −1 (two pixels per frame), and it does not have sub-pixel accuracy, as opposed to the manually-tracked velocity sets.This causes the data points to be tiered with respect to the y-axis for smaller automatic velocity magnitudes.Interrogation cells with velocities an order of magnitude greater than these small velocities are much closer to a 1:1 comparison.
Figure 9(a) shows some interrogation cells where the automatically-tracked U n is less than 0.5 mm s −1 but the manually-tracked U n is greater than 0.5 mm s −1 .This results from the automatic tracking algorithm reading only one velocity component at the same location as the manual reading, both spatially and temporally.This is due to the minimum velocity threshold.Figure 10 depicts this, where the automatic u 1,n is one pixel per frame and is therefore set to zero.Thus, the automatic tracking algorithm only reads u 2,n (0.47 mm s −1 ) whereas the manually tracked coordinates (red points in figure 10) give both u 1,n and u 2,n with U n = 0.85 mm s −1 .While this is an unfortunate pitfall of the automatic tracking algorithm, it could be mitigated by improving image resolution or increasing the size of the interrogation cells.
The cumulative distribution functions (CDF) are nearly identical, with a few small differences (figure 9(b)).One being the minimum velocity-directly attributable to the minimum velocity threshold applied to the automatic tracking algorithm-and the other being that the automatically-tracked CDF is slightly to the left of the manually-tracked CDF, indicating that the automatic tracking algorithm underpredicts the reference velocity (much like what was determined from the qualitative comparison in figure 8 and the linear regression analysis in figure 9(a)).The agreement between the two CDF curves improves as U n increases.The overall performance of the automatic tracking algorithm is discussed in section 4 below.

Small-scale laboratory experiment.
A regular grid of 32 × 32-pixel interrogation cells is used for the Eulerian framework on the small-scale laboratory experiment data (figure 11).The automatically-tracked bulk U is 0.81 mm s −1 .Its average u 1 is 0.11 mm s −1 with a standard deviation of 0.56 mm s −1 , indicating a bias for rightward (positive) movement but still containing leftward (negative) movements.Its average u 2 is 0.63 mm s −1 with a standard deviation of 0.46 mm s −1 , so the overall velocity is much stronger in the vertical direction and almost entirely downward.These magnitudes are expected as the dye plume moved downwards and to the right.Significant heterogeneity of the flow field expressed in the standard deviations is also expected because of pore-to-pore size variability and the tortuous paths around grains.One can see an upwards vector in figure 11 (magnitude of about 0.5 mm s −1 , middle-left portion of the vector field).This reading is not an error because Hilliard et al (2021) found that tortuous paths within porous media may generate locations where the local flow direction opposes the bulk direction.Vectors appear to originate from both pore spaces and within grains.The measurements are performed in a planar two-dimensional manner, but the porous medium is a threedimensional structure, and the injected dye moves as such.These three-dimensional movements, however, are not discernable in the planar two-dimensional images, but are picked up because of the finite thickness of the laser light sheet.Vectors that appear to be located within a grain are actually readings from the dye moving around a grain in a threedimensional manner.
Taking advantage of the contrast between the solid and liquid phases when the porous bed is saturated with dye and using methods described in Hilliard et al (2021), a portion of the grain bed architecture in the small-scale laboratory experiments is reconstructed using the VGStudio MAX software package (figure 12).This portion of the porous medium (close to the left glass wall of the flow cell as seen in figure 3(a)) offers the best light illumination and subsequently the highest quality contrast between the solid and liquid phases.Regions farther away from the illuminating laser see deteriorated illumination because the THV does not have a perfectly uniform RI from grain to grain, causing slight RI mismatches between the solid and liquid phases.Having the liquid saturated with fluorescent dye also hinders illumination across the flow cell due to the absorption of the light by dye molecules.Because of these unfortunate experimental barriers, contrast quality lessens across the images.Reducing the overall reconstructed volume size is done as both a preprocessing step (cropping the raw images) and a postprocessing step (cropping the model in VGStudio MAX).The sub-volume is considered to have properties representative of the entire porous bed by converging quickly to constant parameters as the sub-volume size was increased from a few pores up to several hundred.
The artificial, spherical calibration grains are reconstructed in the software, and both have a radius of 0.59 mm, exactly what they were sketched with.This model reconstruction has a void fraction of 39.6% (defined as interstitial volume over total volume; note that the artificial calibration grains were taken out of the model for this calculation), similar to what is reported by Hilliard et al (2021)-42.8%.The 3.2% difference is apparent from time spent stirring and packing the THV into the thin flow cell so that the porous medium would be tightly-knit, replicating what is used in the large-scale laboratory experiment.

Large-scale laboratory experiment.
In the large-scale laboratory experiment, the Lagrangian framework is advantageous because the dye plume appears more as a coherent structure with a uniform front rather than a multi-front plume Because the dye front in the large-scale laboratory experimental data has a significant breadth (figure 2(c)), it is difficult to choose a single representative point when manually tracking the dye front to provide a reference velocity field.Therefore, three manually-tracked data sets are obtained, with respect to the top flow line in figure 13, to quantify the variability in the manual tracking method for the large-scale laboratory experiment.Results give a variability of 20% between the three manually-tracked data sets, indicating that the reference velocity field has some inherent uncertainty when the dye front is rather large.This variability, along with lesser image resolution, apparently causes the larger bias and random errors when compared to the synthetic and small-scale laboratory experiment errors.

Method accuracy & limitations
Experimentally mapping porous media flow fields using RIMcoupled PIV techniques can give precise results (Li et al 2017, Carrel et al 2018, Hilliard et al 2021), so a maximum average bias error of 5% (from the large-scale laboratory experiment) from the proposed methodology may seem large and unacceptable.However, PIV techniques are not well suited to investigate porous media flows due to pore and grain surface contamination from seeding particles.Furthermore, validations for measured porous media flow fields are lacking in physical subsurface model applications.The study done by Reeder et al (2018) obtained subsurface flow lines of a dune structure in an open-channel flow by computational means, and while they did perform two separate modeling approaches to provide cross checks of the results, they were unable to perform experimental validation of the subsurface velocity field.Janssen et al (2012) performed an extensive experimental and multiphysics computational fluid dynamics study of coupled surface-subsurface flows, where their experimental approach to track the pore-scale flow fields depended on tracking the center of mass of the dye plumes, which were visualized by measuring the dye concentrations.The resulting subsurface flow paths are validated with simulated particle tracks, but only qualitatively and depict comparisons similar to those shown in figure 8 of this study.Therefore, while the errors presented here may seem large, the validation is a first-of-its-kind analysis, and the errors are acceptable within the water resources research community with interest in subsurface flows.
Since the synthetic data is numerically generated with specified velocities, the performance evaluation of this methodology on the synthetic data is done with the ground truth.However, the tracking performance of the small-scale and large-scale laboratory experiments is evaluated with manually tracked results.While the synthetic data analysis (section 3.1) showed that manual tracking can be done to provide an adequate reference velocity field, it still has some errors.In hindsight, established velocimetry methods such as PIV or PTV could have been used to evaluate the real porous media velocity fields presented in this study to provide a more accurate ground truth, though the seeding problems discussed in section 1 could be challenging to overcome.
All work presented here was done studying steady-state flows in a planar two-dimensional fashion, but most porous media flows are three-dimensional and can be unsteady.While tomographic flow measurement techniques are emerging, research is needed to evaluate the feasibility of the proposed methodology when paired with volumetric techniques.If it were usable with three-dimensional measurement techniques, the entire flow field in the small-scale laboratory experiments could be mapped, not just in the downwelling region where out-of-plane dye diffusion is negligible.In addition, time dependent flows could cause increased errors in the methodology due to accelerations of the dye plume front, but these errors could be mitigated by an increased data recording rate or using a regular grid of interrogation cells that are dynamically sized.Another limitation is image resolution.The automatic tracking algorithm cannot give sub-pixel accuracy, and therefore the minimum discernable velocity is directly related to the resolution of the camera.When the magnitude of the dye front velocity is near the minimum velocity threshold, the error increases considerably.Even with these limitations, the methodology presented here is novel and has a broad range of applications in the study of subsurface flows.

Conclusion
The proposed approach allows for quantitative measurements of grain bed architecture and porous media flows.The bed architecture is reconstructed with image analysis of the contrast between the presence (porous space) and absence (solid) of fluorescent dye, and flow fields are obtained by tracking the advection-dominated movement of the fluorescent dye plume front.
The flow fields reconstructed by the proposed methodology are compared with the known flow field for the synthetic data and with the manually tracked data for the laboratory experiments, since manually tracking a dye front provides an accurate representation of the true flow field as shown by manually tracking the synthetic data.Our methodology is applicable in both the Lagrangian and Eulerian framework, with the choice depending on the nature of the flow and the balance between their merits and limitations.The synthetic data analysis indicates that the Lagrangian framework is favorable when the dye plume can be treated as a single object with a single front; the bias error between the true and measured velocity is less than 0.1% with a random error of less than 4%.However, the Eulerian framework is preferred for porescale studies, such as the small-scale laboratory experiment, because the dye plume forms multiple tortuous fingers which do not allow the treatment of the plume as a single object with a single front.Applying the Eulerian framework to the smallscale laboratory experiment gives bias and random errors of 0.5% and 33%, respectively.
Whereas the large-scale laboratory experiment data is taken from a similar porous medium as the small-scale laboratory experiment, the larger-scale analysis and lesser image resolution results in a fluorescent dye plume with a single front rather than a group of individual dye fingers.Thus, the Lagrangian framework is applied to this data set, and provides a flow field along the flow path of the dye with bias and random errors of 5% and 44%, respectively.
Recent developments have enabled researchers to use machine learning methodologies to analyze porous media flow fields with neural networks (Santos et al 2020, Wang et al 2021, Alhubail et al 2022, Kashefi and Mukerji 2023).To obtain accurate flow models from neural networks, training of the network must be done with some sort of ground truth data, such as computational models or experimental data.This dye tracking methodology could be used to train neural networks, the results of which could then be used to analyze the porous media flow fields.Kashefi and Mukerji (2023) presented a point-based deep learning scheme to predict porous media flows.They successfully predicted velocity and pressure fields within a numerically generated synthetic porous medium by training their network with computational results.More specifically, they used no more than 12% of their entire point cloud that spatially defines the pore-scale flow paths to train the network (i.e.12% of their inquiry points were used as 'sensors').This approach of using sparse data to train neural networks is commonly referred to as weakly supervised learning.Since the dye front tracking methodology presented here has exhibited localized errors, the weakly supervised learning approach combined with experimentally obtained dye front tracking data could be robust for developing neural networks that predict porous media flows, and is a possible future direction for this research.
The conceptual idea behind this methodology-tracking the advection-dominated movement of a dye front-is quite fitting for porous media flows, and while some limitations do persist, the presented methodology offers a useful alternative to traditional optically-based velocimetry methods when the use of seeding particles is unacceptable.

Figure 1 .
Figure 1.Progression of a fluorescent dye plume moving through a porous medium.

Figure 2 .
Figure 2. Instantaneous images from the (a) synthetic data set, (b) small-scale laboratory experiments, and (c) large-scale laboratory experiments.The coordinate axes in (a) apply to all data sets presented in this study.

Figure 3 .
Figure 3.The small-scale laboratory experiments are done in a thin flow cell (a) (25.2 cm tall, 19.5 cm wide, 1.9 cm thick) filled with THV grains.A ruler is taped to the front glass wall to maintain approximately 1 mm of head differential between the left and right sides of the flow cell.A sketch of the experimental setup is depicted in (b).

Figure 4 .
Figure 4.The large-scale laboratory experiments are done in the (a) Aquatic Imaging Flume at the University of Idaho Center for Ecohydraulics Research.A sketch of the experimental setup is depicted in (b).The laser light sheet originated from beneath the flume and shone through a pane glass window into the sediment bed.

Figure 5 .
Figure 5.A profile view with dimensions of a model salmon redd constructed in the AIF.The entire model is made of THV grains.

Figure 6 .
Figure 6.The bounding box (red-line and dashed black-line rectangles) envelops the dye object (blue and gray shapes) in each frame.The changes in the bounding box positions and dimensions are used to calculate velocities between successive frames.

Figure 7 .
Figure 7. Velocity fields obtained by automatically tracking the synthetic data using the (a) Eulerian and (b) Lagrangian framework.The gray lines in (a) represent the regular grid superimposed on the data for the Eulerian framework.

Figure 8 .
Figure 8.Comparison of the manually-tracked (left color bar) and automatically-tracked (right color bar) velocity fields for the small-scale laboratory experiments for 4% of the manually tracked points.Each location shows velocity readings for both the manual (reference) and automatic tracking methods, comparing both magnitude and direction.

Figure 9 .
Figure 9. Comparing manually-tracked and automatically-tracked velocity magnitudes for the small-scale laboratory experimental data set via (a) linear regression and (b) cumulative distribution functions.Each data point in (a) represents one interrogation cell.

Figure 10 .
Figure 10.Tracking a dye front within a 32 × 32-pixel interrogation cell from the small-scale laboratory experiment to validate the automatic tracking algorithm.The blue lines are the bounding box sides, and the red points are where the manual readings are located.

Figure 11 .
Figure 11.Velocity field of the small-scale laboratory experiment data set acquired by automatically tracking the dye front using the Eulerian framework.Each vector represents a successful velocity reading of an interrogation cell.

Figure 12 .
Figure 12.A computationally reconstructed sub-volume of the porous medium used in the small-scale laboratory experiments (∼15 × 24 × 14 mm).

Figure 13 .
Figure 13.Smoothed velocity field of the large-scale laboratory experiment acquired by automatically tracking the dye front using the Lagrangian framework.
1269-81 Janssen F, Cardenas M B, Sawyer A H, Dammrich T, Krietsch J and De Beer D 2012 A comparative experimental and multiphysics computational fluid dynamics study of coupled surface-subsurface flow in bed forms Water Resour.Res.48 1-16 Kashefi A and Mukerji T 2023 Prediction of fluid flow in porous media by sparse observations and physics-informed PointNet Neural Netw.167 80-91 Kristoffersen A S, Erga S R, Hamre B and Frette Ø 2014 Testing fluorescence lifetime standards using two-photon excitation and time-domain instrumentation: rhodamine B, coumarin 6 and lucifer yellow J.Fluorescence 24 1015-24 Li Y, Kazemifar F, Blois G and Christensen K T 2017 Micro-PIV measurements of multiphase flow of water and liquid CO 2 in 2D heterogeneous porous micromodels Water Resour.Res.53 6178-96 Malik N A, Dracos T and Papantoniou D A 1993 Experiments in fluids particle tracking velocimetry in three-dimensional flows part II: particle tracking Exp.Fluids 15 279-94 MATLAB 1993 Image Processing Toolbox ™ Moreto J R, Reeder W J, Budwig R, Tonina D and Liu X 2022 Experimentally mapping water surface elevation, velocity, and pressure fields of an open channel flow around a stalk Geophys.Res.Lett.49 e2021GL096835 Ohmi K and Li H-Y 2000 Particle-tracking velocimetry with new algorithms Meas.Sci.Technol.11 603-16 Otsu N 1979 A threshold selection method from gray-level histograms IEEE Trans.Syst.Man Cybern.9 62-66 Quick A M, Reeder W J, Farrell T B, Tonina D, Feris K P and Benner S G 2016 Controls on nitrous oxide emissions from the hyporheic zones of streams Environ.Sci.Technol.50 11491-500 Reeder W J, Quick A M, Farrell T B, Benner S G, Feris K P and Tonina D 2018 Spatial and temporal dynamics of dissolved oxygen concentrations and bioactivity in the hyporheic zone Water Resour.Res.54 2112-28 Rubol S, Tonina D, Vincent L, Sohm J A, Basham W, Budwig R, Savalia P, Kanso E, Capone D G and Nealson K H 2018 Seeing through porous media: an experimental study for unveiling interstitial flows Hydrol.Process.32 402-7 Santos J E, Xu D, Jo H, Landry C J, Prodanović M and Pyrcz M J 2020 PoreFlow-Net: a 3D convolutional neural network to predict fluid flow through porous media Adv.Water Resour.138 103539 Stuart T A 1953 Water currents through permeable gravels and their significance to spawning salmonids, etc Nature 172 407-8 Su L K and Dahm W J A 1996 Scalar imaging velocimetry measurements of the velocity gradient tensor field in turbulent flows.II.Experimental results Phys.Fluids 8 1883-906 ToninaD and Buffington J M 2009a Hyporheic exchange in mountain rivers I: mechanics and environmental effects Geogr.Compass 3 1063-86 Tonina D and Buffington J M 2009b A three-dimensional model for analyzing the effects of salmon redds on hyporheic exchange and egg pocket habitat Can.J. Fish.Aquat.Sci.66 2157-73 Voermans J J, Ghisalberti M and Ivey G N 2017 The variation of flow and turbulence across the sediment-water interface J. Fluid Mech.824 413-37 Wang K, Chen Y, Mehana M, Lubbers N, Bennett K C, Kang Q, Viswanathan H S and Germann T C 2021 A physics-informed and hierarchically regularized data-driven model for predicting fluid flow through porous media J. Comput.Phys.443 110526 Weitzman J S, Samuel L C, Craig A E, Zeller R B, Monismith S G and Koseff J R 2014 On the use of refractive-index-matched hydrogel for fluid velocity measurement within and around geometrically complex solid obstructions Exp.Fluids 55 1862 Wirner F, Scholz C and Bechinger C 2014 Geometrical interpretation of long-time tails of first-passage time distributions in porous media with stagnant parts Phys.Rev. E 90 013025 (Some figures may appear in colour only in the online journal) Oncorhynchus tshawytscha) embryos in the Waitaki river, New Zealand N.Z.J. Mar.Freshw.Res.27437-44 DeVries P 2012 Salmonid influences on rivers: a geomorphic fish tail Geomorphology 157-158 66-74 Ershov D et al 2022 TrackMate 7: integrating state-of-the-art segmentation algorithms into tracking pipelines Nat.Methods 19 829-32 Freixa A, Rubol S, Carles-Brangarí A, Fernàndez-Garcia D, Butturini A, Sanchez-Vila X and Romaní A M 2016 The effects of sediment depth and oxygen concentration on the use of organic matter: an experimental study using an infiltration sediment tank Sci.Total Environ.54020-31 Gendrich C P and Koochesfahani M M 1996 A spatial correlation technique for estimating velocity fields using molecular tagging velocimetry (MTV) Exp.Fluids 22 67-77 Hilliard B, Budwig R, Skifton R S, Durgesh V, Reeder W J, Bhattarai B, Martin B T, Xing T and Tonina D 2023 Data for measuring porous media velocity fields and grain bed architecture with a quantitative PLIF-based technique, HydroShare (available at: www.hydroshare.org/resource/a79d513a08064ecd85f781bb9dfb642d)Hilliard B, Reeder W J, Skifton R S, Budwig R, Basham W and Tonina D 2021 A biologically friendly, low-cost, and scalable method to map permeable media architecture and interstitial flow Geophys.Res.Lett.48e2020GL090462 Hu H and Koochesfahani M M 2006 Molecular tagging velocimetryand thermometry and its application to the wake of a heated circular cylinder Meas.Sci.Technol.17