Effects of microarrangement of solid particles on PCE migration and its remediation in porous media

Groundwater can be stored abundantly in granula-composed aquifers with high permeability. The microstructure of granular materials has important effect on the permeability of aquifers and the contaminant migration and remediation in aquifers is also influenced by the characteristics of porous media. In this study, two different microscale arrangements of sand particles are compared to reveal the effects of microstructure on the contaminant migration and remediation. With the help of fractal theory, the mathematical expressions of permeability and entry pressure are conducted to delineate granular materials with regular triangle arrangement (RTA) and square pitch arrangement (SPA) at microscale. Using a sequential Gaussian simulation (SGS) method, a synthetic heterogeneous site contaminated by perchloroethylene (PCE) is then used to investigate the migration and remediation affected by the two different microscale arrangements. PCE is released from an underground storage tank into the aquifer and the surfactant is used to clean up the subsurface contamination. Results suggest that RTA can not only cause more groundwater contamination, but also make remediation become more difficult. The PCE remediation efficiency of 60.01–99.78 % with a mean of 92.52 and 65.53–99.74 % with a mean of 95.83 % is achieved for 200 individual heterogeneous realizations based on the RTA and SPA, respectively, indicating that the cleanup of PCE in aquifer with SPA is significantly easier. This study leads to a new understanding of the microstructures of porous media and demonstrates how microscale arrangements control contaminant migration in aquifers, which is helpful to design successful remediation scheme for underground storage tank spill.

1 Introduction

Groundwater is an essential natural resource for water supply to domestic, agricultural, and industrial activities as well as ecosystem health (Boswinkel, 2000; Valipour, 2012, 2015; Yannopoulos et al., 2015; Valipour and Singh, 2016). Unfortunately, with the rapid development of economic activities such as mining, agriculture, landfills and industrial activities (Bakshevskaia and Pozdniakov, 2016; Cui et al., 2016; H. Liu et al., 2016; An et al., 2016; Shen et al., 2017), more and more contaminants released from human activities are contaminating the precious groundwater resource and subsurface environment (Dawson and Roberts, 1997; Liu, 2005; Hadley and Newell, 2014; Carroll et al., 2015; Essaid et al., 2015; Huang et al., 2015; Y. Liu et al., 2016; Schaefer et al., 2016; Weathers et al., 2016). Among the contaminants detected in groundwater, dense nonaqueous phase liquids (DNAPLs) such as perchloroethylene (PCE) and other polycyclic aromatic hydrocarbons (PAHs), are highly toxic and carcinogenic (Dawson and Roberts, 1997; Hadley and Newell, 2014). When DNAPLs are released into aquifer from underground storage tanks, they will infiltrate through the entire aquifer and form residual ganglia and pools of DNAPLs due to their large densities, high interfacial tension, and low solubility. The residual ganglia and pools of DNAPLs can serve as long-term sources of groundwater contamination which are harmful to the subsurface environment and human beings (Bob et al., 2008; Liang and Lai, 2008; Liang and Hsieh, 2015). Consequently, it is very important to explore DNAPL migration in aquifers and associated remediation of groundwater contamination.

When DNAPLs migrate in aquifers on a macroscopic scale, the transport properties such as permeability, diffusivity and dispersivity are closely related to the aquifer's microstructures and can affect DNAPL behavior (Yu and Li, 2004; Yu, 2005; Yun et al., 2005; Feng and Yu, 2007; Yu et al., 2009). Therefore, characterizing the effect of microstructures on macroscopic properties is a key point of research on heterogeneity of porous media (Mishra et al., 2016). In the classical Kozeny–Carman equation, the permeability K is related to porosity n, surface area S and the Kozeny constant c, where c is affected by the porosity, solid particles and microgeometric structures (Bear, 1972; Yu et al., 2009). According to fractal theory, natural porous media can be treated as fractal objects (Pfeifer and Avnir, 1983; Katz and Thompson, 1985; Krohn, 1988). For example, the tortuosity of flow paths in porous media is deeply studied by various proposed fractal models (Yu and Cheng, 2002; Yu et al., 2009; Cai et al., 2010), indicating the effectiveness of fractal methods. Based on fractal concepts, mathematic models are proposed to depict the permeability and invasion of fluids in some special porous media (Yu and Cheng, 2002; Yu et al., 2009; Cai et al., 2010). Furthermore, fractal method is also used to explore the effect of microstructure of biological media on associated thermal conductivity while this kind of material has a complex randomly distributed vascular-tree structure at microscale (Li and Yu, 2013).

In this study, we focus on the effect of microarrangement of sand particles on macroscopic DNAPL migration and associated remediation for underground storage tank spills. With the help of fractal theory, the microstructures of two different microscale arrangements of sand particles are explored. Afterwards, the mathematical relationships between porosity, permeability, and entry pressure are derived for regular triangle arrangement (RTA) and square pitch microscale arrangement (SPA). An idealized heterogeneous contaminated site is generated using a sequential Gaussian simulation (SGS) method. An underground storage tank releases PCE into heterogeneous aquifer composed of granular material. After a long-term migration, PCE contamination is alleviated using a surfactant remediation method. A multicomponent, multiphase model simulator, UTCHEM, is then used to simulate the entire process of DNAPL migration and remediation. Effects of arrangements of sand particles on migration and remediation of DNAPLs are comparatively analyzed based on the simulations to reveal how the microstructure of porous media controls the contaminant migration and remediation on a macroscopic scale.

2 Methodology

2.1 Fractal models of two different microscale arrangements of sand particles

The porous media can be treated as the bundle of tortuous capillary tubes, and the relationship between the diameter and the length of the capillary tube is as follows (Yu and Cheng, 2002):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{L}_{\mathrm{t}}\left(\mathit{\lambda }\right)={\mathit{\lambda }}^{\mathrm{1}-D_{\mathrm{t}}}{L}_{\mathrm{s}}^{D_{\mathrm{t}}},\end{array}$$

where L_s is the straight length between the tortuous flow path's end point, λ is the diameter of the capillary tube, and D_t is the fractal dimension of tortuosity for porous media, $\mathrm{1}<D_{\mathrm{t}}<\mathrm{2}$ (Yu and Cheng, 2002).

If an infinitesimal element consisting of a bundle of tortuous capillary tubes from porous media is selected, the total number of capillary tubes in infinitesimal element can be calculated by the power–law relation:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }N(L\ge \mathit{\lambda })={\left({\displaystyle \frac{{\mathit{\lambda }}_{max}}{\mathit{\lambda }}}\right)}^{D_{\mathrm{f}}},\end{array}$$

where D_f is the fractal dimension for pore areas in porous media, $\mathrm{1}<D_{\mathrm{f}}<\mathrm{2}$ (Yu and Cheng, 2002); λ_max is the maximum diameter of capillary tubes.

Afterward, the derivative of Eq. (2) can be achieved:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }-\mathrm{d}N(L\ge \mathit{\lambda })=D_{\mathrm{f}}{\mathit{\lambda }}_{max}^{D_{\mathrm{f}}}{\mathit{\lambda }}^{-(D_{\mathrm{f}}+\mathrm{1})}{\mathrm{d}}_{\mathit{\lambda }}.\end{array}$$

The total number of capillary tubes in an infinitesimal element can be derived from Eq. (3):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{N}_{\mathrm{t}}(L\ge {\mathit{\lambda }}_{min})={\left({\displaystyle \frac{{\mathit{\lambda }}_{max}}{{\mathit{\lambda }}_{min}}}\right)}^{D_{\mathrm{f}}},\end{array}$$

where λ_min is the minimum diameter of capillary tubes.

Dividing Eq. (3) by Eq. (4) can achieve the following:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }-{\displaystyle \frac{{\mathrm{d}}_{N(L\ge \mathit{\lambda })}}{N_{\mathrm{t}}}}=D_{\mathrm{f}}{\mathit{\lambda }}_{min}^{D_{\mathrm{f}}}{\mathit{\lambda }}^{-(D_{\mathrm{f}}+\mathrm{1})}{\mathrm{d}}_{\mathit{\lambda }}=f\left(\mathit{\lambda }\right){\mathrm{d}}_{\mathit{\lambda }},\end{array}$$

where f(λ) is the probability density function, $f\left(\mathit{\lambda }\right)=D_{\mathrm{f}}{\mathit{\lambda }}_{min}^{D_{\mathrm{f}}}{\mathit{\lambda }}^{-(D_{\mathrm{f}}+\mathrm{1})}$.

The probability density function satisfies the following relationship:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}f\left(\mathit{\lambda }\right){\mathrm{d}}_{\mathit{\lambda }}=\mathrm{1}-{\left({\displaystyle \frac{{\mathit{\lambda }}_{min}}{{\mathit{\lambda }}_{max}}}\right)}^{D_{\mathrm{f}}}.\end{array}$$

Considering $(\frac{{\mathit{\lambda }}_{min}}{{\mathit{\lambda }}_{max}}{)}^{D_{\mathrm{f}}}=\mathrm{0}$, the above Eq. (6) becomes the following:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}f\left(\mathit{\lambda }\right){\mathrm{d}}_{\mathit{\lambda }}=\underset{{\mathit{\lambda }}_{min}}{\overset{{\mathit{\lambda }}_{max}}{\int }}f\left(\mathit{\lambda }\right){\mathrm{d}}_{\mathit{\lambda }}=\mathrm{1}-{\left({\displaystyle \frac{{\mathit{\lambda }}_{min}}{{\mathit{\lambda }}_{max}}}\right)}^{D_{\mathrm{f}}}=\mathrm{1}.\end{array}$$

With regard to fluid flow in capillary tubes, the flow rate Q can be calculated by the Hagen–Poiseulle equation:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }Q={\displaystyle \frac{\mathit{\pi }{\mathrm{r}}^{\mathrm{4}}\mathrm{\Delta }P}{\mathrm{8}\mathit{\mu }{L}_{\mathrm{s}}}}={\displaystyle \frac{\mathit{\pi }{\left(\frac{\mathit{\lambda }}{\mathrm{2}}\right)}^{\mathrm{4}}\mathrm{\Delta }P}{\mathrm{8}\mathit{\mu }{L}_{\mathrm{s}}}}={\displaystyle \frac{\mathit{\pi }{\mathit{\lambda }}^{\mathrm{4}}\mathrm{\Delta }P}{\mathrm{128}\mathit{\mu }{L}_{\mathrm{s}}}},\end{array}$$

where μ is fluid's viscosity, and ΔP is the pressure gradient across the capillary tube.

The differentiation of flow rate of capillary tubes is as follows (Yu and Cheng, 2002):

$$\begin{array}{ll}{\displaystyle }{\mathrm{d}}_q& {\displaystyle }=[-{\mathrm{d}}_{N(L\ge \mathit{\lambda })}]{\displaystyle \frac{\mathit{\pi }{\mathit{\lambda }}^{\mathrm{4}}\mathrm{\Delta }P}{\mathrm{128}\mathit{\mu }{L}_{\mathrm{t}}\left(\mathit{\lambda }\right)}}\\ {\displaystyle }& {\displaystyle }=D_{\mathrm{f}}{\mathit{\lambda }}_{max}^{D_{\mathrm{f}}}{\mathit{\lambda }}^{-(D_{\mathrm{f}}+\mathrm{1})}{\mathrm{d}}_{\mathit{\lambda }}\cdot {\displaystyle \frac{\mathit{\pi }{\mathit{\lambda }}^{\mathrm{4}}\mathrm{\Delta }P}{\mathrm{128}\mathit{\mu }{L}_{\mathrm{t}}\left(\mathit{\lambda }\right)}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathit{\pi }}{\mathrm{128}}}{\displaystyle \frac{\mathrm{\Delta }P}{\mathit{\mu }}}{\displaystyle \frac{D_{\mathrm{f}}{\mathit{\lambda }}_{max}^{D_{\mathrm{f}}}}{L_{\mathrm{t}}\left(\mathit{\lambda }\right)}}{\mathit{\lambda }}^{\mathrm{3}-D_{\mathrm{f}}}{\mathrm{d}}_{\mathit{\lambda }}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathit{\pi }}{\mathrm{128}}}{\displaystyle \frac{\mathrm{\Delta }P}{\mathit{\mu }}}{\displaystyle \frac{D_{\mathrm{f}}{\mathit{\lambda }}_{max}^{D_{\mathrm{f}}}}{{\mathit{\lambda }}^{\mathrm{1}-D_{\mathrm{t}}}{L}_{\mathrm{s}}^{D_{\mathrm{t}}}}}{\mathit{\lambda }}^{\mathrm{3}-D_{\mathrm{f}}}{\mathrm{d}}_{\mathit{\lambda }}\\ \text{(9)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathit{\pi }}{\mathrm{128}}}{\displaystyle \frac{\mathrm{\Delta }P}{\mathit{\mu }}}{\displaystyle \frac{D_{\mathrm{f}}{\mathit{\lambda }}_{max}^{D_{\mathrm{f}}}}{L_{\mathrm{s}}^{D_{\mathrm{t}}}}}{\mathit{\lambda }}^{\mathrm{2}+D_{\mathrm{t}}-D_{\mathrm{f}}}{\mathrm{d}}_{\mathit{\lambda }}.\end{array}$$

Integrating the individual flow rate from λ_min to λ_max can achieve the total flow rate (Yu and Cheng, 2002):

$$\begin{array}{ll}{\displaystyle }Q=& {\displaystyle }\int {\mathrm{d}}_q=\underset{{\mathit{\lambda }}_{min}}{\overset{{\mathit{\lambda }}_{max}}{\int }}{\displaystyle \frac{\mathit{\pi }}{\mathrm{128}}}{\displaystyle \frac{\mathrm{\Delta }P}{\mathit{\mu }}}{\displaystyle \frac{D_{\mathrm{f}}{\mathit{\lambda }}_{max}^{D_{\mathrm{f}}}}{L_{\mathrm{s}}^{D_{\mathrm{t}}}}}{\mathit{\lambda }}^{\mathrm{2}+D_{\mathrm{t}}-D_{\mathrm{f}}}{\mathrm{d}}_{\mathit{\lambda }}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathit{\pi }}{\mathrm{128}}}{\displaystyle \frac{\mathrm{\Delta }P}{\mathit{\mu }}}{\displaystyle \frac{D_{\mathrm{f}}}{\mathrm{3}-D_{\mathrm{f}}+D_{\mathrm{t}}}}{\displaystyle \frac{\mathrm{1}}{L_{\mathrm{s}}^{D_{\mathrm{t}}}}}{\mathit{\lambda }}_{max}^{D_{\mathrm{f}}}\\ {\displaystyle }& {\displaystyle }\cdot \left({\mathit{\lambda }}_{max}^{\mathrm{3}-D_{\mathrm{f}}+D_{\mathrm{t}}}-{\mathit{\lambda }}_{min}^{\mathrm{3}-D_{\mathrm{f}}+D_{\mathrm{t}}}\right)\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathit{\pi }}{\mathrm{128}}}{\displaystyle \frac{\mathrm{\Delta }P}{\mathit{\mu }}}{\displaystyle \frac{D_{\mathrm{f}}}{\mathrm{3}-D_{\mathrm{f}}+D_{\mathrm{t}}}}{\displaystyle \frac{\mathrm{1}}{L_{\mathrm{s}}^{D_{\mathrm{t}}}}}{\mathit{\lambda }}_{max}^{\mathrm{3}+D_{\mathrm{t}}}\\ \text{(10)}& {\displaystyle }& {\displaystyle }\cdot \left[\mathrm{1}-{\left({\displaystyle \frac{{\mathit{\lambda }}_{min}}{{\mathit{\lambda }}_{max}}}\left){}^{D_{\mathrm{f}}}\right({\displaystyle \frac{{\mathit{\lambda }}_{min}}{{\mathit{\lambda }}_{max}}}\right)}^{\mathrm{3}+D_{\mathrm{t}}-\mathrm{2}{D}_{\mathrm{f}}}\right].\end{array}$$

Due to $\mathrm{1}<D_{\mathrm{t}}<\mathrm{2}$ and $\mathrm{1}<D_{\mathrm{f}}<\mathrm{2}$, $\mathrm{3}+D_{\mathrm{t}}-\mathrm{2}{D}_{\mathrm{f}}>\mathrm{0}$. Simultaneously, $(\frac{{\mathit{\lambda }}_{min}}{{\mathit{\lambda }}_{max}}{)}^{D_{\mathrm{f}}}\cong \mathrm{0}$, $\mathrm{0}<(\frac{{\mathit{\lambda }}_{min}}{{\mathit{\lambda }}_{max}}{)}^{\mathrm{3}+D_{\mathrm{t}}-D_{\mathrm{f}}}<\mathrm{1}$. Therefore, Eq. (10) can be simplified as follows:

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }Q=\int {\mathrm{d}}_q={\displaystyle \frac{\mathit{\pi }}{\mathrm{128}}}{\displaystyle \frac{\mathrm{\Delta }P}{\mathit{\mu }}}{\displaystyle \frac{D_{\mathrm{f}}}{\mathrm{3}-D_{\mathrm{f}}+D_{\mathrm{t}}}}{\displaystyle \frac{\mathrm{1}}{L_{\mathrm{0}}^{D_{\mathrm{t}}}}}{\mathit{\lambda }}_{max}^{\mathrm{3}+D_{\mathrm{t}}}.\end{array}$$

Substituting Darcy's law, $Q=\frac{kA\mathrm{\Delta }P}{\mathit{\mu }{L}_{\mathrm{0}}}$, in Eq. (11) will obtain the permeability of porous media:

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }k={\displaystyle \frac{\mathit{\pi }}{\mathrm{128}}}{\displaystyle \frac{D_{\mathrm{f}}}{\mathrm{3}+D_{\mathrm{t}}-D_{\mathrm{f}}}}{\displaystyle \frac{L_{\mathrm{0}}^{\mathrm{1}-D_{\mathrm{t}}}}{A}}{\mathit{\lambda }}_{max}^{\mathrm{3}+D_{\mathrm{t}}}.\end{array}$$

To obtain the fractal dimension of tortuosity D_t, the expression of tortuosity (τ) can be obtained from Eq. (1):

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }\mathit{\tau }={\displaystyle \frac{L_{\mathrm{t}}\left(\mathit{\lambda }\right)}{L_{\mathrm{s}}}}={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{1}-D_{\mathrm{t}}}{L}_{\mathrm{s}}^{D_{\mathrm{t}}}}{L_{\mathrm{s}}}}={\left({\displaystyle \frac{L_{\mathrm{s}}}{\mathit{\lambda }}}\right)}^{D_{\mathrm{t}}-\mathrm{1}}.\end{array}$$

Then D_t is given by the following (Yu and Li, 2004):

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }{D}_{\mathrm{t}}=\mathrm{1}+{\displaystyle \frac{\ln \mathit{\tau }}{\ln \left(\frac{L_{\mathrm{s}}}{\mathit{\lambda }}\right)}}.\end{array}$$

Figure 1 Two different microscale arrangements of solid particles: (a) RTA and (b) SPA.

RTA and SPA are shown in Fig. 1. An equilateral triangle and a square are selected from the two microstructures as unit cells (Fig. 1a and b). The unit cell of the equilateral triangle is composed of three solid particles with the pore among them, while the unit cell of square is composed of four solid particles. For the unit cell of RTA in Fig. 1a, corresponding porosity is given by the following:

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }n={\displaystyle \frac{A_{\mathrm{a}}-\mathit{\pi }{R}_{\mathrm{v}}^{\mathrm{2}}/\mathrm{2}}{A_{\mathrm{a}}}},\end{array}$$

where n is porosity, A_a is the total area of equilateral triangle, and R_v is the average radius of solid particles. The total area of equilateral triangle can be achieved:

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{A}_{\mathrm{a}}={\displaystyle \frac{\mathit{\pi }{R}_{\mathrm{v}}^{\mathrm{2}}}{\mathrm{2}(\mathrm{1}-n)}}.\end{array}$$

The side length of the equilateral triangle in Fig. 1a can be calculated as follows:

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{L}_{\mathrm{a}}=R_{\mathrm{v}}\sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\sqrt{\mathrm{3}}(\mathrm{1}-n)}}},\end{array}$$

where L_a is the side length.

Figure 2 Three kinds of Sierpinski gasket [30]: (a) L_a=2; (b) L_a=3; and (c) L_a=5.

The area of irregular pore among solid particles is given by the following:

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle }{A}_{\text{ap}}=A_{\mathrm{a}}-{\displaystyle \frac{\mathit{\pi }{R}_{\mathrm{v}}^{\mathrm{2}}}{\mathrm{2}}}={\displaystyle \frac{\mathit{\pi }{R}_{\mathrm{v}}^{\mathrm{2}}n}{\mathrm{2}(\mathrm{1}-n)}},\end{array}$$

where A_ap is the area of pore in the unit cell.

Approximate the pore in the equilateral triangle as a circle, then the maximum diameter of pore can be obtained:

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }{\mathit{\lambda }}_{max,a}=R_{\mathrm{v}}\sqrt{{\displaystyle \frac{\mathrm{2}n}{\mathrm{1}-n}}},\end{array}$$

where ${\mathit{\lambda }}_{max,a}$ is the diameter of the capillary tube in equilateral triangle. The fluid does not only pass the central pore of the unit cell, but also flow through the gap between adjacent particles. The gap length and the average diameter of the capillary tube perpendicular to the plane of equilateral triangle are calculated as follows:

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{L}_{\mathrm{a}}=L_{\mathrm{a}}-\mathrm{2}{R}_{\mathrm{v}}=R_{\mathrm{v}}\left(\sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\sqrt{\mathrm{3}}(\mathrm{1}-n)}}}-\mathrm{2}\right),\text{(21)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{a}}={\displaystyle \frac{{\mathit{\lambda }}_{max,\mathrm{a}}+\mathrm{\Delta }{L}_{\mathrm{a}}}{\mathrm{2}}}={\displaystyle \frac{R_{\mathrm{v}}}{\mathrm{2}}}\left(\sqrt{{\displaystyle \frac{\mathrm{2}n}{\mathrm{1}-n}}}+\sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\sqrt{\mathrm{3}}(\mathrm{1}-n)}}}-\mathrm{2}\right),\end{array}$$

where ΔL_a is the gap length between solid particles; λ_a is the average diameter of capillary tubes in the equilateral triangle.

Generally, the tortuosity of flow path in porous media is the ratio of the length of tortuous flow path to the straight length of flow path along the flow direction (Taiwo et al., 2016):

$$\begin{array}{}\text{(22)}& {\displaystyle }{\displaystyle }\mathit{\tau }={\displaystyle \frac{L_{\mathrm{t}}}{L_{\mathrm{s}}}},\end{array}$$

where L_t is the length of tortuous flow path, and L_s is the straight length of flow path along the flow direction.

For the flow path shown in Fig. 1a, L_t and L_s are respectively as follows:

$$\begin{array}{}\text{(23)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{t}}=(h_o-R_{\mathrm{v}})+{\displaystyle \frac{\mathit{\pi }{R}_{\mathrm{v}}}{\mathrm{2}}}=R_{\mathrm{v}}\left(\sqrt{{\displaystyle \frac{\sqrt{\mathrm{3}}\mathit{\pi }}{\mathrm{2}(\mathrm{1}-n)}}}+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}-\mathrm{1}\right)\text{(24)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{s}}=h_o=R_{\mathrm{v}}\sqrt{{\displaystyle \frac{\sqrt{\mathrm{3}}\mathit{\pi }}{\mathrm{2}(\mathrm{1}-n)}}},\end{array}$$

where h_o is the altitude of the equilateral triangle, $h_o=\frac{\sqrt{\mathrm{3}}}{\mathrm{2}}{L}_{\mathrm{a}}=R_{\mathrm{v}}\sqrt{\frac{\sqrt{\mathrm{3}}\mathit{\pi }}{\mathrm{2}(\mathrm{1}-n)}}$.

Consequently, the tortuosity of RTA is yielded:

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle }\mathit{\tau }={\displaystyle \frac{L_{\mathrm{t}}}{L_{\mathrm{s}}}}=\mathrm{1}+{\displaystyle \frac{\frac{\mathit{\pi }}{\mathrm{2}}-\mathrm{1}}{\sqrt{\frac{\sqrt{\mathrm{3}}\mathit{\pi }}{\mathrm{2}(\mathrm{1}-n)}}}}.\end{array}$$

The D_f is determined using Sierpinski gasket (Fig. 2) in fractal theory (Yu and Cheng, 2002). The shaded area represents solids of porous media and the white area represents pore. The pore area fractal dimensions in Fig. 2a–c are 0.000, 1.000 and 1.594, respectively ($\mathrm{1}=L_{\mathrm{a}}^{D_{\mathrm{f}}}={\mathrm{2}}^{D_{\mathrm{f}}}$, $\mathrm{3}=L_{\mathrm{a}}^{D_{\mathrm{f}}}={\mathrm{3}}^{D_{\mathrm{f}}}$, $\mathrm{13}=L_{\mathrm{a}}^{D_{\mathrm{f}}}={\mathrm{5}}^{D_{\mathrm{f}}}$). Based on the Sierpinski gasket, the dimensionless pore area in RTA (Fig. 1a) is approximated as follows:

$$\begin{array}{}\text{(26)}& {\displaystyle }{\displaystyle }{A}_{\text{apd}}=(L_{\mathrm{a}}^{+}{)}^{D_{\mathrm{f}}},\end{array}$$

where A_apd is the dimensionless pore area of RTA; $L_{\mathrm{a}}^{+}=L_{\mathrm{a}}/{\mathit{\lambda }}_{min}$. Equation (26) can be solved to achieve D_f:

$$\begin{array}{}\text{(27)}& {\displaystyle }{\displaystyle }{D}_{\mathrm{f}}={\displaystyle \frac{\ln A_{\text{apd}}}{\ln L_{\mathrm{a}}^{+}}}.\end{array}$$

The porosity equals to the ratio of the dimensionless pore area of RTA (A_apd) to the dimensionless total area of RTA ($A_{\mathrm{a}}^{+}$):

$$\begin{array}{}\text{(28)}& {\displaystyle }{\displaystyle }n={\displaystyle \frac{A_{\text{apd}}}{A_{\mathrm{a}}^{+}}},\end{array}$$

where $A_{\mathrm{a}}^{+}=\frac{A_{\mathrm{a}}}{\mathit{\pi }{\mathit{\lambda }}_{min}^{\mathrm{2}}/\mathrm{4}}=\frac{\frac{\mathit{\pi }{R}_{\mathrm{v}}^{\mathrm{2}}}{\mathrm{2}(\mathrm{1}-n)}}{\mathit{\pi }\frac{{\mathit{\lambda }}_{min}^{\mathrm{2}}}{\mathrm{4}}}=\frac{\mathrm{2}{R}_{\mathrm{v}}^{\mathrm{2}}}{{\mathit{\lambda }}_{min}^{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{1}-n}=\frac{(d^{+}{)}^{\mathrm{2}}}{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{1}-n}$; $d^{+}=\frac{\mathrm{2}{R}_{\mathrm{v}}}{{\mathit{\lambda }}_{min}}$, $L_{\mathrm{a}}^{+}=\sqrt{A_{\mathrm{a}}^{+}}$.

From Eq. (28), the dimensionless pore area of RTA (A_apd) is given by the following:

$$\begin{array}{}\text{(29)}& {\displaystyle }{\displaystyle }{A}_{\text{apd}}=n\cdot A_{\mathrm{a}}^{+}.\end{array}$$

The dimensionless total area of RTA ($A_{\mathrm{a}}^{+}$) can be written as follows:

$$\begin{array}{}\text{(30)}& {\displaystyle }{\displaystyle }{A}_{\mathrm{a}}^{+}=(L_{\mathrm{a}}^{+}{)}^{\mathrm{2}}.\end{array}$$

Afterward, $L_{\mathrm{a}}^{+}$ is calculated as follows:

$$\begin{array}{}\text{(31)}& {\displaystyle }{\displaystyle }{L}_{\mathrm{a}}^{+}=\sqrt{A_{\mathrm{a}}^{+}}=\sqrt{{\displaystyle \frac{(d^{+}{)}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-n}}}=d^{+}\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-n)}}}.\end{array}$$

Substituting Eqs. (29) and (31) into Eq. (27) will derive D_f of RTA:

$$\begin{array}{ll}{\displaystyle }{D}_{\mathrm{f}}& {\displaystyle }={\displaystyle \frac{\ln A_{\text{apd}}}{\ln L_{\mathrm{a}}^{+}}}={\displaystyle \frac{\ln (n\cdot A_{\mathrm{a}}^{+})}{\ln \left(\sqrt{A_{\mathrm{a}}^{+}}\right)}}=\mathrm{2}+{\displaystyle \frac{\ln \left(n\right)}{\ln \left(\sqrt{A_{\mathrm{a}}^{+}}\right)}}\\ \text{(32)}& {\displaystyle }& {\displaystyle }=\mathrm{2}+{\displaystyle \frac{\ln \left(n\right)}{\ln \left(d^{+}\sqrt{\frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-n)}}\right)}}.\end{array}$$

For the unit cell of square shown in Fig. 1b, the porosity is as follows:

$$\begin{array}{}\text{(33)}& {\displaystyle }{\displaystyle }n={\displaystyle \frac{A_{\mathrm{b}}-\mathit{\pi }{R}_{\mathrm{v}}^{\mathrm{2}}}{A_{\mathrm{b}}}},\end{array}$$

where A_b is the total area of the square. Equation (33) can also be expressed as the area of unit cell:

$$\begin{array}{}\text{(34)}& {\displaystyle }{\displaystyle }{A}_{\mathrm{b}}={\displaystyle \frac{\mathit{\pi }{R}_{\mathrm{v}}^{\mathrm{2}}}{\mathrm{1}-n}}.\end{array}$$

Again, the side length of the square is as follows:

$$\begin{array}{}\text{(35)}& {\displaystyle }{\displaystyle }{L}_{\mathrm{b}}=\sqrt{A_{\mathrm{b}}}=R_{\mathrm{v}}\sqrt{{\displaystyle \frac{\mathit{\pi }}{\mathrm{1}-n}}}.\end{array}$$

Consequently, the area of an irregular pore in the square is given by the following:

$$\begin{array}{}\text{(36)}& {\displaystyle }{\displaystyle }{A}_{\text{bp}}=A_{\mathrm{b}}-\mathit{\pi }{R}_{\mathrm{v}}^{\mathrm{2}}={\displaystyle \frac{n\mathit{\pi }{R}_{\mathrm{v}}^{\mathrm{2}}}{\mathrm{1}-n}},\end{array}$$

where A_bp is the area of pore in the square.

Using the following equation, the pore as a circle and obtain corresponding maximum diameter can be approximated:

$$\begin{array}{}\text{(37)}& {\displaystyle }{\displaystyle }{\mathit{\lambda }}_{max,b}=\mathrm{2}{R}_{\mathrm{v}}\sqrt{{\displaystyle \frac{n}{\mathrm{1}-n}}},\end{array}$$

where ${\mathit{\lambda }}_{max,b}$ is the maximum diameter of the capillary tube perpendicular to the plane of the square. Similarly, fluid flows through the central pore in the square and the gap between adjacent particles. As a result, the gap and average diameter of the capillary tube are expressed as follows:

$$\begin{array}{}\text{(38)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{L}_{\mathrm{b}}=L_{\mathrm{b}}-\mathrm{2}{R}_{\mathrm{v}}=R_{\mathrm{v}}\left(\sqrt{{\displaystyle \frac{\mathit{\pi }}{\mathrm{1}-n}}}-\mathrm{2}\right),\text{(39)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{b}}={\displaystyle \frac{{\mathit{\lambda }}_{max,b}+\mathrm{\Delta }{L}_{\mathrm{b}}}{\mathrm{2}}}={\displaystyle \frac{R_{\mathrm{v}}}{\mathrm{2}}}\left(\mathrm{2}\sqrt{{\displaystyle \frac{n}{\mathrm{1}-n}}}+\sqrt{{\displaystyle \frac{\mathit{\pi }}{\mathrm{1}-n}}}-\mathrm{2}\right),\end{array}$$

where ΔL_b is the gap length between the adjacent two solid particles, and λ_b is the average diameter of the capillary tube.

For the tortuous flow path in Fig. 1b, L_t and L_s are respectively given by the following:

$$\begin{array}{}\text{(40)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{t}}=\mathrm{\Delta }{L}_{\mathrm{b}}+\mathit{\pi }{R}_{\mathrm{v}}=R_{\mathrm{v}}\left(\sqrt{{\displaystyle \frac{\mathit{\pi }}{\mathrm{1}-n}}}-\mathrm{2}+\mathit{\pi }\right),\text{(41)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{s}}=L_{\mathrm{b}}=R_{\mathrm{v}}\sqrt{{\displaystyle \frac{\mathit{\pi }}{\mathrm{1}-n}}}.\end{array}$$

Afterward, the tortuosity of SPA yields the following:

$$\begin{array}{}\text{(42)}& {\displaystyle }{\displaystyle }\mathit{\tau }={\displaystyle \frac{L_{\mathrm{t}}}{L_{\mathrm{s}}}}=\mathrm{1}+{\displaystyle \frac{\mathit{\pi }-\mathrm{2}}{\sqrt{\frac{\mathit{\pi }}{\mathrm{1}-n}}}},\end{array}$$

The procedure of deriving D_f of SPA is similar to the procedure of calculating D_f of RTA. Similarly, the D_f and porosity of SPA (Fig. 1b) are given by the following:

$$\begin{array}{}\text{(43)}& {\displaystyle }& {\displaystyle }{D}_{\mathrm{f}}={\displaystyle \frac{\ln A_{\text{bpd}}}{\ln L_{\mathrm{b}}^{+}}},\text{(44)}& {\displaystyle }& {\displaystyle }n={\displaystyle \frac{A_{\text{bpd}}}{A_{\mathrm{b}}^{+}}},\end{array}$$

where A_bpd is the dimensionless pore area of SPA, $L_{\mathrm{b}}^{+}=L_{\mathrm{b}}/{\mathit{\lambda }}_{min}$, $A_{\mathrm{b}}^{+}$ is the dimensionless total area of SPA, and $A_{\mathrm{b}}^{+}=\frac{A_{\mathrm{b}}}{\mathit{\pi }{\mathit{\lambda }}_{min}^{\mathrm{2}}/\mathrm{4}}=\frac{\frac{\mathit{\pi }{R}_{\mathrm{v}}^{\mathrm{2}}}{\mathrm{1}-n}}{\mathit{\pi }\frac{{\mathit{\lambda }}_{min}^{\mathrm{2}}}{\mathrm{4}}}=\frac{\mathrm{4}{R}_{\mathrm{v}}^{\mathrm{2}}}{{\mathit{\lambda }}_{min}^{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{1}-n}=(d^{+}{)}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{1}-n}$.

The dimensionless pore area of SPA (A_bpd) can be yielded from Eq. (44):

$$\begin{array}{}\text{(45)}& {\displaystyle }{\displaystyle }{A}_{\text{bpd}}=n\cdot A_{\mathrm{b}}^{+}.\end{array}$$

$L_{\mathrm{b}}^{+}$ can be calculated as follows:

$$\begin{array}{}\text{(46)}& {\displaystyle }{\displaystyle }{L}_{\mathrm{b}}^{+}=\sqrt{A_{\mathrm{b}}^{+}}=\sqrt{(d^{+}{)}^{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-n}}}=d^{+}\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-n}}}.\end{array}$$

Substituting Eqs. (45) and (46) into Eq. (43), D_f of SPA can be derived:

$$\begin{array}{ll}{\displaystyle }{D}_{\mathrm{f}}& {\displaystyle }={\displaystyle \frac{\ln A_{\text{bpd}}}{\ln L_{\mathrm{b}}^{+}}}={\displaystyle \frac{\ln \left(n\cdot A_{\mathrm{b}}^{+}\right)}{\ln \left(\sqrt{A_{\mathrm{b}}^{+}}\right)}}\\ \text{(47)}& {\displaystyle }& {\displaystyle }=\mathrm{2}+{\displaystyle \frac{\ln \left(n\right)}{\ln \left(\sqrt{A_{\mathrm{b}}^{+}}\right)}}=\mathrm{2}+{\displaystyle \frac{\ln \left(n\right)}{\ln \left(d^{+}\sqrt{\frac{\mathrm{1}}{\mathrm{1}-n}}\right)}}.\end{array}$$

The entry pressure of a tortuous capillary tube (P_c) is defined by Young–Laplace equation as follows (Ahn and Seferis, 1991):

$$\begin{array}{}\text{(48)}& {\displaystyle }{\displaystyle }{P}_{\mathrm{c}}={\displaystyle \frac{\mathit{\omega }}{\mathit{\lambda }}}{\displaystyle \frac{\mathrm{1}-n}{n}},\end{array}$$

where P_c is the entry pressure, λ is the diameter of the capillary tube, ω equals to Fσcosθ in which θ is the contact angle between fluid and solid, σ is the surface tension of the wetting fluid, and F is the form factor depending on the capillary tube alignment and the flow direction.

2.2 Dealing with the heterogeneity of porous media

In this study, SGS is used to generate random realization of a heterogeneous porosity field. SGS is a stochastic simulation method combining sequential principles and a Gaussian method. It assumes variable fit to a Gaussian random field. The Gaussian distribution function is constructed at each simulated spatial location based on the characteristics of variation function and afterward randomly selects a value as the variable at the location. In the SGS method, observation data are transformed to Gaussian distributions or normal distributions. Based on current sample data, the conditional probability distribution of points to be simulated is calculated by the SGS method and then simulation is performed based on a semivariogram model. Each simulated value, together with measured data and previous simulation data, becomes the conditional data set for the next step. As simulation proceeds, the conditional data set increases. Previous research suggests 50–400 realizations are required to obtain a statistically stable mean realization (Eggleston et al., 1996; Hu et al., 2007).

2.3 Modeling PCE migration and its remediation

The DNAPL migration and remediation are modeled using a multicomponent, multiphase, and multicomposition simulator named UTCHEM (University of Texas Chemical Compositional Simulator) (Delshad et al., 1996). As an extension to Delshad's work, UTCHEM was developed by the University of Texas as a comprehensive and practical tool. In numerous applications, UTCHEM has proved to be particularly useful in simulation of contaminant migrations and has been a popular multiphase-flow, multiconstituent, reactive transport model used widely in groundwater simulations. UTCHEM accounts for chemical, physical and biological reactions; complex non-equilibrium sorption; decay and geochemical reactions; and surfactant-enhanced solubilization and mobilization of DNAPLs. Moreover, heterogeneous properties of porous media are also considered. As a result, UTCHEM has been adapted for a variety of environmental applications such as surfactant-enhanced aquifer remediation (SEAR). In this study, DNAPL migration and remediation for cleaning up DNAPL contamination in idealized heterogeneous sites are simulated by UTCHEM.

Figure 3 (a) Two-dimensional view of contaminated domain and (b) locations of injection and extraction wells.

3 Application to a synthetic heterogeneous PCE contaminated site

3.1 Site description

The idealized domain of synthetic application is a two-dimensional confined aquifer (Fig. 3). The length, width and depth of aquifer are 101, 25 and 25 m, respectively. Idealized aquifer is discretized into 101 grids horizontally and 25 layers vertically (Fig. 3b). The spacing of each grid is uniformly 1 m along x and z axes, and the longitudinal and transverse dispersivities are set as to 1.0 and 0.1 m, respectively. Horizontal and vertical correlation length values are each 5 m. The top and bottom borders of aquifer are defined as no-flow boundaries, while the left and right borders are defined as constant potential boundaries to create a groundwater flow from left to right under a low hydraulic gradient of 0.005 m m⁻¹ (Liu et al., 2003; Liu, 2005; Qin et al., 2007). The porous media of idealized aquifer is assumed to be heterogeneous and mixed by different grades of sands.

Figure 4 (a) The individual porosity field generated by the sequential Gaussian simulation (SGS) method; (b) the frequency of an individual porosity field; (c) the individual permeability field of RTA obtained from an individual porosity field; (d) the frequency of individual permeability field for RTA; (e) the individual permeability field of SPA obtained from an individual porosity field; (f) the frequency of individual permeability field for SPA; (g) the obtained individual entry-pressure field of RTA; (h) the frequency of individual entry-pressure field of RTA; (i) the obtained individual entry-pressure field of SPA; and (j) the frequency of individual entry pressure of SPA.

The porosity of aquifer is assumed to be spatially and uniformly distributed with an average value of 0.220 and SD (standard deviation) of 0.060. In this study, porosity follows normal distribution and its SD represents the enhanced geological heterogeneity. A total of 200 realizations of the porosity field are generated using SGS. One of the 200 realizations of heterogeneous field is shown in Fig. 4a. Simultaneously, statistical assessment is undertaken on the individual realization of the porosity field, and corresponding histograms are shown in Fig. 4b. We find that the frequency of the individual realization of the porosity field is close to normal distribution, which conforms to the situation that most characteristics of natural aquifer can be expressed as a normal distribution (Montgomery et al., 1987). Based on the heterogeneous porosity field, the fractal dimension of tortuosity D_t, the fractal dimension for pore areas D_f and the diameter of capillary tubes in porous media, permeability is obtained by Eq. (12). Figure 4c shows the individual heterogeneous permeability field selected from the 200 realizations of RTA, and the result of the associated frequency analysis is shown in Fig. 4d. The permeability field fits the lognormal distribution well, which has been presented by many studies which show that the parameter of aquifer penetrability follows lognormal distribution (Montgomery et al., 1987; Veneziano and Tabaei, 2004). Compared to the histogram of the porosity field in Fig. 4b, the shape of permeability is similar. The individual heterogeneous permeability field of SPA is shown in Fig. 4e. Corresponding frequency analysis of SPA reveals that the permeability field follows lognormal distribution, while some difference appears compared with RTA (Fig. 4f). The average permeability of individual realization of RTA is 2.012 × 10⁻¹² m², and the average permeability of individual realization of SPA is 1.618 × 10⁻¹² m². For 200 realizations, the average permeability of RTA and SPA are 2.120 × 10⁻¹² and 1.706 × 10⁻¹² m², indicating the permeability of RTA is slightly bigger than SPA.

The average pore diameters of two different microscale arrangements of particles are derived using corresponding fractal models. In detail, the average diameter of RTA is calculated by Eq. (21) and the average diameter of SPA is calculated by Eq. (39). Consequently, the entry pressure of the two kinds of microscale arrangements can be obtained by Eq. (48). The individual entry-pressure fields of two microscale arrangements and associated frequency analysis are shown in Fig. 4g–j. From the frequency of entry pressure in Fig. 4h and j, the entry pressures of both RTA and SPA are the lognormal distributions. However, the average entry pressure of individual realization of RTA is 1.980 kPa, while the average entry pressure of SPA is 1.481 kPa. For 200 realizations of the entry-pressure field, the average entry pressure of RTA is 1.922 kPa and the average entry pressure of SPA is 1.442 kPa. The differences of average entry pressure between RTA and SPA imply the microstructure of aquifer has an effect on the macroscopic characteristics.

The purpose of this study is to explore the effects of microstructure of aquifer on DNAPL migration and remediation. A PCE spill event (the leaking of underground storage tank) occurs on the top of the aquifer and a surfactant remediation is designed to clean up the contaminated aquifer. The total duration of 300 days is divided into four stages: (1) 300 m³ PCE is released from underground storage tank into aquifer at the top layer of spill position shown in Fig. 3a during 0–30 days, (2) PCE migrates in aquifer freely during 30–100 days, (3) surfactant is injected into aquifer during 100–150 days, and (4) water is flushed during 150–300 days. In the first stage, PCE is released as a point pollution source in the center grid block at the top layer of the aquifer, in which spill is at a constant rate of 10 m³ day⁻¹. After PCE coming into the heterogeneous aquifer, PCE is migrating freely under the effects of gravity and the natural hydraulic gradient condition. The PCE not only migrates downward through the aquifer, but can also be trapped by capillary forces as residual ganglia and globules. During the long-term PCE migration period, PCE is contaminating groundwater and expanding plume. To clean up the contaminated aquifer, 4 % surfactant solution is injected into aquifer through the two injection wells (Fig. 3b) at a constant rate of 80 m³ day⁻¹, and, simultaneously, contaminated groundwater is extracted through production well at a constant rate of 160 m³ day⁻¹. Surfactant can reduce the interfacial tension between the DNAPL and aqueous phase to promote solubilization and mobilization of DNAPL. After surfactant injection, the contaminated aquifer is flushed by water over a long period of 150 days. Based on the distributions of porosity, permeability and entry pressure of two microscale arrangements, the entire PCE migration and remediation process is simulated by a multicomponent, multiphase model simulator, UTCHEM (Delshad et al., 1996). The parameters used in the simulation are listed in Table 1. Simulation results of two different microscale arrangements are compared to reveal the effect of microstructure on the DNAPL migration and remediation.

Table 1 Parameters used in the simulation.

Figure 5 Simulated PCE saturation for individual realization of RTA over the entire migration and remediation periods (0 ∼ 300 day).

3.2 Results and discussion

3.2.1 PCE migration and its remediation based on single realizations

The simulation results of PCE migration for individual realization of the porosity field for RTA are shown in Fig. 5a–f. When PCE is released into an aquifer at the top layer of spill position, PCE almost infiltrates vertically under the effect of gravity (Fig. 5a). Due to the heterogeneity of the aquifer, some preferential flow appears and the PCE plume becomes irregular (Fig. 5b). After 30 days, PCE plume almost touches the bottom of aquifer (Fig. 5c). When the PCE leakage is stopped, PCE migrates continuously in aquifer for 70 days (Fig. 5d–f). The released PCE is migrating downward and entrapped by capillary forces as residual ganglia and globules. The heterogeneity of the aquifer makes PCE migrate along a preferential pathway. When the PCE plume touches the zones of low permeability and high entry pressure, it will bypass these zones and migrate continuously, which leads to an increasing variability in PCE distribution. After the PCE plume reaches the bottom of aquifer, PCE begin accumulate and form a contaminant pool at the bottom. At t = 100 days, a PCE pool is formed at the bottom of aquifer, moving toward the right boundary.

Figure 6 Simulated PCE saturation for individual realization of SPA over the entire migration and remediation periods (0 ∼ 300 day).

Figure 6a–f show the simulated PCE saturation for individual realization of porous media for SPA during migration period. Under the effects of gravity and the natural hydraulic gradient, PCE is migrating and the contaminant plume becomes larger and larger. The heterogeneity of the aquifer significantly changes the migration paths and leads to irregular morphology of the PCE plume (Fig. 6a–c). However, due to the different microarrangements of the aquifer, the entry-pressure distribution is also different, which leads to some differences. After the PCE injection, the simulated PCE saturation in Fig. 6d–f indicates that further trapping and spreading of the PCE occurs during this period. Compared with the simulation results of RTA in Fig. 5, the PCE plume slightly seems similar in Fig. 6. Moreover, PCE infiltrates more quickly in porous media of RTA in Fig. 5. After 70 days, the PCE plume has touched the bottom for RTA (Fig. 5e), while the PCE plume based on SPA still keeps a significant distance from the bottom (Fig. 6e).

To clean up the DNAPL, 4 % surfactant solution is injected through two injection wells at a constant rate of 80 m³ day⁻¹ over 50 days. Afterwards, water flushing is applied during 150–300 days. The locations of the injection and production wells are presented in Fig. 3b. The production well is rightly installed at the location of the PCE spill position and two injection wells are located 39 m to the left and right of the production well. Figure 5g–l show the PCE remediation results of individual realization for RTA. During the early remediation period, the effect of cleaning up DNAPL is not yet apparent (Fig. 5g–i). When the water flushing begins, the surfactant solution circulates throughout the contaminated aquifer (Fig. 5j–l). At t = 200 days, 237.01 m³ PCE is removed from contaminated aquifer, occupying 79.00 % of the total released PCE (Fig. 5j). As time goes on, 268.30 m³ PCE is removed from the aquifer and remediation efficiency reaches 89.43 %.

The same surfactant remediation is also conducted for individual realization of SPA. Compared with the remediation for RTA, the remediation effect is more apparent for SPA (Fig. 6g–l). As the remediation progresses, more DNAPL is removed and less DNAPL remains at the bottom of aquifer. At t = 200 days, 267.68 m³ PCE is removed from the contaminated aquifer and the corresponding remediation efficiency rises to 89.23 %. At t = 300 days, 285.32 m³ PCE is cleaned up and remediation efficiency reaches 95.11 %. From the results of remediation, it is obvious that microstructure has an effect on remediation for macroscopic-scale aquifer. Results suggest contaminated aquifer of RTA is hard to clean up by surfactant remediation while SPA can improve DNAPL remediation efficiency.

Figure 7 (a) PCE volume in aquifer vs. time, RTA represents RTA and SPA represents SPA; (b) changes in GTP as a function of time; (c) cumulative DNAPL removal as a function of time; (d) variation of GTP value as a function of cumulative DNAPL removal percent; (e) the change of the center of PCE plume during the entire periods of migration and remediation; (f) the change of the depth of PCE plume center during the entire periods; (g) variation of second PCE plume moment on the horizontal axis; and (h) variation of second PCE plume moment in vertical axis.

3.2.2 PCE migration and SGS realizations

PCE migration and remediation processes are simulated for 200 realizations of the porosity field for porous media of RTA and SPA. The variations of contaminant mass, the ganglia-to-pool (GTP) ratio and moments of PCE plume vs. time are presented in Fig. 7a–h. During 0–30 days, the PCE in aquifer increases linearly at a constant rate of 10 m³ day⁻¹ (Fig. 7a), which corresponds to the contaminant spill stage. Afterward, PCE volume keeps constant during the second stage (30–100 days), while PCE volume in aquifer is reduced when surfactant is injected into aquifer. After surfactant insertion and water flushing of the contaminated aquifer, most DNAPL is cleaned up. The residual DNAPL mass remained in aquifer of 0.67–119.89 m³ with a mean of 22.42 and 0.79–103.33 m³ with a mean of 12.51 m³ are achieved for 200 heterogeneous realizations based on the RTA and SPA, respectively. The average remediation efficiency of SPA is undoubtedly higher than RTA, indicating the aquifer of SPA is easier to clean up. PCE plume architectures are quantified by measuring the GTP ratio in Fig. 7b. Over entire periods, curves of GTP value show obvious oscillations. Surfactant has the ability to promote solubilization, and mobilization of DNAPL can reduce GTP value. As a result, when surfactant is injected at t = 100 day, the GTP value reduces quickly. When surfactant injection is ended and water flushing begins, the GTP value increases with steep flank slope. Finally, GTP values reach 0.10–0.41 with a mean of 0.21 and 0.15–0.42 with a mean of 0.28 for 200 heterogeneous realizations based on the RTA and SPA, respectively.

Figure 7c shows cumulative PCE removal from the contaminated aquifer vs. flushing time for RTA and SPA. During the surfactant injection period, 100–150 days, the DNAPL removal is not apparent. However, DNAPL is removed effectively and quickly during the water-flushing period. Through long-term remediation, the removal of PCE from the contaminated aquifer reached 179.89–298.98 m³ with a mean of 277.29 and 196.45–298.87 m³ with a mean of 287.21 m³ for 200 realizations based on RTA and SPA, respectively. Average remediation efficiency of SPA (95.83 %) is noticeably higher than average remediation efficiency of RTA (92.52 %).

Figure 7d shows the GTP value as a function of cumulative PCE removal for the contaminated aquifer. The GTP remains at a relatively low level before 30 % of the DNAPL is removed from the aquifer. When 40 % of the total 300 m³ PCE is removed, GTP values are increasing and corresponding curves appear as a wave crest because the high-saturation zones of the PCE plume are dissolved and turned into ganglia. After the wave crest, the GTP values decline quickly with steep flank slope due to PCE ganglia removal through water flushing. Finally, GTP values increase at the end of the remediation process for 200 realizations, indicating that most of PCE is removed and most of residual PCE turns into ganglia.

For the center of the PCE plume on the horizontal axis, associated variations vs. time are similar for 200 realizations based on RTA and SPA (Fig. 7e). Significantly, the PCE plume vertical-infiltration rate in aquifer of RTA is slightly faster than PCE infiltration in the aquifer of SPA for 200 realizations (Fig. 7f). Simultaneously, the second PCE plume moments in the horizontal direction of RTA are different from the second PCE plume moments in the horizontal direction of SPA (Fig. 7g). After PCE migration under natural conditions at t = 100 days, the second PCE plume moments in the horizontal direction are 10.61–40.50 m² with a mean of 21.51 and 10.99–36.38 m² with a mean of 20.75 m² for 200 realizations based on RTA and SPA, respectively. At t = 300 day, the second PCE plume moments in the horizontal direction change to 0.81–34.88 m² with a mean of 5.79 and 1.03–24.57 m² with a mean of 4.64 m² for RTA and SPA, respectively. The horizontal second moment of RTA is always larger than the horizontal second moment of SPA, indicating that the PCE plume in the aquifer of RTA is wider than the PCE plume in the aquifer of SPA, and RTA can cause more groundwater contamination. Similarly, the second moments in vertical direction for RTA are larger than the second moments in vertical direction for SPA.

This study takes an important step toward exploring how microscale arrangements control contaminant migration on a small-aquifer scale. Results are essential to the macroscopic aquifer composed of porous media without large heterogeneity, such as sandy aquifers containing rich groundwater resources. However, there are many problems associated with the upscaling of aquifers in real-world conditions (Dagan et al., 2013; Pacheco, 2013; Pacheco et al., 2015). Due to the large heterogeneity of natural aquifers, research results may be very different and can not be extrapolated to complex regional aquifer on a large scale. However, the finding in this study is absolutely applicable for natural aquifers with similar heterogeneities. If the heterogeneity and anisotropy of natural aquifers are very different, the effect of the microscale arrangements on the macroscopic contaminant migration and remediation will be different. Although real-world conditions are complex, the new findings achieved from this research are very significant for understanding the effect of microscale arrangement on contaminant behaviors on an aquifer scale. The upscaling problem of the results obtained on the simulation scale (100 m × 25 m × 25 m) is the basis and the upscaling problem with more complex heterogeneity conditions is needed to be further investigated. Various research on the upscaling problem is done through experiment and simulation (Wu et al., 2017a–d). Based on this research, the microstructure of porous media is developed and the contaminates migration in porous media is explored using fractal methods in this study, implying the experimental results are very significant for real-world problems on an aquifer scale. Our next procedure involves applying these models in a real-world aquifer with complex heterogeneity conditions and modifying our models and method according to realistic conditions.

4 Conclusions

The microstructure of aquifers has an important effect on macroscopic-scale characteristics of contaminant migration and remediation. In this study, we focus on the DNAPL migration and remediation in heterogeneous aquifers composed of granular porous media with RTA and SPA. The microscale models of RTA and SPA are developed to obtain the mathematical expressions of permeability and entry pressure using fractal methods. A total of 200 realizations of the porosity field are generated using the SGS method, and PCE is released from an underground storage tank into a heterogeneous aquifer. To clean up contamination caused by the underground storage tank spill, a surfactant remediation technique is used to remove contaminants in the aquifer. The entire process of DNAPL migration and remediation is simulated by a multicomponent, multiphase model simulator, UTCHEM. Results suggest RTA not only cause more groundwater contamination than RTA, but also the contaminated aquifer of RTA is harder to clean up compared with SPA. The second PCE plume moments in the horizontal direction are 10.61–40.50 m² with a mean of 21.51 and 10.98–36.38 m² with a mean of 20.75 m² for 200 realizations based on RTA and SPA, respectively, after long-term migration at t = 100 days. Furthermore, the second PCE plume moments in the horizontal direction at t = 300 day are 0.807–34.88 m² with a mean of 5.79 and 1.025–24.57 m² with a mean of 4.64 m² for RTA and SPA, respectively, after long-term remediation. Simultaneously, the residual DNAPL mass remaining in the aquifer is 0.67–119.89 m³ with a mean of 22.42 for RTA and 0.79–103.33 m³ with a mean of 12.51 m³ for SPA, indicating that the remediation efficiency of SPA (65.53–99.74 % with a mean of 95.83 %) is mostly higher than the remediation efficiency of RTA (60.01–99.78 % with a mean of 92.52 %). This study reveals that the microstructure of an aquifer has an important effect on contaminant movement and associated remediation efficiency on a macroscopic scale, which is very essential and significant for dealing with accidental underground storage tank spills and for identifying subsurface contaminant sources in the future.