A comparison of solution schemes for modeling reactive transport with nonlinear sorption Anshuman Singh 1 , Xuan Huang 1 , Richelle Allen-King 2 , Alan J Rabideau 1 1 Department of Civil, Structural and Environmental Engineering, 2 Department of Geology, University at Buffalo, NY 14260. AGU Fall Meeting 2011, Paper Number: H33G-1402 (c) Split-operator with non equilibrium sorption using fast kinetics: This is a two- step solution scheme where solution method for first step, i.e. advection and diffusion, is similar to (a). In the second step, reaction is assumed as first order non equilibrium sorption (Equation IV a & b). First order reaction rate constant based on value of Damköhler number(D a ) >10 results in fast reaction similar to equilibrium conditions. The Isotherm The non-linear dual mode Polanyipartion isotherm was selected for the numerical experiments which is represented by the equation (V). Where, q is the sorbed phase mass (μg/Kg) and the parameters for the isotherm for the numerical experiment were selected as Q 0 =1081μg/Kg; A= -0.1456; B=1.876 from S W = 150000 μg/L; K P =0.1060 L/Kg (from unpublished laboratory data for soil from Borden site for perchloroethylene (PCE)) The considerations for time and space descritization for the numerical experiments were based on a range of Courant number (Cr) and mesh Péclet number (Pe). The true solution for the ADR equation was assumed at low Cr value of 0.005. Péclet number was kept less than 2 for numerical accuracy. The results are presented in Figure I to III. Nonlinear sorption isotherms can be easily incorporated into existing transport codes via the simple 2-step split-operator scheme, with suitable attention to the numerical time step. For the test problem, the three solution schemes produce equivalent results for similar Courant number and mesh Péclet number. However, the computational times for the fast kinetic solution scheme are larger compared to the other two algorithms. For the test problem, the solutions obtained with Cr < 0.1 – 0.5 approached the “true” solution in all cases, with the precise threshold dependent on the mesh Pe. Additional numerical experiments (not shown) indicate that the “lagged retardation factor” approach can be much less accurate for highly nonlinear isotherms. Strategic Environmental Research and Development Program (SERDP). Department of Civil, Structural and Environmental Engineering and Department of Geology, University at Buffalo. A flexible algorithm has been developed to model contaminant transport influenced by “dual mode” and other complex nonlinear sorption processes. The algorithm is designed to support straightforward modification of existing advective-dispersive-reactive models for groundwater transport. In the new time- stepping algorithm, the advective-dispersive equation is solved assuming no sorption, followed by the solution of one simple non-linear algebraic mass balance equation at each spatial node. Implementation is demonstrated using the public domain MOUSER software. Results of numerical experiments are presented to illustrate the efficiency and accuracy of the algorithm as a function of time step, including comparisons with commonly used “fast kinetic” and “lagged retardation factor” numerical procedures. The MOUSER Model Models implemented in MOUSER are applicable to any physical system described by the one-dimensional form of the advective-dispersive-reactive equation (ADRE) and the specified boundary conditions (equation I) where C and q represent the aqueous and sorbed concentrations. (a)Split operator algorithm with equilibrium sorption: This is a two step process where the advection and diffusion portion of equation I is solved to obtain an intermediate concentration C* at each node using finite element method implemented in MOUSER. The second step involves equilibrating the aqueous and sorbed phases for each node, i, using the following algebraic equation representing mass conservation of each node. (b) Lagged Retardation factor: The lagged retardation factor algorithm in MOUSER is based on a one step solution of ADRE using retardation factor for reaction isotherm based on the concentrations from the previous time step. Model, Algorithms and Isotherm Prologue (V) 10 log 0 C K Q q P B C S A W + = − Results of Numerical Experiment 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0 5 10 15 20 25 30 PCE concentration(mg/L) Distance(m) Cr=0.1 Cr=0.5 Cr=1.0 True solution Pe=0.62 Figure I: Comparison of ADRE solution using split-operator algorithm with equilibrium sorption 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 5 10 15 20 25 30 PCE concentration(mg/L) Distance(m) Cr=0.1 Cr=0.5 Cr=1.0 True solution Pe=0.31 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0 5 10 15 20 25 30 PCE concentration(mg/L) Distance(m) True solution Cr=0.1 Cr=0.5 Cr=1.0 Pe=0.62 Figure II: Comparison of ADRE solution using lagged retardation factors. 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0 5 10 15 20 25 30 PCE concentration(mg/L) Distance(m) True solution Cr=0.1 Cr=0.5 Cr=1.0 Pe=0.31 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 5 10 15 20 25 30 PCE concentration(mg/L) Distance(m) Cr=1.0 Cr=0.5 Cr=0.1 True solution Pe=0.62 Figure III: Comparison of ADRE solution using non-equilibrium fast kinetic solution scheme 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0 5 10 15 20 25 30 PCE concentration(mg/L) Distance(m) Cr=1.0 Cr=0.5 Cr=0.1 True solution Pe=0.31 Conclusions Acknowledgements (I) Reaction Diffusion and Advection 2 2 t q n x C D x C v t C b ∂ ∂ − ∂ ∂ + ∂ ∂ − = ∂ ∂ ρ (II) . . * 1 1 t i b i t i b t i q nC q nC ρ ρ + = + + + (III) 1 2 2 x C D x C v t C C q n b ∂ ∂ + ∂ ∂ − = ∂ ∂ ∂ ∂ + ρ b) IV ( ) ( and a); IV ( ) ( i eq i i eq b i q q t q q q n t C i i − − = ∂ ∂ − − = ∂ ∂ α ρ α Retardation Factor