Pore-Scale Hydrodynamics in a Progressively Bioclogged Three-Dimensional Porous Medium: 3-D Particle Tracking Experiments and Stochastic Transport Modeling.

Biofilms are ubiquitous bacterial communities that grow in various porous media including soils, trickling, and sand filters. In these environments, they play a central role in services ranging from degradation of pollutants to water purification. Biofilms dynamically change the pore structure of the medium through selective clogging of pores, a process known as bioclogging. This affects how solutes are transported and spread through the porous matrix, but the temporal changes to transport behavior during bioclogging are not well understood. To address this uncertainty, we experimentally study the hydrodynamic changes of a transparent 3-D porous medium as it experiences progressive bioclogging. Statistical analyses of the system's hydrodynamics at four time points of bioclogging (0, 24, 36, and 48 h in the exponential growth phase) reveal exponential increases in both average and variance of the flow velocity, as well as its correlation length. Measurements for spreading, as mean-squared displacements, are found to be non-Fickian and more intensely superdiffusive with progressive bioclogging, indicating the formation of preferential flow pathways and stagnation zones. A gamma distribution describes well the Lagrangian velocity distributions and provides parameters that quantify changes to the flow, which evolves from a parallel pore arrangement under unclogged conditions, toward a more serial arrangement with increasing clogging. Exponentially evolving hydrodynamic metrics agree with an exponential bacterial growth phase and are used to parameterize a correlated continuous time random walk model with a stochastic velocity relaxation. The model accurately reproduces transport observations and can be used to resolve transport behavior at intermediate time points within the exponential growth phase considered.

Introduction

Porous media flows are strongly affected by the ubiquitous structural heterogeneities of porous networks. Heterogeneity stems from wide distributions of pore sizes and length scales, which induce non‐Gaussian or anomalous transport. This phenomenon is responsible for incomplete mixing or enhanced spreading, persistently spanning scales from the pore to the field (Berkowitz et al., 2000; Dentz et al., 2011; Edery et al., 2014; Gouze et al., 2008; Le Borgne et al., 2011a; Le Borgne & Gouze, 2008). At the pore scale, anomalous transport exhibits many different characteristics such as non‐Gaussian velocity distributions (Bijeljic et al., 2013; Matyka et al., 2016), high temporal correlation of Lagrangian velocities forming a spatial Markov process (Le Borgne et al., 2011b), intermittency of velocities along trajectories (de Anna et al., 2013), and superdiffusive spreading (Holzner et al., 2015; Kang et al., 2014). The intensity of the anomalous transport is related to the heterogeneity of the porous medium and has been investigated in media of different complexities, ranging from simple beadpacks to fractured sandstones and carbonates (Bijeljic et al., 2013; Meyer & Bijeljic, 2016; Morales et al., 2017; Siena et al., 2014).

The distributions of pore sizes and length scales present in porous media can have a dynamic component driven by biological or physicochemical processes, such as the development of bacterial biofilms (Seymour et al., 2004, 2007; Stoodley et al., 1994), mineral dissolution and precipitation (Daccord & Lenormand, 1987; Linga et al., 2017; Menke et al., 2015; Noiriel et al., 2013), particle filtration in particle‐laden flows (Bianchi et al., 2018; Shen & Ni, 2017), or gas exchange in multiphase flows (Datta et al., 2013; Jiménez‐Martínez et al., 2016; Kazemifar et al., 2016; Klump et al., 2007). The resulting structural changes (e.g., porosity and permeability changes) influence the pore‐scale hydrodynamics (e.g., pore‐scale velocity distributions) and alter the system's anomalous transport intensity. In this work, we focus on the clogging of pores by bacterial biofilms. These sessile bacterial communities are of interest due to their ubiquitous presence in natural or industrial systems (Costerton et al., 1999; Hall‐Stoodley et al., 2004). Their relevance regarding porous media flows is in applications such as bioremediation (Bouwer & Zehnder, 1993), microbial enhanced oil recovery (Head et al., 2003; Lazar et al., 2007), and water treatment (Gujer & Boller, 1986; Leverenz et al., 2009).

Biofilm development influences these different applications on a wide range of scales (Battin et al., 2007; Ginn et al., 2002; Morgenroth & Milferstedt, 2009). Therefore, different experimental approaches have been used to study the impact of these systems. At the field scale, the influence of biofilm development is generally monitored with geophysical techniques (Atekwana & Atekwana, 2010) or with the use of tracer tests (Li et al., 2010). At the Darcy scale, several authors considered the influence of biofilm development on macroscale properties (e.g., porosity and permeability) using soil columns (Cunningham et al., 1991; Lappan & Fogler, 1996; Vandevivere & Baveye, 1992; Taylor & Jaffé, 1990; Thullner et al., 2002). At the pore scale, the use of optical visualization techniques revealed how biofilm development altered the properties of two‐dimensional pore networks (Dupin & McCarty, 1999; Kim & Fogler, 2000; Stoodley et al., 1994).

Recently, Coyte et al. (2017), showed that the growth of biofilm patches locally reduced the size of specific pores, which significantly affected the flow field. Biofilms were experimentally (Durham et al., 2012; Nadell et al., 2017; Seki et al., 2006) and numerically (Bottero et al., 2013; Pintelon et al., 2009) observed to induce the formation of preferential flow pathways and stagnation zones in two‐dimensional (2‐D) systems. Numerically, Graf von der Schulenburg et al. (2009) combined lattice Boltzmann simulations with an individual‐based biofilm model to investigate the growth of biofilms in a 3‐D porous media. By comparing 2‐D and 3‐D results, the authors underlined the necessity of considering three‐dimensional systems, as the results in 2‐D or 3‐D differed substantially. For instance, they showed that the variance of velocity distributions increased more slowly in 3‐D than in 2‐D and the formation of preferential flow pathways was strongly delayed in 3‐D, warranting experimental data to corroborate the findings. Ultimately, biofilms were shown to cause an increase of non‐Fickian transport dynamics (Knecht et al., 2011; Kone et al., 2014; Seymour et al., 2004, 2007).

The pioneering work by Seymour et al. (2004, 2007) provided direct access to pore‐scale velocity distributions in a progressively bioclogged porous medium. Additionally, this work revealed a transition from Gaussian to non‐Gaussian transport dynamics due to biofilm development in a homogeneous beadpack. The authors used magnetic resonance microscopy to quantify the evolution of propagators (i.e., displacement probability density functions). The transition of the propagators was then modeled qualitatively for different degrees of bioclogging using a conceptual continuous time random walk based on a power law Lévy wait time and a Gaussian jump length distributions. However, this approach does not honor the temporal correlation of the intermittent pore‐scale flow, and thus does not account for the spatial Markovian nature of porous media flows which was known to be critical for accurate transport modeling (de Anna et al., 2013). The importance of this temporal correlation is expected to increase with the heterogeneity induced by increasing bioclogging (Le Borgne et al., 2011b).

The goal of this study is to experimentally quantify the influence of biofilm growth on pore‐scale hydrodynamics to better understand and model flow and transport processes in bioclogged systems. Lagrangian data from 3‐D Particle tracking velocimetry are used to determine velocity distributions, mean square displacement and particle displacement distributions of flow particles at different bioclogging stages. Then, the hydrodynamic quantities are statistically analyzed and their temporal evolution linked to expected trends of biofilm growth. Lastly, these quantities are used to parameterize a correlated continuous time random walk model that accurately captures transport observations and quantifies flow transition from uniform to preferential.

Materials and Methods

Porous Media and Working Fluid

The porous medium used in this study consisted of Nafion pellets (NR50 1100 EW, Ion Power, Munich, Germany) of diameter , with similar physicochemical properties to that of sand grains (Downie et al., 2012; Leis et al., 2005). The refractive index of the Nafion pellets was matched with a working fluid made up of 11% (w/v) glucose aqueous solution. The optimal glucose concentration was defined using a protocol similar the one presented in Downie et al. (2012). To ensure unlimited bacterial growth, the working solution was saturated with oxygen, nutrients (nitrogen, phosphorus), and electron acceptors (NaNO3, K2HPO4, and NaH2PO4⋅2(H2O) resulting in an influent solution with C:N:P molecular ratio of 1,000:1:1. In order to optimize the refractive index matching, ∼20 g of the Nafion pellets used were heated up in 250 mL of the glucose solution to 65°C for 1 h and cooled in a fresh solution overnight. This cycle was repeated three times.

Biofilm Growth and Hydrodynamic Conditions

The bacterial inoculum used in this study was isolated from the Chriesbach River (Dübendorf, Switzerland; Desmond et al., 2018). Frozen bacterial stock contained in 2 mL Eppendorf tubes was added to 100 mL of the growth medium and incubated for 20–24 h at 30°C while stirring at 200 rpm until a midlogarithmic growth phase was reached (OD600 0.52 ± 0.096). Subsequently, 10 mL of the inoculum were added to 90 mL of fresh growth medium. This incubation procedure was repeated three times for the bacteria to adapt to the synthetic carbon source of the growth medium. An effective exponential growth rate of 0.097 ± 0.012 1/h was measured by considering the OD600 (0.0505 ± 0.011) initially and once midlogarithmic phase was reached. The Nafion grains prepared following the protocol described above were added to the inoculum for the last incubation cycle to allow uniform initial bacterial attachment on the surface of the grains. A custom‐built PMMA flow cell (inner dimensions: 38 × 38 × 16 mm3) with point inlet and outlet openings at the bottom and top of the cell, respectively, was wet packed with the inoculated Nafion grains.

A biofilm was allowed to grow in the fully water saturated packed flow cell for 48 h under a constant volumetric flow rate (10 mL/min, corresponding to an initial average residence time of ∼40 s), which was set with a peristaltic pump (Ismatec, Glattbrugg, Switzerland). Note that the flux was kept constant. At this flow rate, the initial (prior to biofilm growth) Darcy velocity q was of 0.274 mm/s and the average interstitial pore velocity v was measured to be ≈1 mm/s, which yields an effective porosity of 0.27. This results in an initial Reynolds number of , where ν is the kinematic viscosity of the fluid. At the peak of observed biofilm growth, an observed threefold increase in v only raised Re to ≈1.5, which is still within Darcy's law validity range ( ). In this system, the Péclet number, defined as the ratio between advective, v, and diffusive transport rates, , was determined to be , where is the molecular diffusion coefficient of water.

3‐D‐PTV Measurements

In this work, three‐dimensional particle tracking velocimetry (3‐D‐PTV) measurements were performed in a refractive index matched (RIM; Budwig, 1994) porous medium to characterize the evolution of pore‐scale hydrodynamics during biofilm growth. The open source 3‐D‐PTV code employed here allows to quantitatively characterize different kinds of flows by providing position, velocity and acceleration of tracer particles along trajectories (Gülan et al., 2012; Hoyer et al., 2005; Lüthi et al., 2005; Michalec et al., 2015; Schmidt et al., 2016). 3‐D‐PTV was first used to study flows in porous media by Moroni and Cushman (2001). Holzner et al. (2015) and Morales et al. (2017) used 3‐D‐PTV to study the intermittent nature of flows in heterogeneous 3‐D synthetic sandstone packings. Shen and Ni (2017) used 3‐D‐PTV to study the dynamics of filtration and clogging by particulate matter. Here 3‐D‐PTV was used to perform measurements over the course of the biofilm growth, providing a quantitative description of the influence of biofilm growth on flow and transport.

3‐D‐PTV measurements were performed at time points T = 0, 24, 36, and 48 h of biofilm culturing (Carrel et al., 2018a). To perform the 3‐D‐PTV measurements, a syringe connected to the flow cell was mounted on a syringe pump (Lambda Vit‐Fit, Lambda, Baar, Switzerland). Fluorescent tracer particles (Red Polyethylene Microspheres, Cospheric, Santa Barbara, CA, USA) of diameter 70 μm, density g/cm3, and neutral buoyancy (Stokes number ), were seeded to the flow and tracked in time and 3‐D space with 3‐D‐PTV. The Péclet number of the tracers can be considered to be infinite due to the size of the particles and their corresponding negligible diffusivity. The tracer particle concentration used was of 0.02 g/L (volume fraction concentration ≈ 0.002%), to which 0.2 ppm of surfactant (Tween, Cospheric, Santa Barbara, USA) was added to prevent particle agglomeration during the 3‐D‐PTV measurements. We note that the tracer particles are too large to enter the submicrometer sized inner biofilm channels (Carrel et al., 2017; Stoodley et al., 1994) and that the particle to Nafion pellet diameter ratio is of ∼35. Thus, the effective porosity reduction of the porous medium due to interception of the tracer particles can be considered to be insignificant (Sakthivadivel & Einstein, 1970; Yao et al., 1971). Additionally, due to the extremely low volume fraction concentration, particle‐particle interactions can be neglected.

Once the flow fluid was seeded, the particles were then illuminated with a 100 W pulsed Nd‐YLF laser (Darwin Duo, Quantronix, Hamden, USA). Images were acquired with a Photron Fastcam SA5 at 50 Hz with 1 Mega pixel resolution. An image splitter was used to mimic a four camera system that permits stereoscopic viewing of the sample. A redundant four views system allows highly accurate particle recognition and tracking (Lüthi et al., 2005; Maas et al., 1993; Malik et al., 1993). For particle recognition, a preprocessing step consisting of a running image substraction with a lag of 50 images and a high pass filter were applied to improve the signal to noise ratio of individual video frames. Subsequent rejoining of fragmented trajectories was achieved with a 6D position‐velocity search that reconnects trajectory fragments when their velocity and acceleration vectors match with a small tolerance (Xu, 2008). A Savitzky‐Golay filter was then implemented in order to smooth trajectories (Saha et al., 2014). The final subpixel accuracy of the particle positions along trajectories was estimated to be of about 50 μm (Holzner et al., 2015). Short or immobile trajectories with total displacement shorter than 1 mm or 0.5 s, respectively, were discarded from the analysis. This yielded more than high quality trajectories to analyze in the bioclogged porous medium for each time point (4,465, 3,604, 1,953, and 2,660 for the time points T = 0, 24, 36, and 48 h, respectively). The velocity and Lagrangian acceleration was calculated using centered finite differences of position and velocity, respectively, along particle trajectories. Given that pressure transducers were not available in the experiment we estimated permeability indirectly as follows. We calculated the Laplacian of velocity after interpolating the velocity data onto an Eulerian grid with a spacing of 0.2 mm. This allowed estimating the pressure difference between inlet and outlet by integrating the local pressure gradient along streamlines, where ρ is the density of the fluid and D/Dt denotes the material derivative. The permeability was then estimated from Darcy's law as , where L is the length of flow cell. In the course of the experiment, the growing biofilm increasingly attenuated the signal emitted by the tracer particles. The experiment was stopped after 48 h of biofilm growth, as tracer particles could be detected well only up to 50% of the flow cell depth. In order to compensate for the increasing light attenuation and to image particle tracers over the whole flow cell depth, images were acquired from both front and back sides of the flow cell. Over the course of the experiment, the average trajectory length remained relatively constant at around 25% of the total flow cell length (26, 28, 33, and 22% for the time points T = 0, 24, 36, and 48 h, respectively). With the growth of biofilm the transparacy of the media decreased which led to a decrease in the number of reconstructed particle tracks. The RIM quality was not affected significantly. We assessed the RIM by placing a pattern of straight lines behind the test sample and by calculating the scaled root mean square error of linear fit regressions. Before biofilm growth, RIM was excellent and comparable to results in Downie et al. (2012). During bioclogging, the RIM quality decreased only slightly, i.e., from about 50 μm at T = 0 h to about 70 μm at T = 36 h. This decrease is not substantial and on the same order of the accuracy of the 3‐D determination of tracer particle locations. It therefore did not significantly affect the particle tracking. Finally, to account for entry effects at a point inlet, data at the first advective time scale τ of each trajectory were discarded. Here , where λ is the length scale of the system and is the average velocity (see Table 1).

Table 1

Average Velocity Magnitude , Velocity Magnitude Variance , Correlation Length Scale λ of the t‐Lagrangian Velocity Magnitudes, Scale Parameter v

	T = 0 (h)	T = 24 (h)	T = 36 (h)	T = 48 (h)	c_i(1/h)	Ri2
〈υ〉 (mm/s)	1.3	1.4	2.2	2.9	0.015	0.85
συ2 (m2/s2)	9.82E‐07	1.35E‐06	3.16E‐06	9.91E‐06	0.045	0.86
σln(υ)2	0.651	0.825	0.975	1.686	0.017	0.81
λ (mm)	1.38	2.1	2.1	2.7	0.013	0.92
v0	0.518	0.616	0.705	1.093	0.013	0.79
α	1.929	1.619	1.418	0.914	0.012	0.95
ϕ	0.27	0.24	0.19	0.12	0.012	0.89
ϕ/ϕ(t=0)	1	0.89	0.70	0.44	0.012	0.89
k (cm/s)	1.74	1.40	1.03	0.50	0.014	0.92
k/k(t=0)	1	0.80	0.60	0.29	0.014	0.92
Δs (mm)	0.0138	0.021	0.021	0.027

Note. c is the exponential coefficient and R is the coefficient of determination for the exponential change in time of the listed hydrodynamic parameters.

Results

Experimental Pore‐Scale Characterization of the Progressive Bioclogging

Figure 1 top row shows photographs of the flow cell during biofilm culturing (top) after 0, 36, and 48 h. The photographs illustrate how the biofilm progressively develops in the flow cell. Note that the transparency of the flow cell decreases with biofilm growth. Figure 1 bottom row shows equivalent projections of 3‐D‐trajectories (more than 106 data points for each case) color coded with the logarithm of the norm of the velocity vector.

Figure 1

Photographs illustrating (top) progressive changes in the porous media with increasing bioclogging of the flow cell and (bottom) particle trajectories obtained by 3‐D‐PTV for the time points T = 0, 36, and 48 h. The bright spot in the top right corner of the photographs is a reflection by a Teflon‐coated plastic screw that is used to close the opening of the flow cell. The trajectories are color coded with the logarithm of the norm of the velocity vector.

The trajectories display intermittent Lagrangian velocities inherent to fluid flows in porous media (de Anna et al., 2013; Holzner et al., 2015; Kang et al., 2014; Morales et al., 2017), accelerating at pore throats and decelerating at pore bodies (Holzner et al., 2015). The flow cell is operated with a constant flow rate, so that biofilm growth in the system increases the average longitudinal velocity, as illustrated by the increasingly darker color coding of the trajectories. Note that the photographs and the trajectories of Figure 1 are only qualitatively aligned for visual comparison. As larger amounts of biofilm form in the flow cell, tracer particles are confined to fewer channels in the domain. Mass conservation corroborates such channel formation as it is observed in the ∼3 fold increase in longitudinal velocity (see Table 1). Consequently, regions not sampled by the particles can be considered to be stagnation zones occupied by biofilm or newly formed stagnant zones, and thus do not contribute to carrying the flow. With increasing bioclogging, more of these stagnation zones appear. It is visible qualitatively that these stagnant zones in Figure 1 (trajectories for T = 48 h) correlate with the opaque bioclogged regions in Figure 1 (photograph for T = 48 h).

Lagrangian Velocity Probability Density Functions

In order to quantitatively investigate the influence of progressive bioclogging on porous media flow, we first consider probability density functions (PDFs) of isochronic Lagrangian velocities (t‐Lagrangian) along the trajectories (Moroni & Cushman, 2001; Moroni et al., 2007; Holzner et al., 2015; Morales et al., 2017). Under ergodic conditions and for incompressible flow, the t‐Lagrangian velocity PDF is equivalent to the Eulerian velocity PDF (Dentz et al., 2016). The PDFs of t‐Lagrangian velocities are obtained by where N is the number of trajectories, is the initial point of the ith trajectory, T is its duration. denotes the longitudinal or transverse components of the particle velocity, is its magnitude, and is the sampling interval.

The PDFs of the t‐Lagrangian velocities in longitudinal and transverse directions as well as for the velocity magnitude for the four time points during biofilm growth are illustrated in Figures 2a and 2b. As listed in Table 1, upon biofilm growth and subsequent effective porosity reduction, the average velocity in the system increases due to fluid mass conservation.

Figure 2

Probability density functions of the t‐Lagrangian (a) longitudinal velocity component and (b) magnitude obtained at different time points. (c) The autocorrelation function of the velocity magnitude, indicating the corresponding correlation lengths (characteristic length of the system) in the legend. (d) The temporal evolution of the average velocity , the variance of the natural logarithm of the velocity , the correlation length λ, and of the characteristic velocity v 0, as normalized by their value at the time point T = 0 h. The dashed line shows an exponential fit with growth rate c = 0.015 h−1. The evolution of α (with dotted line representing the corresponding exponential fit with h−1) and of the effective porosity are illustrated on the same y axis. The evolution of the normalized velocity variance is shown on the secondary y axis (with the dashed‐dotted line showing the corresponding fit,  = 0.045 h−1).

As Figures 2a and 2b show, the PDFs of the longitudinal velocity component and of the velocity magnitude exhibit a substantial increase in variance and tailing heftiness. The variance increase of the velocity magnitude (see in Table 1) reflects the formation of preferential flow paths of high velocities, as the variance is dominated by the contribution of velocities much higher than the mean . Therefore, we also report the variance of the natural logarithm of the velocity (see in Table 1), which increase nearly three fold and give insight on the increase of width toward low velocities. Higher probability of finding low velocities is a result of the formation of zones of almost stagnant flow.

Velocity Autocorrelation Functions

The spatial autocorrelation length of Lagrangian velocities has previously been identified to reflect the spatial Markovian nature of flows in porous media (Le Borgne et al., 2008). We therefore determine the Lagrangian velocity magnitude autocorrelations during progressive bioclogging. An increase in spatial correlation is indicative of persistent fast velocities along trajectories and suggests the formation of preferential flow paths. The autocorrelation functions of the s‐Lagrangian velocity magnitude (Dentz et al., 2016) sampled equidistantly along a trajectory are defined as where the velocity covariance along a single trajectory is defined by where is the velocity magnitude along the trajectory which starts at a, measured in distance s traveled along the trajectory, and L is the length of the trajectory. The velocity variance along the trajectory is given by . Finally, the correlation length, λ, is obtained by integration of .

Figure 2c shows the autocorrelation functions of the s‐Lagrangian velocity magnitudes obtained from the 3‐D‐PTV measurements for all time points and the corresponding correlation lengths obtained by integration for the different data sets. As these results show (see Figure 2c and Table 1), the correlation length measured along trajectories increases from approximately to , where 2.5 mm is the average Nafion pellet diameter. This implies that with increasing growth of biofilm, the pore structure of the medium becomes more spatially correlated.

Approximation of the Velocity Magnitude With a Gamma Distribution

Hereinafter, we consider the t‐Lagrangian velocity magnitude for all time points presented in Figure 2b. We find that the PDF of the t‐Lagrangian velocity magnitude can be approximated by the following gamma distribution:

This function captures the increasing heterogeneity in the system upon biofilm growth in the evolving function parameters α, the exponent, (probability weight at low velocities, refer to Table 1 for values for each time point) and , the scaling parameter (a characteristic velocity reflecting the increase in the average velocity). The reader is refered to Table 1 for parameter values at each time point.

Figure 3 shows the 3‐D‐PTV velocity magnitudes for all time points and the corresponding gamma distributions with shape and scale parameters estimated by maximum likelihood estimation (MLE). As listed in Table 1, the exponent α which represents the connectivity of the pore network (Holzner et al., 2015) decreases substantially. Holzner et al. (2015) presented a theoretical consideration that the maximal and minimal values of α were 2 and −2 for completely serial or parallel pore arrangements, respectively. The decrease of α observed in Table 1 from 1.929 to 0.914 therefore suggests that the structure of the pore network evolves from a predominantly parallel pore arrangement toward a more serial pore arrangement. The characteristic velocity v0 increases substantially, similar to the average velocity increase upon biofilm growth and the corresponding reduction of the void space. This finding is in agreement with other studies at the column and catchment scale, which have also observed gamma distributed velocity and residence times (McGuire et al., 2005; Wen & Li, 2017).

Figure 3

Probability density functions of the velocity magnitudes obtained with 3‐D‐PTV, the corresponding MLE gamma distributions and the PDFs of the velocity obtained from the CTRW model for all time points.

Temporal Evolution of the Hydrodynamic Quantities

The temporal evolution of various hydrodynamic characteristics of the system is presented in Figure 2d. We expect that the measured changes of these quantities mirror the exponential increase of biomass in the system over time. The growth of biomass equates to porosity reduction. This structural change, combined with mass conservation of a system at constant flux, explains the exponential increase in average pore velocity given that where q is the fixed Darcy flux and is the effective porosity. A similar exponential evolution in the velocity variance, correlation length scale, and gamma function parameters (v 0 and α) is observed and justified by the same logic. Table 1 summarizes the exponential coefficients associated with each hydrodynamic quantity following the expression , where X 0 is the value of the considered quantity at the time point T = 0 h and c is the corresponding exponential constant. The structural changes that lead to increased flow heterogeneity also affect the average flow properties through a gradual reduction of permeability, which, however, decreases more sharply than porosity (see Table 1 and compare relative values of k/k(t = 0) and ). This indicates a gradual reduction of pore connectivity due to bioclogging (Carrel et al., 2018b).

Correlated Continuous Time Random Walk

To investigate the influence of the formation of preferential flow paths and stagnation zones on transport processes, we employ a correlated continuous time random walk (CTRW) model that simulates pore‐scale particle motion (de Anna et al., 2013; Holzner et al., 2015; Le Borgne et al., 2011b). CTRW models provide an efficient modeling approach for transport in velocity fields that vary over a characteristic length scale (Dentz et al., 2016). As discussed in the previous section, s‐Lagrangian particle velocities vary over a correlation scale λ, or in other words particle velocities persist over the characteristic length scale λ, which changes according to the growth of the biofilm. We assume that the evolution of λ is much slower than the typical particle transport scales. Thus, particle motion is modeled by the following recursion relation:

Here and a space increment. The velocity series is modeled by a stochastic relaxation process (Dentz et al., 2016) to account for the correlation of velocities inherent to the spatial Markovian nature of flows in porous media (de Anna et al., 2013; Le Borgne et al., 2011b, 2008). The steady state PDF of the spatial Markov chain is obtained from through flux weighting as (Dentz et al., 2016)

Particles can maintain or change their velocity at turning points after n steps according to a Bernoulli process with persistence probability or randomly sample from the velocity distribution by according to . The single quotation mark indicates the normalization of the spatial increments (see Table 1) and of λ by the reference length scale λ (for T = 0, 24, 36, and 48). The persistence of the velocity p over the spatial increments is based on an exponential decorrelation over twice the computed correlation length in order to account for incomplete mixing at pore throats. In the following, particle transport is simulated by equation (5) for the different time points and the corresponding correlation lengths λ. The PDFs of the t‐Lagrangian velocity magnitude obtained from the correlated CTRW in Figure 3 show good agreement with the experimental data and with the fitted gamma distributions.

Displacement Statistics

The displacement s(t) at a given time t is given in terms of the CTRW model by , where . We study the mean and centered mean squared displacements obtained from the 3‐D‐PTV data, and interpret them in the light of the correlated CTRW model discussed above. This gives some insight into the heterogeneity of the porous medium and the role of biofilm growth for particle dispersion. The angular brackets denote the average over all particles.

Figure 4a shows the evolution of the first and second centered displacement moments of the 3‐D‐PTV experimental data and of the CTRW for all time points. The mean displacement is consistently well captured by the CTRW model. All experimental mean‐squared displacements (MSDs) show a transition from a ballistic regime ( ) to a superdiffusive regime with a temporal scaling greater than the expected asymptotic Fickian regime ( ). Here a transition between ballistic and superdiffusive regimes occurs at about one advective time scale. Both regimes as well as the transition are captured by the CTRW model. Note the discontinuity and interruption of the experimental data after 5 to 10 advective time scales τ due to the finite size of the experimental field of view and of the trajectory lengths. In Figure 4, the insets (b) and (c) show the evolution of the mean‐squared displacements obtained from the CTRW (b) and a magnification of the superdiffusive regime (c). The vertical shift of the MSDs with progressive bioclogging is due to the increasing variance of the velocity magnitude (see Figure 4b or refer to Table 1). The exponent of these MSDs slightly increases with bioclogging (see Figure 4c), suggesting that the formation of preferential flow pathways and stagnation zones upon biofilm growth leads to an increase of anomalous transport.

Figure 4

(a) Experimental (3‐D‐PTV) and numerical (CTRW) mean m and mean‐squared displacements obtained for all time points. The spatial moments presented here are dimensionless, as they are rescaled by the corresponding length scales (λ or ). The time is rescaled by the advective time scale τ. (b) CTRW mean‐squared displacements for all the different time points. (c) Magnification of Figure 4b illustrating a moderate increase in superdiffusive behavior with increasing biofilm growth. The continuous lines in Figure 4a or 4b and 4c, respectively, indicate Fickian scaling ( ).

Conditional Displacement PDFs

In order to further investigate how well the presented model captures transport processes, we consider the conditional displacement PDFs along particle trajectories (also known as propagators), defined as from s to for a given time lag . Figure 5 illustrates the 3‐D‐PTV and the CTRW propagators obtained for all time points in the porous medium considered. These propagators were computed for normalized time lags corresponding to 1, 2.5, and 5 advective time scales. The semilog propagators in the inset show that the transport regime remains superdiffusive for all the time lags considered for all time points, and that a Gaussian displacement distribution is not reached for the time lags defined. Interestingly, the skewness (i.e., third spatial moment) of the propagators increases with progressive bioclogging, as the hetereogeneity of the flow and anomalous transport in the system increase as well. The CTRW model performs very well for the clean porous medium (T = 0 h). For later time points, however, the increased skewness of the propagators is less well captured by the model.

Figure 5

Conditional displacement probability density functions for the 3‐D‐PTV and the CTRW data for all the different time points. The semilogarithmic insets show that the propagators never reach a fickian regime where they would follow a Gaussian distribution.

Discussion

In this study, 3‐D‐PTV measurements were performed during biofilm growth in a porous medium to quantify how biofilm growth affects pore‐scale hydrodynamics and the resulting transport behavior. The growth of biofilm induced an increase of anomalous transport in the considered system. A correlated continuous time random walk model parameterized by the gamma distributed pore‐scale velocity magnitudes accurately captures the mean and mean‐squared displacement changes, and qualitatively reproduces displacement probability density functions with lesser skewness.

The average velocity and further hydrodynamic quantities follow an exponential evolution with biofilm growth, which for all quantities except the variance can be described by an exponential time scale of approximately 0.015 h−1. We conjuncture that this behavior is due to the exponential growth of the bacteria in the considered system. During incubation, the effective growth coefficient of the bacteria was measured to be of approximately 0.09 h−1, which is about one order of magnitude larger than the exponential coefficients observed for the hydrodynamic quantities. However, this difference could be attributed to the fundamental differences in growth kinetics of planktonic and sessile bacteria. Additionally, detachment of bacteria might also have induced a lower exponential growth coefficient. It is logical to expect that the effective porosity, and consequently the average velocity, will increase at a rate that is imposed by the increase of biomass in the porous medium. Similar obtained rates for geometrical parameters of the medium (correlation length λ and connectivity parameter α) are less trivial to justify physically. We nevertheless propose a similar logic to explain their temporal evolution, given that these characteristics depend on average pore‐scale geometry (i.e., length scale and pore arrangement) and therefore are proportional to the amount of biomass produced. Differently, the rate of change in velocity variance is found to be about three times higher with biomass growth. The large increase in variance can be justified by the patchy growth of biomass throughout the porespace. This induces blockage to flow in certain pores (which governs transitions in slow velocities) and develops preferential paths in others (which governs the high velocities). The inclusion of much slower and much faster observed velocities is captured in the increase of distribution variance, which is a quantity independent of mean velocity.

In the literature, the flow velocity distribution is often discussed separately for two distinct regimes, i.e., the slow flow zones associated with stagnant zones and the fast velocities or preferential pathways carrying most of the fluid transport. The question on a mechanistic understanding between pore‐scale organization of the media and the resulting flow velocity distribution for both regimes has been addressed since many years and still remains a challenge. Saffman (1959) used Poiseuille's law as an approximation of the flow in individual pores to relate velocity and pore diameter, which gives a map of pore diameter to pore velocity. This approach was used in Holzner et al. (2015) to obtain the distribution of pore velocities according to a Gamma function, which implies a power law behavior for small velocities and exponential decay of fast velocities that is consistent with their experiments and other previous studies in heterogeneous media (see references in Holzner et al., 2015). de Anna et al. (2017) recently succeeded in deriving a direct relationship between the pore size distribution and the velocity distribution for the slow velocities. Namely, they showed that the velocity distribution follows a power law where the power is half of the power given by the pore size distribution. The fast velocities are given by the conductive backbone and its distribution varies from medium to medium (e.g., Matyka et al., 2016). Alim et al. (2017) presented an analytical model, where the local correlations between adjacent pores, which determine the distribution of flows propagated from one pore downstream, predict the high tail of the flow distribution. They showed that, starting from a highly ordered bead pack and increasing the level of disorder, the fast velocities transitioned from Gaussian toward exponential distribution. In summary, at the pore scale, the distribution of pore sizes controls the low tail (de Anna et al., 2017) and correlations of fluxes between adjacent pores control the high tail of the velocity distribution. In our work, the growing biofilm reduces pore spaces in slow velocity zones which will increase the probability of small radii and consequently the low tail of the velocity PDF. Given that low velocity zones become less conductive, more and more flux is carried in faster velocity regions, which presumably increases correlations between subsequent pores along preferential pathways. This nonlocal effect seems thus at the root of the intensification of the high tail of the velocity PDF we observed here in the bioclogged media.

As the porespace is reduced upon biofilm growth, the velocities in narrowing pores increase on average. However, we also observe an increase of the probability of low velocities. This increase can be attributed to the narrowing pore size as according to Poiseuille's law for flow through a pipe (an idealization of a single pore throat). In such a case, the maximum velocity is proportional to the pore diameter squared for a constant pressure gradient. The increase of probability of low velocities due to the narrowing pore size suggests the formation of almost stagnant zones, which marginally contribute to the total flow. Note however that locally, neither the mass flux nor the pressure gradient are fixed and therefore, the structural changes caused by the biofilm growth induce an increase of the variance of the velocity distributions, whose PDFs become increasingly non‐Gaussian. Bijeljic et al. (2013), Datta et al. (2013), Siena et al. (2014), Matyka et al. (2016), and Meyer and Bijeljic (2016) made similar observations for porous media with different pore network characteristics (effective porosity, tortuosity, pore size distribution, connectivity, etc.), showing that the non‐Gaussianity of the velocity PDFs was caused by the structural heterogeneity of the porous medium considered or by different degrees of saturation for partially saturated porous media. Additionally, several authors have also reported increasing exponents of the superdiffusive regime of the MSDs for porous media or fracture networks exhibiting increasing degrees of heterogeneity (Kang et al., 2017; Meyer & Bijeljic, 2016).

Note that porosity‐permeability relationships, which have been frequently studied in the literature (Cunningham et al., 1991; Lappan & Fogler, 1996; Taylor & Jaffé, 1990; Thullner et al., 2002; Vandevivere & Baveye, 1992), refer to the relation of the average flow velocity to the average pore volume. As shown above, the observed (anomalous) transport behavior can be attributed to the low end of the velocity PDF and cannot be understood by considering the mean velocity alone. Thus, in order to predict the influence of biofilm growth on flow and transport behavior it would be desirable to determine the relation between flow velocity and pore size distribution, and between pore size distribution and biofilm growth. Velocity models based on Poiseuille's law (de Anna et al., 2017; Dentz et al., 2018; Holzner et al., 2015; Saffman, 1959) relate the behavior of the velocity PDF at low velocities to the distribution of small pore sizes and provide possible starting points for such an analysis. Poiseulle's law would imply that uniform biofilm growth leads to a rescaling of the low end of the velocity PDF. Figure 3, however, shows that the PDF at low velocities does not follow a simple scaling for growing biofilm. An analysis of the relation of biofilm growth and the evolution of the velocity PDF along these lines is beyond the scope of this paper and postponed to future work.

The results presented here were obtained for four different time points during the early stages of the bioclogging process. Due to the observed exponential evolution of the hydrodynamic quantities in the bioclogged system, the CTRW presented here could be parameterized for any time point whithin the duration of the experiment (i.e., the bacterial exponential growth phase). Exponential growth of biofilms in porous media has been observed previously (Bottero et al., 2013; Graf von der Schulenburg et al., 2009; Hassanpourfard et al., 2015; Pintelon et al., 2012), and studied both experimentally and numerically. The work by Bottero et al. (2013) and Hassanpourfard et al. (2015) found that the exponential growth phase is finite, and that biofilm growth asymptotically approaches a steady state (in terms of biofilm surface coverage or effective porosity, medium permeability, etc.). Bottero et al. (2013) pointed out that steady state biofilm growth is the result of a balance between biomass decay/detachment and growth. The transition from exponential to steady state growth observed by Bottero et al. (2013) and Hassanpourfard et al. (2015) can be described by a logistic model. In the present study, it was not possible to collect data at long enough temporal scales to observe the transition from exponential growth to steady state.

Interestingly, the propagators obtained in Seymour et al. (2004, 2007) exhibited a substantially stronger peak for very small displacements than in the data presented in Figure 5. This could be partially attributed to the Péclet number of our system, which was computed from 3‐D‐PTV data and can be considered to be infinite. Simply, the diffusion experienced by the tracer particles is negligible, meaning that diffusive displacements are not captured in the data presented here. Another aspect lies within differences intrinsic to the experimental approaches used in both studies. Due to the limitations of the 3‐D‐PTV method used in this work, we are only able to sample flow in pores corresponding to the effective porespace and outside the biofilm structure. As a result, completely stagnant or dead end pores, which do not participate in flow and are not filled with biomass are not sampled. Similarly, the tracer particles used are too large to penetrate the biomass itself, thus limiting data collected in very low flow regions. In contrast, the MRM method used by Seymour et al. (2004) samples flow within the biofilm as well as in all porespaces. This suggests that for low Péclet numbers common in porous media, the permeability of the biofilm itself could affect the transport of solutes (Davit et al., 2013; Deng et al., 2013; Dupin et al., 2001; Ezeuko et al., 2011; Pintelon et al., 2012; Thullner & Baveye, 2008). This has been shown for permeable benthic biofilms, which substantially contributed to transport and reaction processes in gravel bed streams (Aubeneau et al., 2015, 2016; Li et al., 2017).

Finally, we examine how structural changes in the pore geometry due to bioclogging affect the medium permeability. Different relationships between porosity and permeability of the media have been proposed (e.g., Seki & Miyazaki, 2001; Thullner, 2010; Vandevivere & Baveye, 1992). These relationships are either empirical or based on growth models of biofilm in pores (either single pores or simplified pore networks), see Thullner (2010) and references therein. The relationships based on growth models can be subdivided into models that assume that biofilms grow homogeneously or in the form of patches. Given that the small pores are clogged first, one might expect to observe a significant change in the flow's heterogeneity, as expressed by the decrease of the shape parameter α, while the permeability is supposedly less affected because highly conductive pores are still open. Changes in permeability would become significant only when the porosity reaches a certain threshold corresponding to clogging of highly conductive zones by biofilm growth. Threshold‐based permeability‐porosity models along these lines have been proposed in Bernabé et al. (2004) for binary mixtures of high‐permeability and low‐permeability materials and in Thullner (2010) for bioclogged media. We note that our measured decrease in the shape parameter α is similar to the relative porosity decrease, while the permeability decreases faster (Table 1). Thus, even though permeability decreases quite sharply, our data are not consistent with threshold‐based models, but rather with previously observed trends for permeability decreases in porous media coated with patchy biofilms (e.g., Graf von der Schulenburg et al., 2009; Thullner, 2010). Addressing the porosity‐permeability relationship more conclusively, however, is beyond the scope of this work.

Conclusion

In this work, we successfully monitored the transition from normal to anomalous transport in a progressively bioclogged three‐dimensional porous medium. We observed an exponential evolution of different hydrodynamic quantities consistent with exponential growth kinetics of the biofilm. A two parameter gamma function provided a remarkable description of the bulk and the high tail of the velocity distribution. The shape parameter α is considered to represent the connectivity of the medium and accounts for a shift from a mainly parallel arrangement of pores toward a more serial one. This too is consistent with the observed formation of stagnation zones and preferential flow pathways. A correlated continuous time random walk model that employs stochastic relaxation of velocities is able to capture accurately the evolution of the first two spatial moments, but is somewhat limited to reproduce the third moment. Altogether, these insights shed light on the processes behind and provide useful modeling tools to capture the effects of growing biofilm forming bacteria on mass transport processes in porous media. We anticipate that the approach introduced here could also be applied to model transport in other porous media where the heterogeneity is either inherent to the original pore structure or develops dynamically as a consequence of physicochemical processes.