Contaminant transport model in transient and unsaturated conditions applied to laboratory column test with tailings

Abstract

Mining is an important economic activity in the modern world. However, despite the generated benefits, mining produces tremendous volumes of tailings, an environmental liability with numerous adverse effects. Researches about contaminant transport in tailings dam are important to assess the degree of contamination and to propose preventive or remedial measures. In geotechnical practice, the flow of solutes is generally characterized by numerical solution of the Richards equation to describe water movement followed by advection-dispersion equation to describe contaminant movement. This study aimed to model and simulate contaminant transport in a laboratory column test, using a new analytical formulation and mathematical codes, through tailings in transient unsaturated conditions. The analytical solution for the Richards equation was used to simulate the variation in the volumetric water content and to determine the transient contaminant plume using the advection-dispersion equation subsequently. The models were used to calibrate experimental data from hydraulic characterization and contamination tests. Finally, the normalized contaminant plume (c_w/c₀) was simulated as a function of time and space. Comparisons with experimental data showed that the analytical formulations adequately expressed the process of contaminant infiltration through the unsaturated porous medium. The formulations offered effectively and are configured as a new approach to solve various contamination problems in transient unsaturated conditions, providing insights into many complex processes that occur in the lab tests and requires far less computational effort compared with current programs to modeling the solute transport using numerical solutions, as the versatile commercial Software HYDRUS.

Keywords
Unsaturated porous medium; Unsaturated flow; Contaminant transport; Tailing

1. Introduction

A wide range of everyday human activities, such as tailings dams, release different types of contaminants to the soil surface. These contaminants pollute the subsurface through uncontrolled spills, leaks, dumps and disposal (Rutsch et al., 2008; Seferou et al., 2013; Adnan et al., 2018; Akbariyeh et al., 2018). Contaminant transport may lead to severe problems with the quality of the groundwater and the porous medium. A deeper understanding of hydrodynamic and hydrogeochemical processes of the unsaturated zone is essential to assess the vulnerability of an aquifer to contamination. In this sense, many studies related to water flow and solute movement in the subsurface have been made to a wide range of conditions involving different scales, types of soils, contaminants and models of resolution (Bertolo, 2001; Rutsch et al., 2008; Mustafa et al., 2016; Sopilniak et al., 2017; Joshi & Gupta, 2018; Akbariyeh et al., 2018; Godoy et al., 2019).

In geotechnical practice, the flow of solutes is generally characterized by numerical solution of the Richards equation to describe water movement followed by advection-dispersion equation to describe contaminant movement (Richards, 1931). When applied, especially, in the field scale, this numerical solution is highly complex, computationally expensive. It inhibits insights that can be obtained from analytical solutions. One of the most versatile commercial software is HYDRUS.

HYDRUS is a Microsoft Windows based modeling environment for simulating water flow, heat, and solute transport in two- and three-dimensional variably saturated and unsaturated media. The program numerically solves the Richards equation for saturated-unsaturated water flow and advection-dispersion equations for heat and solute transport (Šimůnek et al., 2016).

This study aims to model and simulate the contaminant transport in a laboratory column test, using a new analytical formulation and mathematical codes, through tailings in transient unsaturated conditions, providing insights into many complex processes that occur in the lab tests and requires far less computational effort compared with current programs to modeling the solute transport using numerical solutions, as the Software HYDRUS.

This paper innovates by developing a computational code that allows coupling the analytical solution of the Richards Equation proposed by Cavalcante & Zornberg (2017) to solve the infiltration problem taking into account the phenomena of transport of contaminants to unsaturated porous media, solved from the analytical solution proposed by Brenner (1962), described in terms of volumetric water content.

2. Contaminant transport modeling in an unsaturated porous medium

The study of contaminant transport in unsaturated porous media requires understanding the physical and biophysicochemical phenomena associated with water and solute flow in the soil.

In nature, a porous medium may exist in saturated or unsaturated conditions, where the underground water table delimits each zone. Below the water level, the soil is normally saturated. In this zone, the voids in the soil are entirely filled with water, and the pressure is positive. Above the water level, the soil may be unsaturated. In this zone, the pore pressure is negative and determined by the difference between the air and water pressure in the voids. The negative pressure is known as total suction, which consists of the sum of the matric and osmotic suctions.

The unsaturated zone plays a crucial role in contamination processes. This zone is essentially a natural filter, which reduces or attenuates the contamination of microbiological, physical and chemical constituents to nonhazardous levels by acting as a channel through which various liquid or vapor compounds circulate, attenuate and transform as they move from the surface to the unsaturated zone. Furthermore, contaminants can remain in this zone for decades, affecting the plants and animals and contaminating aquifers long after a spill.

When contaminants begin their trajectory on the surface of the terrain and are dragged by the waters that infiltrate the unsaturated zone, they are subjected to a series of transport mechanisms, which define how a specific compound will move in the porous medium. The movement of these compounds depends not only on the flow of the fluid in which these substances are dissolved but also on the soil-contaminant interactions related to the chemical and biological processes to which these substances are subjected.

The analysis of contaminant transport from saturated to unsaturated porous media is relatively simple. For this purpose, reformulating the transport equations in terms of the volumetric water content of the porous medium suffices because varying the required parameters in an unsaturated porous medium reduces the conductive cross-section through which a contaminant can flow (Fityus et al., 1999).

2.1 Flow in unsaturated porous media

Water flow in unsaturated porous media is governed by the Richards equation, which is a mass conservation equation that combines the Darcy-Buckingham law and the continuity equation.

The continuity equation is based on the principle of mass conservation, i.e., on the fact that the difference between the masses of fluid entering and leaving an infinitesimal control volume equals the rate of change in mass storage in that volume. Considering that water, at the tension levels of this study, is incompressible, the following formulation is derived:

$\frac{\partial {\theta }_w}{\partial t}=-\frac{\partial v_z}{\partial z}$ (1)

where θ_w = θ_w(z,t) is the volumetric water content [L³.L^-3], and v_z = v_z(z,t) is the fluid velocity [L.T^-1] in the z-direction. The fluid velocity can be defined using the Darcy-Buckingham law (Buckingham, 1907; Narasimhan, 2004), which is the unsaturated version of Darcy's law, as follows:

$v_z=-k_z\left(\psi \right)\frac{\partial \varphi }{\partial z}$ (2)

where k_z(ψ) is the hydraulic conductivity expressed in terms of suction [L.T^-1] in the z-direction; ψ = ψ(z,t) is the total suction of water [M.L^-1.T^-2]; and ϕ= ϕ(z,t) is the hydraulic head of the fluid [L], which is defined as follows:

$\varphi =-z+\frac{\psi }{{\rho }_{w}g}$ (3)

where z is the elevation head [L]; ρ_w is the density of water [M.L^-3]; and g is the gravitational acceleration [L.T^-2]. By combining Equations 2 and 3, the following expression is derived:

$v_z=-k_z\left(\psi \right)\left(\frac{1}{{\rho }_{w}g}\frac{\partial \psi }{\partial z}-1\right)$ (4)

Now, by combining Equations 1 and 4, the Richards equation for one-dimensional, unsaturated, transient flow in the z-direction is obtained as follows:

$\frac{\partial {\theta }_w}{\partial t}=\frac{\partial }{\partial z}\left[k_z\left(\psi \right)\left(\frac{1}{{\rho }_{w}g}\frac{\partial \psi }{\partial z}-1\right)\right]$ (5)

Cavalcante & Zornberg (2017) rewrote the Richards equation as a modified version of the Fokker-Planck equation to enable an analytical solution to the transient flow problem in unsaturated conditions. The modified version proposed by Cavalcante & Zornberg (2017) is as follows:

$\frac{\partial {\theta }_w}{\partial t}=\frac{\partial }{\partial z}\left(D_z\left({\theta }_w\right)\frac{\partial {\theta }_w}{\partial z}\right)-a_s\left({\theta }_w\right)\frac{\partial {\theta }_w}{\partial z}$ (6)

where D_z is the unsaturated diffusivity of water in the z-direction [L².T^-1], which is expressed as follows:

$D_z\left({\theta }_w\right)=\frac{k_z\left({\theta }_w\right)}{{\rho }_{w}g}\frac{\partial \psi }{\partial {\theta }_w}$ (7)

and a_s is the unsaturated advective flow [L.T^-1], as given by the following equation:

$a_s\left({\theta }_w\right)=\frac{\partial k_z\left({\theta }_w\right)}{\partial {\theta }_w}$ (8)

Using this representation of the Richards equation and the specific hydraulic functions (constitutive models), Cavalcante & Zornberg (2017) analytically solved the unsaturated flow problem for transient flow conditions. The adopted hydraulic functions involve a logarithmic relationship between the suction and the volumetric water content (soil water retention curve) [ψ(θ_w)] and a linear relationship between the unsaturated hydraulic conductivity and the volumetric water content [k_z(θ_w)]. Therefore, Cavalcante & Zornberg (2017) proposed the following equation to determine ψ(θ_w):

$\psi \left({\theta }_w\right)=\frac{1}{\delta }\ln \left(\frac{{\theta }_w-{\theta }_r}{{\theta }_s-{\theta }_r}\right)$ (9)

where θ_s is the saturated volumetric water content [L³.L^-3]; θ_r is the residual volumetric water content [L³.L^-3]; and δ is a fitting parameter [M^-1.L¹.T²]. In contrast, Cavalcante & Zornberg (2017) proposed the following equation to determine k_z(θ_w):

$k_z\left({\theta }_w\right)=k_s\left(\frac{{\theta }_w-{\theta }_r}{{\theta }_s-{\theta }_r}\right)$ (10)

where k_s is the saturated hydraulic conductivity [L.T^-1].

Cavalcante & Zornberg (2017) presented analytical solutions for the Richards equation, with various initial and boundary conditions, using the characteristic curve and the unsaturated hydraulic conductivity function defined by Equation 9 and Equation 10. The constitutive model by Cavalcante & Zornberg (2017) is the most convenient when compared with other constitutive models because only one fitting parameter is needed; also, this model makes it possible to linearize the Richards equation and, therefore, to solve the problem analytically.

This study used a solution for a soil column of finite length, with a constant infiltration rate imposed on the upper boundary. Specifically, in this case, the initial condition is described by the uniform initial water content as follows:

${\theta }_w\left(z\mathrm{,0}\right)={\theta }_i$ (11)

where θ_i is the initial volumetric water content [L³.L^-3].

A Neumann boundary condition is adopted, which involves a constant infiltration rate imposed on the upper boundary of the domain as follows:

${\left.\left({\overline{D}}_z\frac{\partial {\theta }_w}{\partial z}-k_z\right)\right|}_{z=0}=v_0$ (12)

where v₀ is the infiltration rate [L.T^-1]. The maximum infiltration rate that can be physically imposed corresponds approximately to the saturated hydraulic conductivity of the porous medium (k_s). Specifically, the maximum infiltration rate that can be imposed is as follows:

$v_{\mathrm{0,}\max }=\frac{{\theta }_s{k}_s}{\left({\theta }_s-{\theta }_r\right)}$ (13)

For a soil column of finite length L, the lower boundary condition adopted in this study is described as follows:

$\frac{\partial {\theta }_w}{\partial z}\left(L,t\right)=0$ (14)

This lower boundary condition implies that the water content and, therefore, the suction reach constant values at a given depth. Thus, at a given depth, the hydraulic gradient in the z-direction is equal to one.

The analytical solution that corresponds to this initial condition and the boundary conditions is as follows (Cavalcante & Zornberg, 2017):

${\theta }_w\left(z,t\right)={\theta }_i+\left(\frac{v_0}{k_s}\left({\theta }_s-{\theta }_r\right)-{\theta }_i\right)D\left(z,t\right)$ (15)

where:

$D\left(z,t\right)=1-{\displaystyle \sum _{m=1}^{\infty }\frac{\frac{2{\overline{a}}_{s}L}{{\overline{D}}_z}{\beta }_m\left[{\beta }_m\cos \left(\frac{{\beta }_{m}z}{L}\right)+\frac{{\overline{a}}_{s}L}{2{\overline{D}}_z}\sin \left(\frac{{\beta }_{m}z}{L}\right)\right]\exp \left(\frac{{\overline{a}}_{s}z}{2{\overline{D}}_z}-\frac{{\overline{a}}_s^{2}t}{4{\overline{D}}_z}-\frac{{\beta }_m^2{\overline{D}}_{z}t}{L^2}\right)}{\left[{\beta }_m^2+\frac{{\overline{a}}_{s}L}{{\overline{D}}_z}+{\left(\frac{{\overline{a}}_{s}L}{2{\overline{D}}_z}\right)}^2\right]\left[{\beta }_m^2+{\left(\frac{{\overline{a}}_{s}L}{2{\overline{D}}_z}\right)}^2\right]}}$ (16)

and β_m, are the eigenvalues that correspond to the positive roots of the following equation:

${\beta }_m\cot \left({\beta }_m\right)-\frac{{\beta }_m^2{\overline{D}}_z}{{\overline{a}}_{s}L}+\frac{{\overline{a}}_{s}L}{4{\overline{D}}_z}=0$ (17)

The results can be calculated accurately by considering only the first four terms of the series described by Equation 17. In this case, Equation 17 can be approximated as follows:

$$\begin{array}{l}D\left(z,t\right)=\frac{1}{2}\text{erfc}\left({{\rm Z}}_{-1}\right)+\sqrt{\frac{{\overline{a}}_s^{2}t}{\pi {\overline{D}}_z}}\exp \left(-\frac{{\left(z-{\overline{a}}_{s}t\right)}^2}{4{\overline{D}}_{z}t}\right)-\\ \text{}\text{}-\frac{1}{2}\left(-1+\frac{{\overline{a}}_{s}z}{{\overline{D}}_z}+\frac{{\overline{a}}_s^{2}t}{{\overline{D}}_z}\right)\exp \left(\frac{{\overline{a}}_{s}z}{{\overline{D}}_z}\right)\text{erfc}\left({{\rm Z}}_{+1}\right)+\\ +\sqrt{4\frac{{\overline{a}}_s^{2}t}{\pi {\overline{D}}_z}}\left[1+\frac{{\overline{a}}_s}{4{\overline{D}}_z}\left(2L-z+{\overline{a}}_{s}t\right)\right]\\ \exp \left(\frac{{\overline{a}}_{s}L}{{\overline{D}}_z}-\frac{1}{4{\overline{D}}_{z}t}{\left(2L-z+{\overline{a}}_{s}t\right)}^2\right)-\\ -\frac{{\overline{a}}_s}{{\overline{D}}_z}\left[2L-z+\frac{3{\overline{a}}_{s}t}{2}+\frac{{\overline{a}}_s}{4{\overline{D}}_z}{\left(2L-z+{\overline{a}}_{s}t\right)}^2\right]\\ \exp \left(-\frac{{\overline{a}}_{s}L}{{\overline{D}}_z}\right)erfc\left(\frac{2L-z+{\overline{a}}_{s}t}{2\sqrt{{\overline{D}}_{z}t}}\right)\end{array}$$ (18)

$Z_{\pm 1}=\frac{z\pm {\overline{a}}_{s}t}{2\sqrt{{\overline{D}}_{z}t}}$ (19)

where erfc(Z) is the complementary error function, which is defined as follows:

$erfc\left(Z\right)=1-\frac{2}{\pi }{\displaystyle \underset{0}{\overset{Z}{\int }}ex}p\left(-t^2\right)dt$ (20)

2.2 Contaminant transport

The movement of a solute through a saturated porous medium primarily results from two simultaneous mechanisms: advection and hydrodynamic dispersion. The hydrodynamic dispersion corresponds to the combined action of the mechanical dispersion and molecular diffusion mechanisms. Besides that, there are interactions between the porous medium and the solute, defined as sorption, limiting the contaminant transport. All these mechanisms are described as functions of the volumetric water content to analyze the effect of this variation on the contaminant plume.

Solutes are transported in a porous medium by the flow of the percolating fluid without changing their concentration in the solution, which is a process known as advection. The mass flow by advection is as follows:

$J_a={\theta }_w{v}_p{c}_w$ (21)

where J_a = J_a(z,t) is the mass flow by advection of the contaminant per unit area and per unit time [M.L^-2.T^-1]; c_w = c_w(z,t) is the contaminant concentration [M.L^-3]; θ_w = θ_w(z,t) is the volumetric water content [L³.L^-3] obtained by Cavalcante & Zornberg (2017); and v_p = v_p (θ_w) is the percolation rate [L.T^-1]. The percolation rate is the effective flow velocity that corresponds only to the pores through which the solute actually percolates and is defined as follows (Hillel, 2003; Bear & Cheng, 2010):

$v_p=\frac{k_{z}i}{{\theta }_w}$ (22)

where i is the hydraulic gradient [L.L^-1].

As contaminant molecules are transported, they tend to deviate from the main trajectory, some at a higher velocity than others, thus causing dilution of the solution. The components resulting from the dilution are the result of mechanical dispersion; i.e., they are a consequence of velocity heterogeneity inside the system and molecular diffusion due to different concentration gradients (Freeze & Cherry, 1979). Therefore, the sum of those terms can be expressed by the hydrodynamic dispersion coefficient as follows:

$D_h=D_m+D^{*}$ (23)

where D_h = D_h(z,t) is the hydrodynamic dispersion coefficient [L².T^-1]; D_m = D_m(z,t) is the mechanical dispersion coefficient [L².T^-1]; and D *= D*(z,t) is the molecular diffusion coefficient [L².T^-1].

Mechanical dispersion is predominantly a function of both the heterogeneities in hydraulic conductivity and porosity (intrinsic properties of the porous media) and of the fluid flow. The mechanical dispersion coefficient, D_m, is given by:

$D_m={\alpha }_d{v}_p$ (24)

where α_d is the longitudinal dispersivity coefficient [L].

The mass flow by mechanical dispersion in unsaturated porous media is described by Fick's first law as follows:

$J_m=-{\theta }_w{D}_m\frac{\partial c_w}{\partial z}$ (25)

where J_m = J_m(z,t) is the mass flow by mechanical dispersion of the contaminant per unit area per unit time [M.L^-2.T^-1].

The molecular diffusion mechanism in porous media will differ from that in free water because of the effects of porosity and tortuosity. It is caused by the random movement of contaminant molecules in the liquid phase of the porous medium. The rate of movement is determined by the concentration gradient in the medium. The dissolved solute moves from an area of higher concentration to an area of lower concentration. This phenomenon occurs independently of the flow velocity, and the molecular diffusion coefficient, D*, is defined as follows:

$D^{*}=\tau D_0$ (26)

where D₀ is the molecular diffusion coefficient in the aqueous solution [L².T^-1] and τ = τ(z,t) is the tortuosity factor [L.L^-1], which is defined as follows (Viotti et al., 2005):

$\tau ={\frac{\left({\theta }_s-\theta w\right)}{{\theta }_s^2}}^{10/3}$ (27)

The mass flow by molecular diffusion in unsaturated porous media is described by Fick's first law as follows:

$J^{*}=-{\theta }_w{D}^{*}\frac{\partial c_w}{\partial z}$ (28)

where J* = J*(z,t) is the mass flow by molecular diffusion of the contaminant per unit area and per unit time [M.L^-2.T^-1].

During contaminant transport through unsaturated porous media, interactions occur between the contaminants and the soil particles that delay or accelerate the contamination process. Bear & Cheng (2010) state that most processes involve contaminant mass transfer from the liquid phase to the solid phase (sorption). In some cases, the reverse (desorption) process occurs. Sorption is geochemically quantified using the retardation factor, R (Fetter, 1999).

According to Van Genuchten & Dalton (1986) and Bear & Cheng (2010), the retardation phenomenon in an unsaturated porous medium can be formulated using a linear sorption isotherm, assuming that the percolation rate, v_p, is constant in space and time whenever solute concentration levels are low. In this case, the retardation factor in an unsaturated porous medium is calculated as follows:

$R=1+\frac{{\rho }_d{K}_d}{{\theta }_w}$ (29)

where R = R(z,t) is the retardation factor [dimensionless], ρ_d is the dry density of the soil [M. L^-3], and K_d is the equilibrium distribution coefficient [L³.M^-1].

The equilibrium distribution coefficient measures the amount of chemical substance adsorbed onto soil per amount of water, and it is obtained using the batch equilibrium adsorption test. The linear sorption isotherm is appropriate for cases in which the sorption potential proportionally increases with the concentration and is defined as follows:

$c_s=K_d{c}_w$ (30)

where c_s = c_s(z,t) is the adsorbed contaminant concentration [M.M^-1].

It is important to highlight that most of the contaminants of concern coming from tailings would be metals speciated as cations and cations complexes (e.g., the nickel and cobalt), and their interaction with the porous media is complicated and, under different geochemical conditions, unable to be expressed with reversible and linear sorption represented by an equilibrium distribution coefficient. However, sometimes, this most straightforward form of adsorption linear isotherm is very appropriate at low contaminant concentrations and, for this reason, is used in many studies (Srinivasan & Mercer, 1988; Barone et al. 1992; Mallants et al., 2011).

Column and batch techniques are either commonly used in the laboratory to determine equilibrium distribution coefficients. It should be noted that aspects of transport kinetics, adsorption kinetics, and porous media surface area are not replicated well in batch tests because the results are determined in the equilibrium sorption behavior of a specific compound-sediment combination. Batch tests, despite their limitations, are used in geotechnical practice due to relatively fast and straightforward. However, column experiments always remain limited in their transferability to real-world conditions because of experimental restrictions (such as those imposed by limitations of scale), which means that processes that might co-occur in nature cannot be fully reproduced in a laboratory. Moreover, chemical processes cannot be separated completely from physical processes in column tests (Porro et al., 2000; Banzhaf & Hebig, 2016).

It is important to highlight that when the sorption isotherms are not linear, it is possible to manipulate mathematically and linearize them. Despite the inherent bias of this methodology, as expressed by Foo & Hameed (2010), linearization remains a confident option in the literature, applied in over 95% of the liquid-phase adsorption systems. Therefore, the next challenge in the adsorption field is the identification and clarification of both isotherm models in various adsorption systems.

Models of contaminant transport in physical equilibrium (the advection-dispersion equation – ADE) are based on the classic description of uniform solute flow and transport. The matrix of the porous medium consists of solid particles or impermeable aggregates separated by pores or fractures, wherein solute flow and transport occur (Šimůnek & Van Genuchten, 2008). This formulation is based on the mass conservation equation (or continuity equation) for only one direction. The z-axis of a representative elementary volume in the ADE is as follows:

$\frac{\partial c_w}{\partial t}=-\frac{v_p}{R}\frac{\partial c_w}{\partial z}+\frac{D_h}{R}\frac{{\partial }^2{c}_w}{\partial z^2}$ (31)

where z is the spatial coordinate [L], and t is the time [T].

Brenner (1962) analytically solved the ADE for the following boundary conditions:

Initial condition

$c_w\left(z\mathrm{,0}\right)=0$ (32)

Upper boundary condition

${\left.\left(v_p{c}_w-D_h\frac{\partial c_w}{\partial z}\right)\right|}_{z=0}=\left\{\begin{array}{l}{v}_p{c}_{0}0<t<t_0\\ 0t>t_0\end{array}\right.$ (33)

where c₀ is the initial concentration of the applied solute and is constant in the upper boundary [M.L^-3], and t₀ is the time of application of the displacing solution.

Lower boundary condition

$\frac{\partial c_w}{\partial z}\left(L,t\right)=0$ (34)

where L is the study sample length [L].

Under these hypotheses, the analytical solution by Brenner (1962) is as follows:

$c_w\left(z,t\right)=\left\{\begin{array}{l}A\left(z,t\right)0<t<t_0\\ A\left(z,t\right)-A\left(z,t-t_0\right)t>t_0\end{array}\right.$ (35)

where:

$$\begin{array}{l}A\left(z,t\right)=\frac{1}{2}erfc\left(\frac{Rz-v_{p}t}{2\sqrt{D_{h}Rt}}\right)+\sqrt{\frac{v_p{}^{2}t}{\pi D_{h}R}}\exp \left[-\frac{{\left(Rz-v_{p}t\right)}^2}{4D_{h}Rt}\right]-\frac{1}{2}\left(1+\frac{v_{p}z}{D_h}+\frac{v_p{}^{2}t}{D_{h}R}\right)\exp \left(\frac{v_{p}z}{D_h}\right)\\ .erfc\left(\frac{Rz+v_{p}t}{2\sqrt{D_{h}Rt}}\right)+2\sqrt{\frac{v_p{}^{2}t}{\pi D_{h}R}}\left[1+\frac{v_p}{4D_h}\left(2L-z+\frac{v_{p}t}{R}\right)\right]\exp \left[\frac{v_{p}L}{D_h}-\frac{R}{4D_{h}t}{\left(2L-z+\frac{v_{p}t}{R}\right)}^2\right]\\ -\frac{v_p}{D_h}\left[2L-z+\frac{3v_{p}t}{2R}+\frac{v_p}{4D_h}{\left(2L-z+\frac{v_{p}t}{R}\right)}^2\right]\exp \left(\frac{v_{p}L}{D_h}\right)erfc\left[\frac{R\left(2L-z\right)+v_{p}t}{2\sqrt{D_{h}Rt}}\right]\end{array}$$ (36)

Another way of presenting the analytical solution by Brenner (1962) is in its dimensionless form as follows:

$C\left(Z,T\right)=\left\{\begin{array}{l}B\left(Z,T\right)0<T<T_0\\ B\left(Z,T\right)-B\left(Z,T-T_0\right)T>T_0\end{array}\right.$ (37)

where:

$$\begin{array}{l}B\left(Z,T\right)=\frac{1}{2}erfc\left(\frac{RZ-T}{2\sqrt{\frac{RT}{P}}}\right)+\sqrt{\frac{PT}{\pi R}}\exp \left[-\frac{P\text{}\text{}\text{}{\left(RZ-T\right)}^2}{4RT}\right]-\frac{1}{2}\left[1+P\left(Z+\frac{T}{R}\right)\right].\exp \left(PZ\right)\\ erfc\left(\frac{RZ+T}{2\sqrt{\frac{RT}{P}}}\right)+2\sqrt{\frac{PT}{\pi R}}\left[1+\frac{P}{4}\left(2-Z+\frac{T}{R}\right)\right]\exp \left[P-\frac{PR}{4T}{\left(2-Z+\frac{T}{R}\right)}^2\right]\\ -P\left[2-Z+\frac{3T}{2R}+\frac{P}{4}{\left(2-Z+\frac{T}{R}\right)}^2\right]\exp \left(P\right)erfc\left[\frac{R\left(2-Z\right)+T}{2\sqrt{\frac{RT}{P}}}\right]\end{array}$$ (38)

where C = C(Z,T) is the normalized concentration (C = c_w /c₀); T is the normalized time (T = v_pt / L); Z is the normalized spatial coordinate (Z = z / L); and P is the Peclet number (P = v_p L / D_h).

After applying the consolidation theory induced by micro-collapses, it is necessary for its validation by comparing values obtained in the field. Thus, the consistency of the method was through the Guelph Permeameter. The depth and place of the test were the same that the undisturbed samples. In the end, some comparisons of the alternative and classical theories were executed.

3. Materials and methods

3.1 Materials

Data from laboratory tests conducted by Rodríguez-Pacheco (2002) were used in this study, e.g., the soil water retention curve, the sorption isotherm, and the breakthrough curve of deformed samples of nickel tailings. The samples were collected from the Cuban nickel industry located in the northeast of the province of Holguín in Cuba, where ore is extracted from lateritic nickel ore deposits in the municipalities of Mayari and Moa. Ore is extracted from those deposits using the open-pit mining method (Rodríguez-Pacheco, 2002). Nickel and cobalt concentrates are extracted using the metallurgical process of ammoniacal carbonate leaching (ACL) and of sulfuric acid leaching (SAL). In this paper, the transport of contaminants through the tailings generated by the SAL process was studied. These beneficiation processes generate large volumes of tailings, which are mixed, diluted in water, and transported in pipes as a viscous liquid (pulp) to tailings dams (Rodríguez-Pacheco, 2002; Sosa, 2016).

The experiments and details on the composition of SAL Tailing are described in Rodríguez-Pacheco (2002). In synthesis, the predominant mineral is hematite (about 70%), followed by aluminum (about 10%) and quartz (3%). The SAL residue has little organic matter content, on average 0.6%. The SAL Tailing has a high specific weight (Table 1) with very fine granulometry.

Table 1
Characteristics and physical properties of the column used in the contaminant transport tests (Rodríguez-Pacheco, 2002).

Column characteristics						Physical indices					Velocity
L (cm)	D (cm)	V (cm3)	M (g)	Vp (cm3)	ρd (g/cm3)		ρ (g/cm3)	ρs (g/cm3)	n	θ	vp (cm/h)
5	1.6	10.55	15.52	6.06	1.56		2.15	3.88	0.61	0.61	1.2

• where, with respect to column characteristics, L is the column length [L], D is the diameter [L], V is the total volume [L³], M is the mass of tailing [M], and _p is the pore volume [L³].

Rodríguez-Pacheco (2002) used two suction control methods to collect experimental data for the soil water retention curve to cover the most extensive possible range of suction values. Thus, the psychrometric method was chosen for suction values ranging from 100 to 10 000 kPa, in small, cylindrical samples (15 mm in diameter and 12 mm in height) compacted with controlled moisture and an initial dry density to complement the suction-controlled oedometer test (axis translation technique) for suction values ranging from 10 to 900 kPa, which were difficult to measure accurately using the psychrometric method.

The soil water retention curve of the experimental data is shown in Figure 1 for drying and wetting trajectories determined using the psychrometric and suction-controlled oedometric techniques for remolded samples with an initial void ratio of 1.75. For the tests performed by Rodríguez-Pacheco (2002), Figure 1 shows that the two suction measurement techniques overlap well. However, more tests are needed to reach such a conclusion. This finding indicates that, in this study, the osmotic suction corresponds to a small portion of the total suction because the psychrometric technique measures the total suction. The oedometric technique with axis translation measures only the matric suction.

Figure 1
Soil water retention curve for the drying and wetting trajectories of SAL tailings (Rodríguez-Pacheco, 2002).

Figure 1 also shows the phenomenon of hysteresis, defined by the different drying and wetting paths. In general, more water is retained by the system during drying than is adsorbed by the system at the same magnitude of suction during wetting. The possible causes of the phenomenon are the contact angle between a fluid and a solid, entrapment of air, and different spatial connectivity of pores during drying and wetting (Gennes et al., 2004).

Rodríguez-Pacheco (2002) performed a laboratory batch test to construct the Ni²⁺ solute sorption isotherm. This test was performed at a controlled temperature of 22 ± 2 °C, with an initial pH of 6.9. The metal was dissolved in 0.01 mM KNO₃ with an electrolyte support of pH = 5.5. This solution is the same as that used by Rodríguez-Pacheco (2002) in column tests. Rodríguez-Pacheco (2002) used a Ni salt, (NO₃)₂6H₂O, to prepare the solution of the metal. In Figure 2, the experimental Ni²⁺ sorption isotherm data reveal a linear behavior.

Figure 2
Ni²⁺ sorption isotherm from SAL tailings (Rodríguez-Pacheco, 2002).

In results showed by Rodríguez-Pacheco (2002), nickel sorption had different behaviors when tested in SAL and ACL tailings. In ACL, sorption showed a non-linear isotherm. In SAL, a linear isotherm is well adjusted (Figure 2). It is important to note that the equilibrium distribution coefficient depends on the solute concentration. Also, according to Godoy et al. (2019), the sorption phenomenon is strongly correlated to the cation exchange capacity and significantly correlated to mesoporosity and microporosity of the porous medium.

Rodríguez-Pacheco (2002) conducted laboratory column tests with a stainless-steel column. He placed a polyvinyl chloride (PVC) membrane, which is inert and resistant to high pressures, to construct the breakthrough curve of Ni²⁺. The tailings samples were packed in the columns in compressed layers and subjected to vibrations until reaching a uniform density of 1.55 g/cm³. The columns were saturated with an electrolyte solution of 0.1 mM KNO₃ for 24 hours to fully saturate the pores, which would avoid any cavities, which would favor preferential flow. The same electrolyte solution was passed through the column until steady-state conditions were reached, which was achieved when the flow rate, pH, and electrical conductivity of the solution entering and leaving the column were equal. The tailings columns were loaded with solutes in a steady-state and a constant downward flow. Effluents were collected at varying time intervals. Table 1 outlines the different parameters of the column used in the contaminant transport tests.

The column test was performed by continuously injecting 91 pore volumes of solution for the sorption process and 127 pore volumes of solute-free solution for the desorption process. Table 2 outlines the main conditions and characteristics of the column test of Ni²⁺ with sorption and desorption processes in SAL tailings used by Rodríguez-Pacheco (2002).

Table 2
Characteristics and concentrations in the column test (Rodríguez-Pacheco, 2002).

Contaminant	Vps	Vpd	v(cm/h)	C0 (mg/L)	Cw (mg/L)
Ni2+	91	127	1.2	534	5.72

• where _ps is the pore volume of the solution injected in sorption and _pd is the pore volume of the contaminant-free solution injected in desorption.

Figure 3 presents the experimental data of the breakthrough curve of Ni²⁺ constructed from the column test with SAL tailings. This figure shows that 14 pore volumes must be passed to reach a normalized concentration equal to or close to 1 (c_w /c₀ ≈ 1).

Figure 3
The breakthrough curve of Ni²⁺ from SAL tailings (Rodríguez-Pacheco, 2002).

3.2 Method

Figure 4 presents a flowchart of the method used to implement the model for contaminant transport in unsaturated porous media.

Figure 4
Flowchart of the method of analysis.

The first stage included calibration of the solute flow and the experimental transport data through analytical formulations. Then, the strategy proposed by Bear & Cheng (2010) was used to fit the experimental data. The authors recommend using the least-squares minimization method, which consists of minimizing the error, E, expressed as the squared sum of the difference between a sample of actual values and a sample of estimated values, and is described by the following equation:

$E={\displaystyle \sum _{i=1}^n{\left[{\left(c/c_0\right)}_m-{\left(c/c_0\right)}_c\right]}^2}$ (39)

where (c/c₀)_m is the difference between the experimental normalized concentration and(c/c₀)_c is the theoretically calculated normalized concentration.

The second stage included contaminant transport simulations for unsaturated conditions. Modeling was performed using analytical formulations and mathematical codes to simultaneously analyze the solute flow and the effect of varying the volumetric water content on the contaminant plume using the data obtained from the calibration process. To simulate the contaminant infiltration process and, therefore, the variation in the volumetric water content in unsaturated conditions, the analytical solution by Cavalcante & Zornberg (2017) was used. This solution presents a distribution of the volumetric water content within the porous medium that varies in space and time. In turn, the analytical solution of the ADE for contaminant transport by Brenner (1962), rewritten in terms of the volumetric water content, was used to simulate the advance of the contaminant plume.

3.3 Model calibration

The calibrations of the experimental data for the hydraulic characterization tests and the contamination of SAL tailings are presented in this section.

Figure 5 shows the experimental soil water retention curve for drying and wetting trajectories with an initial void ratio of 1.75. These experimental data were fitted using the constitutive model by Cavalcante & Zornberg (2017). The fitted parameters and their error values are presented in Table 3. The analysis of the error values shows that the constitutive model by Cavalcante & Zornberg (2017) fitted the experimental data of the soil water retention curve very well for wetting and drying trajectories.

Figure 5
Constitutive Model of Cavalcante & Zornberg (2017) applied for drying and wetting trajectories of the SAL tailings.

Table 3
Parameters of the soil water retention curve fitted using the Constitutive Model of Cavalcante & Zornberg (2017).

Trajectory	e	θs	θr	δ	E
Drying	1.75	0.636	0.049	0.00069	0.0162
Wetting	1.75	0.635	0.036	0.00201	0.0199

• where e is the void ratio, θs is the saturated volumetric water content, θr is the residual volumetric water content, δ is the fitting parameter, and E is the adjust error.

Figure 6 shows the Ni⁺² solute sorption isotherm as a linear curve in the experimental data. Therefore, these data were fitted using a linear regression isotherm. The value of the equilibrium distribution coefficient, K_d = 1.31 L/kg, was calculated from this fitting adjustment.

Figure 6
Ni²⁺ sorption isotherm fitted using the linear model for the SAL tailings.

In Figure 7, the experimental elution (sorption – desorption) curve of Ni²⁺ is shown. These experimental data were model using the analytical solution for the ADE by Brenner (1962). The parameters R and P, calibrated using the ADE, and the error value is presented in Table 4. Therefore, the analytical model by Brenner (1962) reliably represents the experimental elution curve of Ni⁺².

Figure 7
Breakthrough curve of Ni²⁺ modeled by Brenner (1962).

Table 4
Calibrated parameters of the ADE model – Brenner (1962).

Solute	R	P	Dh (m2/s)	αd (m)	E
Ni2+	4.22	8	2.08×10-8	0.00625	0.39

• where R is the retardion factor, P is the Peclet number, _h is hydrodynamic dispersion coefficient, _d is longitudinal dispersivity coefficient and E is the adjust error.

The retardation factor of the studied contaminant is greater than 4 (Table 4). Therefore, the tailings have a good capacity to sorb the studied solute. The value of the Peclet number (P = 8) shows that advection and mechanical dispersion are the predominant mechanisms in the column test and that molecular diffusion processes are negligible. The value of the error (E = 0.39) shows that the model by Brenner (1962) reliably fits the experimental data of the breakthrough curve.

4. Results and discussion

The results and analyses of the modeling of the contaminant transport in SAL tailings in unsaturated conditions are presented in this section. In the simulations, a continuous source of contamination is placed in the upper limit of the SAL tailings of the column under analysis. Also, the initial volumetric water content of the SAL tailings is considered to be already equal to the residual volumetric water content (θ_r), with an initial concentration equal to zero; i.e., before the unsaturated flow phenomenon occurs, the tailings profile already contains moisture, and the initial concentration of nickel is zero. This model also assumes that the tailings column is a finite medium.

The parameters used for the one-dimensional modeling of contaminant transport in unsaturated conditions, which describe both the hydraulic and contaminant transport parameters, were calculated as follows:

• The saturated volumetric water content (θ_s = 0.635), residual volumetric water content (θ_r = 0.036) and fitting hydraulic (δ = 0.00201) parameters were calculated by calibrating the experimental data of the soil water retention curve of the SAL tailings using the constitutive model by Cavalcante & Zornberg (2017).

• The hydraulic gradient (i = 1.1) and real average saturated hydraulic conductivity (k_s = 5.26×10^-6 m/s) of the tailings were obtained from the doctoral thesis of Rodríguez-Pacheco (2002).

• The column length (L = 0.5 m) and the dry density of the tailings (ρ_d = 1560 kg/m³) were calculated from the column test.

• In the column length, z = 0 indicates the top of the sample, where the contaminant transport starts.

• The Ni²⁺ equilibrium distribution coefficient (K_d = 1.31×10^- m³/kg) was calculated by calibrating the experimental data from the batch test using a linear model.

• The values of the molecular diffusion coefficient of Ni²⁺ in an aqueous solution (D₀ = 6.79×10^-10 m²/s) are outlined in Fetter (1999).

• The dispersivity coefficient (α_d = 0.00625 m) of Ni²⁺ was calculated by calibrating the experimental data from the column test using the analytical solution of the ADE model by Brenner (1962).

Figure 8 shows the time history of volumetric water content (θ_w) for different locations. The curves of the volumetric water content start at the value of the residual volumetric water content (θ_r = 0.036) in unsaturated conditions and end at the value of the saturated volumetric water content (θ_s = 0.635) when all voids are filled with the contaminant. This figure shows a transition zone characterized by the marked increase in the volumetric water content with infiltration time.

Figure 8
Predicted time history of volumetric water content at different locations.

The volumetric water content is the ratio of the volume of water to the unit volume of soil. This parameter has an important role for contaminant transport over space and time. As can be observed in Figure 8, the volumetric water content is higher near the top of the specimen, and lower at greater depths. It can also be seen that in the first few minutes, the increasing of the volumetric water content is more pronounced.

Figure 9 illustrates the time history of matric suction at different locations. The data in this figure show that, as expected, the matric suction is maximal when the volumetric water content is minimal (Figure 8) and equal to the residual volumetric water content (θ_r) in unsaturated conditions. Minimal matric suction is also observed when the volumetric water content (Figure 8) is maximal and equal to the saturated volumetric water content (θ_s) in saturated conditions. This figure also shows a transition zone characterized by a marked decrease in matric suction as the volumetric water content increases (Figure 8).

Figure 9
Predicted time history of matric suction at different locations.

Matric suction is the free energy change in a unit volume of water when isothermally transferred from the soil water state to the free water state (Zhang & Lu, 2019). In Figure 9, one can observe that the matric suction is lower near the top of the specimen, and higher at greater depths. It is correct because the top of the specimen became saturated firstly (Figure 8). It can also be seen that in the first few minutes, the decreasing of the matric suction is more pronounced.

Figure 10 shows the time history of unsaturated hydraulic conductivity at different locations. The data in this figure show that, as expected, the hydraulic conductivity is minimal when the volumetric water content (Figure 8) is minimal and equal to θ_r in unsaturated conditions. The data also show that the hydraulic conductivity is maximal when the volumetric water content (Figure 8) is maximal and equal to θ_s in saturated conditions. This figure also shows a transition zone characterized by a marked increase in hydraulic conductivity as the volumetric water content increases (Figure 8).

Figure 10
Predicted time history of unsaturated hydraulic conductivity at different locations.

Unsaturated hydraulic conductivity refers to a measure of water-retaining ability of the soil when the pores are not saturated with water. As can be observed in Figure 10, the unsaturated hydraulic conductivity is higher near the top of the specimen, and lower at greater depths. It is correct because the top of the specimen became saturated firstly (Figure 9). It can also be seen that in the first few minutes, the increasing of the unsaturated hydraulic conductivity is more pronounced. Combining Figures 9 and 10, one can confirm that the matric suction influences the behavior of unsaturated soils in terms of permeability.

Figures 1111b, and 11c illustrate the time history, at different locations, of the degree of saturation, the contaminant percolation rate, and the tortuosity factor, respectively. The data in these figures show that, as expected, the degree of saturation, the percolation rate, and the tortuosity factor are practically zero when the volumetric water content (Figure 8) is minimal and equal to θ_r in unsaturated conditions. Conversely, when the volumetric water content (Figure 8) is maximal and equal to θ_s in saturated conditions, the degree of saturation, percolation rate, and tortuosity factor are maximal. The curves have similar shapes because these physical quantities are directly proportional to the volumetric water content (Figure 8).

Figure 11
Predicted time history, at different locations, of: (a) Degree of saturation, (b) Contaminant percolation rate, and (c) Tortuosity factor.

As can be observed in Figure 11, the degree of saturation, the percolation rate and the tortuosity factor are higher near the top of the specimen, and lower at greater depths. It is correct because the top of the specimen became saturated firstly (Figure 9). It can also be seen that in the first few minutes, the increasing of the unsaturated hydraulic conductivity is more pronounced. At near the top of the specimen, the degree of saturation is reached after 180 min. The tortuosity is an intrinsic property of a porous material usually defined as the ratio of actual flow path length to the straight distance between the ends of the flow path (Bear, 1988). Therefore, the tortuosity factor is expected to change during the transport of the contaminant.

Figure 12 shows the time history of the hydrodynamic dispersion coefficient of Ni²⁺ at different locations. This figure shows that the hydrodynamic dispersion coefficient is minimal when the volumetric water content (Figure 8) is minimal and equal to θ_r in unsaturated conditions. In these conditions, the hydrodynamic dispersion coefficient equals the molecular diffusion coefficient because the percolation rate is zero; therefore, contaminants are only transported at a microscopic scale. Conversely, the hydrodynamic dispersion coefficient is maximal when the volumetric water content (Figure 8) is equal to θ_s in saturated conditions. In these conditions, the hydrodynamic dispersion coefficient equals the mechanical dispersion coefficient; therefore, contaminants are predominantly transported at a macroscopic scale. This figure also shows a transition zone characterized by a marked increase in the hydrodynamic dispersion coefficient as the volumetric water content increases (Figure 8).

Figure 12
Predicted time history of hydrodynamic dispersion coefficient of Ni²⁺ at different locations.

The hydrodynamic dispersion coefficient is one of the most important parameters in the prediction of contaminant concentrations in soil using advection–dispersion models. It is a measure for describing the mixing processes of solutes in porous media. As can be observed in Figure 12, the hydrodynamic dispersion coefficient is higher near the top of the specimen, and lower at greater depths. It can also be seen that in the first few minutes, the increasing of the hydrodynamic dispersion coefficient is more pronounced.

Figure 13 illustrates the time history of the retardation factor of Ni²⁺ at different locations. A comparison of Figures 13 and 8 shows that the retardation factor is maximal when the volumetric water content is minimal and equal to the residual volumetric water content θ_r in unsaturated conditions. Therefore, contaminants are transported only at a microscopic scale by the mechanism of molecular diffusion because the percolation rate is virtually zero. Also, the retardation factor is minimal when the volumetric water content (Figure 8) is maximal and equal to the saturated volumetric water content θ_s; thus, the percolation rate is also maximal in saturated conditions. Therefore, contaminants are predominantly transported at a macroscopic scale due to advection and mechanical dispersion mechanisms. Figure 13 also shows a transition zone characterized by the marked increase in the retardation factor when the volumetric water content increases (Figure 8).

Figure 13
Predicted time history of retardation factor of Ni²⁺ at different locations.

The retardation factor is a measure of the ability of the ground to restrain the migration of the contaminant. It shows how many times the migration of the substance is subjected to adsorption slower than the actual speed of water flow in the pore spaces (Mikołajków, 2003). In Figure 9, one can observe that the retardation factor is lower near the top of the specimen, and higher at greater depths. It is correct because the top of the specimen became saturated firstly (Figure 8) and the unsaturated hydraulic conductivity is greater (Figure 10). It can also be seen that in the first few minutes, the decreasing of the retardation factor is more pronounced.

Figure 14 shows a time history of the normalized concentration of Ni²⁺ at different locations. An analysis of Figure 14 shows that, as expected, as the depth increases, the elution curve shifts to the right, thereby increasing the time for the normalized concentration to approach one in the ascending section of the curve (sorption process) as well as the descending section of the curve (desorption process). Figure 14 also shows that as the infiltration time increases, the normalized concentration increases proportionally until it reaches its maximum in the sorption process. As the infiltration time increases, the normalized concentration decreases until reaching its minimum in the desorption process.

Figure 14
Predicted time history of normalized concentration of Ni²⁺ at different locations.

As can be observed in Figure 14, the normalized concentration of Ni²⁺ is higher near the top of the specimen, and lower at greater depths, during the sorption process. It can also be seen that it becomes lower near the top of the specimen, and higher at greater depths, during the desorption process. For this experiment, after 1500 min, the desorption process begins.

5. Conclusions

Based on the experimental soil water retention curve of SAL tailings, the osmotic suction of these tailings is negligible. This is important because the model by Cavalcante & Zornberg (2017) adopted in this study disregards osmotic suction effects and considers the gradient of osmotic suction as approximately zero, since it is related to the concentration of salts in the soil and its variations over space, in general, are not significant. Although for the data of this research, the results obtained by the two types of equipment resulted in a small contribution of osmotic suction, it is important to emphasize that for contaminated soils, this is not always observed.

The constitutive models proposed by Cavalcante & Zornberg (2017) have only a single fitting parameter in contrast to other constitutive models, such as the model by Van Genuchten (1980), which makes it easier to understand and mathematically handle the model because, in addition to being a single parameter, this parameter has a physical meaning. Except for the hydraulic fitting parameter (δ), all other parameters (θ_r, θ_s, and k_s) of the constitutive model can be directly calculated from geotechnical tests. These constitutive models reliably fit the experimental data of the soil water retention curve of the SAL tailings for the wetting and drying trajectories.

The analytical solution of the Richards equation by Cavalcante & Zornberg (2017) for unsaturated and transient flow conditions represents the variation in the volumetric water content in one-dimensional space and time. Controlling the volumetric water content was crucial in modeling the contaminant infiltration process and, therefore, the advance of the contaminant plume in space and time in unsaturated conditions in SAL tailings. Furthermore, this relation allows us to use the various analytical solutions for contaminant transport in saturated conditions for unsaturated conditions by replacing the porosity with the volumetric water content.

The unsaturated hydraulic conductivity function was derived using constitutive models proposed by Cavalcante & Zornberg (2017), based on the fitting parameter δ, obtained from the soil water retention curve and the saturated hydraulic conductivity. These parameters contributed to a reliable representation of the hydraulic conductivity of the SAL tailings and, therefore, to the determination of the percolation rate, which is used to determine the contaminant mass flow by advection and hydrodynamic dispersion.

The calibration of the experimental data of the breakthrough curve of Ni²⁺ showed that the retardation factor is greater than 4. Hence, the SAL tailings have a good capacity of sorption for the studied solute, which is environmentally beneficial due to the ability of the SAL tailings to retain and decrease the advance of this metal through its porous matrix.

The analytical solution by Cavalcante & Zornberg (2017) adequately represented the contaminant infiltration process and, therefore, the variation in the volumetric water content in unsaturated conditions in space and time, and the analytical solution by Brenner (1962) adequately represented the contaminant plume in space and time through the SAL tailings, with the initial and boundary conditions established in our problem. Controlling the volumetric water content was crucial for analyzing its effect on the contaminant transport mechanisms and contaminant plume.

Based on analytical solutions, the present model has simplifications with specific initial and boundary conditions and a linear sorption isotherm (mathematically expressed in terms of the equilibrium distribution coefficient). Despite its limitations, it adequately represented the contamination transport of nickel in SAL tailing in Cuba.

The model advantages are that it allows validating numerical approaches, providing insights of many complex processes that occur in the field and in the lab and field and requires far less computational effort compared with current programs to modeling the solute transport using numerical solutions, as the versatile commercial Software HYDRUS 2D/3D (Šimůnek et al., 2016).

As expected, in the process of simulating contaminant transport for SAL tailings in unsaturated conditions, the volumetric water content is minimal and equal to the residual volumetric water content, which implies that the hydraulic conductivity and, consequently, the percolation rate are virtually zero. Furthermore, the retardation factor is maximal. Therefore, in these conditions, contaminants are transported only at a microscopic scale by the mechanism of molecular diffusion because the percolation rate is zero. As the infiltration progresses, the volumetric water content increases until reaching its maximum value, which is equal to the saturated volumetric water content; therefore, the hydraulic conductivity and the percolation rate also increase until reaching their maximum values, and the retardation factor decreases until reaching its minimum value. This occurs when the tailings are in saturated conditions. Therefore, contaminants are predominantly transported at a macroscopic scale by the mechanisms of advection and mechanical dispersion.

List of symbols

ACL Ammoniacal Carbonate Leaching

ADE Advection-Dispersion Equation

SAL Sulfuric Acid Leaching

Ni²⁺ Nickel (II)

(c/c₀)_c Theoretically calculated normalized concentration.

(c/c₀)_m Experimental normalized concentration

a_s Unsaturated advective flow

C Normalized concentration

c₀ Initial concentration of the applied solute

c_s Adsorbed contaminant concentration

c_w Contaminant concentration

D * Molecular diffusion coefficient

D₀ Molecular diffusion coefficient in the aqueous solution

D_h Hydrodynamic dispersion coefficient

D_m Mechanical dispersion coefficient

Dz Unsaturated diffusivity of water in the z-direction

e Void Ratio

E Adjust Error

g Gravitational acceleration

i Hydraulic gradient

J* Mass flow by molecular diffusion of the contaminant per unit area and per unit time

J_a Mass flow by advection of the contaminant per unit area and per unit time

J_m Mass flow by mechanical dispersion of the contaminant per unit area per unit time

K_d Equilibrium distribution coefficient

k_s Saturated hydraulic conductivity

k_z Hydraulic conductivity expressed in terms of suction in the z-direction

L Sample length

P Peclet number

R Retardation factor

T Normalized time

t Time

t₀ Time of application of the displacing solution

v₀ Infiltration rate

v_0,max Maximum infiltration rate

v_p Percolation rate

v_z Fluid velocity in the z-direction

Z Normalized spatial coordinate

z Spatial coordinate

α_d Longitudinal dispersivity coefficient

β_m, eigenvalues

δ Fitting parameter, Cavalcante & Zornberg (2017) coefficient

ϕ Hydraulic head

θ_i Initial volumetric water content.

θ_r Residual volumetric water content

θ_s Saturated volumetric water content

θ_w Volumetric water content

ρ_d Dry density of the soil

ρ_w Density of water

τ Tortuosity factor

ψ Total suction of water

Acknowledgements

This study was financed in part by the Coordination for the Improvement of Higher Education Personnel – Brasil (CAPES) – Finance Code 001. The authors also acknowledge the support of the National Council for Scientific and Technological Development (CNPq Grant 304721/2017-4, 435962/2018-3, 140923/2020-9 and 305484/2020-6), the Foundation for Research Support of the Federal District (FAPDF) (Projects 0193.002014/2017-68 and 0193.001563/2017), and the University of Brasília.

• Erratum

  The article entitled “Contaminant transport model in transient and unsaturated conditions applied to laboratory column test with tailings” (DOI https://doi.org/10.28927/SR.2022.076021), from Soils and Rocks 45(2):e2022076021 (2022), was published with a typo that caused misspelling of the first author’s name.
  On page 1, where it reads:
  Eliu James Carbaja¹ 0000-0003-2551-8845
  It should read:
  Eliu James Carbajal¹ 0000-0003-2551-8845
  On pages 1-15, where the text reads:
  Carbaja et al.
  It should read:
  Carbajal et al.
  The publisher apologizes for the errors.

References

• Adnan, S., Iqbal, J., Maltamo, M., & Valbuena, R. (2018). GIS-based DRASTIC model for groundwater vulnerability and pollution risk assessment in the Peshawar District, Pakistan. Arabian Journal of Geosciences, 11(16), 458. http://dx.doi.org/10.1007/s12517-018-3795-9
  » http://dx.doi.org/10.1007/s12517-018-3795-9
• Akbariyeh, S., Bartelt-Hunt, S., Snow, D., Li, X., Tang, Z., & Li, Y. (2018). Three-dimensional modeling of nitrate-N transport in vadose zone: roles of soil heterogeneity and groundwater flux. Journal of Contaminant Hydrology, 211, 15-25. http://dx.doi.org/10.1016/j.jconhyd.2018.02.005
  » http://dx.doi.org/10.1016/j.jconhyd.2018.02.005
• Banzhaf, S., & Hebig, K.H. (2016). Use of column experiments to investigate the fate of organic micropollutants: a review. Hydrology and Earth System Sciences, 20(9), 3719-3737. http://dx.doi.org/10.5194/hess-20-3719-2016
  » http://dx.doi.org/10.5194/hess-20-3719-2016
• Barone, F.S., Rowe, R.K., & Quigley, R.M. (1992). A laboratory estimation of diffusion and adsorption coefficients for several volatile organics in a natural clayey soil. Journal of Contaminant Hydrology, 10(3), 225-250. http://dx.doi.org/10.1016/0169-7722(92)90062-J
  » http://dx.doi.org/10.1016/0169-7722(92)90062-J
• Bear, J. (1988). Dynamics of fluids in porous media (761 p.). Courier Corporation.
• Bear, J., & Cheng, A. (2010). Modeling ground water flow and contaminant transport. New York: Springer. https://link.springer.com/book/10.1007%2F978-1-4020-6682-5
  » https://link.springer.com/book/10.1007%2F978-1-4020-6682-5
• Bertolo, R.A. (2001). Hidrodinâmica and hidrogeoquímica da Zona não saturada do Aqüífero Adamantina em Urânia-SP [Thesis Doctoral, University of São Paulo]. University of São Paulo (in Portuguese).
• Brenner, H. (1962). The difussion model of longitudinal mixing in beds of finite length. Numerical values. Chemical Engineering Science, 17, 229-243. http://dx.doi.org/10.1016/0009-2509(62)85002-7
  » http://dx.doi.org/10.1016/0009-2509(62)85002-7
• Buckingham, E. (1907). Studies on the movement of soil moisture (USDA Bureau of Soils, Bull. no. 38). Washington, DC: Government Printing Office.
• Cavalcante, A.L.B., & Zornberg, J.G. (2017). Efficient approach to solving transient unsaturated flow problems. I: analytical solutions. International Journal of Geomechanics, 17(7), 04017013. http://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0000875
  » http://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0000875
• Fetter, C.W. (1999). Contaminant hydrogeology (500 p.). New Jersey: Prentice Hall.
• Fityus, S.G., Smith, D.W., & Booker, J.R. (1999). Contaminant transport through an unsaturated soil liner beneath a landfill. Canadian Geotechnical Journal, 36(2), 330-354. http://dx.doi.org/10.1139/t98-112
  » http://dx.doi.org/10.1139/t98-112
• Foo, K.Y., & Hameed, B.H. (2010). Insights into the modeling of adsorption isotherm systems. Chemical Engineering Journal, 156(1), 2-10. http://dx.doi.org/10.1016/j.cej.2009.09.013
  » http://dx.doi.org/10.1016/j.cej.2009.09.013
• Freeze, R. A. & Cherry, J.A. (1979). Groundwater (604 p.). Pretice Hall.
• Gennes, P.-G., Brochard-Wyart, F., & Quéré, D. (2004). Capillarity and wetting phenomena: drops, bubbles, pearls, waves (308 p.). New York, USA: Springer Verlag.
• Godoy, V.A., Zuquette, L.V., & Gómez-Hernández, J.J. (2019). Spatial variability of hydraulic conductivity and solute transport parameters and their spatial correlations to soil properties. Geoderma, 339, 59-69. http://dx.doi.org/10.1016/j.geoderma.2018.12.015
  » http://dx.doi.org/10.1016/j.geoderma.2018.12.015
• Hillel, D. (2003). Introduction to environmental soil physics (498 p.). Academic press.
• Joshi, P., & Gupta, P.K. (2018). Assessing groundwater resource vulnerability by coupling GIS-Based DRASTIC and solute transport model in Ajmer District, Rajasthan. Journal of the Geological Society of India, 92(1), 101-106. http://dx.doi.org/10.1007/s12594-018-0958y
  » http://dx.doi.org/10.1007/s12594-018-0958y
• Mallants, D., Van Genuchten, M.T., Jacques, D., & Seetharam, S. (2011). Leaching of Contaminats to Groundwater. In F.A. Swartjes (Ed.), Dealing with contaminated sites: from theory to practical application (pp. 787–850). Dordrecht, Netherlands: Springer Verlag. https://doi.org/10.1007/978-90-481-9757-6_18
  » https://doi.org/10.1007/978-90-481-9757-6_18
• Mikołajków, J. (2003). Laboratory methods of estimating the retardation factor of migrating mineral nitrogen compounds in shallow groundwater. Geological Quarterly, 47, 91-96.
• Mustafa, S., Bahar, A., Aziz, Z.A., & Suratman, S. (2016). Modelling contaminant transport for pumping wells in riverbank filtration systems. Journal of Environmental Management, 165, 159-166. http://dx.doi.org/10.1016/j.jenvman.2015.09.026
  » http://dx.doi.org/10.1016/j.jenvman.2015.09.026
• Narasimhan, T.N. (2004). Darcy’s law and unsaturated flow. Vadose Zone Journal, 3(4), 1059. http://dx.doi.org/10.2113/3.4.1059
  » http://dx.doi.org/10.2113/3.4.1059
• Porro, I., Newman, M.E., & Dunnivant, F.M. (2000). Comparison of batch and column methods for determining strontium distribution coefficients for unsaturated transport in basalt. Environmental Science & Technology, 34(9), 1679-1686. http://dx.doi.org/10.1021/es9901361
  » http://dx.doi.org/10.1021/es9901361
• Richards, L.A. (1931). Capillary conduction of liquids through porous mediums. Physics, 1(5), 318-333. http://dx.doi.org/10.1063/1.1745010
  » http://dx.doi.org/10.1063/1.1745010
• Rodríguez-Pacheco, R.L. (2002). Estudio experimental de flujo y transporte de cromo, níquel y manganeso en residuos de la zona minera de Moa (Cuba): influencia del comportamiento hidromecánico [Thesis Doctoral, Universidad Politecnica de Cataluña]. Universidad Politecnica de Cataluña (in Spanish). https://www.tdx.cat/handle/10803/6223
  » https://www.tdx.cat/handle/10803/6223
• Rutsch, M., Rieckermann, J., Cullmann, J., Ellis, J.B., Vollertsen, J., & Krebs, P. (2008). Towards a better understanding of sewer exfiltration. Water Research, 42(10-11), 2385-2394. http://dx.doi.org/10.1016/j.watres.2008.01.019
  » http://dx.doi.org/10.1016/j.watres.2008.01.019
• Seferou, P., Soupios, P., Kourgialas, N.N., Dokou, Z., Karatzas, G.P., Candasayar, E., Papadopoulos, N., Dimitriou, V., Sarris, A., & Sauter, M. (2013). Olive-oil mill wastewater transport under unsaturated and saturated laboratory conditions using the geoelectrical resistivity tomography method and the FEFLOW model. Hydrogeology Journal, 21(6), 1219-1234. http://dx.doi.org/10.1007/s10040-013-0996-x
  » http://dx.doi.org/10.1007/s10040-013-0996-x
• Šimůnek, J., & Van Genuchten, M.T. (2008). Modeling nonequilibrium flow and transport processes using HYDRUS. Vadose Zone Journal, 7(2), 782-797. http://dx.doi.org/10.2136/vzj2007.0074
  » http://dx.doi.org/10.2136/vzj2007.0074
• Šimůnek, J., Van Genuchten, M.T., & Šejna, M. (2016). Recent developments and applications of the HYDRUS computer software packages. Vadose Zone Journal, 15(7), 1-25. http://dx.doi.org/10.2136/vzj2016.04.0033
  » http://dx.doi.org/10.2136/vzj2016.04.0033
• Sopilniak, A., Elkayam, R., Rossin, A.V., & Lev, O. (2017). Emerging organic pollutants in the vadose zone of a soil aquifer treatment system: pore water extraction using positive displacement. Chemosphere, 190, 383-392. http://dx.doi.org/10.1016/j.chemosphere.2017.10.010
  » http://dx.doi.org/10.1016/j.chemosphere.2017.10.010
• Sosa, E.R. (2016). Caracterização e aproveitamento dos rejeitos oriundos de processos hidrometalúrgicos do níquel e cobalto com um enfoque geoambiental [Thesis Doctoral em Geotecnia. Publicação: G.TD-123/16]. Department of Civil and Environmental Engineering, University of Brasília. (in Portuguese). https://repositorio.unb.br/handle/10482/22164
  » https://repositorio.unb.br/handle/10482/22164
• Srinivasan, P., & Mercer, J.W. (1988). Simulation of biodegradation and sorption processes in ground water. Ground Water, 26(4), 475-487. http://dx.doi.org/10.1111/j.1745-6584.1988.tb00414.x
  » http://dx.doi.org/10.1111/j.1745-6584.1988.tb00414.x
• Van Genuchten, M.T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44, 892-898. http://dx.doi.org/10.2136/sssaj1980.03615995004400050002x
  » http://dx.doi.org/10.2136/sssaj1980.03615995004400050002x
• Van Genuchten, M.T., & Dalton, F.N. (1986). Models for simulating salt movement in aggregated field soils. Geoderma, 38(1), 165-183. http://dx.doi.org/10.1016/0016-7061(86)90013-3
  » http://dx.doi.org/10.1016/0016-7061(86)90013-3
• Viotti, P., Papini, M.P., Stracqualursi, N., & Gamba, C. (2005). Contaminant transport in an unsaturated soil: laboratory tests and numerical simulation model as procedure for parameters evaluation. Ecological Modelling, 182(2), 131-148. http://dx.doi.org/10.1016/j.ecolmodel.2004.07.014
  » http://dx.doi.org/10.1016/j.ecolmodel.2004.07.014
• Zhang, C., & Lu, N. (2019). Unitary definition of matric suction. Journal of Geotechnical and Geoenvironmental Engineering, 145(2), 02818004. http://dx.doi.org/10.1061/(ASCE)GT.1943-5606.0002004
  » http://dx.doi.org/10.1061/(ASCE)GT.1943-5606.0002004

Publication Dates

• Publication in this collection
  13 May 2022
• Date of issue
  2022

History

• Received
  20 Oct 2021
• Accepted
  21 Feb 2022