Numerical and experimental study on solute transport through physical aquifer model

The hydrogeological properties of the subsurface zone are responsible for maintaining the quality of groundwater resources (Al Kuisi et al. 2009). Major processes such as advection, diffusion, dispersion, and sorption are accountable for deteriorating groundwater quality. The migration of dissolved contaminants in porous media is mainly controlled by advective transport through the fractures combined with the diffusion of solutes into the relatively immobile pore water in the soil matrix between fractures. The diffusion process, influencing both reactive and non-reactive tracers, can attenuate solute migration relative to the velocity of flow in the fractures and is often referred to as matrix diffusion (Chiogna et al. 2010). The attenuation of the solute particles due to chemical mixing depends on the nature of the sediment, the contaminant, and the geochemical environment (Belhachemi & Addoun 2011).

Analytical solutions of one-dimensional solute transport problems, subject to different initial and boundary conditions, infinite and the semi-finite domain have been reported by many researchers (Ogata & Banks 1961; Neville et al. 2000; Gao et al. 2010; Joshi et al. 2012; Leij et al. 2012; Sharma et al. 2015; Masciopinto & Passarella 2018). Solutions of two-, three-dimensional deterministic ADE's have been investigated in numerous studies and are still being actively studied (Freeze & Cherry 1979; Latinopoulos et al. 1988; Batu 1989; Leij et al. 1991, Batu 1993; Leij et al. 1993; Serrano 1995). A two-dimensional solute transport model for variably saturated porous media with different chemical ions has been investigated by Šimůnek & Suarez (1994). Early arrival and long tailing of breakthrough curves are observed for solute transport through heterogeneous porous media, which shows anomalous transport behaviour (Levy & Berkowitz 2003). Thus, the advection-dispersion equation with a constant dispersion model is less adequate for simulating anomalous transport through heterogeneous porous media. Hence, a time-dependent dispersion model is used to capture the anomalous transport behaviour of a solute through saturated porous media (Huang et al. 1995; Gao et al. 2010). Marinoschi et al. (1999) explored a technique for the analytical solution of one-, two- and three-dimensional solute transport problems with time-dependent dispersion coefficients subjected to different boundary conditions. An analytical solution of a spatially variable coefficient advection-diffusion equation has been solved by Zoppou & Knight (1999) in up to three dimensions with the assumption that the velocity component is proportional to the distance and that the diffusion coefficient is proportional to the square of the corresponding velocity component. James & Jawitz (2007) have discussed a two-dimensional reactive transport ADE model using a splitting technique where advective, dispersive, and reactive parts of the equation were solved separately. An explicit finite-volume Godunov method was used to approximate the advective part, while a hybridized mixed finite element method was used to solve for the dispersive step and the backward Euler method for the reactive component.

Most of the researchers assumed that the dispersion coefficient is time-independent (Corey et al. 1970; Fenske 1979; Selvadurai 2004). Fried (1975) identified the scale effect of the dispersion coefficient due to temporal and spatial variability of the dispersion coefficient. Both field and laboratory-scale experiments show that in subsurface transport problems the dispersivity parameters can be considered as time-dependent (Fenske 1979). Stochastic analyses have shown that the dispersion may depend on travel time and may increase until it reaches an asymptotic value (Gelhar et al. 1979; Gelhar & Axness 1983). Analytical solutions of the solute transport equation considering time-dependent dispersion coefficients have been developed by Pickens & Grisak (1981), Yates (1992), and Barry & Sposito (1989). Yates (1990) obtained one-dimensional analytical solutions for uniform flow with boundary conditions of both constant concentration and constant flux type. Basha & El-Habel (1993) developed an analytical solution of the one-dimensional transport equation with time-dependent dispersion coefficients for an infinite domain. Aral & Liao (1996) derived a general analytical solution of the two-dimensional solute transport equation with time-dependent dispersion coefficients for an infinite domain aquifer.

The method of moment analysis is an important tool for computing mass recoveries in tracer experiments, travel velocities of the plume, and the description of the shape of the plume in terms of dispersivity, skewness, and kurtosis (Aris 1956). Thus, the method of the moment was employed for the study of many natural gradient field tracer tests (Freyberg 1986). The relationship between the spatial and temporal moments and the properties of an evolving solute plume is based on work by Aris (1956). Goltz & Roberts (1987) developed Aris's method of moments, which can be used to analyze the temporal and spatial moment behavior of concentration distributions obtained using both equilibrium and physical non-equilibrium solute transport models.

The spatial moments of the bromide tracer distribution were used to calculate the tracer mass, velocity, and dispersivity during the large-scale tracer test at the Canadian Air Force Base in Borden, Ontario (Freyberg 1986; Farrell et al. 1994). Knorr et al. (2016) performed a saturated dual-porosity soil column experiment and proposed a single fissure dispersion model, which includes the diffusive mass exchange between mobile and immobile water regions. Xu et al. (2019) reviewed the recent advancement in the context of experimental studies of dilution and reactive mixing in saturated porous media at steady-state conditions. Apiratikul (2020) studied the application of an analytical solution for the advection-dispersion reaction model to predict the three-dimensional plot of breakthrough curves and mass transfer for sorption of heavy metal ions using a fixed-bed column.

In this study, an attempt has been made to describe the two-dimensional advection-dispersion transport equation with time-dependent dispersion coefficients. A two-dimensional numerical model is developed using the implicit finite difference method. Also, experiments were conducted in the lab for solute transport through porous media consisting in the rectangular tank, and results of observed breakthrough curves were predicted for both conservative and non-conservative solutes.

The Forward-time central-space (FTCS) method is a finite difference method used for numerically solving the heat equation and similar parabolic partial differential equations. It is a first-order method in time, explicit in time, and is conditionally stable when applied to the heat equation. When used as a method for advection equations, or more generally hyperbolic partial differential equations, it is unstable unless a retardation term is included (Dehghan 2004). The FTCS method is based on the central difference in space and the forward Euler method in time, giving first-order convergence in time and second-order convergence in space. Transport phenomenon is central for understanding many procedures in numerous sciences and is described by the partial differential equation.

The Gauss-Seidel iteration method has been used to solve the set of linear simultaneous equations, and the convergence was found to be satisfactory for almost all cases. Both the Peclet number and Courant number are kept less than one to reduce the numerical error. Therefore, other iterative techniques, which may provide faster convergence, were not explored. The numerical model has been validated with an analytical solution Aral & Liao (1996), as shown in Figure 1. Following input parameters, pore water velocity v = 0.1 m/day, longitudinal dispersion coefficient D_L = 0.035 m²/day, transverse dispersion coefficient D_T = 0.0035 m²/day, and injected plume size of 4.6 m × 5.6 m were used during simulation. The numerical results of breakthrough curves were predicted at y = 7 m and flow direction in the x-axis. Value of asymptotic coefficient C_A = 0 indicates a constant dispersion model, and C_A = 50 day indicates an asymptotic dispersion model. Hence, the numerical model has validated both cases of constant and asymptotic time-dependent dispersion models.

Figure 2(a) shows that the simulated results for the mean travel distance match very well with the observed values. As expected, for ideal transport, the mean travel distance increases linearly with time. Similarly, the simulated results of the second longitudinal and transverse moments match well with the experimental data (Figure 2(b)). Figure 2(c) represents the simulation of transverse variance of experimental field data.

Figure 3(a) and 3(b) show the behavior of mean travel distance and second spatial moment; that is, the variance of solute in the presence of different values of constant dispersion coefficient. In the presence of a smaller value of dispersion coefficient (D_L = 0.035 and 0.175 m²/day), the mean travel distance of the solute plume increases linearly with an increase in travel time. However, in the case of the higher value of dispersion coefficient (D_L = 0.35 m²/day), the mean travel distance of a solute plume increase non-linearly during higher transport time. The behavior of the second spatial moment is increasing non-linearly with the increase in transport time, as shown in Figure 3(b).

Figure 4 shows the variation of the longitudinal dispersion coefficient in the presence of different values of asymptotic constants. The behavior of the longitudinal dispersion coefficient increases non-linearly during short transport time while its value remains constant for higher transport time. It is also shown that the magnitude of the dispersion coefficient is higher for smaller values of asymptotic constant. However, a larger value of the asymptotic constant leads to reduction of the magnitude of the dispersion coefficient.

Figure 5 shows the variation of the longitudinal dispersion coefficient in the presence of different values of linear constants. The behavior of the longitudinal dispersion coefficient increases linearly with an increase in time in the presence of a linear constant. It is also shown that the magnitude of the dispersion coefficient is higher for smaller values of linear constant. However, a more considerable value of the linear constant leads to reduction of the magnitude of the dispersion coefficient.

Figure 6(a) shows the variation of mean travel distance in the presence of different values of asymptotic constants. The behaviour of mean travel distance increases linearly with an increase in transport time. It is also seen that the magnitude of mean travel distance remains the same in the presence of different values of the asymptotic constant. It is also shown that the magnitude of the dispersion coefficient is higher for smaller values of asymptotic constant. However, the larger value of the asymptotic constant leads to reduction of the magnitude of the dispersion coefficient. Figure 6(b) shows the behaviour of second spatial moments; that is, the variance of solute plume for different values of asymptotic constants. The higher value of asymptotic constants leads to reduce the variance; that is, spreading of the solute plume.

Figure 7(a) shows the variation of mean travel distance in the presence of different values of constants for linear time-dependent dispersion coefficients. The behaviour of mean travel distance increases linearly with an increase in transport time. It is also seen that the magnitude of mean travel distance remains the same in the presence of different values of asymptotic constant during early transport time. However, the magnitude of mean travel distance decreases non-linearly for the smaller value of constant during large transport time. But, the behaviour of variance; that is, second spatial moments, is decreasing with the increase of linear constant, as shown in Figure 7(b). Behaviour of variance is non-linear during early transport time but increases at higher transport time.

Figure 8(a) and 8(b) show the behavior of mean travel distance and variance of solute plume for different values of retardation coefficients. The value of the retardation coefficient equal to 1 represents the non-reactive solute plume, and the result shows that the magnitude of mean travel distance increases linearly with an increase in transport time, as shown in Figure 9(a). Values of retardation coefficients equal to 2, 5, and 10 show that solute plume is reactive. It means some of the restive solutes are adsorbed with a solid particle of soil media. Hence, higher values of retardation coefficients lead to reduction in the magnitude of mean travel distance with the increase of transport time. Similarly, the behavior of variance of the solute plume is shown in Figure 8(b). In this case, also, the magnitude of variance is decreasing for the higher value of retardation coefficients.

This section describes the behavior of movement of the solute plume in the form of a contour plot in two-dimensional porous media. Figure 9(a) shows the initial location of the solute plume in the two-dimensional porous media. The value of transport parameters; that is, pore water velocity, = 0.1 m/day and longitudinal dispersivity = 0.05 m, and transverse dispersivity = 0.005 m, is used during the simulation of the results. Figure 9(b) and 9(c) show the contour plots of the solute plume at transport times of 200 days and 600 days, respectively. These results show that the solute plume is moving in the flow direction, and it is also spreading in both longitudinal and transverse directions.

Higher values of longitudinal and transverse dispersivity equal to 0.5 m and 0.05 m are used for predicting the results of the movement of the solute plume. Results of various contour plots of the solute plume at transport time of 200 days and 600 days are shown in Figure 10(a) and 10(b). These results show that the spreading of the solute plume is more in case of the higher value of dispersivity as compared to a smaller value of dispersivity.

Figure 11(a)–11(c) show the results of contour plots for reactive solute plume for different values of retardation coefficient. Results show that the size of the solute plume is reduced in the presence of a higher value of retardation coefficient.

The experimental setup consists of a 5 m × 1 m × 1 m dimensioning iron box, provided with a plexiglass material sheet of a thickness of 2 mm surrounding the lateral boundary (Figure 12). The bottom surface is provided with waterproof plywood sealed with multipurpose sealant to prevent leakage during the experiment. The iron box is painted with galvanized paint to prevent rusting and weathering effects. The inlet and outlet ends of the tank model were provided with hollow PVC pipes of diameter 2.54 cm and are sealed with a sealant material around the pipe peripherals. For investigating the tracer transport, the length of each pipe under investigation was prepared at 0.2 m, and perforations were provided at the pipe surface by the drilling machine for a full depth. The basement of the pipe was sealed to restrict the flow, and the pipe was covered with thin, clean white clothes for making the entrance of water and tracer particles into the pipes. In the homogeneous sandy porous medium, these pipes were inserted at regular intervals of 0.2 m along the transverse direction with a centrally located pipe at a distance of 0.5 m throughout the length of the porous media. The longitudinal interval between the pipes was kept at 0.5 m. The pipes were also inserted vertically at a depth interval of 0.05 m and 0.1 m from the topsoil surface to monitor the flow depth-wise.

The porous material was obtained from excavating the soil available at the site-location to a depth of about 5 m from the top surface and was filled up to the depth of 30 cm from the bottom of the tank. The soil was excavated to conduct the experimental studies on a porous material having representative physiochemical properties of agricultural fields. This type of study will be useful to understand the problem of the transport of fertilizers and salts in the subsurface so that the risk of groundwater quality can be minimized. The physical properties of the soil material were obtained by performing experiments in the geotechnical lab. During the experiments, the test material was sieved using standard mechanical sieves to obtain the particle size distribution of the soil material. The graphical plot on a semi-logarithmic scale, as shown in Figure 13, represents the variation of grain sizes present in the soil sample, in mm, with respect to the percentage finer by weight. The material characteristic properties were studied by determining the value of the coefficient of uniformity (C_U) and coefficient of curvature (C_C) of the test sample. Table 2 shows the characteristic values of the coefficients along with the median grain sizes. The soil was classified using the Indian standard soil classification system, and the material was found to be poorly graded sandy soil, falling under the SP category.

The layering of sand in the tank model was done by compacting the sand in three equivalent layers at optimum moisture content to a depth of about 2 feet (0.6 m) from the invert level. The packaging of the sandy layers is done at a maximum dry density so as to minimize the voids present in the porous medium. The bulk density and dry density of the sample as determined by the sand replacement method were found to be 1.85 g/cm³ and 1.54 g/cm³, respectively. The pore water velocity is estimated to be 0.86 cm/min. The hydraulic conductivity of 1.12 × 10⁻³ cm/min was measured by using the falling head method. The porosity of the medium was determined to be 0.39.

The numerical model has been used to simulate observed breakthrough curves for non-reactive as chloride and reactive solute as fluoride through porous media. Experimental data of observed breakthrough curves have been predicted at 100 cm, 200 cm, 300 cm, 400 cm, and 500 cm down gradient direction in the flow direction.

Experimental data of observed breakthrough curves are shown for discharge flux rate equal to 1.248 cm/min and pore water velocity equal to 3.2 cm/min. During experiments, constant flux was maintained by keeping the constant head at the inlet and outlet of the tank and it was verified by maintaining a constant flow at the outlet ensuring no leakage in the tank. All the breakthrough curves are simulated well predicted in the flow direction and are shown in Figure 14(a)–14(e). From the above analysis, it is shown that the spreading and peak of breakthrough curves are higher for smaller distances, and the magnitude of the peak of breakthrough curves is reduced at a large travel distance in the flow direction. Estimation of the correlation coefficient, root mean square error, and NSE during simulation of observed breakthrough curves at a different distance in the flow direction is shown in Table 3.

Figure 15(a) shows the simulation of observed breakthrough curves for fluoride predicted at 500 cm down gradient distance in the flow direction. Input parameters; that is, pore water velocity v = 3.2 cm, longitudinal hydrodynamic dispersion coefficient D_L = 8 cm²/min, and transverse dispersion coefficient D_T = 0.8 cm²/min, are used during simulation. Afterward, different values of retardation coefficient, R = 1, 1.0824a, and 1.257, are used. It is seen that the value of the retardation coefficient equal to 1.0824 gives the best simulation of the observed breakthrough curve. This value of retardation coefficient equal to 1.0824 is used to simulate the observed breakthrough curves at 100 cm, 200 cm, 300 cm, and 400 cm down gradient distance in the flow direction. Estimation of the correlation coefficient, root mean square error, and NSE during simulation of observed breakthrough curves at a different distance in the flow direction is shown in Table 4.

Based on the above study, the following conclusions are mentioned below:

• 1.

  Experimental data of spatial moments for bromide and chloride chemicals (Freyberg 1986) are simulated well using an asymptotic time-dependent dispersion model. The higher value of dispersion coefficient leads to a decrease in the value of mean travel distance and increases the value of variance during large transport time. The behaviour of variance is non-linear during small transport time in the presence of a higher value of dispersion coefficient.

• 2.

  Variation of longitudinal dispersion coefficients has been predicted in the presence of linear and asymptotic constants. The value of the dispersion coefficient decreases with an increase in the values of both linear and asymptotic constants. However, the behaviour of the dispersion coefficient is non-linear during early transport time, and its magnitude remains constant for large transport time. It is also shown that magnitude of mean travel distance and variance of the solute plume is decreasing with an increase in the value retardation coefficient. It means that the sorption of reactive solute retards the movement of the solute plume in the flow direction.

• 3.

  Experimentally observed breakthrough curves of chloride and fluoride solutes were predicted and simulated using a numerical model. Thus, estimated parameters such as a dispersion coefficient equal to 8 cm²/min and retardation factor equal to 1.0824 are found through numerical simulation of observed data of breakthrough curves. The movement of conservative chloride solute within the porous medium is faster than reactive fluoride due to sorption with soil media. Finally, this study shall be useful to understand the problem of the transport of fertilizers and salts in the subsurface so that the risk to groundwater quality can be minimized.

The authors declared that there is no conflict of interest.

All relevant data are included in the paper or its Supplementary Information.