This work, for the first time, developed an unconditionally stable, implicit numerical solution of a one-sided spatial fractional advection-dispersion equation (s-FADE) on a bounded domain to model contaminant transport through porous media during a leaching process. We used the zero-value Robin and Neumann conditions with a fractional-order derivative as the inlet and outlet boundaries, respectively, and a constant initial concentration of the contaminant in the whole domain as an initial condition. By conducting numerical analyses, we found that there is a good agreement between the results of the numerical solution and of its analytical counterpart, validating the accuracy of the numerical solution. In the case of the α variations meaningfully altered the heterogeneity degree of the porous medium, these variations significantly affected the results of the numerical solution, especially at intermediate times. Contrary to the α parameter, the numerical solution was insensitive to the spatial node number. The noisy observational data led to the maximum uncertainty in the estimated fractional dispersion coefficient. The application of the developed numerical solution to simulate sodium chloride (NaCl) leaching curves (LCs) at a laboratory scale revealed that the numerical solution described the leaching process in the homogeneous and heterogeneous sands better than COMSOL Multiphysics®, a commercial finite element software package. According to the fitting results of the developed numerical solution to the experimental data, the mean α values were 1.876 and 1.379 for the homogeneous and heterogeneous sands, respectively. The α values correctly represented the heterogeneity degrees of the used sands and the deviation rate of the transport process from the Fickian transport. Overall, the proposed numerical solution is a reliable numerical model to simulate the leaching process in homogeneous and heterogeneous porous media with bounded domains.
