The impact of nanoparticle aggregation on their size exclusion during transport in porous media: One- and three-dimensional modelling investigations

Greater particle mobility in subsurface environments due to larger size, known as size exclusion, has been responsible for colloid-facilitated transport of groundwater contaminants. Although size exclusion is not expected for primary engineered nanoparticles (NP), they can grow in size due to aggregation, thereby undergoing size exclusion. To investigate this hypothesis, an accurate population balance modelling approach and other colloid transport theories, have been incorporated into a three-dimensional transport model, MT3D-USGS. Results show that incorporating aggregation into the transport model improves the predictivity of current theoretical and empirical approaches to NP deposition in porous media. Considering an artificial size-variable acceleration factor in the model, NP breakthrough curves display an earlier arrival when aggregation is included than without. Disregarding the acceleration factor, aggregation enhances NP mobility at regions close to the injection point at a field scale and causes their retention at greater distances through alteration of their diffusivities, secondary interaction-energy minima, and settling behaviour. This results in a change of residual concentration profiles from exponential for non-aggregating dispersions to non-monotonic for aggregating dispersions. Overall, aggregation, hitherto believed to hinder the migration of NP in subsurface porous media, may under certain physicochemical conditions enhance their mobilities and deliver them to further distances.


Other equations of the model.
Detachment can be considered using an alternative approach in which the attached particle concentration can also be tracked but the Kd parameter used to consider acceleration factor in Eqs. (3) and (4) is not included anymore: 1,2 where Sk is the mass concentration of deposited phase particles of size class k [MM -3 ] and is detachment rate coefficient for size class k [T −1 ]. These equations can be solved using the MT3D code when a non-equilibrium sorption model is considered. In doing so, model parameters for each size class can be converted as: = βk/ and = βk ⁄(ρb ). In MT3D notation 3 , βk is the first-order mass transfer rate between the mobile and retained phases for species k. It should be noted that detachment can be excluded from this model formulation in the MT3D code by assigning a large value to . 4,5 In scenarios including size exclusion, a variable Kd was assumed to vary with size according to the following linear equation: where , , and are Kd values for a given size class k, for size class with smallest particles, and for size class with largest aggregates, respectively; ak is the radius of size class k and amax is the radius of largest aggregate size class.
The collision kernel for environmental colloids is commonly given as the sum of three mechanisms: perikinetic collisions (Brownian), orthokinetic collisions (shear-induced aggregation under fluid motion), and differential settling (collection of smaller aggregates by the larger ones during sedimentation). 6 Expressing collision frequencies based on the volume (or mass) of aggregates as a representative variable, 7,8 using fractal dimension relationships and considering permeability drag effects 9 the following relationships yield: where μ is the dynamic viscosity of the suspending fluid, is the Boltzmann constant (1.381 × 10 −23 ), is the temperature (°), vi and vj are representative solid volumes of each aggregate in size classes i and j, Df is the fractal dimension of aggregates, G is the volume-averaged fluid velocity gradient or shear rate (herein it was assumed as zero), v0 is the volume of individual primary particles or smallest size class, g is the acceleration due to gravity, and Ui and Uj are the settling velocities of agglomerates at respective size classes, and are the drag coefficient correction factor of aggregates at respective size classes, defined as the ratio of drag force exerted on a permeable aggregate to drag force exerted on an impervious aggregate with the same size, 10,11 and and are the fluid collection efficiency of aggregates at respective size, defined as the ratio of flow through an aggregate to total flow approaching the aggregate. 11,12 The superposition of the three collision frequencies gives the total rate of collisions, ( , ): The drag coefficient correction factor, Ω, is given as: 10,11,13,14 where is the non-dimensional permeability of the porous aggregate given as: where is permeability [L 2 ] as described later.
The fluid collection efficiency, , can be determined from the Brinkman equation: 11,12 = 1 − − 3 (S10) where = 3 3 (1 − tanh( ) ) (S11) = − 1 ( 5 + 6 3 − tanh( ) (3 5 + 6 3 )) (S12) where 0 is the radius of primary particles and is porosity given as: 12,17,18 The representative volume of aggregates in each size class, , can be calculated as: The mass of each aggregate can then be calculated as: 18 = 0 (S19) where ρ0 is the density of the particles and mk is the mass of each aggregate in size class k. The conversion between the mass concentration Ck, and the number concentration nk, is then performed via:

Calculation of attachment rate coefficient
Two approaches are used to calculate the attachment rate coefficient, : First using a combination of colloid filtration theory (CFT) 19,20 and Derjaguin, Landau, Verwey, and Overbeek where lA and lB are the lower and the upper integration limits that will be discussed later, and is the mean value of the dimensionless total interaction energy determined from extended DLVO.
The total interaction energy [Kg m 2 S -2 ], , is calculated as a sum of three interaction energies, i.e., van der Waals attraction, , electrostatic repulsion, , and Born repulsion, : where h is the separation distance [L] between the colloid and the solid-water interface.
The sphere-plate van der Waals interaction energy can be calculated using the expression of Gregory: 25 where 123 is the combined Hamaker constant [ML 2 T −2 ] for the particle-water-surface system, and is a characteristic wavelength assumed equal to 100 nm. 21 The combined Hamaker constant for hydroxyapatite NP interacting with water-quartz interface is 2.94×10 -21 J. 26 The electrostatic repulsion energy for a sphere-plate interaction is given as: 27 where 0 is the permittivity of a vacuum which is taken as 8.854× 10 −12 C/(V·m), is the dielectric constant of the medium [-] taken as 78.5, 1 and 2 are the surface potentials (mV) for the porous media grain and the interacting particle, respectively. Although 1 and 2 are usually assumed approximately equal to the zeta potential of the two surfaces, 28 in a more accurate way they the can be determined from the following equation: 29-31 where is the zeta potential (V), ℎ is the distance between the surface of the charged particles and the slipping plane-usually taken as 5 , is the inverse of the Debye length (nm) given as: where is Avogadro's number ( where is the collision diameter assumed as 0.26 nm in order to achieve a primary minimum depth at 0.157 nm, a commonly accepted distance of closest approach, following Bradford and Torkzaban. 21 The values of lA and lB in Eq. (S23) are determined as described in detail by Bradford and Torkzaban. 21 In brief, five types of the interaction energy profile are considered, ranging from favourable to unfavourable interaction condition evaluated using 1 , 2 , and standing for primary minimum, secondary minimum, and the energy barrier in the total interaction energy profile as shown in Fig. S11. Based on these classes and ignoring the detachment, the values of lA and lB can be determined from Table S3. 21 In the alternative approach, i.e., calculation of the attachment rate coefficient according to the ANN-based empirical correlations, 22 the empirical matrices of neural network were employed in the MATLAB code developed in the present study. The code used these networks to calculate both attachment and detachment rate coefficients for each size class using 20 experimental characteristics listed in Table S1. It should be note that these empirical networks were already exerted in an MS Excel spreadsheet and were made available online as electronic supporting information of Babakhani et al. 22 All Fortran and MATLAB codes developed in the present study are available upon reasonable request to the author. :   Table S1. Parameters used in 1-D simulations following previous studies. 22           (c,f). Deposition is considered using a size-variable Katt calculated using CFT-DLVO. No acceleration factor has been included here.