Correlation Equation for Predicting the Single-Collector Contact

Jun 9, 2015 - ABSTRACT: The single-collector contact efficiency (η0) for physicochemical colloid filtration under horizontal flow in saturated porous...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/Langmuir

Correlation Equation for Predicting the Single-Collector Contact Efficiency of Colloids in a Horizontal Flow Jing Li,*,† Xiaohu Xie,‡ and Subhasis Ghoshal† †

Department of Civil Engineering, McGill University, Montreal, Canada H3A 2K6 School of Computer Science, McGill University, Montreal, Canada H3A 0E9

‡

S Supporting Information *

ABSTRACT: The single-collector contact efficiency (η0) for physicochemical colloid filtration under horizontal flow in saturated porous media was calculated using trajectory analysis in three dimensions. Past studies have developed correlation equations for colloids with densities close to that of water, such as bacteria and latex particles. A new correlation equation was developed for predicting η0 based on a large number of trajectory simulations to account for higher-density particles representative of metal colloids. The correlation equation was developed by assuming Brownian diffusion, interception, and gravitational sedimentation contributed to η0 in an additive manner. Numerical simulations for colloid trajectory analysis used for calculating η0 were based on horizontal flow around a collector under the action of van der Waals attractive forces, gravity, and hydrodynamic forces as well as Brownian motion. The derived correlation equation shows excellent agreement with existing correlation equations for particles with density close to that of water. However, the correlation equation presented in this study shows that η0 of high-density colloids, such as metal particles, transported under horizontal flow deviates from that predicted by existing correlations for colloids larger than 4 μm and under low approach velocities. Simulations of trajectory paths show that a significantly reduced contact of high-density colloids larger than 4 μm in size with a collector is due to gravity forces causing trajectory paths to deviate away from the underside of collectors. The new correlation equation is suitable for predicting the single collector efficiency of large particles (several hundred nanometers to several micrometers) and with a large amount of density transport in the horizontal flow mode but is unsuitable for particles with a quite small size (several to tens of nanometers) and for the particle with a large amount of density flow in the vertical flow mode. The trajectory analysis was conducted for particles under favorable deposition conditions.

1. INTRODUCTION Advancing the understanding of the transport and deposition of engineered nanoparticles in saturated porous media is essential to the assessment of their fate and transport in subsurface media, their removal from sand filters, and their mobility from injection wells during remediation operations.1−5 The colloid filtration theory (CFT) is widely used in modeling the fate and transport of nano as well as submicrometer particles and biocolloids.6−11 In CFT, the probability of contact between colloids and a collector is quantified by the single-collector contact efficiency (η0), and the probability of successful attachment of the colloids in contact with the collector is quantified by the attachment efficiency or sticking efficiency (α). The value of η0 is calculated from correlation equations that account for particle, collector, and fluid properties.12−15 The value of α is usually obtained from the experimental transport data and η0. Significant attempts, however, have been devoted to developing correlations to predict α.16−20 With the CFT equation, the colloid deposition rate coefficient in saturated porous media is usually given by eq 1 (under the condition that the potential for © 2015 American Chemical Society

colloids to migrate along the solid−water interface before finding a retention location was not considered21) ⎡ 3(1 − f )U ⎤ k=⎢ α ⎥η0 ⎣ 2dsf ⎦

(1)

where k is the first-order rate constant, f is the porosity of media, ds is the diameter of granular porous media, α is the attachment efficiency, and U is the approach velocity. Certain modified forms of eq 1 have also been reported.22 The contributions of Brownian diffusion and various forces acting on colloids during its transport around collectors, such as hydraulic and gravity forces, van der Waals forces, and electrostatic double-layer interaction forces for the calculation of the single-collector contact efficiency (η0), have been outlined in several studies.12−15,22−32 The contribution of Brownian motion to η0 is dominant for submicrometer Received: March 20, 2015 Revised: May 27, 2015 Published: June 9, 2015 7210

DOI: 10.1021/acs.langmuir.5b01034 Langmuir 2015, 31, 7210−7219

Article

Langmuir

Figure 1. Happel’s sphere-in-cell porous media model in horizontal flow origination. (a) Components of gravity force on a particle near a 3-D collector. (b) Symmetrical trajectory (green) of colloids on top and bottom of the collector for light particles. (c) Digressed trajectory (red) of colloids around the collector under significant gravity forces caused by heavy particles.

particles.15,22 There are two common approaches used for the estimation of η0: (1) Lagrangian trajectory analysis of the colloid that predicts the path of colloid around model collectors such as the Happel sphere in cell,12,13,22,30,32 hemisphere in cell,28,31 spheres in series, unit bed porous media geometries, or a constricted tube23,24,26,32 and (2) the numerical solution of the convective−diffusion equation describing the concentration of colloids around similar collectors (Eulerian approach).14,15,27,33 A summary of studies that have presented different correlation equations for η0 are presented in Table S1. Several studies have provided correlation equations to facilitate the calculation of η0.12,13,15,28,30 Two popular correlation equations for η0 are the equations developed by Rajagopalan and Tien13 (referred henceforth as the RT equation) and by Tufenkji and Elimelech (referred henceforth as the TE equation).14 The RT equation was developed based on the trajectory analysis of colloids around a Happel’s spherein-cell collector model for downward flow, and force and torque balance equations for gravity, surface force, and hydraulic drag on the colloid were solved. The RT equation does not account for Brownian diffusion in their colloid trajectory simulation, but rather a Brownian diffusion term was directly superimposed with terms describing interception and sedimentation. A recent study done by Nelson and Ginn22 suggested that the simple superimposition of the Brownian diffusion term into the η0 correlation equation could lead to overestimation. Several subsequent studies have incorporated Brownian diffusion directly into equations for calculating η0 in colloid transport around a collector by representing the contributions of Brownian diffusion as a Brownian force or by directly incorporating Brownian motion paths.15,22,26,28,34 The TE equation was also developed on the basis of Happel’s sphere in cell, but the deposition rate was calculated on the basis of the numerical simulation of the convective−diffusive equation (Eulerian approach). The TE equation was developed by including the effects of the van der Waals force and hydrodynamic interactions in the Brownian motion of particles, which were not included in the development of the RT equation. The RT and TE equations have been derived for twodimensional analyses of colloid transport around collectors in a vertical flow. The equations need to be modified for horizontal flow conditions, particularly when the colloids have a high density. Such conditions are relevant to the transport of nanoscale zero-valent iron (NZVI) particles or other engineered metal particles in subsurface porous media. The direct injection of NZVI for the remediation of sites contaminated by chlorinated organic compounds, metalloids,

and heavy metals has been suggested.35−38 In such scenarios, NZVI particles with densities in the range of 6 to 7 g/cm3 travel horizontally following injection from wells. For vertical flow, the gravity force is aligned with the flow direction, and thus the trajectory of the particle can be calculated on the basis of two dimensions and gravity is decomposed in radial and tangent directions (as shown in Figure S1). For horizontal flow, the gravity force acts perpendicular to the flow direction and induces a downward trajectory of the particle, causing it to deviate from the flow direction and the horizontal plane at which the particle was approaching the collector. During vertical flow, the gravity force does not cause such a deviation from the vertical plane when colloids transport around collector, resulting in trajectory patterns which can be described in two dimensions. The trajectory of the colloid during horizontal flow around the collector can be accounted for by resolving the gravity force in three dimensions (as shown in Figure 1a). In addition to accounting for horizontal flow and the effects of density for metal nanoparticles, we also consider the effects of solution viscosity. The presence of polymers or polyelectrolytes, commonly used to coat engineered nanoparticles or added to achieve shear thinning of the aqueous medium to promote the colloidal stability and transport of particles, leads to higher solution viscosity because of the free polymer fraction present in solution.39−42 The changes in viscosity can directly influence the rate of deposition of colloids. Herein, we report modified equations based on the RT model in 3-D for predicting η0 for particle transport around a collector under horizontal flow conditions. Furthermore, the correlation equations were improved from the RT equation by explicitly incorporating Brownian diffusion in trajectory simulations in this study, as done in the study by Nelson and Ginn,22 but with autoadapted time steps achieved by using the Runge−Kutta−Fehlberg method rather than a fixed time step used in a correlation equation developed by Nelson and Ginn (referred henceforth as the NG equation).30 The model was used to estimate η0 within a wide range of parameter values for the particle density (0.95 to 7 g/cm3), viscosity of the suspension (1 × 10−3−26.3 × 10−3 kg/m/s), flow velocity (4 × 10−7−2 × 10−3 m/s), and collector diameter (0.05−1.2 mm) (as shown in Table 1).

2. NUMERICAL SIMULATION FOR η0 AND VERIFICATION EXPERIMENTS 2.1. Geometry of the Collector and Forces Acting on the Particle. η0 is calculated in the RT model as13 7211

DOI: 10.1021/acs.langmuir.5b01034 Langmuir 2015, 31, 7210−7219

Article

Langmuir

The stream function in the liquid shell was obtained from Rajagopalan and Tien’s study13 and is depicted in detail in section II of the SI. As a result, the velocity field of the flow in the liquid shell is

Table 1. Summary of Parameter Values Used in Numerical Calculations parameters

range

particle diameter, dp collector (grain) diameter, ds particle density, ρp dispersant viscosity υ, approach velocity, U Hamaker constant, A temperature, T porosity, f

0.01−10 μm 0.05−1.2 mm 0.95−7.00 g/cm3 1.0 × 10−3−26.3 × 10−3kg/m/s 4 × 10−7−2 × 10−3 m/s 1 × 10−20−4 × 10−20 J 298 K 0.26−0.48

η0 =

I πb UC0 2

⎛ 1 ∂ψ ⎞ ⎛ 1 ∂ψ ⎞ ⎟e + ⎜ ⎟e V = −⎜ ⎝ r sin θ r ∂θ ⎠ r ⎝ r sin θ ∂r ⎠ θ

where ψ is the stream function (additional details in eq S1). All equation terms are explained in the Nomenclature section at the end of the article. The forces and torques exerted on the colloid due to gravity, surface interactions, and hydraulic drag are described and listed in detail in Table S2. The primary difference between the forces working on particles in 2-D and 3-D collectors was the gravity forces, as expressed in the following equations:

(2)

where I is the rate at which a particle collides onto the collector, b is the summation of the collector radius and the thickness of the liquid shell, which is dependent on the radius and porosity of the collector (as described in the Equation S1), U is the approach velocity of fluid, and C0 is the number concentration of colloid particles in fluid. In the RT model, the trajectories of the suspended particles across the collector were aligned with the direction of the fluid flow and were symmetrical around the particle. Thus, the trajectory can be expressed in 2-D with polar coordinates (r, θ).13 On the basis of θs, the angle below which all trajectories can contact the collector, the value of η0 can be determined as

η0 = sin 2 θs

(4)

f 2G =

⎛4⎞ 3 ⎛4⎞ 3 ⎜ ⎟πa (ρ − ρ )g = ⎜ ⎟πa (ρ − ρ )g f f ⎝3⎠ p p ⎝3⎠ p p

[− cos θ × er + sin θ × eθ)

f 3G =

(5)

⎛4⎞ 3 ⎛4⎞ 3 ⎜ ⎟πa (ρ − ρ )g = ⎜ ⎟πa (ρ − ρ )g f f ⎝3⎠ p p ⎝3⎠ p p

[− sin θ × cos ϕ × er − cos θ × cos ϕ × eθ − sin φ × eϕ) (6)

From the balances of forces and torques, the velocities of particles in the liquid shell around the 3-D collector in the r, θ, and φ directions were obtained as ⎛ ∂u ⎞ 4 m⎜ ⎟ = − 6πvapur f rt − 6πvapAC y 2 f rm − πap3(ρp − ρf ) ⎝ ∂t ⎠r 3 3 ⎧ ⎫ 2Aap α ⎪ ⎪ ⎨ 2 ⎬ g cos φ sin θ − ⎪ 2⎪ 3 (2 a ) δ δ + ⎩ ⎭ p 2 2 ⎫ ⎧ ⎤⎡ ⎡ 2ξ ξ ⎪ εapκ (ξc + ξp ) ⎪ e−κδ ⎤ c p ⎬×⎢ 2 ⎥ +⎨ − e−κδ ⎥⎢ 2 ⎪ ⎪ ⎥⎦⎣ (1 − e−2κδ ⎦ 2 ⎭ ⎢⎣ (ξc + ξp ) ⎩ (7)

(3)

where θs is the terminal angle of the integration. In contrast, for the 3-D collector, in addition to r and θ coordinates used in the 2-D collector, a third polar coordinate, φ, is employed (as shown in Figure 1). Because randomized Brownian motion was included in the trajectory simulations which can result in deviations from the trajectory path, calculated on the basis of flow, gravity, surface forces, and hydrodynamic drag, we did not compute the limiting φ (similar to the aforementioned θs) to determine the initial area in which all trajectories would contact the collector. Instead, in our analyses, a sufficiently large number of particles are considered to be uniformly distributed in the flow cross-section (shown in Figure 1b,c). This is systematically performed by dividing the flow cross-section into 20 equal and concentric areas as well as calculating the trajectories starting from the central section moving toward the circumferential regions, and the number of particle contacts with the collector is simultaneously recorded. Once a concentric area with no trajectory contacting the collector is encountered, the calculations are skipped for the regions in the outer layers. The trajectories were calculated for 6000 particles, and it was verified in several simulations that there was minimal change in the number of contacts of particle and collector by increasing the particle number to 10 000. The total calculated contact events divided by the total particle number (6000) approaching the collector yielded η0. 2.2. Particle Velocity. Even though the particle trajectory is considered in three dimensions in this study, the fluid velocity was considered in two dimensions, as described by the RT model. This is because the stream function is not altered due to the constant volume flux through the collector surface.43 Thus, the stream function can still be depicted in the 2-D radial coordinate system (r, θ).

⎛ ∂u ⎞ m⎜ ⎟ = − 6πvapuθ f θt + 6πvap 2w1f θr + 6πvap ⎝ ∂t ⎠θ [BC yf1mθ + DC y 2 f 2mθ ] −

4 πvap3(ρp − ρf )g cos θ cos ϕ 3

⎛ ∂u ⎞ 4 m⎜ ⎟ = 6πvapuϕf ϕt + πvap3(ρp − ρf )g sin ϕ ⎝ ∂t ⎠ϕ 3

(8)

(9)

Equations 7−9, which present particle velocities in threedimensional coordinates on the forces in the corresponding coordinate, can be simplified as ⎞⎡ ∂ur U⎛ 1 = t⎜ + 1 + δ +⎟⎢− AC+(1 + δ +)2 f rm + NG sin θ cos ϕ ∂t f r ⎝ NR ⎠⎢⎣ −

NLoas(δ +)2 ⎤ ⎥ 4(as + ap) ⎥⎦

(10)

∂uθ U + = [BC s2 + DC+(1 + δ +)s3 + NG cos θ cos φ] ∂t rs1 (11)

∂uφ ∂t 7212

=

U NG sin φ rs1

(12) DOI: 10.1021/acs.langmuir.5b01034 Langmuir 2015, 31, 7210−7219

Article

Langmuir

Figure 2. Predictions of η0 from the TE, RT, and NG equations as well as the 3-D numerical simulated η0 with particle diameter under the following conditions: (a) ρp = 1.05 g/cm3, υ = 10−3 kg/m/s, U = 0.508 cm/min, ds = 340 μm; (b) dp = 2000 nm, ds = 340 μm, f = 0.37, U = 0.508 cm/min, υ = 10−3 kg/m/s; (c) dp = 2000 nm, U = 0.508 cm/min, f = 0.37, υ = 10−3 kg/m/s, ρp = 1.05 g/cm3; (d) dp = 200 nm, ds = 340 μm, f = 0.37, U = 0.508 cm/min, ρp = 1.05 g/cm3.

In the above equations, NR, NG, and NLo are dimensionless numbers that characterize the various forces acting on the particle. The wall-effect correction factors (s 1−s 3) are incorporated into the calculations. 2.3. Brownian Diffusion. Brownian motion is simulated by the method reported by Nelson and Ginn.22 Briefly, the method can be described as follows. pn is the three-dimensional displacement vector at time tn, and Δt is the time step. The displacement with time is calculated by22

∑ biki i=1

tn + 1 = tn + Δt

(13) (14)

(15)

Δp = [Δx , Δy , Δz]T

(16)

(19)

(20)

where ηD is the contribution of Brownian motion to the contact, ηI is the contribution of interception to the collision, and ηG is the contribution of gravitational sedimentation to the collision. Although previous studies13,14,44 have suggested that all three terms were power functions of multiple dimensionless parameters, Nelson and Ginn have suggested that the gravitational sedimentation term was not a simple power function of NG.30 The correlation of these three terms is determined by using the methodology developed by Tufenkji and Elimelech.14 Additional details are provided in section S4 in the SI. 2.5. Verification Experiments. A set of column transport experiments were conducted to verify the particle settling mechanism in horizontal flow revealed in the above mathematical simulations with a type of iron oxide (with a

where pn+1 is the ultimate displacement of the particle at time tn+1 and Δp is the three-dimensional displacement vector generated by Brownian motion, which can be expressed in Cartesian coordinates as Δx = nx 2D∞Δt

Δz = nz 2D∞Δt

η0 = ηD + ηI + ηG

where pn̅ +1 is the three-dimensional partial displacement vector, which is derived from the forces (London van der Waals force, gravity, and hydrodynamic forces) at time tn+1. Brownian motion is not used in the calculations in the above steps. After p̅n+1 is computed, the displacement rendered by Brownian motion is added to the above calculations: pn + 1 = pn̅ + 1 + Δp

(18)

where nx, ny, and nz are random variables of normal (0, 1). D∞ is the diffusion coefficient of the particle whose value is the same as the bulk diffusion coefficient. The time step Δt was adaptively controlled in the Runge−Kutta−Fehlberg method, by which the range of Δt is bounded to 10−7 to 10−13 s. 2.4. Regression Analysis to Obtain the Correlation Equation. The superposition of the three deposition mechanisms of diffusion, interception, and sedimentation to yield the correlation equation of η0 is assumed in this study, as in the cases of developing RT and TE equations

4

pn̅ + 1 = pn +

Δy = ny 2D∞Δt

(17) 7213

DOI: 10.1021/acs.langmuir.5b01034 Langmuir 2015, 31, 7210−7219

Article

Langmuir diameter of 3 μm and a density of 2.25 g/m3) in two different sizes of representative collectors (sand, diameters are 1140 and 340 μm).

shown in Figures S2 and S3. The reason for the deviation is discussed in the following section. The effect of collector diameter on η0 is shown in Figure 2c. The numerically simulated η0 on the collector decreased with the increase in collector diameter in the TE, RT, and NG equations for a particle with a low density of 1.05 g/cm3. However, by increasing the particle density to 6.58 g/cm3 with an identical size of 2000 nm, the change in η0 with collector diameter showed distinct trends from different prediction methods (Figure S4). In particular, 3-D simulations and the NG equation show an approximately constant η0 with collector diameter, whereas the TE and RT equations exhibit a significant increase in η0 with the increase in collector diameter. Figure 2d showed the effect of the viscosity of the suspension on η0. For a small particle diameter such as 200 nm, no noticeable difference in η0 between the 3-D simulation and the RT, TE, and NG equations was observed for the whole range of viscosity of 1 × 10−3−26.3 × 10−3 kg/m/s at a particle density of 1.05 g/cm3. The influence of flow velocity on η0 was investigated, and the results are shown in Figure S5. The trends of predicted η0 from TE, RT, and NG equations as well as from our 3-D simulation for horizontal flow around a collector showed satisfactory agreement for particles with a small density (1.05 g/cm3), which provides validation of the numerical simulations presented in this study. The comparison of the trend of η0 with the porosity change between different equations and the 3-D simulations is shown in Figure S6. Briefly, the η0 values remain constant in TE, NG, and our 3-D simulations within the whole porosity range of 0.28−0.48 because the gravitational settling term, which is dominant in the η0 value for a 2000 nm particle diameter and a 6.58 g/cm3 of particle density, is not porosity-dependent. The η0 values determined by 3-D simulations are very close to TE equation predictions and are higher than those obtained from the NG equation. In the NG equation, the authors considered a lower particle concentration at the collector boundary than in the bulk aqueous phase, and this leads to the differences in the predicted η0 values.22 Only the RT equation showed a trend of decreasing η0 with the increase in porosity. This is due to the change in values of the porosity-dependent parameter, As (eq S3), in the gravity settling term as listed in Table S1. 3.2. New Correlation Equation for Colloid Transport. From the regression analysis, the overall correlation equation for predicting the single collector contact efficiency in the 3-D model can be expressed as

3. RESULTS AND DISCUSSION 3.1. Comparison of the Numerical Results with the RT, TE, and NG Equations. The validation of our model was conducted by comparing several simulations for various parameter values with those predicted by RT, TE, and NG equations. The results from the RT, TE, and NG equations here were obtained by using the reported correlation equations derived on the basis of conditions of vertical flow in 2-D for a Happel sphere collector.13,30,45 The trajectory simulation of particles in the 2-D Happel sphere collector was also conducted by using our method (without the decomposition of gravity in three dimensions). The simulated results (Supporting Information II) showed that for latex particles the 2-D simulation results were significantly consistent with those predicted by the TE equation, revealing that the method developed here was appropriate. The results shown in Figure 2 present a comparison of results with a series of parameters: particle diameter, particle density, collector diameter, and viscosity of the suspension. In Figure 2a, no significant difference in η0 values existed between the TE prediction and our 3-D simulation, suggesting that our 3-D model methodologies are reliable. However, there was a significant difference displayed between our 3-D model simulation or the TE equation and the RT and NG results. The 3-D model and the TE equation predict increasingly higher η0 values than those predicted by the RT equation for particle diameter larger than 2000 nm at a particle density of 1.05 g/ cm3 (density similar to that of latex particles or bacteria). This can be explained by the difference in the ηG expressions in the RT and TE equations: The coefficient value in the ηG expression in the former was 2.4 × 10−3, but it was 0.22 in the latter, indicating a lower estimation of ηG in the RT equation compared to that in the TE equation. η0 values obtained from the NG equation were even lower than those obtained from TE equation over the whole range of particle size, as shown in Figure 2a, primarily because of the difference in ηG expressions and the factor of γ2 for all three terms in the NG equation, which is attributable to the consideration of the boundary condition for developing the NG equation as well as the smaller diffusion contribution caused by the combined effect of diffusion and gravity for large particles.30 The influence of particle density on η0 as predicted by RT, TE, and NG equations as well as 3-D simulation is shown in Figure 2b. Overall there was good agreement between the 3-D numerical simulations with the RT model at the particle diameter of 2000 nm within the fixed density range (0.95−7.00 g/cm3). The NG equation, however, predicts lower η0 values compared to predictions from other equations or methods, as the boundary condition was considered for developing the NG equation, leading to a low particle concentration around the collector surface and thus a smaller collector efficiency. Additionally, given that the NG equation fits a narrow particular particle density range (from 0.95 to 1.8 g/cm3), for particles whose density is larger than 1.8 g/cm3 the equation might not be suitable for predicting η0 and is shown as a much lower η0 prediction. With changing particle size, the numerically calculated η0 showed a different dependence on particle density from those predicted by the TE, RT, and NG equations, as

η0 = 14.16As1/3NR −0.04434NPE−0.8538 + 0.621As NR1.6548NLO 0.105 +

⎡ 0.3507NR −0.07⎢⎣

15.38

⎤

⎦ 15.38 + 0.78NG1.2 ⎥

[1 + exp(2.6683 − 6.796NG)]

(21)

and thus ηD = 14.16As1/3NR −0.04434NPE−0.8538

(22)

ηI = 0.621As NR1.6548NLO 0.105

(23)

ηG =

⎡ 0.3507NR −0.07⎢⎣

⎤ 15.38 ⎦ 15.38 + 0.78NG1.2 ⎥

[1 + exp(2.6683 − 6.796NG)]

(24)

From eq 21, it can be observed that major differences between our 3-D correlation and the RT, NG, and TE 7214

DOI: 10.1021/acs.langmuir.5b01034 Langmuir 2015, 31, 7210−7219

Article

Langmuir

equations can be attributed to the different expressions of the gravitational term: ηG appears to increase with the increase in NG in the latter three equations but not in eq 21. The 3-D model simulated a very different particle contact pattern for particles with a density of 6.58 g/cm3 compared to those with a density of 1.05 g/cm3. As shown in Figure S7, the particles with a density of 1.05 g/cm3 and a diameter of 4000 nm are deposited uniformly on the top and botton of the collector, whereas when the density increased to 6.58 g/cm3, particles are deposited mainly on the top part of the collector. Particles did not make contact with the collector from the bottom side, as a result of the significant settling of particles out of the liquid shell trajectory space around the collector, caused by their considerable gravity (as shown in the schematic in Figure 1c), verifying the result shown in Figure 3 that the η0 values obtained from 3-D simulation are lower than those predicted by the TE and RT equations. Furthermore, the result that η0 values remain nearly constant in the 4000 to 10 000 nm particle diameter range, as shown in Figure 3, can be ascribed to the likelihood that with the increase in particle size in this range the population of particles which can be transported onto the top side of the collector would decrease because of gravity settling. On the other hand, with the increase in particle size, the contact probability of particles on the top side of the collector would be increased due to the higher gravity forces. It is likely that these two effects balance out to yield a constant η0. The gravity settling term in eq 21 is not directly power-lawdependent on NG, unlike in the RT and TE equations. The comparison of η0 values with the particle density predicted by eq 21, numerical calculations, and the TE, RT, and NG equations is shown in Figure S8. For larger particles such as 4000 nm, the increasing trend of η0 with particle density was slightly diminished in our 3-D simulations and eq 21, compared to the TE, RT, and NG equations, when the particle denstiy exceeded 4 g/cm3. This result supports our hypothesis that when the particle density increases to a certain value, gravity forces would result in a very significant gravitational settling on particles in the horizontal flow mode, leading to lower contact in the bottom hemisphere of the collector and yielding smaller η0 than those from TE and RT equations, under identical conditions. However, the NG equation gives a smaller η0 than the prediction from eq 21 over the whole range of particle density The effect of the approach velocity on η0 is presented in Figure 4. For a particle diameter of 4000 nm, η0 simulated from our 3-D model showed a distinct difference from those predicted in TE, RT, and NG equations at low approach velocities: the 3-D model yielded a lower η0 than those predicted in the TE and RT equations. This is attributable to the relatively significant gravity settling of particles under lower hydraulic drag forces in our 3-D model in this low flow velocity range, yielding less possibility for particle contact. Nelson and Ginn have already pointed out that η0 values predicted from the TE and RT equations commonly exceed unity in the very low approach velocity range.30 However, our 3-D model did not yield η0 values in excess of unity for identical parameter values. It should be noted that in the very low velocity range, η0 values increased then decreased with the increase in the approach velocity from our 3-D model (shown in Figure 4). This is because particles could not be transported onto the top side of the collector at very low velocity, but with the increase in velocity, the proportion of particles transported on the top side increased, leading to the increase in the contact efficiency. The

equations are the following: (1) The gravitational sedimentation term, ηG (eq 24), derived on the basis of the 3-D collector in horizontal flow is not a simple power function of the dimensionless parameters as in the ηG expressed in TE and RT equations. However, the gravitational term in eq 24 is similar to the expression of ηG in the NG equation.30 (2) The diffusion term in eq 22 differs from the ηD expressions in the RT and TE equations: ηD ≈ N−0.04434 in this study is different from ηD ≈ R N−0.081 in the TE equation, and NR is not contained in the ηD R expression in the RT equation. This difference indicates that the dependence of η0 on the ratio of particle radius to collector radius is more significant in the TE equation and is likely due to the hydrodynamic retardation effects considered in developing the TE equation. (3) The expression of ηI in eq 23 also differs from those in the TE equation (ηI = 0.55 AsN1.675 N0.125 R A ) and the 15/8 1/8 RT equation (ηI = 0.72 AsNR NLo ). Because NLo = (3/4) NA, ηI expressions in the TE and RT equations are very close to each other. Overall, the key feature of eq 21 is the ηG expression, which can be used to more precisely predict the single-collector contact efficiency of particles with high gravity forces under horizontal flow. The comparison of the change in η0 with the particle diameter for higher particle density (ρp = 6.58 g/cm3), obtained from numerical simulation, eq 21, and TE, RT, and NG equations, is shown in Figure 3. η0 values predicted from eq 21

Figure 3. Comparison of predictions of single-collector contact efficiency (η0) based on rigorous numerical simulations of trajectories of particles, the new correlation, and the predictions from TE, RT, and NG equations under the following conditions: ρp = 6.58 g/cm3, υ = 10−3 kg/m/s, U = 0.508 cm/min, and ds = 340 μm.

match those numerical simulation data fairly well. However, unlike the results shown in Figure 2a, according to eq 21, the dependence of η0 on particle diameter at high particle density (ρp = 6.58 g/cm3) significantly diverged from those obtained from TE, RT, and NG equations. The predicted η0 values from the 3-D simulation match closely with the RT and TE equations up to a particle diameter of 3000 nm, but the η0 values remain constant at around 0.3 with increasing particle diameter in the 3-D simulation and eq 21. In contrast, the η0 values increase with the increase in particle diameter for the TE, NG, and RT equations, even though the predicted η0 values from the NG equation are much lower than those from the other two equations over the whole range of particle density. The discrepancy between eq 21 and the NG, TE, and RT 7215

DOI: 10.1021/acs.langmuir.5b01034 Langmuir 2015, 31, 7210−7219

Article

Langmuir

gravitational settling term in the equation as well as the smaller ηD caused by the combined effect of Brownian diffusion and gravity, showing η0 values that are always lower than the predictions from our 3-D model as well as the TE and RT equations. In addition, the NG equation limits the particle density range from 0.95 to 1.8 g/cm3, and thus it might be not suitable for predicting η0 for particles with larger densities and in turn presenting lower η0. The essential reason for the smaller η0 predictions in the NG equation is that the authors believed the existence of the particle concentration gradient under the boundary condition, yielding a lower particle concentration in the shell and thus leading to a lower contact possibility.30 Furthermore, as considered in the combined effect of Brownian diffusion and gravity, the diffusion term in the NG equation gives a lower ηD prediction, particularly for large particles. Additionally, the sedimentation term of the NG equation is a positively nonpower functional of the gravity group, showing a smaller increasing trend in η0 with ηG.30 The interception and gravity sedimentation predominantly determine the η0 value for particles with a 4000 nm diameter and a 6.58 g/cm3 density. With increases in the collector size, these relatively heavy particles would tend to deposit more in the top of the collector during vertical flow (as expected in the RT and TE equation) due to the effect of gravity, thus resulting in the increased η0. However, the sedimentation term in the NG equation is insignificantly dependent on NR (a function of N−0.05 ), and thus R showed an approximately constant η0 value with the increase in collector diameter, but for horizontal flow around a 3-D collector, with an increase in collector size, a higher density of colloids would be subject to increased difficulty to move to the top part of the collector. On the other hand, if particles do move to the top part of the collector, then the higher gravity forces would be favorable for particles to deposit on the top side of the collector. These two factors lead to a nearly constant η0. However, at the identical particle diameter and smaller particle density of 1.05 g/cm3, η0 values decreased with the increase in the collector diameter in our 3-D simulation, simliar to those from the RT, TE, and NG equations (shown in Figure 2c). This is because the primary contribution for particles to collide with the collector was not gravity but surface forces for light particles. With the increase in the collector size, the surface forces decreased because they are inversely correlated with the distance between the centers of the particle and the collector, resulting in a decrease in η0. 3.3. Verification of Gravity Settling in Horizontal Flow in Column Transport Experiments. When the transport experiments with iron oxide particles were complete, significant particle settling on the bottom side of the column was observed, particularly in the case of coarse sand (with a diameter of 1140 μm), as shown in Figure S9a. A newly published paper also observed the significant gravity settling in a horizontal column for clay particles.46 But in the case of fine sand (with a diameter of 340 μm), the settling of particles onto the bottom side of the column was not obvious, as shown in Figure S9b. This verifies the above mathematical simulations that heavy particles would be subject to more difficulty in depositing on the bottom side of the large collector and would exhibit an accumulative effect to settle on the bottom side of the horizontally placed column. The 3-D mathematical simulation of the particle trajectory in horizontal flow mode for calculating η0 and the subsequent established correlation equation in this study were conducted for particles under favorable deposition conditions. The 3-D

Figure 4. Comparison of predictions of single-collector contact efficiency (η0) based on rigorous numerical simulations of trajectories of particles, the new correlation, and the predictions from TE, RT, and NG equations under the following conditions: dp = 4000 nm, υ = 10−3 kg/m/s, ρp = 6.58 g/cm3, and ds = 340 μm.

NG equation extends the CFT applicability to a lower velocity range up to 1 × 10−7 m/s and gives lower η0 values because of the different boundary conditions.30 The comparison of η0 values as a function of collector diameter predicted by eq 21, numerical calculations, and the TE, RT, and NG equations is shown in Figure 5. At a particle

Figure 5. Comparison of predictions of single-collector contact efficiency based on rigorous numerical simulations of trajectories of particles, the new correlation, and the predictions from the TE, RT, and NG equations under the following conditions: dp = 4000 nm, υ = 10−3 kg/m/s, U = 0.508 cm/min, ρp = 6.58 g/cm3.

diameter of 4000 nm and a particle density of 6.58 g/cm3, η0 values remained constant from our numerical simulations and the NG equation prediction, unlike the trends of increasing η0 with the increase in the single-collector contact efficiency in the predictions from the TE and RT equations. Although η0 also approximately remained constant from the NG equation, it was smaller than that from our 3-D model predictions. As discussed in Figure 2a,b, the lower prediction of η0 in the NG equation is caused by the factor of γ2 for all three terms and the expected smaller gravitational settling contribution resulting from the 7216

DOI: 10.1021/acs.langmuir.5b01034 Langmuir 2015, 31, 7210−7219

Langmuir

■

trajectory simulation (with the direct combination of Brownian motion into trajectory path simulation) was confirmed to be appropriate from the comparison of simulations of light particles (with the density close to that of water) in 2-D and 3-D. For heavy particles, the 3-D simulations revealed that their trajectory deviated from those of light particles and the heavy particles performed differently in the top and bottom hemispheres of the collector. η0 calculated from our 3-D model can actually reveal what occurred when the particle (particularly for the engineered metal particles) approached the collector in natural flow mode, but it is obvious that gravity settling causes the particle distribution in the direction perpendicular to the flow direction to be nonuniform, unlike that in the vertical flow mode. Thus, the usually used conventional convection−dispersion equation for particle deposition calculations might be not suitable for large engineered metal particles in natural flow mode. The column transport experimental results also verified this: a smaller singlecollector efficiency for the particle is expected in the horizontal column but less transport is obtained, due to gravity settling in the macroporous media.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: (514) 398-6842. Fax: (514) 398-7361. Present Address

(J.L.) Qihe Green-Spring Environmental Technology Company, Shandong Province, China. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The research was funded by the Natural Sciences and Engineering Research Council of Canada, by Golder Associates Ltd., and by the Fonds de Recherche du Québec - Nature et Technologies. J.L. was supported in part by a McGill Engineering doctoral award.

■

NOMENCLATURE

Symbols

A ap as b C0

Hamaker constant radius of a particle radius of the spherical collector radius of the Happel cell the number concentration of particles in the suspension D∞ diffusion coefficient in an infinite medium, D∞ = kBT/(6πυap) ds diameter of a spherical collector, dc = 2as dp diameter of a particle, dp = 2ap er, eθ, eφ unit vectors in the r, θ and φ directions f tr, ftθ, f mr , fm1θ, f m2θ, frθ drag correction factors in Table S2, functions of δ f ji force vector on the particle, where the superscript specifies the source of the force g magnitude of the gravity vector, g = 9.81 g/ s2 r t m m gφ, gφ, g1φ, g2φ torque correction factors in Table S2, functions of δ k first-order rate constant kB Boltzmann constant, 1.3805 × 10−23 J/K NG gravity group, NG = 2ap2(ρp − ρf)g/9υU NLo London group, NLo = A/9πυap2U NPE Peclet number, NPE = Uds/D∞ NR relative size group, NR = ap/as p defined in eq S2, p = as/b r radial coordinate r* dimensionless radial coordinate, r* = r/as s1 s1 = (ftθgrφ − frθgtφ)/grφ s2 s2 = (frθgm1φ − y+fm1θgrφ)/grφ s3 s3 = (frθgm2φ − y+f m2θgrφ)/grφ t time variable T absolute temperature U approach velocity V liquid flow field y y = (r − as)

4. CONCLUSIONS The correlation equation and the framework for trajectory analyses suggested here will improve the calculations of η0 and colloid deposition rate coefficients during the horizontal flow of colloids. The correlation equation is particularly effective in calculating η0 for a higher density (e.g., metal colloids) of 4 μm and higher than that in the TE, RT, or NG equation. The flow of metal colloids in the subsurface may often be in the horizontal direction, such as when magnetic (aggregated) nanoscale iron particles are injected through wells during remediation or other aggregated metal nanoparticles are released into aquifers. Given the increasing trends in the production of nanomaterials and the use of metal nanoparticles, a greater understanding is needed for the transport behavior of such colloids in groundwater environments where horizontal flow around collector soil grains is very common. This study shows that when accounting for the effects of gravity forces for high density, larger colloids becomes important for calculating the rate and location of contact of these colloids with the collector surface, and these effects are not adequately captured by commonly used single-collector contact efficiency correlation equations, which can cause erroneous predictions of η0 during horizontal flow as demonstrated in this study.

■

Article

ASSOCIATED CONTENT

S Supporting Information *

Summary of literature reporting correlation equations for single-collector contact efficiency, equations for the stream function for fluid flow around a collector, schematic showing gravity forces on a particle on a 2-D collector, figures showing additional simulation results for the estimation of the contact efficiency of particle density, viscosity, approach velocity, and collector diameter, and a figure showing particle trajectories contacting the collector for high-density and low-density colloids. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/ acs.langmuir.5b01034.

Greek Letters

a(δ, ap, λe) retardation correction, a function of δ, ap, and λe δ surface-to-surface separation between the collector and the particle δ+ dimensionless parameter, δ+ = δ/ap f porosity of the porous media 7217

DOI: 10.1021/acs.langmuir.5b01034 Langmuir 2015, 31, 7210−7219

Article

Langmuir η0 ηD, ηI, ηG θ θs φ κ λe υ ρf ρp ψ ur, uθ, uϕ

■

deposition during injection of a polyelectrolyte-modified Fe0 nanoparticle at high particle concentration in saturated sand. J. Contam. Hydrol. 2010, 118, 152−164. (17) Phenrat, T.; Song, J. E.; Cisneros, C. M.; Schoenfelder, D. P.; Tilton, R. D.; Lowry, G. V. Estimating Attachment of Nano- and Submicrometer-particles Coated with Organic Macromolecules in Porous Media: Development of an Empirical Model. Environ. Sci. Technol. 2010, 44, 4531−4538. (18) Bai, R.; Tien, C. Particle Deposition under Unfavorable Surface Interactions. J. Colloid Interface Sci. 1999, 218, 488−499. (19) Bai, R.; Tien, C. A New Correlation for the Initial Filter Coefficient under Unfavorable Surface Interactions. J. Colloid Interface Sci. 1996, 179, 631−634. (20) Elimelech, M. Predicting collision efficiencies of colloidal particles in porous media. Water Res. 1992, 26, 1−8. (21) Li, X.; Lin, C.-L.; Miller, J. D.; Johnson, W. P. Role of Grain-toGrain Contacts on Profiles of Retained Colloids in Porous Media in the Presence of an Energy Barrier to Deposition. Environ. Sci. Technol. 2006, 40, 3769−3774. (22) Nelson, K. E.; Ginn, T. R. Colloid Filtration Theory and the Happel Sphere-in-Cell Model Revisited with Direct Numerical Simulation of Colloids. Langmuir 2005, 21, 2173−2184. (23) Paraskeva, C.; Burganos, V.; Payatakes, A. Three-dimensional trajectory analysis of particle deposition in constricted tubes. Chemical Eng. Commun. 1991, 108, 23−48. (24) Cushing, R. S.; Lawler, D. F. Depth filtration: Fundamental investigation through three-dimensional trajectory analysis. Environ. Sci. Technol. 1998, 32, 3793−3801. (25) Tobiason, J. E.; O’melia, C. R. Physicochemical aspects of particle removal in depth filtration. J. - Am. Water Works Assoc. 1988, 54−64. (26) Johnson, W.; Li, X.; Yal, G. Colloid retention in porous media: Mechanistic confirmation of wedging and retention in zones of flow stagnation. Environ. Sci. Technol. 2007, 41, 1279−1287. (27) Ryan, J. N.; Elimelech, M. Colloid mobilization and transport in groundwater. Colloids Surf., A 1996, 107, 1−56. (28) Ma, H.; Pedel, J.; Fife, P.; Johnson, W. P. Hemispheres-in-cell geometry to predict colloid deposition in porous media. Environ. Sci. Technol. 2009, 43, 8573−8579. (29) Petosa, A. R.; Jaisi, D. P.; Quevedo, I. R.; Elimelech, M.; Tufenkji, N. Aggregation and deposition of engineered nanomaterials in aquatic environments: Role of physicochemical interactions. Environ. Sci. Technol. 2010, 44, 6532−6549. (30) Nelson, K. E.; Ginn, T. R. New collector efficiency equation for colloid filtration in both natural and engineered flow conditions. Water Resour. Res. 2011, 47, W05543. (31) Ma, H.; Pazmino, E. F.; Johnson, W. P. Gravitational Settling Effects on Unit Cell Predictions of Colloidal Retention in Porous Media in the Absence of Energy Barriers. Environ. Sci. Technol. 2011, 45, 8306−8312. (32) Wei, Y.-T.; Wu, S.-c. Development of a trajectory model for predicting attachment of submicrometer particles in porous media: Stabilized NZVI as a case study. Environ. Sci. Technol. 2010, 44, 8996− 9002. (33) Li, Z.; Zhang, D.; Li, X. Tracking colloid transport in porous media using discrete flow fields and sensitivity of simulated colloid deposition to space discretization. Environ. Sci. Technol. 2010, 44, 1274−1280. (34) Li, Z.; Zhang, D.; Li, X. Tracking colloid transport in real pore structures: Comparisons with correlation equations and experimental observations. Water Resour. Res. 2012, 48, W05533. (35) Kanel, S. R.; Greneche, J.-M.; Choi, H. Arsenic (V) removal from groundwater using nano scale zero-valent iron as a colloidal reactive barrier material. Environ. Sci. Technol. 2006, 40, 2045−2050. (36) Bennett, P.; He, F.; Zhao, D.; Aiken, B.; Feldman, L. In situ testing of metallic iron nanoparticle mobility and reactivity in a shallow granular aquifer. J. Contam. Hydrol. 2010, 116, 35−46. (37) Zhang, W. x.; Elliott, D. W. Applications of iron nanoparticles for groundwater remediation. Remediation 2006, 16, 7−21.

the single collector contact efficiency collection efficiencies due to diffusion, sedimentation, and interception, respectively angular coordinate the terminal angle of the integration angular coordinate Debye−Huckel reciprocal length wavelength of electron oscillation, λe = 100 nm viscosity of suspension density of the liquid density of the particle stream function velocity vectors in different directions

REFERENCES

(1) Tufenkji, N.; Ryan, J. N.; Elimelech, M. Peer Reviewed: The Promise of Bank Filtration. Environ. Sci. Technol. 2002, 36, 422A− 428A. (2) Quevedo, I. R.; Olsson, A. L. J.; Tufenkji, N. Deposition Kinetics of Quantum Dots and Polystyrene Latex Nanoparticles onto Alumina: Role of Water Chemistry and Particle Coating. Environ. Sci. Technol. 2013, 47, 2212−2220. (3) Heidmann, I. In Metal Oxide Nanoparticle Transport in Porous Media−An Analysis about (Un)Certainties in Environmental Research; Journal of Physics: Conference Series IOP Publishing: 2013; p 012042. (4) Raychoudhury, T.; Tufenkji, N.; Ghoshal, S. Aggregation and deposition kinetics of carboxymethyl cellulose-modified zero-valent iron nanoparticles in porous media. Water Res. 2012, 46, 1735−1744. (5) Ouyang, Y.; Shinde, D.; Mansell, R.; Harris, W. Colloid-enhanced transport of chemicals in subsurface environments: A review. Crit. Rev. Environ. Sci. Technol. 1996, 26, 189−204. (6) Tufenkji, N. Application of a dual deposition mode model to evaluate transport of Escherichia coli D21 in porous media. Water Resour. Res. 2006, 42, W12S11. (7) He, F.; Zhang, M.; Qian, T.; Zhao, D. Transport of carboxymethyl cellulose stabilized iron nanoparticles in porous media: Column experiments and modeling. J. Colloid Interface Sci. 2009, 334, 96−102. (8) Saleh, N.; Kim, H.-J.; Phenrat, T.; Matyjaszewski, K.; Tilton, R. D.; Lowry, G. V. Ionic Strength and Composition Affect the Mobility of Surface-Modified Fe0 Nanoparticles in Water-Saturated Sand Columns. Environ. Sci. Technol. 2008, 42, 3349−3355. (9) Ryan, J. N.; Elimelech, M.; Baeseman, J. L.; Magelky, R. D. Silicacoated titania and zirconia colloids for subsurface transport field experiments. Environ. Sci. Technol. 2000, 34, 2000−2005. (10) Bradford, S. A.; Yates, S. R.; Bettahar, M.; Simunek, J. Physical factors affecting the transport and fate of colloids in saturated porous media. Water Resour. Res. 2002, 38, 63−1−63−12. (11) Bolster, C. H.; Walker, S. L.; Cook, K. L. Comparison of and Transport in Saturated Porous Media. J. Environ. Quality 2006, 35, 1018−1025. (12) Yao, K.-M.; Habibian, M. T.; O’Melia, C. R. Water and waste water filtration. Concepts and applications. Environ. Sci. Technol. 1971, 5, 1105−1112. (13) Rajagopalan, R.; Tien, C. Trajectory analysis of deep-bed filtration with the sphere-in-cell porous media model. AIChE J. 1976, 22, 523−533. (14) Tufenkji, N.; Elimelech, M. Correlation equation for predicting single-collector efficiency in physicochemical filtration in saturated porous media. Environ. Sci. Technol. 2004, 38, 529−536. (15) Long, W.; Hilpert, M. A Correlation for the Collector Efficiency of Brownian Particles in Clean-Bed Filtration in Sphere Packings by a Lattice-Boltzmann Method. Environ. Sci. Technol. 2009, 43, 4419− 4424. (16) Phenrat, T.; Kim, H.-J.; Fagerlund, F.; Illangasekare, T.; Lowry, G. V. Empirical correlations to estimate agglomerate size and 7218

DOI: 10.1021/acs.langmuir.5b01034 Langmuir 2015, 31, 7210−7219

Article

Langmuir (38) Tratnyek, P. G.; Johnson, R. L. Nanotechnologies for environmental cleanup. Nano Today 2006, 1, 44−48. (39) Vecchia, E. D.; Luna, M.; Sethi, R. Transport in porous media of highly concentrated iron micro-and nanoparticles in the presence of xanthan gum. Environ. Sci. Technol. 2009, 43, 8942−8947. (40) Johnson, R. L.; Nurmi, J. T.; O’Brien Johnson, G. S.; Fan, D.; O’Brien Johnson, R. L.; Shi, Z.; Salter-Blanc, A. J.; Tratnyek, P. G.; Lowry, G. V. Field-Scale Transport and Transformation of Carboxymethylcellulose-Stabilized Nano Zero-Valent Iron. Environ. Sci. Technol. 2013, 47, 1573−1580. (41) Truex, M. J.; Vermeul, V. R.; Mendoza, D. P.; Fritz, B. G.; Mackley, R. D.; Oostrom, M.; Wietsma, T. W.; Macbeth, T. Injection of Zero-Valent Iron into an Unconfined Aquifer Using Shear-Thinning Fluids. Groundwater Monit. Rem. 2011, 31, 50−58. (42) Liu, Y.; Phenrat, T.; Lowry, G. V. Effect of TCE concentration and dissolved groundwater solutes on NZVI-promoted TCE dechlorination and H2 evolution. Environ. Sci. Technol. 2007, 41, 7881−7887. (43) Batchelor, G. K. An Introduction to Fluid Dynamics; Cambridge University Press: 2000. (44) Prieve, D. C.; Ruckenstein, E. Effect of London forces upon the rate of deposition of Brownian particles. AIChE J. 1974, 20, 1178− 1187. (45) Tufenkji, N.; Elimelech, M. Correlation Equation for Predicting Single-Collector Efficiency in Physicochemical Filtration in Saturated Porous Media. Environ. Sci. Technol. 2003, 38, 529−536. (46) Chrysikopoulos, C. V.; Syngouna, V. I. Effect of Gravity on Colloid Transport through Water-Saturated Columns Packed with Glass Beads: Modeling and Experiments. Environ. Sci. Technol. 2014, 48, 6805−6813.

7219

DOI: 10.1021/acs.langmuir.5b01034 Langmuir 2015, 31, 7210−7219