MNM1D: A Numerical Code for Colloid Transport in Porous Media: Implementation and Validation

Problem statement: Understanding the mechanisms that control the trans port and fate of colloidal particles in subsurface environments is a crucial issue faced by several researchers in the last years. In many cases, natural colloids have been sh own to play a major role in the spreading of strongly sorbing contaminants, while manufactured m icro-and nanoparticles, which are nowadays widely spread in the subsurface, can be toxic thems elve . On the other hand, in recent years studies have been addressed to the use of highly reactive m icro-and nanoparticle suspensions for the remediation of contaminated aquifers. Provide the s et of partial-differential model equations and its numerical solution for the colloid transport under transient hydrochemical conditions, that have been previously shown to be extremely important in micro -and nanoparticle transport in porous media. Approach: This study presented a novel colloid transport mode l, called MNM1D (Micro-and Nanoparticle transport Model in porous media in 1D geometry), able to simulate the colloid behavior in porous media in the presence of both constant and t r sient hydrochemical parameters (namely ionic strength). The model accounts for attachment and de tachment phenomena, that can be modeled with one or two linear and/or langmuirian interaction sites. The governing equations were solved using a finite differences approach, herein presented and discusse d in d tails. Results: Both qualitative and quantitative comparisons with results of well-established colloi d transport models, based both on analytical and numerical solutions of the colloid transport equati on, were performed. The MNM1D results were found to be in good agreement with these solutions. Conclusion: The shown good agreement between MNM1D and the other models indicated that this code can r epresent in the future a useful tool for the simula tion of colloidal transport in groundwater under transient hydrochemical conditions.


INTRODUCTION
Micro-and nanoparticles, both natural (clays, oxides, bacteria) and anthropogenic (coatings, paints, chemical reactants) are widely spread in the subsurface [1,2] . In many cases, nanoparticles can be toxic themselves, or can act as vehicles of contaminants sorbed on their surface [3][4][5][6] , due to their high surface area and charge. Besides, in recent years studies have been addressed to the use of highly reactive micro-and nanoparticle suspensions for the remediation of contaminated aquifers [7] . As millimetric zerovalent iron is known to degrade a large number of organic compounds [8] , in the last years studies focused on the use of colloidal suspensions of nanoscale iron particles, that are known to have a specific area (and reaction kinetics) orders of magnitude higher than millimetric iron [9][10][11] . Knowledge of the main phenomena controlling the spreading and fate of both natural and synthesized colloidal particles in groundwater is then of a great importance, as they can represent a risk as well as a resource for human health.
The present study illustrates the MNM1D code (Micro-and Nanoparticle transport Model in porous media), an innovative colloid transport model developed by the authors of this study [12] , able to simulate colloid deposition and release in saturated porous media during transients of ionic strength. The dependence of attachment/detachment coefficients on salt concentration is explicitly embedded into the model. The model validation through both analytical (Stanmod) and numerical (Hydrus1D) solutions is reported.

MATERIALS AND METHODS
Model equations: Transport of colloids is governed by advection-dispersion phenomena and by particle-particle and particle-soil physico-chemical interactions [6] . In some cases, other important interactions can play a significant role, e.g., acid-base interactions, steric repulsions (in the presence of longchain polymers adsorbed on the particle surface) and magnetic interactions (as in the case of iron nanoparticle suspensions used for groundwater remediation) [7,11,13] .
Colloidal particle transport in saturated porous systems is usually modeled using a modified advectiondispersion equation. It describes the dual-phase, nonequilibrium interactions between particles in the liquid (water) and solid (grains) phase. Colloid deposition on the grain surface is generally referred to as attachment, colloid release as detachment. Classical filtration theory CFT (also called clean-bed filtration theory) considers first-order attachment kinetics, while detachment is supposed to be negligible [14] . However, this theory was shown not to be appropriate in many cases: non-linear kinetic attachment behaviors, namely blocking or ripening effects, are often observed and cannot be described using the CFT [3,[15][16][17][18][19][20] . Besides, one of the most relevant factors affecting particle mobility was shown to be the solution chemistry [21][22][23][24][25] . In this sense, advanced theories in modeling colloid release have been proposed by Grolimund et al. [26,27] .
To date, several 1D forms of the model have been proposed, taking into account the various phenomena described above. However, all these models can be described by a general form of the modified advectiondispersion equation, which includes a general exchange term for the liquid-solid phase transfer: The generic non-equilibrium exchange term f(c) describes the mass transfer between solid and liquid phase and can be designed in several forms, according to the mechanisms that have to be modeled. In particular, if a linear exchange term is hypothesized, f(c) becomes: Where: On the other hand, when blocking phenomena occur, i.e., when the colloid concentration in the solid phase is limited to a fixed value s max , f(c) is expressed as: If the collector grains, i.e. the porous matrix, is homogeneous, only one interaction site type is included in the model. On the contrary, if the solid matrix is supposed to be heterogeneous for the affinity to colloidal particles, two or more interaction site types can be modeled, leading to: Finally, if the influence of hydrochemical parameters, i.e., ionic strength, is to be included into the model, the partial differential equations for colloid transport are to be coupled with the advection-dispersion equation for a conservative tracer, that allows simulating transients of salt concentration. Attachment/detachment coefficients can then be tied to the tracer concentration. Empirical functions, based on laboratory results and theoretical considerations, have been proposed by the authors of this study in order to state a functional relationship between ionic strength and attachment, detachment and blocking coefficients [12] .
The MNM1D code here described solves the colloid transport problem in transient ionic strength conditions, considering two interaction sites (site S 1 , with blocking and site S 2 , with a linear attachment model). Consequently, Eq. 1 becomes: Where: c t = The conservative salt concentration [M L −3 ] s 1 = The particle concentration on solid phase [L −3 ] for blocking sites of type 1 (reversible attachment with blocking) s max,1 (c t ) = The maximum particle concentration for sites 1 and depends on the salt concentration k a,1 (c t ) and k d,1 (c t ) = The coefficients for sites 1 s 2 = The particle concentration on solid phase [L -3 ] for linear sites of type 2 k a,2 (c t ) and k d,2 (c t ) = The corresponding coefficients

Model implementation:
The micro-and nanoparticle transport model in porous media (MNM1D) code: The model Eq. 5, coupled with the functional relationships reported in [12] , was implemented in a Matlab environment using a finitedifferences model. The Euler implicit scheme was employed for time derivatives: Where: k = The index for the time step j = The index for the space discretization nodes Both upwind and central-in-space schemes were implemented for space derivatives, giving the user the possibility of choosing the most suitable solution technique for the transport problem to be simulated: where, γ = 0.5 for the central-in-space scheme and γ = 1 for the upwind scheme.
As the problem is non linear, due to the presence of the blocking function in the exchange term for interaction site of type 1, the model implements a Picard iterative scheme. Controls were introduced to avoid numerical oscillations.
The model can be used both for direct problems and for inverse solution. In the latter case, the interpolation procedure involves a large-scale fitting, based on the interior-reflective Newton method described by Coleman and Li [28] . An application of the MNM1D code for the solution of the inverse problem was presented by the researchers [12] , where it was applied for fitting a set of laboratory tests run in transient ionic strength conditions. In that case, both interaction sites were considered to be active, but the code allows to inactivate one of them when the solid matrix can be considered homogeneous.
To the authors' knowledge, no solution is available to date for the simulation of colloid transport with timeand space-varying attachment/detachment coefficients. Consequently, the MNM1D code was compared to well-established analytical and numerical models that provide results for problems with constant coefficients. First, a simplified form, with one linear interaction site (i.e., site of type 2) and constant attachment/detachment coefficients, was compared with results provided by the analytical solution of Stanmod [29] . Secondly, solutions for the complete model Eq. 5 with constant coefficients (i.e., for a two-site model, constant ionic strength) were compared with results provided by Hydrus-1D [30] .

RESULTS
The MNM1D code for one-site colloid transport problems: Comparison with Stanmod and Hydrus results: The Stanmod software [29] , based on analytical solutions of differential equations describing solute transport in porous media, was used for the validation of the following simplified formulation of the model Eq. 5: In particular, the CXTFIT code [29] was used, that gives an analytical solution for non-equilibrium transport problems in 1D geometry. The CXTFIT model is here applied for the resolution of the following non-equilibrium solute transport model: Where: α [T −1 ] = The first order kinetic rate coefficient K d [L 3 M −1 ] = The distribution (equilibrium) coefficient for linear adsorption With respect to the original CXTFIT equation, hypotheses of no equilibrium exchange sites, no solute decay in liquid and solid phase and no solute production in liquid and solid phase have been done. Re-arranging Eq. 9, the Stanmod parameters α and K d can be expressed in term of the attachment/detachment parameters of the colloid transport problem: These relations are to be used when comparing MNM and Stanmod results.
Hydrus-1D is a finite-element model able to simulate flow, heat and reactive solute transport in 1D geometry, in variably saturated porous media [30] . For the one-site model (8) the "Two Kinetic Sites Model" Hydrus-1D solution is to be applied, which is designed for virus, colloid and bacteria transport. With respect to the complete solution, in this application no generation and inactivation neither in liquid nor in solid phase were considered. As in this model two interaction sites are available, the second one was disabled and no blocking function was considered for the first one.
When using the MNM1D code, the coefficients k a,2 and k d,2 were set equal to zero and s max equal to infinite. As for Stanmod, the "Deterministic Non-equilibrium CDE" solver was used, choosing the model parameters according to [29] : Where: L = The column length v d = The darcyan velocity The following parameters were adopted for the simulations: Effective porosity n = 0.4, darcyan velocity v d = 8·10 −5 m sec −1 , hydrodynamic dispersion D = 1·10 −7 m 2 sec −1 , bulk density ρ b =1.5·10 6 g m −3 . Breakthrough curves were then calculated at an outlet distance of L = 0.1 m. Although the MNM1D code use a finite-difference scheme, while Hydrus-1D is a finiteelement-based solution, the spatial and temporal discretization was chosen consistently in both cases. In particular, a time step of ∆t = 2 s and a space discretization of ∆x = 1·10 −3 m were adopted. A central-in-space scheme was used for the MNM1D code, while a Galerkin Finite Element space weighting scheme was employed in Hydrus-1D. In both cases, a stepwise injection in x = 0 m of a solution containing colloids for a total injection time of 3600 s at constant concentration c 0 was simulated. Breakthrough curves are reported for a total time of 6000 s. Table 1 reports the transport parameters and R square values for a set of simulations run with MNM1D, Hydrus-1D and Stanmod consistent with the model Eq. 8. R square coefficients were calculated according to:  Figure 1 reports the simulated curves for two sets of attachment/detachment coefficients (sets A and D in Table 1) given by the three codes. It can be qualitatively seen that the results are in good agreement, for both cases. The corresponding R square coefficients in Table 1 confirm the consistency of the results. A sensitivity analysis was then lead on the influence of the hydrodynamic dispersion coefficient, using the attachment/detachment parameters of case D. Dispersion coefficients ranging between 1·10 −5 m 2 sec −1 and 1·10 −8 m 2 sec −1 were considered. The R square coefficients reported in Table 1 show that, in all cases, the results given by the three codes are always in good agreement, although quite smaller values of R square are found for high dispersivities.    In particular, R square is higher when comparing the MNM1D and Stanmod results, which confirms the reliability of MNM1D code, being Stanmod based on an analytic solution. For a qualitative comparison, Fig. 2 reports the breakthrough curves for cases B, C, D and E obtained with MNM1D and Stanmod.
The MNM1D code for two-site colloid transport problems: Comparison with Hydrus results: The Hydrus-1D software was used for the validation of the two-site attachment/detachment model (5) in constant ionic strength conditions (i.e., attachment, detachment and blocking coefficients are constant during the whole simulation). For two-site, non-linear colloid transport problems, as the ones described in Eq. 5, with constant attachment/detachment coefficients, the "Two Kinetic Sites Model" Hydrus-1D solution is to be applied. Consistently with the laboratory results previously presented by the researchers [12] , a number of sets of attachment/detachment coefficients were identified, each one corresponding to the "typical" behavior of negatively-charged latex colloidal particles flowing through a saturated sand column under ionic strength from 1-300 mM (i.e., respectively, from conditions where a low removal of colloids occur, to conditions that lead to a high colloid removal). The coefficients for six different ionic strength conditions are reported in Table 2. R square coefficients between the MNM1D and Hydrus-1D curves are also reported.
Good agreement was found between the simulated breakthrough curves obtained with the two codes, for all ionic strength investigated. R square values, reported in Table 2, are always very high. However, the coefficient tend to decrease with increasing ionic strength (i.e., the more the colloid behavior differs from the behavior of a tracer, the more the difference between MNM1D and Hydrus-1D results is pronounced). Figure 3 reports the curves for 1 and 100 mM ionic strength. Zooms on the descending part of the curve show how the MNM1D code gives more smoothed results, that have been shown to be more similar to the analytical solution given by Stanmod, when solving the 1-site colloid transport equation.

CONCLUSION
The research presents the numerical implementation of a finite-differences code for the simulation of colloid transport in porous media in transient ionic strength conditions, named MNM1D. Because all other available codes are implemented for constant attachment/detachment coefficients (i.e., They do not allow the simulation of transients in hydrochemical parameters), the MNM1D code was validated for constant ionic strength, comparing results with well-established software, namely Hydrus-1D and Stanmod. Comparisons were performed for both onesite (Hydrus-1D and Stanmod) and two-sites (Hydrus-1D) models. Good agreement was found for all attachment/detachment conditions and for a wide range of hydrodynamic dispersion parameters.
The MNM1D freeware software can be downloaded at the following website: www.polito.it/groundwater/software.