Bimolecular Reactive Transport Experiments and Simulations in Porous Media

Abstract

For reactive transport process in porous media, limited mixing and non-Fickian behavior are difficult to understand and predict. To explore the effects of anomalous diffusion and limited mixing, the column-based experiments of bimolecular reactive migration were performed and simulated by the CTRW-FEM model (continuous time random walk-finite element method). Simulated parameters were calibrated and the correlation coefficients between modeled and observed BTCs (breakthrough curves) were greater than 0.9, indicating that CTRW-FEM can solve over-prediction and tailing problems effectively. Porous media with coarser particle size show enhanced mixing and the non-Fickian behavior is not affected by particle size. β (a parameter of CTRW-FEM) and Da (Damköhler number) of CTRW-FEM under different Pe (Péclet number) values showed logarithmic linear relationship. Model sensitivity analysis of the CTRW-FEM model show that the peak concentration is most sensitive to the average pore velocity and the arriving peak time of peak concentration is most sensitive to β. These findings provide a theoretical basis for handling mixing and non-Fickian behavior patterns under actual environmental conditions.

1. Introduction

In porous media, reactive pollutants like heavy metals, organic compounds and nitrogen may undergo anomalous diffusion [1,2,3,4] and limited mixing [5], because transport patterns are affected by boundary conditions, media characteristics and reactive processes collectively. Such behaviors can lead to the complicated concentration fluctuation in space and time, which make prediction difficult.

The traditional advection-dispersion-reaction equation (ADRE) which is commonly used to model reactive transport processes often overestimates the peak concentration compared with practical experiments’ measurements of reactive transport and fails to portray the effects of anomalous diffusion [6,7,8]. With consideration of boundedness of ADRE, Eulerian continuum models and Lagrangian particle tracking approaches have been proposed to solve over-prediction and tailing problems [9].

Based on the ADRE equation, Eulerian continuum models are developed to upscale the transport equation from pore scale to larger scales and they can be implemented by various approaches such as volume-averaging [10,11], effective upscale parameter [12,13] techniques and the spatial Markov model (SMM) [14]. However, it is still a challenge to across the scales of larger transport and smaller reaction in heterogeneous media [15].

As the common ones of Lagrangian particle tracking approaches, the continuous time random walk (CTRW), the time domain random walk TDRW and discrete-time random walk (RW) have been applied to model the solute transport behavior [16,17]. Among the different approaches, the CTRW-PT (particle tracking) is considered the effective method to describe Fickian and non-Fickian behavior in the reactive-diffusion process [18,19,20]. Nevertheless, the accuracy of PT is closely related to the number of setting particles, and generally the quantity of particles needed is much larger than the number available.

In addition to the PT coupled with the CTRW model, the partial differential governing equations of CTRW (CTRW-PDE) has been widely used in non-Fickian diffusion studies of non-reactive and first-order degradation solutions [5,21]. As for nonlinear bimolecular reaction, the finite element method was applied to solve numerical solution of CTRW partial differential governing equations (CTRW-FEM) for the first time [22]; then, the simulations were developed to compare with the results of experimental measurement to evidence the method [23].

For a simple reaction process, representative bimolecular reactions in the form of $\mathrm{A}+\mathrm{B}\to \mathrm{C}$ [24,25,26,27] were chosen to explore reactive transport patterns which were affected by chemical reactions and fluid dynamics at pore scale. In previous researches, Raje and Kapoor [25] and Gramling, Harvey and Meigs [24] conducted one-dimensional homogeneous glass column experiments to explore mixing effects; moreover, Anna, et al. [28] conducted homogeneous porous medium equipment by soft lithography and analyzed the reaction rate and mixing process, and based on the work of Raje and Kapoor [25], Qian, Zhan, Zhang, Sun and Liu [27] replaced glass with quartz sand and conducted 100 cm column experiments.

It’s well-known that the particle size is an essential determinant in solute transport for its influence on seepage and reactive process. However, the above-mentioned experiments rarely explore the effects of variable grain size on the limited-mixing and non-Fickian behavior. Meanwhile, needs for the practical application of CTRW-FEM for modeling bimolecular reactions arise. In this study, bimolecular reactive transport experiments were conducted in porous media under variable media particle sizes. First, CTRW-FEM was applied to model the bimolecular reactive column experiments. Then, the concentration fluctuations, especially for over-prediction and tailing based on the analytical results of experiments and modeling, were further discussed. Finally, analysis of parameters and model sensitivity were conducted to further application of CTRW-FEM.

2. Materials and Methods

2.1. Laboratory Experiments

The typical bimolecular reaction between 1,2-naphthoquinone-4 sulfonic acid (NQS) and aniline (AN) was chosen as the object of this study, and the product of this reaction is 1,2-naphthoquinone-4-aminobenzene (NQAB). During the reaction, the color of the solution changes obviously. Specifically, AN is colorless and NQS is yellow, while reaction product (NQAB) is red. The reaction in this study is instantaneous and irreversible.

As shown in Figure 1, the pore medium was a plexiglass column that was 30 cm high and 4 cm in diameter. The injection port and sampling port were located above and below the soil column, respectively. Soil columns were loaded depending on the density of media in

Table 1. Bimolecular reactive migration in porous media.

Flow Rate (mL/min)	Media	Density (g/cm3)	V (cm/s)	Porosity	Da	Pe
40.063	fine	1.64	0.131	0.406	4.42	29.48
	medium	1.70	0.134	0.396	49.16	44.98
	coarse	1.76	0.137	0.389	196.65	61.31
60.094	fine	1.64	0.200	0.406	4.42	108.00
	medium	1.70	0.202	0.396	49.16	151.03
	coarse	1.76	0.205	0.389	196.65	208.78
80.126	fine	1.64	0.261	0.406	4.42	224.40
	medium	1.70	0.268	0.396	49.16	307.35
	coarse	1.76	0.273	0.389	196.65	422.00

and purified water was filled to remove gas. Each group of experiments consists of two parts. Firstly, AN was injected until the concentration of AN in the column reached the stable initial value, and then NQS was injected to carry the reactive transport experiment. The initial concentrations of AN and NQS were 0.2 mmol/L, and buffer solutions (9.3 mmol/L KH₂PO₄, 10.7 mmol/L NaHPO₄, and 13.6 mmol/L NaCl) were added to maintain pH at 7 during the whole experiments. Samples were collected every 30 seconds and the absorbances at 276 nm, 370 nm and 476 nm of the liquid were measured by spectrophotometry (752N, Shanghai Precision & Scientific Instrument Co., Ltd., China). Then sample concentrations were calculated using the Beer–Lambert law.

In order to explore the reactive migration process under different media, fine (<0.45 mm), medium (0.5–1 mm) and coarse (1–2 mm) were served at the at flow rates of 40.063, 60.094 and 80.126 mL/min. The flow rates were controlled by a peristaltic pump. Before the experiment, the density and porosity of different media were measured (Table 1).

The pore-scale transport mechanisms (advection, molecular diffusion and reaction) can be quantified by Péclet number (Pe) and Damköhler number (Da) [7]. Pe represents the ratio of the time scale between the convection and molecular diffusion, and Da is the ratio of the time scale between the reaction and molecular diffusion.

$$Pe=Vl/D_d$$

(1)

$$Da=l^2{C}_{o}k/D_d$$

(2)

In Equations (1) and (2), V and $D_d$ refer to Darcy average velocity and molecular diffusion coefficient of a solute in water (0.001 mm²/s for NQAB), respectively; $l$ is characteristic length which can be expressed by the average particle size; $C_o$ and $k$ are the initial concentration and the reaction rate of reactants, respectively. According to the experiments, $C_o$ and $k$ were set as 0.2 mmol/L and 0.438 mM⁻¹ s⁻¹, respectively.

The values of Da for fine, medium, and coarse media were 4.42, 49.16 and 196.65, respectively. For the same media, the value of Da remained equal and greater than 1. These results indicate that the reaction-transport process is affected by the mixing degree [29]. The values of Pe in Table 1 range from 29.48 to 422.0, which expresses mechanical dispersion dominate in solute dispersion compared to molecular diffusion.

2.2. Numerical Simulation

2.2.1. CTRW-FEM Theory

In the theory of CTRW, reactive transport process can be considered by a series of particles undergoing transitions of random jumping length (s) and random wait time (t), and thus CTRW effectively captures the concentration fluctuations caused by heterogeneous media and mixing processes [1,18]. In general, the formulations of partial differential equation (PDE) and particle tracking (PT) are allowed for solutions in CTRW. Here, the CTRW-PDE formulation coupled with nonlinear reactive term which is solved by the finite element method (CTRW-FEM) is applied for bimolecular reactive transport [22]. A brief introduction of CTRW-FEM is sketched to explain the related mechanism.

The CTRW equation is controlled by the diffusion of particles including the probability density functions of p (s) and$\psi \left(t\right)$. The joint distribution of $\psi \left(s,t\right)$ is used in the recursion relation in Equation (3) to explain the possibility of reaching s at time t by different paths for particles.

$$R\left(s,t\right)={\displaystyle \sum }_{s^{\prime }}{{\displaystyle \int }}_0^t\psi \left(s-s^{\prime },t-t^{\prime }\right)\mathrm{R}\left(s^{\prime },t^{\prime }\right)d{t}^{\prime }$$

(3)

Then, the probability distribution function of the particle at site s and time t, $P\left(s,t\right)$, is given in Equation (4); Ψ (t) represents the possibility of a given site at time t in Equation (5).

$$P\left(s,t\right)={{\displaystyle \int }}_0^t\Psi \left(t-t^{\prime }\right)\mathrm{R}\left(s,t^{\prime }\right)d{t}^{\prime }$$

(4)

$$\Psi \left(t\right)=1-{\displaystyle \sum }_{s^{\prime }}{{\displaystyle \int }}_0^t\psi \left(s^{\prime },t^{’}\right)d{t}^{\prime }$$

(5)

Assuming $\psi \left(s,t\right)$= p(s)$\psi \left(t\right)$, the partial differential governing equation of CTRW can be written by Taylor expansion to extend discrete scale to continuous scale in Equations (6)–(9) [15].

$$u\tilde{c}\left(s,u\right)-c\left(s,0\right)=\tilde{M}\left(u\right)\left[-v_{\psi }·\nabla \tilde{c}\left(s,u\right)+D_{\psi }{\nabla }^2\tilde{c}\left(s,u\right)\right]$$

(6)

$$\tilde{M}\left(u\right)=\overline{t}u\frac{\tilde{\psi }\left(u\right)}{1-\tilde{\psi }\left(u\right)}$$

(7)

$$V_{\psi }=\frac{1}{\overline{t}}{{\displaystyle \int }}_{V}p\left(s\right)s ds$$

(8)

$$D_{\psi }=\frac{1}{\overline{t}}\frac{1}{2}{{\displaystyle \int }}_{V}p\left(s\right)ssds$$

(9)

where $u$ is the Laplace variable; $\tilde{c}\left(s,u\right)$ is the Laplace-transformed generalized concentration and $c\left(s,0\right)$ is the initial concentration; $\overline{t}$ is characterized time; $\tilde{M}\left(u\right)$ is the memory function; $V_{\psi }$ and $D_{\psi }$ are effective migration velocity and dispersion coefficient, respectively; $p\left(s\right)$ is the probability density function of random jumping length (s).

For the bimolecular reaction (A+B$\to \mathrm{C}$) with the reaction rate ($K$), solutions of A, B and C are defined as three types of particles. Under the same memory function $M\left(t\right)$, the chemical reaction is assumed to be independent on the flow process and then the reaction term is obtained by Pull change and Taylor expansion, $\tilde{M}\left(u\right)\frac{K{c}_A\left(s,t^{\prime }\right)c_B\left(s,t^{\prime }\right)}{\overline{t}}$. Bimolecular reactive transport equation of CTRW can be written as Equations (10)–(12). The numerical solution process is discussed in Ben-Zvi, Nissan, Scher and Berkowitz [22].

$$\mathrm{u}\tilde{c}\left(s,u\right)-c\left(s,0\right)=\tilde{M}\left(u\right)\left[F\left\{c_i\right\}\left(s,u\right)-\frac{K{c}_A\left(s,u\right)c_B\left(s,u\right)}{\overline{t}}\right]d{t}^{\prime },i=A,B$$

(10)

$$\mathrm{u}\tilde{c}\left(s,u\right)-c\left(s,0\right)=\tilde{M}\left(u\right)\left[F\left\{c_i\right\}\left(s,u\right)+\frac{K{c}_A\left(s,u\right)c_B\left(s,u\right)}{\overline{t}}\right]d{t}^{\prime },i=\mathrm{C}$$

(11)

$$F\left\{c_i\right\}\left(s,u\right)=\left[V_{\psi }·\nabla c_i\left(s,u\right)-D_{\psi }{\nabla }^2{c}_i\left(s,u\right)\right]$$

(12)

Anomalous transport which is affected by reactive transport can be described by $\psi \left(t\right)$ that is the core of the CTRW. Here, the truncated power law (TPL) in Equation (13) is used to depict anomalous transport behavior.

$$\tilde{\psi }\left(u\right)\equiv {\left(1+{\tau }_{2}u{t}_1\right)}^{\beta }\Gamma \left(-\beta ,{\tau }_2{}^{-1}+u{t}_1\right)/\Gamma \left(-\beta ,{\tau }_2{}^{-1}\right), 0<\beta <2$$

(13)

where $t_{1 }$and $t_2$ represent average transfer time and truncation time of particles and ${\tau }_2=t_2/t_1$; β is a dimensionless parameter. As for $\beta$, there is Fick diffusion for β $\ge$ 2 and anomalous diffusion for 0 < β < 2. The smaller $\beta$ is, the more obvious anomalous behavior is [30].

The bimolecular reaction migration process is solved by the finite element method, in which $c\left(s,t\right)$ was divided into the value function $c_J\left(t\right)$ and the shape function$N_J\left(s\right)$ (Equation (14)).

$$c\left(s,t\right)=c_J\left(t\right)N_J\left(s\right)$$

(14)

Then equation was discretized and the weakened formula in Equation (15) was obtained by using the general Galerkin method.

$$A_{IJ}\left(c_J\left(t\right)\right)+B_{IJ}{I}_J\left(t\right)+Q_I\left(t\right)-A_{IJ}{S}_J\left(t\right)=0$$

(15)

where $I_J\left(t\right)=M\left(t\right)⨂c_J\left(t\right)$ is convolution function, $A_{IJ}$, $B_{IJ}$and $Q_I$ are discrete mass, convection and dispersion respectively; $S_J\left(t\right)$ represents discrete source and sink terms. To avoid excessive storage of time step results in convolution term, Prony series is used to replace $M\left(t\right)$. More details can be found in Ben-Zvi, et al. [31].

$$M\left(t\right)\cong a_0\delta \left(t\right)+{\displaystyle \sum }_{p=1}^P{a}_p{e}^{-b_{p}t}=a_0\delta \left(t\right)+{\displaystyle \sum }_{p=1}^P{M}_p\left(t\right)$$

(16)

2.2.2. Modeling Inputs and Simulation

Based on the theoretical analysis of continuous stochastic time random walking, CTRW-FEM MATLAB package is used to simulate the breakthrough curves in laboratory experiments (http://www.weizmann.ac.il/EPS/People/Brian/CTRW/software). In the process of building the model, columns were generalized into a one-dimensional domain (L = 30 cm). Initially, AN was uniformly distributed in the domain at a concentration of 0.2 mmol/L. At the inlet (x = 0), the initial boundary was the Dirichlet boundary (C_NQS= 0.2 mmol/L, C_AN = C_NQAB = 0). At the outlet (x = L), the initial boundary was the Neumann boundary (zero derivatives for all species). The column was uniformly divided into 500 elements. The simulations used a time step of 0.1 s and proceed until 1500 s. The reaction coefficient was set at 0.438 mM⁻¹ s⁻¹. In addition, the average void velocity (V), dispersion coefficient (DL) and (β) were adjusted by trial and error based on the concentration of NQAB. Then, optimal values were obtained by correlation coefficient (r²) and root mean square error (RMSE).

3. Results and Discussion

3.1. Breakthrough Curves

To analyze the concentration fluctuations under different media particle size, the curves of the outlet concentration change with time of product NQAB (40.063, 60.094, and 80.126 mL/min) are shown in Figure 2. The shapes of the breakthrough curves (BTCs) show similar rules. Initially, the column was filled with AN and the mass of product NQAB increased as NQS was injected continually. Subsequently, the mass of product NQAB decreased as AN was consumed. In the order of coarse media, medium media and fine media, the peak value became lower. That demonstrates concentration fluctuations in the reactive migration process are related to the particle size.

Simulation parameters of reactive transport in porous media are displayed in

Table 2. Fitting parameters of CTRW-FEM for reactive migration in porous media.

Flow Rate (mL/min)	Media	V(cm/s)	DL(cm2/s)	β	R2	RMSE
40.063	fine	0.131	0.244	1.65	0.991	0.0076
	medium	0.134	0.274	1.7	0.993	0.0027
	coarse	0.137	0.294	1.8	0.992	0.0024
60.094	fine	0.200	0.374	1.58	0.992	0.0023
	medium	0.202	0.387	1.63	0.996	0.0024
	coarse	0.205	0.395	1.69	0.993	0.0030
80.126	fine	0.261	0.415	1.5	0.907	0.0030
	medium	0.268	0.581	1.53	0.994	0.0019
	coarse	0.273	0.614	1.62	0.999	0.0007

. The correlation coefficients (R²) between the experimental and simulated values are greater than 0.9, and the values of RMSE are less than 0.01. These excellent fitting results indicate that CTRW-FEM can be applied to model 1-D bimolecular reaction process successfully. For the same flow rate, finer media leads to slower velocity (V), smaller dispersion coefficient and larger β.

3.2. Reactant Mixing and Anomalous Diffusion

In porous media, reactive transport mechanisms of limited mixing and anomalous diffusion are often influenced by spatial heterogeneity, especially for preferential flow [32,33,34]. In a completely mixed and homogeneous reactive system, the product concentration predicted by ADRE evolves to a Gauss distribution and the peak concentration is half of the maximum reactant concentration (0.1 Mm).

However, the peak concentration of NQAB in Figure 2 is approximately 0.06 Mm. In the order of coarse media, medium media and fine media, the peak value becomes lower and the true mixing is weaker. The smaller particle size produced a larger porosity and thus the smaller pore velocity (Table 1). This behavior is related to the fluid and reactive mechanism in pores scale (the values of Pe and Da) [7]. In Table 1, the values of Pe and Da increase with particle size, which expresses the convection and reaction is faster than diffusion.

In Figure 3, the typical non-Fickian behavior of late tailing is observed in the semi-logarithmic coordinate curve. For the whole tailing period, it is found that the CTRW-FEM can fit the measured concentration in the early age. However, the CTRW-FEM can’t capture accurately concentration fluctuation in the late age (the concentration of NQAB is less than 0.001 Mm). The difference may be contributed to the measured method used in the study and the spectrophotometry having restricted measurement-precision. For different media, the tailings evaluate similarly. This indicates that non-Fickian behavior in sand columns is not affected by media particles. This might be due to the uniform packing columns and the similar pore structure in the experiments. Therefore, the variation of β there is not caused by anomalous diffusion and related to the limited-mixing.

3.3. Analysis of Parameters

In this paper, CTRW-FEM is proven to simulate reactive transport processes successfully in the presence of limited mixing and anomalous diffusion. However, the relevant fitted parameters remain unclear under multiple conditions. Pe and Da can express the transport mechanism of relative importance of advection, reaction and diffusion by varying the flow rate and particle size. To estimate model parameters effectively, β and DL under different value of Pe and Da are shown in Figure 4 and Figure 5. For the same Da, β decreases with Pe and DL increases with Pe. The double logarithmic graphs in Figure 4 and Figure 5 express the explicit mathematical relationship between β and Pe or DL and Pe. The coefficients of equations are identical for the same media and laws can be used to fit optimal parameters based on Pe. When Pe and Da grow together, β and DL increase. However, the exact mathematical relationship was not found.

3.4. Model Sensitivity Analysis

The effects of input parameters on concentration fluctuation are pivotal to apply CTRW-FEM further. To analyze the simulation response of CTRW-FEM to the changes of input values, a sensitivity analysis of three parameters, namely, the average void velocity (V), dispersion coefficient (DL) and β were performed. According to experimental conditions, the range of V, DL and β were set as 0.1–0.3 cm/s, 0.1–0.7 cm²/s and 1–2, respectively. Ten points were uniformly selected for the local parameter sensitivity analysis, and the relationship between the peak concentration or arrival peak time of product NQAB and the variation in the parameters was analyzed. After running this procedure, the concentration values of the peaks and the time to the peaks were obtained for the different simulations.

Figure 6 shows the sensitivity analysis results for the three parameters. It can be observed that V was the most sensitive parameter for the peak concentration, followed by β and DL. For the arrival peak time, β is the most sensitive parameter, followed by V and DL. Increases of the average pore velocity and dispersion coefficient reduce peak concentration and arrival peak time. This can be explained by the drastic dilution caused by convection. Increases of β enhance peak concentration and arrival time, this can be attributed to the media segregation effects caused by heterogeneity. The fitting functions with different parameter variations are in power form and are described in

Table 3. Fitting functions with different parameter.

Parameter	C/Cmean	T/Tmean
β	$\mathrm{y}=1.001-1.89\times e^{-7.474x}$	$\mathrm{y}=0.6235-2375.2\times e^{-9.11x}$
V	$\mathrm{y}=0.635+0.721\times e^{-0.682x}$	$\mathrm{y}=0.529-3.329\times e^{-1.97x}$
DL	$\mathrm{y}=1.01-0.0032\times e^{-1.168x}$	$\mathrm{y}=0.462-0.65\times e^{-0.189x}$

. The correlation coefficients are greater than 0.99.

4. Conclusions

To explore reactive transport mechanisms in porous media, experiments of bimolecular reactive transport were conducted under various media particle sizes. The CTRW-FEM model was used to simulate the transport process and match the experimental results, and a model sensitivity analysis was conducted for the CTRW-FEM model. The following conclusions are drawn from this research:

(1)

Breakthrough curves of product NQAB in the column-based and numerical experiments are indicative of incomplete mixing and non-Fickian behavior. The correlation coefficients are greater than 0.9. In general, CTRW-FEM can be used to solve over-estimated peak and tail problems in BTCs.

(2)

It was found that coarser particles can lead to enhanced true mixing. Non-Fickian behavior in sand columns for different particle sizes changes slightly.

(3)

β decreases with Pe and DL increases with Pe for the same Da; β and DL increase as Pe and Da grow together. There is a clear linear relationship between lg β and lg Pe or lg Da and lg Pe.

(4)

The results of the sensitivity analysis show that peak concentration is most sensitive to V and the arrival peak time is most sensitive to β.