Study on flow distribution pattern and conductivity of porous media in bioretention cells

ABSTRACT

To evaluate the long-term performance of bioretention cell (BRC), a study was undertaken to assess the flow distribution and conductivity. Despite initial conductivity of the original medium being the common predictor of hydraulic performance, most of the BRCs are affected by conductivity variations during actual operation. This happen due to the fact that microbial behavior plays an important role in the conductivity variations. This linkage may occur when bacteria as inert colloids transports between particles and biodegrades dissolved pollutants, either promoting or retarding flow distribution and conductivity in BRC. Flow distribution was determined by numerical simulation and tracer test, and the correlation between conductivity and flow distribution was revealed by conductivity experiment coupled with flow distribution analysis. Results revealed a non-uniform flow distribution in BRC, and seepage flow in submerged zone was virtually impossible push flow. Conductivity had an inversely proportional relationship with hydraulic efficiency where hydraulic efficiency reached the highest value (0.297) under a low hydraulic conductivity (0.000107 m/s, approximately K/K_ini = 0.79). Primary cause of hydraulic capacity reduction was the initial permeability decrease due to medium structure changes. Results revealed a sharp upward trend followed by a slight decrease, and then, stabilized to a stable infiltration stage. Permeation process of sewage influent was similar to the one of potable water where the permeability reduced to 0.000102 m/s after 450 h and declined continuously. Thus, it is clear that flow distribution and conductivity in bioretention must be estimated more accurately on a microscopic scale.

KEYWORDS:

• Bioretention Cell
• Flow distribution
• Microbial behavior
• Conductivity
• Numerical simulation
• Tracer test

1. Introduction

Hydraulic conductivity (K) is the most important parameter to evaluate operating performance and lifespan of bioretention cell (BRC). The design and maintenance requirements of BRC are directly affected by K value; moreover, K also affects the decontamination capability of BRC. For example, the down flow bioretention (DFBRC) with a lower K can easily lead to more frequent overflows, therefore, decreasing the annual mass of pollutants captured by the system. Microbial community plays an important role in the ecological treatment system [1]; meanwhile, analysis of the experimental data reveals a clear relationship between the spatial structure of particles and the ability of biomass to reduce its hydraulic conductivity [2]. The microbial clogging of ecological treatment system [constructed wetland and BRC] has been an object of intense study in the last few decades. Bergman et al. (2011) [3] suggested a final K between 30% and 70% of the initial values for the infiltration trenches after 2 years and 9 months of operation. Consequently, it is critical that these potential decreases in K are taken into account when designing systems.

At present, the design of BRC is commonly based on the initial hydraulic conductivity [K_ini]. Most of the researchers in system design has confined to the macroscopic K_ini, which is still a research method of ‘black box’. Costumer et al. (2009) [4], undertook a review of 37 biofilters in Australia and found that a large number had measured hydraulic conductivity of around 25–50% of their initial values. According to Wilson et al. (2015) [5], K_ini reduction is attributed to soil compaction during basin construction. Therefore, as the most important design parameter, K_ini is only an initial value based on the packing ratio of the medium, and it could be changed during operation. In addition, even if the media used in the DFBRCs have the similar K_ini, but results still revealed different trend for the permeation performance due to the difference in the flow disorder during operation; the particles compaction degree will be different under various influent loads, which result in the system that cannot achieve the desired efficacy. Therefore, the primary focus of these studies was the macro-performance rather than the process identification. To open ‘black box’, its indispensable approach from the micro scale, the flow distribution must be simulated and analyzed on a porous medium, The formation, decomposition, and exploitation of hydrates are constrained by the state of aqueous interfaces, and the state and behavior of interface water is remarkably different from the well-known behavior of bulk water due to the interface effects [6]. Which considering a spatially varying nature at micro scale, which leads to some inabilities in describing some important interfacial phenomena and properties [7]. Meanwhile, the effect of flow distribution on the conductivity coefficient must be revealed, and thus, the hydraulic conductivity is more accurately estimated during operation a non-uniform flow distribution in BRC and seepage flow in the submerged zone was virtually impossible push flow. The flow distribution simulation in porous mediums in the literature can been seen in Table 1.

Table 1. The flow distribution simulation in porous mediums

Application	Material	Simulation	Porous media -simulation	Reference
Fluid saturations	Sand, oil and water	Lattice Boltzmann Method	Input structure and visual observations	[8]
Permeability of porous materials	Bentheimer sandstone,	Lattice Boltzmann Method	Input structure	[9]
Permeability of filter cakes	Glass beads	Lattice Boltzmann Method	Input structure	[10]
Permeability and tortuosity of porous materials	Glass beads	Lattice Boltzmann Method	Input structure	[11]
Permeability, tortuosity, and packing of porous materials	Lava	Random Walk Simulations	Input structure	[12]
Permeability and tortuosity of porous materials	Glass beads	Random Walk Simulations	Transport properties by tracking ion	[13]
Permeability and tortuosity of porous materials	Sandstone	Random Walk Simulations	Comparison with medical CT	[14]
Bedstructure coordination number and packing distribution	Glass beads	Optimization method	Homogeneity and isotropy	[15]

Numerical simulation has become an important tool for the design and operation of water treatment reactor, the conventional design is based on ideal plug flow reactor (PFR) model, ideal completely stirred reactor (CSTR) model and tank in series (TIS) model [16] with significant differences in their performances. For example, the flow characteristics of the constructed wetland are usually simulated using the ideal PFR model and ideal CSTR model; among these, ideal PFR model was recommended by USEPA and WPCF [17,18], to describe the flow regime in wetlands. However, the tracer concentration would emerge a ‘tail’ after reaching its peak in the actual tracer concentration-time curve. Therefore, the ideal PFR CSTR models cannot completely simulate the actual asymmetric tracer concentration-time curve. The degree of actual reactor flow regime deviating from the classical ideal reactor flow regime is very high.

In particular, the biochemical process can be affected by the reactor flow conditions because the biomass, substrates, and inhibiting compounds can be distributed in different reactor zones; such as dissolved oxygen as an inhibiting compound will inhibit the nitrification reaction when the content is low, and many refractory organics often have their own toxicity and inhibit the growth of microorganisms. This implies that the biochemical process can occur with a different kinetic depending on hydrodynamic conditions [19]. The microbial growth is mainly related to the biomass formation in the filler system, and diffusion and convection provide a powerful condition for nutrient transportation for microbial growth, whose availability therefore can also be influenced by the liquid flow. However, biomass detachment means removal of biofilm pieces, which can be carried away by the flow [20]. Inhibition could also happen due to the accumulation of some inhibitory compounds in specific reactor zones including dead angle, short-circuit flow, and cake ditches. Residence Time Distribution (RTD) is a common analysis method for describing the reactor hydraulic behavior and exploring the existing flow deviations from ideal flow [21]. RTD curve can be used to evaluate the infiltration phenomenon and the degree of mixing rate in bioreactors [22].

The hypothesis of this study was the four-stage theory of ‘decrease-increase-then decline-stability’ of BRC seepage coefficient. Therefore, current work was aimed to provide an overall assessment of the flow distribution and hydraulic conductivity behavior of BRCs. The flow distribution was determined by numerical simulation and tracer test, and the correlation between conductivity and flow distribution was revealed by conductivity experiment coupled with flow distribution analysis results. Current study also aimed to investigate the hydraulic conductivity behavior can be judged more accurately on the microscopic scale, which can provide more accurate control parameters for practical applications.

2. Materials and methods

2.1. Experimental setup

For current study, two DFBRCs were constructed from nonopaque Perspex. The inner dimensions of one reactor were 500 mm (L) × 800 m (H) × 300 mm (W). The inner walls of the reactors were sandblasted to prevent preferential flow along the edge (Figure 1). The inlet PVC tubes (20 mm in diameter) were attached on the top of DFBRC, and the orifice outlets (50 holes, 6 mm) were placed at the bottom of DFBRC and extended into a vertical riser plate. By overflowing the plate, the water level could be kept constant at the bottom of the cell, thereby creating the submerged zone (SZ).

Figure 1. Schematic diagram of experimental unit (seven sampling points: No. 1, No. 2, No. 3, No. 4, No. 5, No. 6, No. 7)

In the filter media, six sampling pipes (No. 1, No. 2, No. 3, No. 4, No. 5, and No. 6) were installed at 100 and 200 mm from the bottom and 125, 250, and 375 mm from the left wall, and one sampling pipe (No. 7) was installed at 50 mm from the bottom in the effluent trough (Figure 1). These pipes were made of 6 mm rubber tubes traversing the entire cell diameter to ensure that the samples were representative of the entire filter.

The design criteria of DFBRCs is according to the Australian FAWB adoption guidelines [23]: a top gravel cover layer (125 mm) to prevent the washout of the overlying soil, a filter layer of washed sand (300 mm, diameter 0.15–1.00 mm), a sand submerged layer (200 mm, diameter 0.25–0.50 mm), following transition layer of coarse sand (50 mm, diameter 0.50–1.00 mm) and the bottom gravel drainage layer (50 mm, diameter 3–5 mm) with an orifice outlet, which was connected to the effluent trough. Consequently, DFBRC had an initial hydraulic conductivity of around 0.135 mm/s with 0.1387 porosity. The reactors were kept without any vegetation in order to assess the influence of compaction (due to hydraulic loading) on hydraulic conductivity. Take great care to waterproof seals at connection points. Use collars on outlet pipes at the point where it traverses the wall. This can be tricky, especially to achieve compaction around the seal. Alternatively it is feasible to use shotcrete to create a large collar extended across the basin surface. A filter fabric can be used around the top of inlet pits and underneath inlets and sediment forebays to prevent preferential flows underneath and down the sides, where the structures are embedded below the filter media surface [23].

2.2. Experimental procedure

Two lab-scale reactors were placed in our laboratory located in Nanjing, China, to ensure that the only inflow water received was by controlled dosing, and these were then flushed 3 times during 3 days with deionized water, to promote natural biofilm development. The design of the experiment was aimed at replicating the real time conditions as much as possible while ensuring the degree of control and repeatability necessary to elucidate the impacts of the experimental parameters.

Lab-scale experiments were conducted over 30 days (from July 1^st, 2016, under summer conditions) in triplicate; two reactor tests (‘Test A’ and ‘Test B’) were designed to emulate the flow distribution and measure the hydraulic conductivity. Test A was performed by tracking the flushing of deionized water dosed with sodium chloride tracers for describing the flow distribution under four hydraulic loading rates (0.5, 1.0, 1.5, and 2.0 m³/m²·d). Contrast to Test A, Test B was continually dosed with synthetic sewage and Test B was packed and manually compressed to the same bulk density as that in Test A. To ensure consistency in the inflow and avoid logistical issues in collecting large swage Test B, synthetic sewage was used during the entire Test B experiment with C/N = 2:1. Synthetic sewage was made from glucose, NH₄Cl, K₂HPO₄, NaHCO₃, FeCl₂ · 4H₂O, humic acid and trace elements to simulate real domestic wastewater, as follows: chemical oxygen demand (COD): 150 mg/L; total suspended solids (TSS): 40 mg/L; total nitrogen (TN): 64 mg/L; total phosphorus (TP): 4 mg/L; Sodium Bicarbonate (NaHCO₃): 90 mg/L; iron (Fe): 10 mg/L; cobalt (Co): 1 mg/L and nickel (Ni): 0.5 mg/L. Subsequently, to measure the change in K with time under 24 different hydraulic loading rates (HLR), from 0.4 m³/m²·d in Test B, the HLR was increased by 0.1 m³/m²·d at intervals of 48 h until HLR reached 2.5 m³/m²·d, immediately, the HLR was raised to 3.0 m³/m²·d after 48 h, and then, the HLR was boosted to 3.5 m³/m²·d after 48 h. The annual mass of pollutants captured by the system was determined by the system treatment capacity, which is the joint action of filler filtration, adsorption and biological reaction under the corresponding HLR, which is subject to the pollutants migration and transformation mechanism in porous media.

2.3. Flow distribution analysis method

2.3.1. Mathematical model

A comprehensive, two-dimensional model was developed for DFBR (Figure 1) and considered using a commercial code Fluent 6.22 (Fluent Inc., USA). Porous media were modeled by adding a momentum source term (S_i) to the standard fluid flow equations (De boer and Didwania, 2004) [24]. Thus, the momentum balance in the porous media could be defined as Equation (1)(1) $\rho \mathrm{v}\cdot \mathrm{\nabla }\mathrm{v}=-\mathrm{\nabla }\mathrm{p}+\mathrm{\nabla }\cdot \left[\mu \left(\mathrm{\nabla }\mathrm{v}+{\left(\mathrm{\nabla }\mathrm{v}\right)}^T\right)\right]+S_i$(1) .

(1) $\rho \mathrm{v}\cdot \mathrm{\nabla }\mathrm{v}=-\mathrm{\nabla }\mathrm{p}+\mathrm{\nabla }\cdot \left[\mu \left(\mathrm{\nabla }\mathrm{v}+{\left(\mathrm{\nabla }\mathrm{v}\right)}^T\right)\right]+S_i$(1)

where ρ is the density of the liquid, v is the vector velocity of the liquid (m/s), p is the static pressure (Pa), and μ is the viscosity (Pa·s). S_i was composed of two parts: a viscous loss term (the first term on the right-hand side of Equation (2)(2) $S_i=-\left(\sum _{j=1}^3{D}_{ij}\mu v_j+\sum _{j=1}^3{C}_{ij}\frac{1}{2}\rho v_{mag}{v}_j\right)$(2) ) and an inertial loss term (second term on the right side of Equation (2)(2) $S_i=-\left(\sum _{j=1}^3{D}_{ij}\mu v_j+\sum _{j=1}^3{C}_{ij}\frac{1}{2}\rho v_{mag}{v}_j\right)$(2) ):

(2) $S_i=-\left(\sum _{j=1}^3{D}_{ij}\mu v_j+\sum _{j=1}^3{C}_{ij}\frac{1}{2}\rho v_{mag}{v}_j\right)$(2)

where S_i is the source term for the i^th (x, y, or z) momentum equation, and D and C are the prescribed matrices. In case of simple homogeneous porous media:

(3) $S_i=-\left(\sum _{j=1}^3\frac{\mu }{\alpha }{v}_i+\sum _{j=1}^3{C}_2\frac{1}{2}\rho v_{mag}{v}_i\right)$(3)

where 1/α is the viscous loss coefficient (m⁻²) and C₂ is the inertial loss coefficient (m⁻¹).

In case where modeling a laminar flow through the porous media in this study which was similar to a packed bed, the viscous loss and inertial loss coefficients in each component direction could be identified as (De boer and Didwania, 2004) [24]:

(4) $1/\alpha =\frac{150}{d_p^2}\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^3}$(4)

(5) $C_2=\frac{3.5}{d_p}\frac{\left(1-\epsilon \right)}{{\epsilon }^3}$(5)

where d_p is the mean particle diameter (mm), and ε is the void fraction.

The computational grid files of the two-dimensional model were generated by Gambit 2.2.30 [24,25]. The velocity inlet was adopted as inlet boundary conditions, and the pressure outlet was adopted as outlet boundary conditions. Interior boundary conditions were used for the filter layer, and fluid boundary conditions were used for the internal unit area, whereas the surface wall was adopted as the surrounding boundary conditions. The simulation conditions and the geometric parameters of the model are presented in Figure 1 and Table 2.

Table 2. Geometric parameters and simulation parameters of the layered bioretention model

Parameters	Cover layer	Distribution layer	Filter layer	Submerged layer	Transition layer	Drainage layer
Height (mm)	25	100	300	200	25	25
Mean diameter (d1, mm)	12	6	2	2	6	15
Void fraction (ε1, %)	44.5	44	30	33	37	45
Viscous resistance coefficient (1/α, m^−2)	3.641 × 10^6	1.534 × 10^7	6.806 × 10^8	4.684 × 10^8	3.265 × 10^7	2.213 × 10^6
Inertial resistance coefficient (C2, m^−1)	1.837 × 10^3	3.834 × 10^3	4.537 × 10^4	3.263 × 10^4	2.684 × 10^3	21.408 × 10^3

2.3.2 Tracer test method

Tracer tests were conducted using sodium chloride (NaCl). The density of NaCl solution was far greater than deionized water. The effluent electrical conductivity (EC) was tested and recorded for RTD calculation. 300 mL of NaCl solution (250 mg/L and 1000 mg/L, respectively) was added into volumetric flasks with different kinds of sand and gravel. Electrical conductivity had no significant change both within 12 and 24 hours and the concentration of NaCl was linear with EC (R² = 0.99878). For each trial, the box was first saturated by deionized water piped from a water storage tank. Background deionized water’s EC was recorded for each trial. Tracer tests were started with a fast pulsing injection of NaCl solution (100 mL). At the interval of 10 minutes, 50 mL of sample was collected from the effluent. EC of each sample was tested by a conductivity analyzer (DDS-307A, Shanghai Precision Instrument Co., Ltd., Shanghai, China). Calibration was made by testing the reference EC at a known NaCl tracer concentration and water temperature. Linear regression was made to guarantee that more than 90% of tracer concentration was captured for RTD curve in each trial. It was measured with NaCl at the interval of 10 minutes, until the conductivity was not displayed. The total amount of collected effluent was the total amount of captured pollutants and obtained the total amount of captured pollutants.

Earlier studies noted that based on the fluid reactor theory, the concentration measured in pulse tracer experiments was equal to the hydraulic residence time distribution density [26]. The measured conductivity was normalized according to Equation (6)(6) $N\left(t\right)=\frac{\left(E\left(t\right)-E_w\right)M_{NaCl}Q}{\left({\lambda }_{N_a}+{\lambda }_{Cl}\right)M}$(6) :

(6) $N\left(t\right)=\frac{\left(E\left(t\right)-E_w\right)M_{NaCl}Q}{\left({\lambda }_{N_a}+{\lambda }_{Cl}\right)M}$(6)

where t is the tracer injection time (h), N(t) is the normalized hydraulic residence time distribution density (h⁻¹), E(t) is the effluent electrical conductivity (s/m), E_w is the background water conductivity (s/m), Q is the inflow rate (m³/h), M_NaCl is the molar mass, 58.44 g/mol, λ_Na and λ_Cl are the conductivities of Na⁺ (5.01 × 10⁻³ (s·m²)/mol) and Cl⁻ (7.63 × 10⁻³ (s·m²)/mol), respectively, and M is the total amount of tracer (g).

The mean residence time (t_m) and the standard deviation for the average time (σ) of the liquid in DFBRC were obtained directly from the particle trajectory model. To compare RTD in different constructions or dissimilar conditions, each parameter must be normalized. Normalized retention time (t_θ) and normalized variance (σ_θ²), which were derived from t_m and σ to estimate the hydraulic performance, were, respectively, normalized by Equation (7)(7) $t_{\theta }=\frac{t_m}{t_n}$(7) and Equation (11)(11) ${\sigma }_{\theta }^2=\frac{{\sigma }^2}{t_n^2}$(11) . However, neither t_θ nor σ_θ² were adequate to compare the variable designs. Then, the assessment of the effects of short circuit flow and mixing ability on hydraulic performance by hydraulic efficiency [λ) was introduced by Persson et al. (1999) [27] as shown in Equation (12)(12) $\lambda =\frac{t_p}{t_n}=t_{\theta }\left(1-{\sigma }_{\theta }^2\right)$(12) . Here, λ is the ratio of the peak hydraulic retention time (t_p) and the nominal residence time (t_n). The hydraulic efficiency of DFBRC was high in the higher λ values.

(7) $t_{\theta }=\frac{t_m}{t_n}$(7)

(8) $t_m=\underset{0}{\overset{\mathrm{\infty }}{\int }}tN\left(t\right)dt/\underset{0}{\overset{\mathrm{\infty }}{\int }}N\left(t\right)dt$(8)

(9) $t_n=\frac{V}{Q}$(9)

(10) ${\sigma }^2=\underset{0}{\overset{\mathrm{\infty }}{\int }}{\left(t-t_m\right)}^{2}N\left(t\right)dt/\underset{0}{\overset{\mathrm{\infty }}{\int }}N\left(t\right)dt$(10)

(11) ${\sigma }_{\theta }^2=\frac{{\sigma }^2}{t_n^2}$(11)

(12) $\lambda =\frac{t_p}{t_n}=t_{\theta }\left(1-{\sigma }_{\theta }^2\right)$(12)

where t_θ is the normalized retention time, t_m is the mean residence time (h) and defined as the gravity position of the residence time distribution curve, t_n is the nominal residence time (h) and defined as the ratio of pore volume and inflow rate, σ² is the variance of the RTD (h²), which represented the dispersion range of the tracer concentration curve relative to the distribution mean, t_p is the peak hydraulic retention time (h) and defined as the effluent time of the tracer maximum concentration, and σ_θ² is the normalized variance that represented the hydraulic distribution divergence.

2.4. Hydraulic conductivity

To measure the change in hydraulic conductivity over the experimental time, conductivity measurements were conducted in one month interval. The hydraulic conductivity of each reactor was calculated by measuring the outflow volume per unit time using the standard constant head method (ASTM International D2434-68, 2006) as per Darcy’s Law.

3. Results and discussion

The hypothesis of this study was the four-stage theory of ‘decrease-increase-then decline-stability’ of BRC seepage coefficient. To evaluate the long-term performance of bioretention cell (BRC), a study was undertaken to assess the flow distribution and conductivity.

3.1. Effect of hydraulic loading rate on flow distribution

The static pressure and the stream function under different hydraulic loading rates were obtained through fluent simulation. As shown in Figure 2, the pressure gradient distribution with the increase of the hydraulic load was gradually from the left to the right side of the reactor. The pressure was mainly distributed in the remaining area of the reactor, which led to an earlier entry into the stable structure at the same time along with the narrower water channel and the larger head loss. The significant reduction in the flow capacity at local or even full regions could be confirmed from the velocity vectors in Figure 3. The results revealed that the flow was smaller and the degree of mixing was weaker. The weak mixing degree not only affected the seepage but also the pollutant reduction.

Figure 2. Static pressure in BRC under: (a) 0.5 m³/m²·d of hydraulic loading rate, (b) 1.0 m³/m²·d of hydraulic loading rate, (c) 1.5 m³/m²·d of hydraulic loading rate, and (d) 2.0 m³/m²·d of hydraulic loading rate

Figure 3. Velocity vectors in BRC under: (a) 0.5 m³/m²·d of hydraulic loading rate, (b) 1.0 m³/m²·d of hydraulic loading rate, (c) 1.5 m³/m²·d of hydraulic loading rate, and (d) 2.0 m³/m²·d of hydraulic loading rate

3.2. Comparative analysis of tracer test and simulation results

The hydraulic loading had a significant impact on the short-circuit flow and dead angle in the vertical flow filtration system. The relationship between the hydraulic load and the hydraulic efficiency had a very important role in understanding pollutant removal [28]. Padilla et al. (1999) and Nützmann et al. [29,30] assumed that there was a positive relation between the degree of saturation and the degree of diffusion, and SZ could improve the flow uniformity. However, both studie [29,30] demonstrated that there was an inverse proportional relationship between the degree of saturation and the degree of diffusion. For this purpose, seven sampling points (Figure 1) in SZ were investigated, and the RTD variety at different hydraulic loading rates and the sampling points at the same hydraulic loading rates were analyzed, respectively. According to Equation (6)(6) $N\left(t\right)=\frac{\left(E\left(t\right)-E_w\right)M_{NaCl}Q}{\left({\lambda }_{N_a}+{\lambda }_{Cl}\right)M}$(6) , effluent electrical conductivity was normalized under different hydraulic loading rates. The RTD distribution curves of seven sampling ports at different hydraulic loading rates were obtained, as shown in Figure 4.

Figure 4. Retention time density distribution curves at different hydraulic loading rates for: (a) sampling point No. 1, (b) sampling point No. 2, (c) sampling point No. 3, (d) sampling point No. 4, (e) sampling point No. 5, (f) sampling point No. 6, and (g) sampling point No. 7

Figure 4. Continued

The results for the variations in the peak appearance times, peak values, peak widths, and tailing length were compared directly as they appeared. According to Darcy Law [31] and Forchheimer Law [32] for the saturated porous media, K would be accelerated under higher HLR, where the HRT also decreased. Therefore, it was clear from the RTD curve that with higher HLR values, the peak widths were very narrow, and the higher the peak values appeared at the early period. The results presented in Figure 4(b-e), f, and g showed that the sampling data were all in accordance with the above rules. However, the peak values of Figure 4(a,b) appeared at the lower HLR, because the sampling point No. 1 and No. 4 at up and down position relation were located in the lower left position of the reactor, while the outlet of the reactor was in the lower right position. The flow was infiltrating along the water head, which led to the worst degree of mixing in the lower left position of the reactor. This trend was in complete agreement with the numerical simulation of flow velocity (Figure 3). The peak value of No. 3 (Figure 4(c)) appeared at the lower HLR but there was no obvious value even if the HLR increased to 2.0 m³/m²·d, which indicated that the partial dead angle existed in the No. 3 region (Figure 1). This might have been due to fewer voids by excessive compaction in the filling process of the filter materials. Low porosity led to a sharp reduction in the flow area and to divert the infiltration direction. This trend was consistent with the results of numerical simulation of pressure values (Figure 2).

It was directly compared with the RTD curve variation of each sampling point at the same HLR (Figure 5). The results revealed that the residence time in the curve showed the trend of significant increase and then, decrease as follows: T_1.0 > T_0.5 > T_1.5 > T_2.0. The hydraulic parameters at different HLR were calculated according to Figure 4 and Equation (7(7) $t_{\theta }=\frac{t_m}{t_n}$(7) –12(12) $\lambda =\frac{t_p}{t_n}=t_{\theta }\left(1-{\sigma }_{\theta }^2\right)$(12) ). These parameters are presented in Table 2. Under different HLR, both normalized retention time (t_θ) and hydraulic efficiency (λ) showed an increased and then decreased trend: t_θ-1.0 > t_θ-1.5 > t_θ-0.5 > t_θ-2.0 and λ_1.0 > λ_1.5 > λ_0.5 > λ_2.0. Data presented in Table 2 revealed that t_θ values were lower than 1 under each hydraulic loading condition, which was an indication of the small ‘dead zone’ in the system. This meant that through molecular diffusion or mechanical dispersion, the tracer entered into the dead zone and it is hard to enter the mainstream channel outflow. It is difficult for the tracer to enter the main infiltrating channel again after flowing into the ‘dead zone’ through molecular diffusion or mechanical dispersion.

Figure 5. Retention time density distribution curves for different hydraulic loading rates : (a) 0.5 m³/m²·d; (b) 1.0 m³/m²·d; (c) 1.5 m³/m²·d; and (d) 2.0 m³/m²·d

During the HLR of 1.0 m³/m²·d, the effective volume ratio for pollutant removal was the largest and the hydraulic efficiency was the highest. Meanwhile, a larger σ_θ² value illustrated that the flow was more close to completely stirred flow of CSTR model because of the larger distribution divergence. However, as shown in Figure 5 and Table 3, in addition to No. 3, the peak appearance times of the other sampling points: No. 1 ≈ No. 2 > No. 4 ≈ No. 5 ≈ No. 6 > No. 7, showed that the flow was close to the plug flow of the PFR model. According to the reactor theory, in the plug flow reactor, the reaction time for the reactant reducing to the same concentration was shorter; that is, the reaction efficiency was higher [33]. It is also revealed that the seepage flow is virtually impossible push flow of PFR model in SZ of DFBRC.

Table 3. Estimated hydraulic parameters based on impulse tracer test data

Hydraulic loading rates (m^3/m^2·d)	tm (h)	tn (h)	tp (h)	σ^2 (h^2)	σθ^2	N (h^−1)	tθ (h)	λ
0.5	8.1	33.65	6.9	168.716	0.149	6.750	0.241	0.205
1.0	6.2	16.82	5.0	55.168	0.195	5.167	0.369	0.297
1.5	3.1	11.28	2.8	67.637	0.098	10.333	0.275	0.248
2.0	1.7	8.41	1.6	4.173	0.059	17.000	0.202	0.190

3.3 Correlation analysis of flow distribution and hydraulic conductivity

Hydraulic conductivity tests were synchronized in two DFBRC with and without sewage influent under different HLR, which shows in Figure 6. The results revealed a non-linear relationship between K and HLR. The reason was the heterogeneous flow distribution increased by the increasing porosity and its non-homogeneous distribution in DFBRC, resulting in an increase of infiltrating flow and exceeding the applicable range of Darcy’s Law. From the mechanics aspect, the study of particle material mechanics found that the particle system had the unique characteristic of the nonlinear response to the external micro action [34].

Figure 6. Correlation between permeability, hydraulic efficiency and hydraulic loads

The results showed that if HLR increased from 0.4 m³/m²·d to 1.0 m³/m²·d, the K value had the similar downward trend with and without sewage influent when hydraulic efficiency (λ) reached the highest at 1.0 m³/m²·d. During this period, the microorganism is in the growth stage, and the complete biofilm is not formed on the particles surface, which makes the microorganism have little contribution to the change of flow field [35]. Therefore, the reason for this trend is that when the pressure increases to a certain value, the original structure between the particle swarm begins to deform and the local cavity collapses. Resulted in the formation of smaller voids with increased water head loss increased, and decrease in K value. When HLR increased to 1.5 m³/m²·d, there had been a fluctuating rises due to the dilatancy principle [36], the larger area of the particles was larger than the area of collapse, flow channel became wide, and K showed a upward trend.

When HLR increased to 2.5 m³/m²·d, there were different trends in the two curves, during infiltration of the deionized water, The K value decreases slightly and then tends to be stable, which is due to the complete collapse of particles. When the head pressure was insufficient to make the structure deformation, the flow path became more stable and the volumes of the individual bridges approached one another within a characteristic time of about 5 min [37], which ensured the smooth flow. During infiltrating the raw water with sewage influent, K decreased gradually, this trend could be explained through the fundamental principles of thermodynamics and biofilm morphology. A closed system tended to become free energy to minimal state, along with the entropy of the system reaching the maximum by thermal diffusion motion, and eventually, forming a uniform distribution equilibrium; therefore, granular structure was ultimately stabilized. At this time, the sewage carrying particulate pollutants and the gradual maturation biofilm had become a restrictive factor to infiltration, resulting in the occurrence of the clogging phenomenon. Formation of thick biofilms may cause pore clogging or the biofilm detachment under the shear force, thus providing a substantial decrease in permeability and a reduction in effluent flow, therefore isolating the pollutant into specific subsurface zones [38]. Thus, understanding biofilm development/detachment and its interaction with fluid flow in porous media becomes important for the successful development and application of many engineering techniques. However, in a porous medium with complex geometry, the local fluid pressure acting normal to the biofilm surface may also play an important role. Local clogging due to biofilm formation may cause an increase in the fluid pressure experienced by the biofilm surface at different locations (Duddu et al. 2009) [39]. The results revealed that when K was low (0.000107 m/s, K/K_ini = 0.79), the hydraulic efficiency was the highest. This trend ensured the contact and retention times; hence, DFBRC had the ability to remove pollutants efficiently.

The accurate control parameters for practical applications are determined as follows: the BRC should be strictly in accordance with the grading requirements of each layer; the HLR should be controlled at about 1.0 m³/m²·d, the K should be controlled at about 0.000107 m/s and K/K_ini = 0.79.

4. Conclusions

Permeability, flow distribution pattern, and hydraulic loading rates and hydraulic conductivity are very important factors in the performance of BRC‒. The flow distribution and the change of microscopic filter structure due to hydraulic loading have an important influence on hydraulic conductivity and thus on the decontamination ability. In this study, the use of a low hydraulic loading rate (approximately, 1.0 m³/m²·d) could improve the degree of diffusion, resulting in a lower hydraulic conductivity (approximately, 0.000107 m/s, K/K_ini = 0.79), which could achieve the desired permeability and treatment effects. On the basis of CFD simulation and experimental verification, the four-stage theory of ‘decrease-increase-then decline-stability’ of BRC seepage coefficient was put forward. This phenomenon reveals the relaxation process of pressure expansion when the packing particles are disturbed by HLR. The appropriate choice of initial hydraulic conductivity can be a key element in BRC design because it can improve the degree of diffusion and, therefore, indirectly increase the pollutant removal by limiting the volume of water bypassing the system through overflow. Hence, it is clear that the flow distribution and conductivity in BRC must be estimated more accurately on a microscopic scale .

Acknowledgements

This research was funded by the National Natural Science Foundation of China (No. 41967043 and No. 52160003) and Gansu Province Natural Science Foundation (20JR5RA461, 20JR10RA145).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was funded by the National Natural Science Foundation of China [No. 41967043 and No. 52160003] and Gansu Province Natural Science Foundation [Natural Science Foundation of Gansu Province 20JR5RA461, 20JR10RA145].