Bimolecular Reactive Transport Experiments and Simulations in Porous Media

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.


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 A + B → 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.

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 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 2 PO 4 , 10.7 mmol/L NaHPO 4 , 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.

= /
(1) In Equations (1) and (2), V and refer to Darcy average velocity and molecular diffusion coefficient of a solute in water (0.001 mm 2 /s for NQAB), respectively; is characteristic length which can be expressed by the average particle size; and are the initial concentration and the reaction rate of reactants, respectively. According to the experiments, and were set as 0.2 mmol/L and 0.438 mM −1 s −1 , 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.   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.
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 2 /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 −1 s −1 , 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.

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 ψ(t). The joint distribution of ψ(s, t) is used in the recursion relation in Equation (3) to explain the possibility of reaching s at time t by different paths for particles.
Then, the probability distribution function of the particle at site s and time t, P(s, t), is given in Equation (4); Ψ (t) represents the possibility of a given site at time t in Equation (5).
where u is the Laplace variable; c(s, u) is the Laplace-transformed generalized concentration and c(s, 0) is the initial concentration; t is characterized time; M(u) is the memory function; V ψ and D ψ are effective migration velocity and dispersion coefficient, respectively; p(s) is the probability density function of random jumping length (s). For the bimolecular reaction (A+B → 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(t), 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, M(u) . 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].
Anomalous transport which is affected by reactive transport can be described by ψ(t) that is the core of the CTRW. Here, the truncated power law (TPL) in Equation (13) is used to depict anomalous transport behavior.
where t 1 and t 2 represent average transfer time and truncation time of particles and τ 2 = t 2 /t 1 ; β is a dimensionless parameter. As for β, there is Fick diffusion for β ≥ 2 and anomalous diffusion for 0 < β < 2. The smaller β is, the more obvious anomalous behavior is [30]. The bimolecular reaction migration process is solved by the finite element method, in which c(s, t) was divided into the value function c J (t) and the shape function N J (s) (Equation (14)).
Then equation was discretized and the weakened formula in Equation (15) was obtained by using the general Galerkin method.
where I J (t) = M(t) c J (t) is convolution function, A IJ , B IJ and Q I are discrete mass, convection and dispersion respectively; S J (t) represents discrete source and sink terms. To avoid excessive storage of time step results in convolution term, Prony series is used to replace M(t). More details can be found in Ben-Zvi, et al. [31].

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 −1 s −1 . 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 2 ) and root mean square error (RMSE).

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.
Water 2020, 12, x FOR PEER REVIEW 6 of 12 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. The correlation coefficients (R 2 ) 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 β.   Simulation parameters of reactive transport in porous media are displayed in Table 2. The correlation coefficients (R 2 ) 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 β.

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.

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.

Analysis of Parameters
In this paper, CTRW-FEM is proven to simulate reactive transport processes successfully presence of limited mixing and anomalous diffusion. However, the relevant fitted parameters r unclear under multiple conditions. Pe and Da can express the transport mechanism of r importance of advection, reaction and diffusion by varying the flow rate and particle size.

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. 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.

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

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 2 /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. The correlation coefficients are greater than 0.99.

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 2 /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. The correlation coefficients are greater than 0.99.

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 2 /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. The correlation coefficients are greater than 0.99.

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.

Conflicts of Interest:
The authors declare no conflicts of interest.