# Spatial moment analysis of multispecies contaminant transport in porous media

## Article information

## Abstract

Spatial moment analysis has been performed on the concentration of the first species in a multispecies solute transport in porous media. Finite difference numerical technique was used in obtaining the solute concentration. A constant continuous source of contaminant was injected at the inlet of the domain. Results suggest that the decaying of solute mass increases as the magnitude of mean fluid velocity increases. The dispersion coefficient is highly time dependent under decaying of solutes with a complex behavior of mixing of solutes. The solute mobility and mixing varies non-linearly with time during its initial period, while the same ceases with higher decay rates of the first species much faster.

**Keywords:**Dispersion, Mobility, Multispecies transport, Porous media, Spatial moments

## 1. Introduction

Movement of contaminants through fractured porous media has received significant attention over the last few decades as the high permeability fracture provides a preferential pathway for the movement of contaminants in the heterogeneous subsurface system. Even though single species transport has been the major focus of research, multispecies transport nevertheless has its own significance. Multispecies contaminant transport is considered important since multiple reactive contaminates occur in many field situations. For example, the nuclear waste can get contaminated with radioactive materials which degrade to produce many daughter products such as tetrachloroethylene (PCE) and trichloroethylene (TCE) and many byproducts [1, 2]. Analytical solutions have been provided by researchers in the past for multispecies contaminant transport in porous media [3–11]. Natarajan and Suresh Kumar [12] provided a finite difference numerical solution for the same. Moreover, Natarajan and Suresh Kumar [13] analysed the effect of non-linear sorption on multispecies solute transport in heterogeneous porous media. Recently, Natarajan [14] investigated the effect of distance dependent and time dependent dispersion coefficient on multispecies contaminant transport in porous media.

Although earlier studies have addressed the techniques that can be adopted to model the transport of multispecies contaminant transport in porous media, studies pertaining to the mobility and spreading characteristics of such solutes is yet to be studied. Spatial moment analysis is a method that is used in analysing the effect velocity, effective dispersion coefficient and dispersivity of the solutes within the system.

Spatial moment analysis was introduced by Aris [15] to study the dispersion of a solute in a fluid flowing through a tube. Later, Marle et al. [16] and Guven et al. [17] used this method to analyse transport of non-reactive solutes for steady horizontal flow in a perfectly stratified aquifer. Freyberg [18] used this method to interpret the data based on advection and dispersion of non-reactive tracers at the Borden site. Valocchi [19] examined the impact of adsorption kinetics upon the spatial moments of the depth averaged contaminant plume using the spatial moment analysis used by Horn [20] and Brenner [21]. Dagan and Cvetkovic [22] used spatial moment analysis to analyse the kinetically sorbing solute plume in a heterogeneous aquifer. Later, the method of spatial moment analysis was extended to heterogeneous porous media by several authors for single species solute transport [23–28]. While spatial moments have been used earlier for solute transport in homogeneous as well as heterogeneous porous media, the objective of the present study is to apply the method of spatial moment analysis for multi-species solute transport in a classical porous domain.

## 2. Governing Equations

The basic equation describing the single species contaminant transport for one dimensional porous media is of the form [29]

Where *c* is the concentration of the contaminant, *x* indicates the dimension, *t* is the time, *D* is the hydrodynamic dispersion coefficient, *v* is the velocity of the fluid, and *f* is the reaction rate.

When first order reaction rate is assumed, Eq. (1) becomes

Where *k* is the first order reaction rate constant. The above equation can be extended to describe the multispecies contaminant transport in porous media. For example, the chlorinated solvent contaminants such as PCE degrades to produce daughter products such as TCE, dichloroethylene (DCE) and vinyl chloride (VC) [5]

[6]

Where *c** _{i}* is the species concentration in the

*i*th generation. Species

*i*, which is produced from species i-1, also reacts to produce species i+1 and this further reacts to produce more species. The generalized form of such a sequential transport system can be described using the following equations.

Where *i* is the index of the species, *k** _{i}* is the first order reaction rates, k

_{0}= 0,

*y*is the yield coefficient and

*n*is the number of species.

Three species solute transport has been considered for this study. The numerical model for the three species transport in porous media is of the form:

The initial and boundary conditions for the above equations are as follows:

Where *C*_{1}*, C*_{2}*, C** _{3}* are the concentration of the first, second and third species;

*k*

_{1}*, k*

_{2}*, k*

*are the first order reaction rates,*

_{3}*y*

*and*

_{1}*y*

*are the stochiometric yield coefficients and*

_{2}*L*is the length of the domain. The parameters used for this model is shown in Table 1.

## 3. Numerical Method and Spatial Moment Analysis

An implicit finite difference numerical technique has been employed to model the partial differential equations. A constant continuous source of parent contaminant is injected at the inlet of the domain. This coupled system of equations is solved numerically using the upwind scheme for advection and second-order central difference finite difference scheme for dispersion. The implicit finite difference code for solving the coupled equations was written in FORTRAN 90 language and the same was validated using analytical solution. The validation plots have been provided in Figs. 1(a) and 1(b).

Having obtained the concentration distribution of the first species contaminant in the porous system at each time level, from the above numerical method, the method of spatial moments as a function of travelling time is calculated. Such spatial moments of point concentration data provide an integrated measure of the concentration field over the entire extent of the domain. In the present model, the contaminant transport parameters are considered by characterizing the three spatial moments of the concentration distribution along the fracture.

The lower order spatial moments have been obtained using a similar approach to Guven et al. [17]. The zeroth moment (M_{0}) is proportional to the total mass of the fluid in the high permeability fracture. The first spatial moment (M_{1}) describes the displacement of the centre of the mass and the second spatial moment (M_{2}) describes the spread of the deviation from the centre of mass. The expressions for the evaluation of the zeroth moment, first moment and second moment are given below.

From these moments, the effective velocity and effective dispersion coefficient can be obtained using the following expressions:

However, the above expressions are valid for pulse sources only [30]. Since the boundary condition in the inlet is assumed to be a constant continuous source, a first derivative of the concentration in the fracture is used to obtain the equivalent pulse in order to use the above expressions (12) to (15). The spatial moment analysis has been carried out only for the first species contaminant in this study.

## 4. Results and Discussion

The results from the numerical model for multispecies solute transport in porous media without sorption has been validated with analytical solution derived with different decay rate coefficients for the three species but same yield coefficient, provided by Slodicka and Balazova [9]. The dataset used for this validation has been provided in Table 2.

### 4.1. Validation of the Numerical Model

The comparison of the results from the numerical model with that of the analytical solution has been provided in Fig. 1(a) below.

It observed from Fig. 1(a) that there is a good agreement between analytical solution and the numerical solution obtained from the present numerical model. Further, the numerical model results were validated with the analytical solution derived with different decay rate coefficients as well as yield coefficients for the three species, provided by Sun et al. [6]. The dataset that was used for this validation is same as the one used for conducting this study (Table 1). Fig. 1(b) illustrates the comparison between the analytical and numerical solution.

It observed from Fig. 1(b) that there is a good agreement between analytical solution and the numerical solution obtained from the present numerical model.

### 4.2. Moments for Different Mean Fluid Velocities

Fig. 2 provides the temporal variation of zeroth moment for different mean fluid velocities varying by an order of magnitude. Zeroth moment (the area under the curve) yields the total mass of solutes present in the physical system at any given time level. For the transport of solutes by advection and dispersion in the absence of solute decay, the temporal variation of zeroth moment will yield a straight line emanating from the zeorth moment (i.e., the ordinate) equal to unity (a line parallel to abscissa). However, it can be observed from Fig. 2 that the magnitude of zeroth moment remains at unity nearly up to 10 d, after which, the zeroth moment starts declining depending upon the magnitude of mean fluid velocity. Thus, it is very evident from Fig. 2 that the solutes do not undergo decay for the first 10 d, after which, the decaying of solutes become very sensitive. It is also observed that the decaying of solutes is very minimal for the fluids with lowest mean fluid velocity, while the same is very high for the fluids with the largest mean fluid velocity. Also, the variation of solute mass after 10 d remains nearly non-linear for the all the cases discussed in Fig. 2. It can thus be concluded from Fig. 2 that the decaying of solute mass increases as the magnitude of mean fluid velocity increases. A further observation indicates that there are actually three distinct regimes of solute mass, particularly at higher mean fluid velocities. When the mean fluid velocity is relatively higher, the zeroth moment follows three distinct regimes as shown in Fig. 2. The first regime pertains to the linear transportation of solutes in the absence of solute decay (up to 10 d); the second regime pertains to the strong non-linear decaying of solute mass (nearly up to 30 d); and the third asymptotic regime again represents the linear transportation of solutes but with the maximum decaying of solutes. Thus, the method of spatial moment analysis clearly provides the details of various time levels up to which the solutes undergo different levels of decay depending upon the magnitude of mean fluid velocity. Such details are highly critical in evaluating the concentration of the contaminants at various downstream points far away from the inlet source.

Fig. 3 provides the temporal distribution of spatial first moment for different mean fluid velocities. The slope of the first spatial moment with respect to time provides the magnitude of the velocity of the solutes. It can be observed from Fig. 3 that the first moment has significant slope with time nearly up to 15 d, after which, the line reaches their respective asymptoticity, which indicates that the velocity of the solute has reached zero as all the solute mass has exhausted due to decaying of solutes. Thus, from Fig. 3, it is the first 15 d that is very critical while assessing the solute velocity. It is also observed that the solute velocity reaches zero value much earlier for the cases with lower mean fluid velocities, while the time taken for the solute velocity to reach zero is larger for the cases with the larger mean fluid velocities. This is similar to the results obtained by Natarajan and Kumar [31] for colloidal transport in fracture matrix coupled system. It can be concluded from Fig. 3 that the transportation of solutes with lower mean fluid velocity experiences the maximum decaying of solutes at the earliest, and subsequently, the velocity of the solute reaches zero value at the earliest. In a practical sense, a maximum decaying of solutes can be expected in a geological unit with a relatively lower hydraulic conductivity for a given hydraulic gradient. For example, transportation through clay will encounter more decaying of solutes than the transportation through sand.

Fig. 4 provides the temporal variation of spatial second moment for the same set of mean fluid velocities as discussed in earlier plots. Half the slope of the spatial second moment with time provides the magnitude of dispersion coefficient. It can be observed from Fig. 4 that the solutes moving with the maximum mean fluid velocity has experienced the maximum dispersion as expected. However, it can be noted that the magnitude of dispersion coefficient has reduced drastically (for example, nearly 1.75 sq.m/d @ 10 d) from its initial value of 4 sq.m/d. Thus, it is interesting to note that the decaying of solutes mitigates the intensity of solute mixing significantly. It can also be noted that the profiles have three distinct regimes: the first regime with steep slope (nearly up to 10 d) indicates a linear variation of dispersion coefficient with time; the second regime with curved profile indicating a non-linear mixing of solutes with time (between 10 and 20 d); and the third asymptotic regime representing the absence of mixing indicating the completing exhaustion of solute mass (after 20 d). This is very much similar to the observation of Suresh Kumar [32] where the effect of mixing reduced with time due to the combined effect of matrix diffusion and decay. Thus, it can be concluded from Fig. 4 that the mixing of solutes is significantly reduced by the decaying of solutes, and that the intensity of reduction in mixing is a function of mean fluid velocity.

### 4.3. Moments for Different Dispersion Coefficients

Fig. 5 provides the temporal variation of zeroth spatial moment for various dispersion coefficients. It can be observed from Fig. 5 that the solute mass does not change under a relatively lower (D = 2 m^{2}/d) dispersion coefficient, while there is a significant loss of solute mass at a later stage (after 10 d) with an increased dispersion coefficient (D = 4 & 6 m^{2}/d). Thus, a higher dispersion coefficient is associated with a higher loss of solute mass. Thus, a larger decaying of solute mass can be expected from a system with a relatively higher dispersion coefficient, i.e., with a relatively heterogeneous system, as against the conventional homogeneous system. It can be concluded from Fig. 5 that the decaying of solute mass is higher in a heterogeneous system than a homogeneous system.

Fig. 6 provides the temporal distribution of first spatial moment for various dispersion coefficients. It can be observed from Fig. 6 that all the profiles are highly non-linear during its early stage (0–10 d) followed by its asymptoticity (after nearly 15 d). The asymptotic nature reflects the absence of solute mobility resulting from complete loss of solutes from decaying, while the initial non-linear variation of profiles reflects the time dependent solute velocity as against the assumed constant solute velocity.

Fig. 7 provides the temporal distribution of second spatial moment for various dispersion coefficients. It can be observed from Fig. 7 that the behavior of mixing of solutes is quite complicated in the sense that there is no monotonous increase or decrease of dispersion coefficients with time for various constant dispersion coefficients considered. For example, the intensity of mixing of solutes for D = 2 m^{2}/d is lesser than that with D = 4 m^{2}/d, and greater than D = 6 m^{2}/d. It can be concluded from Fig. 7 that the dispersion coefficient is highly time dependent under decaying of solutes with a complex behavior of mixing of solutes.

### 4.4. Moments for Different Decay Rate Coefficients

Fig. 8 provides the temporal distribution of zeroth spatial moment for various decay rates of species 1. It can be observed from Fig. 8 that there is a significant loss of solute mass for a relatively lower decay rate of species 1 (k_{1} = 0.2/d), while there is nearly zero loss of solute mass for relatively higher decay rate of species 1 (k_{1} = 0.4 & 0.6/d). Thus, it can be concluded from Fig. 8 that the profiles with a relatively low decay rate of species 1 decays into nothing much faster than the rest of the cases implying that the profiles with relatively high decay rate of species 1 paves way for decaying of other species.

Fig. 9 provides the temporal distribution of first spatial moment for various decay rates of species 1. It can be observed from Fig. 9 that there is a monotonous declining of solute velocity as the decay rate of species 1 increase. The profiles are highly non-linear during its early stage, while they reach asymptoticity at a later stage. It can be concluded from Fig. 9 that the solute velocity varies non-linearly with time during its initial period, while the solute mobility ceases with higher decay rates of species 1 much faster.

Fig. 10 provides the temporal distribution of second spatial moment for various decay rates of species 1. It can be observed from Fig. 10 that there is a monotonous declining of solute mixing as the decay rate of species 1 increase. Similar to Fig. 9, the profiles are highly non-linear during its early stage, while they reach asymptoticity at a later stage. It can be concluded from Fig. 10 that the solute mixing varies non-linearly with time during its initial period, while the solute mixing ceases with higher decay rates of species 1 much faster.

### 4.5. Effective Dispersion Coefficient for Different Fluid Velocities, Dispersion Coefficients and Decay Rate Coefficients

Fig. 11 shows the temporal distribution of effective dispersion coefficients of the solute for different fluid velocities. The dispersion coefficient is 1.5 at the beginning of the time period for all the fluid velocities considered in this study. The dispersion coefficient then increases sharply with increment in fluid velocity at the fifth day and then decreases with time. The dispersion coefficient reaches zero just before the 20^{th} day is reached. The same phenomenon is observed for all the fluid velocities considered. Thus spatial moment analysis helps in determining the effective dispersion coefficient of the solute, which otherwise cannot be known.

Fig. 12 shows the temporal distribution of effective dispersion coefficients of the solute for different dispersion coefficients. The dispersion coefficient is different at the beginning for all the dispersion coefficients. Thereafter, an increment in the dispersion coefficient is observed. The increment is gradual for D = 2 m^{2}/d, but the increment is sharp for other dispersion coefficients. Then, the dispersion coefficients decreases with time. The effective dispersion coefficient reaches zero early for high dispersion coefficients and later for low dispersion coefficients.

Fig. 13 shows the temporal distribution of effective dispersion coefficients of the solute for different decay rate coefficients. The dispersion coefficient is 1.5 initially for all the decay rate coefficients considered in this study. The dispersion coefficient then increases sharply for all the decay rate coefficients and then decreases with time. The dispersion coefficient reaches zero just before the 30^{th} day is reached in all the cases.

## 5. Conclusions

Spatial moment analysis has been carried out on the concentration of the first species of a multispecies contaminant transport in porous media. A constant continuous source was assumed for the first species at the inlet of the domain. The concentration profile of the contaminants was obtained using implicit finite difference technique. The conclusions from this study is summarised as follows:

The decaying of solute mass increases as the magnitude of mean fluid velocity increases.

The transportation of solutes with lower mean fluid velocity experiences the maximum decaying of solutes at the earliest.

The mixing of solutes is significantly reduced by the decaying of solutes and that the intensity of reduction in mixing is a function of mean fluid velocity.

The decaying of solute mass is higher in a heterogeneous system than a homogeneous system.

The dispersion coefficient is highly time dependent under decaying of solutes with a complex behavior of mixing of solutes.

Loss of solute mass is relatively higher with lower decay rate of species 1.

The solute mobility and mixing varies non-linearly with time during its initial period, while the same ceases with higher decay rates of species 1 much faster.

The effective dispersion coefficient increases during the initial time period and then decreases with time duration for different fluid velocities, dispersion coefficients and decay rate coefficients.