Improvement of numerical*simulation of bimolecular reactive transport using time-dependent parameters

LIU Yong （劉詠）, QIAN Jia-zhong （錢家忠）, MA Lei （馬雷）, ZHAO Wei-dong （趙衛(wèi)東）, LUO Qian-kun （駱乾坤）

1. School of Biotechnology and Food Engineering, Hefei University of Technology, Hefei 230009, China, E-mail: Liuyong@hfut.edu.cn

2. School of Resources and Environmental Engineering, Hefei University of Technology, Hefei 230009, China

Improvement of numerical*simulation of bimolecular reactive transport using time-dependent parameters

LIU Yong （劉詠）1, QIAN Jia-zhong （錢家忠）2, MA Lei （馬雷）2, ZHAO Wei-dong （趙衛(wèi)東）2, LUO Qian-kun （駱乾坤）2

1. School of Biotechnology and Food Engineering, Hefei University of Technology, Hefei 230009, China, E-mail: Liuyong@hfut.edu.cn

2. School of Resources and Environmental Engineering, Hefei University of Technology, Hefei 230009, China

(Received September 24, 2012, Revised January 6, 2015)

The overestimation of the contaminant concentration is a main issue in simulating the reactive transport using the common advection-dispersion-reaction equation (ADRE). To solve this problem, a new modeling method is developed. The parameters of the model are identified based on experimental data. The results of the model are compared with previous results in the literature, and the sensitivity of the model is analyzed by examining the model’s response to changes of the model parameters. Main conclusions are as follows: (1) The numerical modeling approach is feasible with an improved simulating accuracy. The predicted values are in agreement with the experimental measurements. The relative errors of the peak concentration between the simulated and experimental values are less than 2.5%. The errors are greatly reduced as compared with the results (about 67.8% as the maximum) based on a traditional ADRE in the literature. (2) There are three parameters (m,0β and D) which can be calibrated on the basis of experimental data in the model. The reactive product concentrations are mainly influenced by the parameters involved in the reactive ratio such as m and0β. The hydrodynamic dispersion coefficient D has almost no influence on the reactive product transport. (3) Our model does not provide better fitting curves for the “early arrival” and the “l(fā)ong tail”. The “early arrival” and the“l(fā)ong tail” are associated with low values of the product AB. Under these conditions, the reaction rate is close to 0, and the model of the ADRE reduces to the advection-dispersion equation (ADE). Further mechanism study is needed in the future.

numerical simulation, bimolecular reactive transport, porous media, ADRE

Introduction

The solute transport, especially, the reactive transport, characterized by a combination of physical and chemical processes, is so complicated that it becomes a challenging issue in the field of hydrogeology in order to scientifically quantify the processes[1-5]. To study the transport mechanism, the bimolecular reaction (A+B?AB) is often chosen as an elementary case of all chemical reactions in porous media. Quantitative analysis and prediction of the reactive transport processes in groundwater systems is typically accomplished by modeling tools. A common model is based on the advection-dispersion-reaction equation (ADRE). However, recent researches show that the common model of the ADRE often has difficulty in simulating the behavior of the reactive transport, because the constant sink/source term in the ADRE can over-predict the pore-scale mixing and over-predict the product concentration[2,6].

For the common model, the equations governing the solute transport for a bimolecular reaction occurring in a saturated porous medium in one-dimensional are expressed as follows:

where A, B are the reactants and C is the product, D is the hydrodynamic dispersion coefficient of the solute, v is the seepage velocity of the flow, and the sink/source terms represent the rate of the transformation of that solute, k is the chemical transformation coefficient, andAc,Bc andCc denote the concentrations of the reactants and the product, respectively. For a second-order reaction in the porous medium, the chemical transformation rate is positive for the product and negative for the reactants, and its magnitude is taken as kcAcB. When k=0, Eqs.(1) and (2) turn into non-reactive transport models. Generally, k is treated as a constant. And Eqs.(1) and (2) are regarded as the traditional (ADRE).

For the ADRE, there are two assumptions in the theoretical framework: (1) the characteristic length associated with the pore space geometry is much smaller than that of the averaging volume, and (2) the solutes undergo instantaneous and complete mixing within a control volume[2]. However, the two assumptions are often not satisfied in the reactive transport processes that involve mixing, and therefore need a different macroscopic description of the underlying processes[7,8].

Laboratory study is useful to verify existing methodologies for the reactive transport modeling. Raje and Kapoor[1]used the data from bimolecular reactive transport column experiments to deal with the inadequacy of the ADRE model. They conducted two groups of experiments under different seepage velocities, initial concentrations and input styles. In their first group of reactive transport experiments, the glass bead-filled column was saturated with aniline (AN) at the concentration of 0.5 M. and 1, 2-naphthoquinone-4-sulfonic acid (NQS) at the concentration of 0.5 M was injected as a step input. The second group experimental column was saturated with NQS at the concentration of 0.25 M, then 0.25 M AN was injected continuously. The product 1, 2-naphthoquinone-4-aminobenzene (NQAB) from the reaction between aniline and 1, 2-naphthoquinone-4-sulfonic acid was measured at the end of the column. The authors reported the breakthrough curves of NQAB at two different seepage velocities, as well as the concentration of the total product as a function of time. The measured produced NQAB was compared to the numerical solution of the traditional ADRE, derived upon the assumption that the concentration of the limiting reactant was instantaneously and completely consumed in the reaction. The authors observed that the numerical solution of the traditional ADRE significantly over-estimated the space-time evolution (in particular, the peak concentration) of NQAB in the system. The peak product concentrations simulated in the two groups were about 68% and 36% more than the observed concentrations, respectively.

To resolve the over-estimation problem, many methods related with different fields were used to simulate the reactive transport[9]. However, some recent studies show that regardless of the above mentioned limitations, a continuous formulation based on the ADRE is still useful as a basis modeling tool for the reactive transport problems[10,11]. Recently, in the context of a transport scenario involving a precipitation reaction, Tartakovsky et al.[12]and Katz et al.[13]observed that a very fine grid discretization could capture some of the features of the investigated process, although an accurate quantitative description was not obtained.

The objective of this paper is to reinterpret the experiments of Raje and Kapoor[1]on the basis of an ADRE formulation but using a kinetic reaction rate suggested in the literature[2]. A finite element method was developed to solve the resulting system of partial differential equations. The results show that the new methods provide a viable means to interpret the experiments and describe the macroscopic chemical behavior of the system. At last, the sensitivity of the model response to changes of the model parameters is analyzed.

1. Theoretical model and method

For facilitating the mathematical discussion, Eqs.(1) and (2) are merged as

where ciis the aqueous concentration of chemical species i( i=A, B, AB), Miis a stoichiometric coefficient (M=-1 for the reactants and M=1 for the reaction product), and r=kcAcBis the space-time dependent rate at which species i is produced (or removed) by the reaction.

In Eq.(3), the chemical transformation coefficient k is treated as a constant and usually calculated by non-reactive transport experiments. This is considered as the reason why over-estimation of the product concentration is caused[2]. To describe an incomplete reaction, we imbue the effects of the incomplete mixing at the pore scale in a time-dependent kinetic reaction term based on the work of Sanchez-Vila et al.[2]and Haggerty et al.[14]. The time-dependent kinetic reaction term can be written as

where t is the time for the reactive transport, the coefficient0β and the exponent m are two parameters need to be calibrated by the time-dependent experimental observations.

Substituting Eq.(4) into Eq.(3) yields the reactive transport equation

where0A and0B can be evaluated analytically using the ADE method provided by Gramling et al.[15].

In the case of the experimental column saturated with A and injected with B substance, the corresponding initial and boundary conditions for A are:

where0C is the initial concentration.The reactive transport equations for B can be obtained in the same manner. Then according to the definition A0=cA+cCand B0=cB+cC, we obtain

Finally, cC(x, t) can be obtained from the following equation

Based on Eqs.(6)-(9), a finite element modeling approach is developed to solve the bimolecular reactive transport problem to simulate the concentration of the reactants and the reaction product at any given point in space and time.

2. Interpretation of the experiments and discussions

The main interest of this study is to assess our proposed modeling method in view of its ability to capture the time-dependent behavior of the system, in terms of the evolution of the breakthrough curves, the model parameter variation and the sensitivity.

Fig.1 Predicted and measured plume profiles at the flow speed of 7.0×10-4m/s

Fig.2 Predicted and measured plume profiles at the flow speed of 9.6×10-4m/s

Table1 Simulated parameters in different flow velocities and initial concentration

We focus on the experiment performed by Raje and Kapoor[1]. Two group concentration breakthrough curves (including the measured and simulated ones) at different seepage velocities were documented for the experiment by Raje and Kapoor[1]. The same experiment is simulated for a comparison of the general modeling strategy and results. The results are listed and discussed below.

2.1 Comparison of peak concentration

To test the validity of our model for the reactive transport of AN and NQS, the data in Raje and Kapoor[1]are used in our model. As mentioned earlier, Raje and Kapoor[1]conducted the test of a bimolecular reaction of AN+NQS?NQAB in porous media under two different seepage velocities and initial concentrations for the reactants (5.0 mM for Run 1, 2.5 mM for Run 2). Our simulations are compared with the experimental and modeling results of Raje and Kapoor[1], as shown in Fig.1 and Fig.2. The experimental conditions for the simulation are listed in Table 1.

Figure 1 shows the simulated and measured plume profiles at the flow speed of 7.00×10–4m/s and 5.0 mM. Fig.1(a) is the breakthrough curve of the reactive product displayed in a rectangular coordinate system, and Fig.1(b) is the breakthrough curve of the reactive product displayed in a semi-logarithmic coordinate system. In the results of Raje and Kapoor[1], the observations could not be properly fitted with the common ADRE. The peak product concentration simulated (7.18×10-5mol/L) is 35.5% higher than the peak experimental concentration (9.73×10-5mol/L). By using the data in literature by Raje and Kapoor[1]including the initial concentration, the seepage velocity, the hydrodynamic dispersion coefficient, and the generally good simulating results with the peak concentration (7.17×10-5mol/L), our model results give a relative error of 0.14% for the peak concentration.

Figure 2 is the simulated and measured plume profiles at the flow speed of 9.60×10-4m/s and 2.5 mM. Figure 2(a) is the breakthrough curve of the reactive product displayed in an arithmetic coordinate system, and Fig.1(b) is the breakthrough curve of the reactive product displayed in a semi-logarithmic coordinate system. As shown in Fig.2, the observations could not be well fitted with the common ADRE[1]. The simulated peak product concentration simulated (1.21×10-4mol/L) is 67.8 % higher than the peak experimental concentration (2.03×10-4mol/L). By using the data in literature by Raje and Kapoor[1]including the initial concentration, the seepage velocity, the hydrodynamic dispersion coefficient, and the generally good simulating results with the peak concentration (1.24×10-4mol/L), our model gives a relative error of 2.47% for the peak concentration.

2.2 Model parameters and sensitivity analysis

Our modeling approach requires the calibration of the three parameters: m and0β involved in the reactive ratio and the hydrodynamic dispersion coefficient D. The D value was used by Raje and Kapoor[1], which was obtained by an independent nonreactive transport experiment. Thus we only need to fit the other two parameters. The fitted parameters are listed in Table 1. The value of m increases with the increase of the seepage velocity, but the value of0β decreases. These results indicate that, when the seepage velocity increases, the value ofdecreases. The smaller the value ofis, the lower the reactive rate will be. If0β value is close to 0, no reaction occurs. This result is reasonable. As compared with Ref.[2], our method has two advantages. First, while Sanchez-Vila et al.[2]gave a quantitative interpretation of the column experiment in the chemical system: CuSO+Na EDTA2?CuEDTA2-+2Na++

To examine the model response to the changes of the model parameters, a sensitivity analysis is carried out. From the original dataset obtained from the literature as explained above, the base simulation is established. From this, the sensitivity analysis is performed, to assess the effect produced by a given variation of the input data on the model output. Spider plots are used for the representation of the results of the sensitivity analysis[16]. The effect produced by a given varia-tion of the input parameters on the model output is shown in Fig.3.

Fig.3 Model sensitivities for different parameters

Figure 3 shows that the model sensitivity changes with model’s parameters. The important effect of the value of m is evident. The second important is the parameter0β. The hydrodynamic dispersion coefficient, D, has a slight influence on the model sensitivity. This is why the dispersion coefficient may be treated as a constant for different reactants and products in the literature[1].

2.3 The “early arrival” and the “l(fā)ong tail”

The “Early arrival” and the “l(fā)ong tail” are often seen in the process of non-reactive transport[17]. For reactive transport, the two phenomena were explored in previous studies, as shown in our review of the literature. The two phenomena remain to be open questions, however, in both non-reactive and reactive transports. Further work is needed to fully disclose the mechanism of these phenomena with proper modeling methods. To clearly identify the phenomena of the“early arrival” and the “l(fā)ong tail”, we suggest to use the breakthrough curves of the reactive product in a semi-logarithmic coordinate system (See Fig.1(b) and Fig.2(b)). Figure 1(b) shows that one sees strong indications of the “early arrival” and the “l(fā)ong tail”. However, the conventional model and our model do not provide a very good fitting for the concentration breakthrough curve. Figure 2(b) shows that the “l(fā)ong tail” is more markedly displayed. The two models also do not provide better fitting curves. These are not a surprising feature. The “early arrival” and the “l(fā)ong tail” are associated with low values of product AB. Under these conditions, the reaction rate is close to 0 (see Eq.(3)), and the transport equations of the product basically reduce to the ADE. The fact that an ADE-based transport model cannot properly reproduce the early arrival and the long tail is a well-known feature that has been extensively discussed in the literature[18]. However, this discrepancy did not influence significantly the ability of the model to capture the remaining salient features of the experiment, in particular, the peak concentration of the product.

3. Conclusions

The reactive solute transport through a porous laboratory column is studied by numerical simulations. The following conclusions are drawn:

(1) A numerical modeling approach, based on the formulation of the reactive transport problem and by conceptualizing the porous system as an equivalent continuum, is developed to interpret the experimental results of Raje and Kapoor[1]. The predicted values are in agreement with the experimental measurements. The relative deviations of the peak concentration between the simulated and experimental values are no more than 2.5%, when an appropriate effective firstorder kinetic term for the reaction rate is considered. Essentially, this term is included in the model as a way to account for incomplete mixing. The errors are greatly reduced as compared with the result (about 67.8% in maximum) based on a traditional ADRE in the literature.

(2) There are three parameters (m,0β and D) which can be calibrated on the basis of experimental data in the model. Results of sensitivity analysis show that reactive product concentration of the reactive transport are mainly influenced by the parameters involved in the reactive ratio such as m and0β. The hydrodynamic dispersion coefficient D has almost no influence on the reactive product transport. This result also implies that the reactive transport is different from non-reactive transport. The mechanism needs to be further studied in the future.

(3) Our model does not provide better fitting curves for the “early arrival” and the “l(fā)ong tail”. this is not a surprising feature. the “early arrival” and the“l(fā)ong tail” are associated with low values of product AB. Under these conditions, the reaction rate is close to 0, and the model of the ADRE reduces to that of the ADE. The assessment and the simulation of the “early arrival” and the “l(fā)ong tail” in the processes of reactive transport is still an open issue. Further developments might include the assessment of the potential of the model to provide further insights on physical processes of reactive transport.

In short, by comparison, we believe that our interpretation can form a basis to promote further researches to assess the potential use of continuum approaches for the quantitative interpretation of reactive transport problems at different scales.

Acknowledgments

This work was supported by the Creative Research Groups of Hefei University of Technology (Grant No. 2009HGCX0233). We would like to thank Dr. Zhang Yong ( DRI, USA) and Dr. Sun Pengtao(UNLV, USA) for their valuable assistance and Dr. Xu Wenqing for her help in manuscript.

[1] RAJE D. S., KAPOOR V. Experimental study of bimolecular reaction kinetics in porous media[J]. Environmental Science and Technology, 2000, 34(7): 1234-1239.

[2] SANCHEZ-VILA X., FERNANDEZ-GARCIA D. and GUADAGNINI A. Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation[J]. Water Resources Research, 2010, 46(12): W12510.

[3] HAN Jiang-bo, ZHOU Zhi-fang and FU Zhi-min et al. Vadose-zone moisture dynamics under radiation boundary conditions during a drying process[J]. Journal of Hydrodynamics, 2014, 26(5): 734-744.

[4] LUO Q., WU J. and SUN X. et al. Optimal design of groundwater remediation systems using a multi-objective fast harmony search algorithm[J]. Hydrogeology Journal, 2012, 20(8): 1497-1510.

[5] ZHANG Y., QIAN J. Z. and PAPELIS C. et al. Improved understanding of bimolecular reactions in deceptively simple homogeneous media: From laboratory experiments to Lagrangian quantification[J]. Water Resources Research, 2014, 50: 1704-1715.

[6] DING D., BENSON D. A. and PASTER A. et al. Modeling bimolecular reactions and transport in porous media via particle tracking[J]. Advances in Water Resources, 2013, 53: 56-65.

[7] JOSE S. C., CIRPKA O. A. Measurement of mixingcontrolled reactive transport in homogeneous porous media and its prediction from conservative tracer test data[J]. Environmental Science and Technology, 2004, 38(7): 2089-2096.

[8] TARTAKOVSKY A. M., TARTAKOVSKY G. D. and SCHEIBE T. D. Effects of incomplete mixing on multicomponent reactive transport[J]. Advances in Water Resources, 2009, 32(11): 1674-1679.

[9] EDERY Y., SCHER H. and BERKOWITZ B. Particle tracking model of bimolecular reactive transport in porous media[J]. Water Resources Research, 2010, 46(7): W07524.

[10] FERNANDEZ-GARCIA D., SANCHEZ-VILA X. and GUADAGNINI A. Reaction rates and effective parameters in stratified aquifers[J]. Advances in Water Resources, 2008, 31(10): 1364-1376.

[11] SANCHEZ-VILA X., FERNANDEZ-GARCIA D. and GUADAGNINI A. Conditional probability density functions of concentrations for mixingcontrolled reactive transport in heterogeneous aquifers[J]. Mathematical Geoscience, 2009, 41(3): 323-351.

[12] TARTAKOVSKY A. M., REDDEN G. and LICHTNER P. C. et al. Mixing-induced precipitation: Experimental study and multiscale numerical analysis[J]. Water Resources Research, 2008, 44(6): W06S04.

[13] KATZ G. E., BERKOWITZ B. and GUADAGNINI A. et al. Experimental and modeling investigation of multicomponent reactive transport in porous media[J]. Journal of Contaminant Hydrology, 2011, 120-121: 27-44.

[14] HAGGERTY R., HARVEY C. F. and FREIHERR VON SCHWERIN C. et al. What controls the apparent timescale of solute mass transfer in aquifers and soils? A comparison of experimental results[J]. Water Resources Research, 2004, 40(1): W01510.

[15] GRAMLING C. M., HARVEY C. F. and MEIGS L. C. Reactive transport in porous media: A comparison of model prediction with laboratory visualization[J]. Environmental Science and Technology, 2002, 36(11): 2508-2514.

[16] ESCHENBACH T. G. Spiderplots versus tornado diagrams for sensitivity analysis[J]. Interfaces, 1992, 22(6): 44-46.

[17] CHEN Zhou, QIAN Jia-zhong and QIN Hua. Experimental study of the non-Darcy flow and solute transport in a channeled single fracture[J]. Journal of Hydrodynamics, 2011, 23(6): 745-751.

[18] BERKOWITZ B., CORTIS A. and DENTZ M. et al. Modeling non-Fickian transport in geological formations as a continuous time random walk[J]. Reviews of Geophysics, 2006, 44: RG2003.

10.1016/S1001-6058(15)60456-5

* Project supported by the National Natural Science Foundation of China (Grant Nos. 41372245, 41272251).

Biography: LIU Yong (1968-), Female, Master, Associate Professor

QIAN Jia-zhong, E-mail: qianjiazhong@hfut.edu.cn