Effects of Cemented Porous Media on Temporal Mixing Behavior of Conservative Solute Transport

Abstract

The cementation of porous media leads to the variation of the pore space and heterogeneity of the porous media. In this study, four porous media (PM1, PM2, PM3, and PM4) with the different radii of solid grains were generated to represent the different cementation degrees of the porous media. The direct simulations of flow and conservative solute transport in PM1–4 were conducted to investigate the influence of the cemented porous media and Peclet number (Pe) on the temporal mixing behavior. Two metrics, scalar dissipation rates (SDR) and dilution index, were employed to quantify the temporal mixing behavior. It was found that the spatial velocity variability of the flow field was enhanced as cementation degree increased. The results of the coefficient of velocity variation ($C{V}_U$) increased from 0.943 to 2.319 for PM1–4. A network consisted of several preferential flow paths was observed in PM1–4. The preferential flow enhanced the mixing of the conservative solute but had a negative influence on the mixing of the solute plume when the cemented solid grains formed several groups, and there were some stagnant regions where the flow was almost immobile. As the Pe increased, for PM1–3, the exponent of the best-fitting power law of the global SDR decreased. At the case of Pe = 400, the slope of the global SDR reduced to around −1.9. In PM4 where the preferential flow was enhanced by the cemented solid grains, the slope of the global SDR increased as the Pe increased. The global SDR results indicated that the temporal mixing behavior followed a Fickian scaling ($SDR\propto p{v}^{-1.5}$) in the early stage (Pv < 0.05), while the mixing behavior turned to be non-Fickian in the late stage. The transition time from the Fickian scaling to the non-Fickian scaling was found to be sensitive to the cementation degree of the porous media.

1. Introduction

The solute transport in porous media plays an important role in subface environment problems (e.g., contaminant degradation [1], nuclear waste disposal [2], and bioremediation [3]). In recent decades, the investigation of dilution and mixing processes in porous media has been focused on by the scholars in various fields [4,5,6,7]. The characterization of the spreading and mixing processes of the conservative solute is instrumental in understanding and assessing reactive solute transport, which is necessary for studying complex chemical biological reaction-transport in groundwater. Although many studies [8,9] have provided new insights into the mechanisms and properties of the mixing behavior at the Darcy scale or the larger field scale, little attention has been focused on porous media at the pore scale. There is always an influence of smaller-scale processes on larger scale behaviors [10]. Therefore, understanding and evaluating mixing processes at pore-scale with suitable methods is very significant [11].

The spatial heterogeneities of the pore space in aquifers lead to complex fluid flow and anomalous solute transport. The occurrence of preferential flow paths and/or cementations, which can be found in several systems, result in the non-Fickian solute transport. The cementation of solid grains is ubiquitous in natural soil and is one of the most common reasons for the variation of pore space and heterogeneity of the porous media. In many laboratory and field scale studies [12,13,14,15], results show that the solute transport followed the non-Fickian transport. By fitting the experimentally measured or directly simulated breakthrough curves of solute plumes with continuous-time random walk (CTRW) models, some studies [16,17,18] indicated a strong link between the non-Fickian transport and heterogeneity of porous media. The varied geometry of pore channels and the cementation of solid grains are very common in natural solids, resulting in the heterogeneity of porous media. Previous studies indicated that the characteristics of pore channels such as connectivity or pore particle distribution could affect the flow field and solute transport processes [19,20]. The evolution of solute mixing in porous media with different degrees of heterogeneity is caused by the changes in pore channel geometry and solid grains due to the cementation of porous media [17]. However, little attention has been focused on this problem until now.

For the heterogeneous system, the flow velocity distribution in the porous media is a fundamental issue for understanding the mixing behavior of solute plume [21]. As the mixing process is highly influenced by the flow field [22,23], molecular diffusion causes the solute molecule movement between streamlines and enhances mixing. Some recent studies characterized topological [24] and kinematic [25] features of the flow field to explore the influence of the flow field on the mixing behavior of the solute plume. In many types of research in contaminant transport, insufficient mixing could restrict mass transformation or transformation rates of contaminants. The conservative solute transport could be considered as a process composed of spreading and mixing. However, the difference between the spreading and mixing needs to be distinguished in detail [26]. Spreading is the change of solute plume shape associated with the deformation and stretching, which are caused by the variability of the flow field. Mixing indicates concentration redistribution and trends to be uniform inside the solute plume. Many approaches have been proposed to quantify and describe the mixing behavior. In homogenous porous media, second central spatial moments can be adopted as a good measure for quantifying mixing. However, it can’t accurately reflect the mixing (or dilution) processes in the heterogeneous system [27]. Kitanidis [28] proposed the dilution index that measures the volume occupied by the solute plume. Rolle et al. [27] proposed a flux-related dilution index based on the concept of dilution index to quantify dilution in steady-state flow with a continuous injection mode. Kapoor and Kitanidis [29] used concentration variance to quantify the mixing. Several scholars quantified the mixing in subsurface problems by the scalar dissipation rate (SDR), which was proposed to describe the turbulent flow and combustion originally. Mixing at pore-scale is a necessary condition for the reaction between different dissolved compounds [28,30]. Thus, quantification of mixing at pore-scale is a key to understanding and predicting the process of contaminant transport or reaction.

The main goal of this study is to investigate the influence of the cemented porous media and Peclet number (Pe) on the mixing process of conservative solute in the different cemented porous media. Four pore-scale domains with the distinct pore channel morphology were considered to represent the different degrees of the cemented porous media and named as PM1, PM2, PM3, and PM4, respectively. The computational fluid dynamics (CFD) simulations of the flow field and solute transport in the porous media were performed and effectively determined the effect of characteristics of pore channel on the flow field and the mixing behavior of conservative solute. Solute mixing processes in the different porous media (e.g., PM1–4) were quantitively analyzed by both the SDR and dilution index.

2. Materials and Methods

2.1. Generation of Cemented Porous Media

Due to the complexity of the natural porous media, it is difficult to characterize the flow and solute transport. Although many numerical simulation methods have been greatly developed in recent years, pore-scale direct simulation of flow and solute transport in three-dimensional porous media is still a challenge. However, two-dimensional pore-scale simulation can also capture most of the physical processes and features [31]. In this paper, a two-dimensional porous medium model with the given porosity was developed to simulate the fluid flow and solute transport.

In certain two-dimensional zones, circular particles with the different radii are generated to represent solid grains. By using Matlab programming, an array that represents solid grains’ random radii is generated first. The random radius distribution of the generated solid grains is assumed to follow the lognormal distribution, which is consistent with the previous study [32]. The number of solid grains n is related to the area of the zone and target porosity of porous medium ($n=S*{\varnothing }_0/(pi*r_{ave}$), where S is the area of the zone, ${\varnothing }_0$ is the target porosity, and $r_{ave}$ is the mean radii). Second, the solid grains with random radius are continuously generated in the zone one by one. To achieve this goal, the location of grain is determined by using a random procedure, which is used to generate the coordinate of the grain in the zone. It should need to be noted that the solid grains must be generated in the zone in order of radius from large to small, which could improve efficiency. Third, the real-time porosity was calculated as a grain was generated in the zone. The porosity of the generated solid grains in the two-dimensional domain is defined as

$$\begin{array}{c}\varnothing =1-{\displaystyle \sum _{i=1}^n}\frac{\pi r_i{}^2-s_i}{S}\end{array}$$

(1)

where n is the number of solid grains, $r_i$ is the radius of $i$th solid grain, $s_i$ is the area of the part of solid grain beyond the boundary, and $S$ is the area of the two-dimensional zone where the solid grains are generated. Until the porosity (∅) reached the target value, the porous medium was obtained. Then, the data of particles’ locations and radii are imported into COMSOL for next work.

In this study, 795 solid grains were randomly distributed in a two-dimensional domain of 70 × 40 mm², as shown in Figure 1. The corresponding porosity of the generated solid grains in this two-dimensional domain was 0.45. It should be mentioned that a checking algorithm during the random distribution process of the solid grains was used to avoid the overlapping solid grains. In addition, to obtain the different degrees of the cementation in the same porous media model, we assumed that a minimum distance between any two solid grains was set to be one-fifth of the mean radius. As can be observed in Figure 1, the solid grains were totally free and non-overlapping. The solid grains’ mean radius and the coefficient of variation (COV) were set to 0.8 mm and 0.2 mm, respectively. The original porous media in Figure 1 was named as PM1. The probability distribution of the radii of the solid grains is shown in Figure 2.

In order to investigate the influence of the cemented solid grains on the solute transport and mixing process, another three solid grains with the different degrees of cementation were generated on the base of PM1 by increasing the radii and referred to as PM2, PM3, and PM4, respectively. All locations of solid grains in PM1–4 were consistent. With the increase of radii, the pore space between solid grains shrunk. As a result, another two porous media (e.g., PM2 and PM3) with the porosities of 0.40 and 0.35 were obtained, respectively. It should be mentioned that the solid grain in both PM2 and PM3 were still set to be totally non-overlapping. As the radius of the solid grains increased, some solid grains tended to contact and began to cement. Once solid grains contacted, the radius of these solid grains stopped increasing, and other radii of the solid grains keep increasing by the same stage value (1/200 of mean radius). As a result, a fraction of the solid grains in PM4 was cemented, and the corresponding porosity was 0.30. Due to the cemented solid grains in PM4, the characteristics of the pore space in PM4 changed most compared to PM1. In order to quantify the degrees of the cementation of PM1–4, the number of cemented groups ($N_c$) and the average number of individual solid grains contained in a cemented group ($N_a$) were introduced as two metrics. The properties of the cementation of PM1–4 are listed in

Table 1. Properties of the cementation of PM1–4.

Porous Media	Porosity	${\mathit{r}}_{\mathit{a}\mathit{v}\mathit{e}}$ (mm)	${\mathit{N}}_{\mathit{c}}$	${\mathit{N}}_{\mathit{a}}$	${\mathit{d}}_{\mathit{a}\mathit{v}\mathit{e}}$ (mm)	${\mathit{N}}_{\mathit{g}}$
PM1	0.45	0.785	0	0	0.36355	1679
PM2	0.40	0.818	0	0	0.29635	1679
PM3	0.35	0.854	0	0	0.22507	1679
PM4	0.30	0.888	154	4.62	0.23562	1112

Note: $r_{ave}$, the average radii of particles; $N_c$, the number of cemented groups; $N_a$, the average number of individual solid grains contained in a cemented group; $d_{ave}$, the average gap between two solid grains where solute plume pass through; $N_g$, the number of gaps.

.

As can be seen in Table 1, for PM1–3, the $d_{ave}$ decreases as the radius of solid grains increase, while the $N_g$ keeps constant due to the fact that the solid grains in PM1–3 are not cemented. When the average radius of PM4 increases by 13.12%, the corresponding porosity decreases and the cementation degree (the ratio of cemented solid grains to all solid grains) becomes 89.13%—associated with the decrease of $N_g$. Although the $N_g$ decreases by 33.77% in PM4, the $d_{ave}$ slightly increases and is larger than that in PM3. This could lead to the spatial velocity variability of the flow field. The spatial velocity variability of the flow fields was quantitively analyzed by introducing a coefficient of velocity variation ($C{V}_U$) in the Section 3.2.

2.2. Flow Field and Solute Transport Models in Porous Media

The fluid in the flow field is assumed to be incompressible and not affected by gravity. The flow field is solved by using the continuity and the Navier–Stokes equations,

$$\begin{array}{c}\nabla ·\mathit{u}=0\end{array}$$

(2)

$$\begin{array}{c}\rho \left(\mathit{u}·\nabla \mathit{u}\right)-\nabla \left(\mu \nabla \mathit{u}\right)=-\nabla p\end{array}$$

(3)

where $\rho$ is the density of the fluid, $\mathit{u}=\left[u,w\right]$ is the velocity vector, $p$ is the fluid pressure, and $\mu$ is the dynamic viscosity of the fluid.

The boundary condition of the flow field is schematically illustrated in Figure 1. In each model, the left and right boundaries are set to be the first-type boundary. The steady-flow is induced from left to right by the given pressure drop. The no-slip boundary condition is applied at the solid–liquid interfaces. There is no flux across the top and bottom boundaries.

Transient solute transport is described by the advection-diffusion equation for conservative non-sorbing solute transport at the pore-scale,

$$\begin{array}{c}\frac{\partial c}{\partial t}=-\nabla ·\left(\mathit{u}c\right)+D_m{\nabla }^{2}c\end{array}$$

(4)

where $c$ is the solute concentration, $t$ is time, and $D_m$ is the molecular diffusion coefficient. The velocity vector in Equation. (4) is from the flow field based on the solution of Equations (2) and (3). A dimensionless parameter related to the time (e.g., pore volume) is introduced as follows,

$$\begin{array}{c}Pv=Qt/A\end{array}$$

(5)

where Q is the flow rate per length, and A is the area of pore space. The porous media is initially saturated with water. At time $t=0$, the solute was injected by the line source pulse with a uniform concentration $c_0$ at the inlet boundary ($x=0$). For the simulations of the solute transport in PM1–4, the line source pulse lasts from $Pv$ = 0 to $Pv$ = 0.01.

2.3. Quantification of Mixing

In this study, the mixing behavior of the conservative solute is quantified by the dilution index and the SDR, respectively. The dilution index, introduced by Kitanidis [28], is capable of quantifying the degree of the mixing and presents a measure of the porous media volume occupied by a solute plume. This metric is expressed as the exponential of the Shannon entropy,

$$\begin{array}{c}E\left(t\right)=EXP\left[-{{\displaystyle \int }}_v^{}p\left(x,t\right)Ln\left(p\left(x,t\right)\right)dV\right]\end{array}$$

(6)

where $P\left(x,t\right)$ is interpretation as the probability distribution of the location at time t, which has units of 1 over volume.

$$\begin{array}{c}p\left(x,t\right)=\frac{c\left(x,t\right)}{{{\displaystyle \int }}_v^{}c\left(x,t\right)dV}\end{array}$$

(7)

Another metric for quantifying the mixing behavior of the conservative solute is the SDR. The SDR is defined as,

$$\begin{array}{c}\chi \left(t\right)={{\displaystyle \int }}_{\mathsf{\Omega }}^{}{D}_m\nabla c\left(\mathit{x},t\right)·\nabla c\left(\mathit{x},t\right)dx\end{array}$$

(8)

The numerical calculation of the concentration gradient in Equation (8) is one of the main difficulties since the concentration gradient can vary sharply over a small distance in the highly heterogeneous porous media. For instance, in the highly heterogeneous flow field, a very fine numerical discretization for both the flow and the concentration field is required to obtain an accurate quantification of concentration gradients, which results in multiplying the amount of computation. To avoid this, Pope SB [33] suggests that the global SDR can be defined by using the integral of the squared concentration instead of the local concentration gradients. Multiplying Equation (4) by c(x, t) and integrating over the spatial domain,

$$\begin{array}{c}\frac{1}{2}\frac{\partial }{\partial t}{{\displaystyle \int }}_{\mathsf{\Omega }}^{}c{\left(x,t\right)}^{2}d\mathsf{\Omega }+\frac{1}{2}{{\displaystyle \int }}_{\mathsf{\Omega }}^{}\nabla ·\left[\mathit{u}c{\left(x,t\right)}^2\right]d\mathsf{\Omega }=\frac{1}{2}{D}_m{{\displaystyle \int }}_{\mathsf{\Omega }}^{}\nabla ·\nabla c{\left(x,t\right)}^{2}d\mathsf{\Omega }-{{\displaystyle \int }}_{\mathsf{\Omega }}^{}{D}_m\nabla c\left(\mathit{x},t\right)·\nabla c\left(\mathit{x},t\right)d\mathsf{\Omega }\end{array}$$

(9)

In this paper, the SDR was only considered within the time from the end of injection to the beginning of the solute through the outlet boundary. In this time period, there was no mass flux out of the domain. Thus, the flow field was divergence-free, and the terms involving a divergence operator in Equation (9) are zero. Then,

$$\begin{array}{c}\chi \left(t\right)=-\frac{1}{2}\frac{d{M}^2\left(t\right)}{dt}=D_m{{\displaystyle \int }}_{\mathsf{\Omega }}^{}\nabla ·\nabla c{\left(x,t\right)}^{2}d\mathsf{\Omega }\end{array}$$

(10)

where $M^2\left(t\right)$ is the concentration second moment and can be expressed as,

$$\begin{array}{c}{M}^2\left(t\right)={{\displaystyle \int }}_{\mathsf{\Omega }}^{}c{\left(\mathit{x},t\right)}^{2}d\mathsf{\Omega }\end{array}$$

(11)

By assuming that the solute transport follows the Fick’s Law, the corresponding one-dimensional analytical solution for the concentration distribution in the absence of reaction is given by

$$\begin{array}{c}{c}_0\left(x,t\right)=\frac{m_0}{\sqrt{4\pi D_{m}t}}EXP\left(-\frac{x^2}{4\pi D_m}\right)\end{array}$$

(12)

By integrating the square of Equation (12) overall domains, the corresponding concentration second moment can be expressed as

$$\begin{array}{c}{{\displaystyle \int }}_{\mathsf{\Omega }}^{}{c}_0{\left(x,t\right)}^{2}d\mathsf{\Omega }=\frac{m_0^2}{\sqrt{8\pi D_{m}t}}=M_0^{}\left(t\right)\end{array}$$

(13)

From Equation (10), the one-dimensional analytical SDR solution can thus be expressed as

$$\begin{array}{c}{\chi }_0\left(t\right)=-\frac{1}{2}\frac{d{M}_0^{}\left(t\right)}{dt}=\frac{1}{8}\frac{m_0^2}{\sqrt{2\pi D_m}}{t}^{-\frac{3}{2}}\end{array}$$

(14)

3. Result and Discussion

3.1. Model Setup

In this study, the void space in the porous media was saturated by the water with standard properties at 20 °C (e.g., $\rho$ = 998.2 kg/m³ and $\mu =1.002\times {10}^{-3}$ Pa∙s). The numerical simulation of conservative solute transport was performed. The corresponding solute molecular diffusion was set to $D_m=1\times {10}^{-9}$ m²/s. The flow field and transient solute transport models based on Equations (2)–(4) were implemented in the COMSOL Multiphysics package (COMSOL Inc., Burlington, MA, USA) using the Galerkin finite-element method. The steady-state flow field was induced by adjusting the pressure drop over the entire porous media to obtain the same Pe. The solved flow field serves as the input for the transient solute transport model. In this work, the total simulation time was 1000 s with the time step was set as 0.5 s. In order to ensure numerical stability and accuracy, the sensitivity analysis for mesh dependency was performed, and the corresponding results showed that when the porous media domain was discretized into ~196,000 triangular elements, the solutions were mesh-independent and the numerical dispersion was negligible.

The Pe is defined by $Pe={\tau }_D/{\tau }_a=\overline{u}\overline{d}/D_m$, which is the ratio of the characteristic diffusion time $({\tau }_D={\overline{d}}^2/D_m)$ to the characteristic advection time (${\tau }_a=\overline{d}/\overline{u}$), where $\overline{u}$ is the mean flow velocity in the porous media and $\overline{d}$ is the average diameter of particles. The Reynolds number is defined as $Re=\overline{u}\overline{d}/\nu$, where $\nu$ is the kinematic viscosity of the fluid. Three different Pe (e.g., Pe = 100, 200, and 400) was considered in this work, and the value was changed by adjusting the pressure drop between the inlet boundary and the outlet boundary.

3.2. Flow Fields

Figure 3 shows the flow fields in PM1–4 with the Pe = 100. The color contrast in PM1 is small and represents a relatively uniform velocity distribution. A network consisting of several preferential flow paths could be found in PM1. The velocity reaches the maximum value in those flow paths. As the radii of solid grains increase and the pore space gradually narrows, the morphology of the pore space become complicated, causing a pronounced flow velocity variability. In PM3 and PM4, the area of the dark blue color regions increases, which represents that the local flow velocity is much smaller than the average flow velocity. The flow velocity in the preferential flow path in PM3 is higher than that in PM1 and PM2. Since the cemented solid grains formed several groups of the solid grains in PM4, the preferential flow paths become more irregular than that in PM1–3. The flow velocity distribution in PM4 is the most non-uniform. In PM4, the preference flow paths, which originally existed in PM1, turn into a continuous dominant path connecting the inlet and outlet boundaries. It can be seen in PM4 that there are some stagnant regions where the fluid flow is almost immobile, especially near the top and bottom boundaries. Thus, the preference flow in PM4 is enhanced.

The $C{V}_U$ proposed by Rolle and Kitanidis [34] was selected as a global measure of spatial velocity variability of the flow field in PM1–4. The $C{V}_U$ was defined as,

$${CV}_U=\frac{{\sigma }_u}{u_{av}}=\frac{\sqrt{{\int }_{V_p}(u_x(x)-u_{x,av}){}^{2}d{V}_p/V_p}}{{\int }_{Vp}{u}_x(x)d{V}_p/V_p}$$

(15)

where $u_x$ the velocity component of the fluid in the flow direction, ${\sigma }_u$ and $u_{av}$ are the variance and average value of $u_x$, respectively, and $V_p$ is the pore space occupied by the fluid.

The properties of flow velocity distribution in PM1–4 were listed in

Table 2. The properties of flow velocity distribution in PM1–4 under the Pe = 100, 200, and 400.

Dimensionless Pe and Re	Porous Media	τa (s)	${\mathit{\sigma }}_{\mathit{u}}$ (-)	${\mathit{u}}_{\mathit{a}\mathit{v}}$ (m/s)	$\mathit{C}{\mathit{V}}_{\mathit{U}}$ (-)	$\overline{\mathit{u}}$ (m/s)	${\mathit{u}}_{\mathit{m}\mathit{a}\mathit{x}}$ (m/s)
Pe = 100 Re = 0.0996	PM1	24.622	$4.953\times {10}^{-5}$	$5.250\times {10}^{-5}$	0.943	$6.373\times {10}^{-5}$	$4.375\times {10}^{-4}$
	PM2	26.776	$5.322\times {10}^{-5}$	$4.847\times {10}^{-5}$	1.098	$6.111\times {10}^{-5}$	$5.766\times {10}^{-4}$
	PM3	29.155	$6.405\times {10}^{-5}$	$4.253\times {10}^{-5}$	1.506	$5.857\times {10}^{-5}$	$9.244\times {10}^{-3}$
	PM4	31.530	$8.296\times {10}^{-5}$	$3.577\times {10}^{-5}$	2.319	$5.632\times {10}^{-5}$	$2.724\times {10}^{-3}$
Pe = 200 Re = 0.1992	PM1	12.311	$9.906\times {10}^{-5}$	$1.050\times {10}^{-4}$	0.943	$1.275\times {10}^{-4}$	$8.749\times {10}^{-4}$
	PM2	13.388	$1.064\times {10}^{-4}$	$9.694\times {10}^{-5}$	1.098	$1.222\times {10}^{-4}$	$1.153\times {10}^{-3}$
	PM3	14.578	$1.281\times {10}^{-4}$	$8.507\times {10}^{-5}$	1.506	$1.171\times {10}^{-4}$	$1.855\times {10}^{-3}$
	PM4	15.765	$1.659\times {10}^{-4}$	$7.155\times {10}^{-5}$	2.319	$1.126\times {10}^{-4}$	$5.449\times {10}^{-3}$
Pe = 400 Re = 0.3985	PM1	6.1555	$1.981\times {10}^{-4}$	$2.100\times {10}^{-4}$	0.943	$2.549\times {10}^{-4}$	$1.749\times {10}^{-3}$
	PM2	6.6940	$2.128\times {10}^{-4}$	$1.939\times {10}^{-4}$	1.098	$2.445\times {10}^{-4}$	$2.306\times {10}^{-3}$
	PM3	7.2887	$2.562\times {10}^{-4}$	$1.701\times {10}^{-4}$	1.506	$2.343\times {10}^{-4}$	$3.696\times {10}^{-3}$
	PM4	7.8825	$3.316\times {10}^{-4}$	$1.431\times {10}^{-4}$	2.319	$2.253\times {10}^{-4}$	$1.090\times {10}^{-2}$

Note: $\overline{u}$, the average velocity; $u_{max}$, the maximum velocity.

. As shown in Table 2, the cemented porous media has a significant influence on the flow field. As the radii of solid grains increase, at the case of Pe = 100, the $u_{max}$ in PM2–PM4 are 31.79%, 110.82%, and 190.24% higher than in PM1, respectively. Combined with the plots of local velocity at the cross-sections, the locations of the high velocity are almost unchanged and where the pore spaces are relatively narrow. The increase of $u_{max}$ is due to the coarctation of pore space. However, as the solid grains become cemented (e.g., PM4), some stagnant regions where the fluid flow is almost immobile appear and result in the $\overline{u}$ and $u_{x,av}$ decrease. The $\overline{u}$ in PM4 is 88.37% of that in PM1. The $u_{x,av}$ in PM4 is only 68.13% of that in PM1. Such opposite trend (the increase of $u_{max}$ and the decrease of $\overline{u}$) indicates that the flow velocity distribution becomes more and more non-uniform during the cementation of the porous media. Meanwhile, the value of ${\sigma }_u$ increases during the cementation of the porous media. This reflects that the preferential flow is enhanced, which is consistent with the feature reflected in Figure 3.

For a given porous media, the $C{V}_U$ is constant under the different Pe while it can be seen that the $C{V}_U$ is considerably higher in PM4 where the calculated $C{V}_U$ is approximately three times than in PM1. As the radii of the porous media increase, the $C{V}_U$ increases. This indicates that the heterogeneity of the flow field is sensitive to the variation of the pore space caused by the cementation of the porous media. This also implies that the $C{V}_U$ is capable of describing the magnitude of the heterogeneity of the flow field in porous media. In this study, the magnitude of the heterogeneity of flow field was defined by the $C{V}_U$. Since the $C{V}_U$ of PM4 is larger than PM1–3, PM1, PM2-PM3, and PM4 is considered as a weakly, mildly and highly heterogeneous porous media, respectively.

3.3. Characteristic of Solute Transport

In order to investigate the influence of the cemented porous media on the temporal solute mixing process, the conservative solute transport simulations were conducted in PM1–4 with the Pe = 100, 200 and 400.

3.3.1. Plume Dilution

The dilution index in PM1–4 was calculated at the Pe = 100, 200 and 400, respectively, as shown in Figure 4. In general, the cemented porous media has a significant influence on the dilution index, while the Pe has little influence on the dilution index for a given porous media. At the first stage (Pv < 0.1) as shown in Figure 4b,c,f–h, the dilution index in PM1 and PM2 are a bit greater than that in PM3 and PM4. As the Pv increases, the dilution index monotonically increases. The value of the dilution index in the highly heterogeneous porous media (e.g., PM4) begins to exceed that in the mildly heterogeneous porous media (e.g., PM2 and PM3). It is worth noting that the dilution index in PM3 is greater than in PM4 during the entire process. Although the porosity of PM4 is the smallest among all porous media, the corresponding peak value of the dilution index in PM3 is the largest (see

Table 3. The peak value of the dilution index in PM1–4 with Pe = 100, Pe = 200, and Pe = 400.

Dimensionless Pe	The Peak Value of Dilution Index
	PM1	PM2	PM3	PM4
Pe = 100	$5.7914\times {10}^{-4}$	$5.0943\times {10}^{-4}$	$6.1390\times {10}^{-4}$	$5.4365\times {10}^{-4}$
Pe = 200	$5.8992\times {10}^{-4}$	$5.9883\times {10}^{-4}$	$6.1525\times {10}^{-4}$	$5.2547\times {10}^{-4}$
Pe = 400	$6.0294\times {10}^{-4}$	$6.0676\times {10}^{-4}$	$6.1612\times {10}^{-4}$	$5.1070\times {10}^{-4}$

). Considering the analysis of the flow field, this indicates that the enhanced preferential flow in PM4 has a negative influence on the dilution of the solute plume. At the same Pv condition, the ratio of the fluid volume that flowed into the porous media to the pore space volume is the same. The solute transports mainly along with the fluid flow. It could be found that, in Figure 5, the concentration distribution for the Pe = 400 is almost same with it for the Pe = 100. Since the dilution index is a measure of porous media volume occupied by a solute plume. Therefore, the change in Pe has little influence on the dilution index for the same porous media.

The zoomed snapshots of the local solute transport in PM1 and PM4 were represented in Figure 5. Due to the increase of the radii of solid grains in PM4, some pore space forms the immobile regions. The local flow velocities in these immobile regions are much smaller than the average fluid flow velocity. For the solute transport via the dominant paths, it is difficult to occupy these regions in a short period of time. This leads to the fact that the dilution index is lower than that in the porous media with high connectivity (e.g., PM2 and PM3). Since the residual solute is captured in these local pores, the dilution index in the highly heterogeneous porous media (e.g., PM4) is relatively high in the late period. Therefore, the enhanced preferential flow has a negative influence on the dilution process in the highly heterogeneous porous media. When the Peclet number was increased, the concentration distribution in certain porous media for the Pe = 400 is almost same with it for Pe = 100. But the concentration gradient of the plume is slightly pronounced.

3.3.2. Global SDR

The global SDR was determined for the conservative solute transport in all porous media. Figure 6 shows the evolution of the global SDR in the porous media with the different pore space morphologies for the Pe = 100, 200, and 400. In Figure 6, the temporal SDR $x\left(t\right)$ is normalized by $x_0\left({\tau }_a\right)$. In all cases, the global SDR in general decreases with the Pv.

In the early stage (Pv < 0.05), the global SDR follows the Fickian scaling $x\left(Pv\right)\propto p{v}^{-1.5}$, as shown in Equation (15). As the Pv increases, the mixing behavior turns to be non-Fickian, which is caused by the fact that the transverse diffusion as new interfaces are created by the shearing of the plume. The transition time from the Fickian scaling to the non-Fickian scaling is important for judging the mixing process. It is found that the transition time is mainly controlled by the degree of the cementation of porous media. However, the Pe has little influence on the transition time. In general, the transition time for the weakly heterogeneous porous media (e.g., PM1) is greater than for highly heterogeneous porous media, as shown in Figure 6. The value of the transition time in PM1 under the case of Pe = 100 is nearly two times than in PM4. This reflects a better mixing process in the weakly heterogeneous porous media at the beginning. The result is consistent with the previous results of the dilution index.

Each SDR in Figure 6 was fitted by a power-law function of Pv. It is found that in all porous media, both the Pe and the cemented porous media have a significant influence on the temporal evolution of the SDR. For a given Pe, the exponent of the best-fitting power law decreases as the porosity of the porous media decreases. However, the exponent in PM4 increases slightly compared to that in PM3. As the Pe increases, for PM1–3, the exponent of the best-fitting power law decreases. In PM4 where the preferential flow was enhanced by the cemented solid grains, the slope of the global SDR increases with the Pe increase. This trend against the cases in PM1–3 indicates that the cementation of porous media has a significant influence on the temporal evolution of solute mixing process, which is consistent with the previous results reflected by the dilution index.

Furthermore, there is a fluctuation of the global SDR around the corresponding power-law fitting line for each porous media. For the highly heterogeneous porous media (e.g., PM4), the temporal evolution of the global SDR fluctuates more than in the other three cases (e.g., PM1–3). While the magnitude of this fluctuation is insensitive to the Pe. This suggests that both the narrowed pore space and the cemented solid grains, which results in the spatial velocity variability in the flow field and the geometry deformation of the front interface between solute and water, have a significant influence on the temporal evolution of the global SDR. Since the streamlines in the highly heterogeneous porous media are deformed and bent, the flow velocity variability is greater than in other cases. The flow velocity variability leads to stretching the plume and smearing out the concentration gradients. In addition, the velocity fluctuations create a potential chance for mixing by generating the disorder and large interfacial areas between the resident and invading water [35]. At the local area, the concentration gradients are smeared out by diffusion and transformed into effective mixing.

In addition, as shown in Figure 6, the exponent of each best-fitting power-law is less than −1.5, reflecting that the temporal mixing is non-Fickian. Therefore, using the one-dimensional analytical SDR will underpredict the true SDR [36], especially at the case of the high Pe or the highly heterogeneous porous media (e.g., PM4). In the non-Fickian regime, the slope of the global SDR reflects the dimensionality of the mixing process [37]. The result shows that the mixing in porous media after transition time is intermediate between the 1D and 2D behaviors.

4. Summary and Conclusions

In this study, in order to simulate the cementation of the porous media, PM1–4 were generated by increasing the radii of solid grains. The corresponding porosity of PM1–4 decreases from 0.45 to 0.30. The flow fields in PM1–4 were solved by using the continuity and the Navier–Stokes equations. The $C{V}_U$ was introduced as a global measure to evaluate the spatial velocity variability of the flow field in PM1–4. It was found that the cemented porous media had a significant influence on the spatial velocity variability of the flow field. A network consisted of several preferential flow paths was found in PM1–4. The $C{V}_U$ increased from 0.943 to 2.319 as the porosity of PM1–4 decreased from 0.45 to 0.30. The preferential flow in PM4 was enhanced since the cemented solid grains formed several groups of the solid grains.

The simulations of the conservative solute transport were performed in PM1–4 by considering the Pe = 100, 200, and 400. The temporal mixing behavior of conservative solute transport in PM1–4 was quantitively analyzed by two metrics, the dilution index and the SDR, respectively. Both the dilution index and the SDR results showed that the cemented porous media had a significant influence on the temporal mixing behavior of conservative solute. The preferential flow enhanced the mixing of the conservative solute but had a negative influence on the mixing of the solute plume when the cemented solid grains formed several groups of the solid grains and there were some stagnant regions where the fluid flow was almost immobile (e.g., PM4).

The global SDR results showed that both the Pe and the cemented porous media had a significant influence on the temporal evolution of the global SDR. For a given Pe, the exponent of the best-fitting power law decreased as the porosity of the porous media decreased. However, the exponent for PM4 increased slightly compared to that for PM3. As the Pe increased, for PM1–3, the exponent of the best-fitting power law decreased. In PM4, where the preferential flow was enhanced by the cemented solid grains, the slope of the global SDR increased as the Pe increased. Furthermore, compared to the scaling of the one-dimensional analytical solution of SDR, the global SDR results indicated that the temporal mixing behavior followed a Fickian scaling in the early stage (Pv < 0.05), while the mixing behavior turned to be non-Fickian in the late stage. In addition, it was found that the cemented porous media had a significant influence on the transition time from the Fickian scaling to the non-Fickian scaling.