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

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


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.

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

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: 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: 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: 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: Now, by combining Equations 1 and 4, the Richards equation for one-dimensional, unsaturated, transient flow in the z-direction is obtained as follows: where v 0 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: For a soil column of finite length L, the lower boundary condition adopted in this study is described as follows: 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): and β m, are the eigenvalues that correspond to the positive roots of the following equation: 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: where erfc(Z) is the complementary error function, which is defined as follows:

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: 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: 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: The mass flow by mechanical dispersion in unsaturated porous media is described by Fick's first law as follows: 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: where D 0 is the molecular diffusion coefficient in the aqueous solution [L 2 .T -1 ] and τ = τ(z,t) is the tortuosity factor [L.L -1 ], which is defined as follows (Viotti et al., 2005): The mass flow by molecular diffusion in unsaturated porous media is described by Fick's first law as follows: 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: 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: 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: 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 ( ) Upper boundary condition where c 0 is the initial concentration of the applied solute and is constant in the upper boundary [M.L -3 ], and t 0 is the time of application of the displacing solution.
Lower boundary condition where L is the study sample length [L]. Under these hypotheses, the analytical solution by Brenner (1962) is as follows: Another way of presenting the analytical solution by Brenner (1962) is in its dimensionless form as follows: where: where C = C(Z,T) is the normalized concentration (C = c w /c 0 ); T is the normalized time (T = v p t / 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.

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.
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 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 2+ 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 3 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 3 ) 2 6H 2 O, to prepare the solution of the metal. In Figure 2, the experimental Ni 2+ sorption isotherm data reveal a linear behavior.
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 2+ . 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 3 . The columns were saturated with an electrolyte solution of 0.1 mM KNO 3 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 2+ with sorption and desorption processes in SAL tailings used by Rodríguez-Pacheco (2002). Figure 3 presents the experimental data of the breakthrough curve of Ni 2+ 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 0 ≈ 1).    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:

Method
where (c/c 0 ) m is the difference between the experimental normalized concentration and(c/c 0 ) 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.

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 6 shows the Ni +2 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.
In Figure 7, the experimental elution (sorptiondesorption) curve of Ni 2+ 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, where V ps is the pore volume of the solution injected in sorption and V pd is the pore volume of the contaminant-free solution injected in desorption.  the analytical model by Brenner (1962) reliably represents the experimental elution curve of Ni +2 . 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.

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 3 ) 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 2+ equilibrium distribution coefficient (K d = 1.31×10m 3 /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 2+ in an aqueous solution (D 0 = 6.79×10 -10 m 2 /s) are outlined in Fetter (1999). • The dispersivity coefficient (α d = 0.00625 m) of Ni 2+ was calculated by calibrating the experimental data from the column test using the analytical solution of the ADE model by Brenner (1962). 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. where R is the retardion factor, P is the Peclet number, D h is hydrodynamic dispersion coefficient, α d is longitudinal dispersivity coefficient and E is the adjust error.  Contaminant transport model in transient and unsaturated conditions applied to laboratory column test with tailings 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.
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).
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).
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 11a, 11b, 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).
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 2+ 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).
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 2+ 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).
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 2+ 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.
As can be observed in Figure 14, the normalized concentration of Ni 2+ 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.

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

Retardation Factor
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 2+ 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.