1. Introduction

The migration of solutes in a porous medium is usually controlled by three mechanisms: Convection, molecular diffusion and mechanical dispersion [1] [2] . During transport of a soluble contaminant from soil surface to groundwater, mechanical dispersion is the major migration mechanism [3] . It is well known that the dispersion of a pollutant in a soil, saturated or not, is influenced both by the intrinsic properties of the soil (such as its nature, porosity, tortuosity and grain size distribution), and by the properties of the fluid flow and transportation (such as viscosity, density and velocity of the fluid) [4] [5] . However, independent of these parameters, it is important in the study of the dispersion process to take into account the nature and origin of the pollutant. Whether organic, inorganic, mineral or heavy metals, domestic, industrial, or pesticides, pollutants involve a large number of phenomena and reactions that may significantly influence their migration in soils [6] [7] .

Major phenomena or reactions encountered during migration of leachates in soil, include: Solubilisation at acidic pH and the precipitation at a basic pH [8] , metal complexation by organic matter [9] [10] , biodegradation of organic matter [11] [12] , adsorption of ions or molecules on the surface of the soil grains [13] and cationic exchange phenomena [14] . Thus, to effectively prevent contamination of groundwater by leachate, it is important to define models which can combine the mechanisms of the fluid flow and transport in soils to the mechanisms related to the pollutants nature.

In this paper, a risk assessment tool was developed based on a model coupling the convection-dispersion to the biodegradation of organic pollutants during leachate migration processes in soils for the prediction of groundwater contamination in the context of developing countries.

2. Background Concepts

2.1. Convection-Dispersion Model

Generally, the change in concentration of a passive pollutant A between two points of a soil column is considered as a one-dimensional transport problem, which is simultaneously subjected to the phenomena of convection and dispersion [15] [16] . This gives:

$\frac{\partial c_A}{\partial t}=D_z\frac{{\partial }^2{c}_A}{\partial z^2}-q_z\frac{\partial c_A}{\partial z}+S\left(c_A\right)$ (2.1)

where: $c_A$ is the concentration of the pollutant $A\left(\text{M}\cdot {\text{L}}^{-3}\right)$ , $D_z$ is the overall coefficient of diffusion-dispersion $\left({\text{L}}^2\cdot {\text{T}}^{-1}\right)$ , is the Darcy velocity, $\left(\text{L}\cdot {\text{T}}^{-1}\right)$ , z and t are the spatial (L) and time (T) variables, and $S\left(c_A\right)$ describes the set of reactions and contributions that take place in the porous medium $\left(\text{M}\cdot {\text{L}}^{-3}\right)$ also called source.

Notice that: $c_A=c_{A_p}+c_{A_s}$ ; with:

$\frac{\partial c_{A_p}}{\partial t}=\epsilon \frac{\partial c_A}{\partial t}$ , change in concentration of A in the soil void; (2.2)

$\frac{\partial c_{A_s}}{\partial t}=\left(1-\epsilon \right)\ast {\rho }_s\ast K_d\frac{\partial c_A}{\partial t}$ , change in concentration of A at soil grain surface.(2.3)

where $\epsilon$ is the porosity of the soil expressed in %, ${\rho }_s$ is the density of the solid matrix $\left(\text{M}\cdot {\text{L}}^{-3}\right)$ , and $K_d$ is the exchange capacity between the fluid and the solid matrix.

2.2. Coupling Biodegradation and Migration Processes

Assuming that the major reaction taking place in the porous medium is a first- order biodegradation reaction of organic matter by oxidation of organic compounds. The source term of (2.1) is written as follows: $S\left(c_A\right)=\frac{\partial c_A}{\partial t}=-{\lambda }_D\cdot c_A$

[17] . Where ${\lambda }_D$ is the degradation rate of the pollutant A. Therefore, Equation (2.1) becomes:

$\frac{\partial c_A}{\partial t}=\frac{D_z}{\epsilon R}\frac{{\partial }^2{c}_A}{\partial z^2}-\frac{v_z}{\epsilon R}\frac{\partial c_A}{\partial z}-\frac{{\lambda }_D}{\epsilon R}\cdot c_A;\text{with}R=1+\frac{\left(1-\epsilon \right)}{\epsilon }\ast {\rho }_s\ast K_d$ (2.4)

where R is the retardation coefficient reflecting the slowing of studied pollutant; it depends on the affinity that has the pollutant to the solid matrix.

The Darcy velocity in Equation (2.4) assumes that the pollutant transport takes place in an empty column. However, the column is filled with a soil of porosity. Therefore, the Darcy velocity is replaced with the interstitial velocity: $u_z=q_z/\epsilon$ , and then the effective dispersion coefficient appears and is expressed as: $D=D_z/\epsilon$ .

For solute transfer in soil, the diffusivity is usually smaller than the dispersion [2] . For a small ratio of the diameter to the length of a soil column (d/L) and for a large fluid velocity, the radial dispersion may be neglected in comparison with

the axial dispersion [18] [19] . Thus, with: $\lambda =\frac{{\lambda }_D}{\epsilon R}$ , the migration mechanisms

of a dissolved biodegradable organic substance in the granular porous medium is expressed by:

$\frac{\partial c_A}{\partial t}=\frac{D}{R}\frac{{\partial }^2{c}_A}{\partial z^2}-\frac{u_z}{R}\frac{\partial c_A}{\partial z}-\lambda \cdot c_A$ (2.5)

Equation (2.5) expresses the coupling “biodegradation - convection dispersion”, taking place during the migration of soluble contaminants in soil. The finite difference method or the finite volume solves this equation.

3. Materials and Methods

Model calibration and verification often require some selected initial values and certain key parameter values. In this study, experimental tests of injections of leachate on soil column were performed according to the protocol for measuring the axial dispersion [4] . Measured data were used to set up inputs to drive the model and provide a basis for adjusting model parameters, in particular dispersion, D.

3.1. Inputs Values: Experimental Tests

The columns of sampled soils were monolithic type. Identification of soil sampling location was made according to ISO 15175:2004 standard, applied to all dumpsites of the municipality of Abomey-Calavi in Benin. The studied soil is a loamy soil. The geotechnical characteristics of the soil of the samples are shown in Table 1.

The characteristics of the leachate used for infiltration tests are presented in Table 2. It was a synthetic leachate obtained after leach tests performed by method in batch (cf. [20] ) on the household waste from the main dumpsites of the municipality of Abomey-Calavi. During the tests, electrical conductivity, BOD₅ and TKN were measured. These data were used to trace the breakthrough curves of these parameters. The biodegradation rate of dissolved organic carbon was determined according to the method presented in [21] and [22] : λ = 0.000037 d⁻¹.

3.2. Model Discretization

The discretization of Equation (2.5) was made here by the finite difference method on uniform mesh. Consider that along a vertical soil column with sufficiently large length L, a liquid with initial concentration $C_0$ flows with interstitial velocity $u_z$ . At a moment $t_0$ , a pollutant of concentration C is injected at the upper end of the column and migrates in the axial direction towards the lower end. It is convenient to define the appropriate boundary conditions (Dirichlet conditions here):

$C\left(0,t\right)=C;C\left(\infty ,t\right)=C_0;C\left(z,0\right)=C_0.$ (3.2)

The differentiation schemes used are:

$\{\begin{array}{l}{\left(\frac{\partial C}{\partial t}\right)}_i^j=\frac{C_i^{j+1}-C_i^j}{\text{Δ}t}\\ {\left(\frac{{\partial }^{2}C}{\partial z^2}\right)}_i^j=\frac{C_{i+1}^j-2C_i^j+C_{i-1}^j}{\text{Δ}{z}^2}\\ {\left(\frac{\partial C}{\partial z}\right)}_i^j=\frac{C_{i+1}^j-C_{i-1}^j}{2\text{Δ}z}\end{array}$ (3.3)

where: L is divided into N intervals with ends or nodes $z_i$ , i ranges from 1 to N + 1; $C_i^j$ is the concentration at the node $z_i=i\text{Δ}z$ at $t=j\text{Δ}t$ . Thus Equation (2.5) is equivalent to the set of:

$\begin{array}{l}{C}_i^{j+1}=\left(\alpha +\beta \right)C_{i-1}^j+\left(1-2\alpha -\gamma \right)C_i^j+\left(\alpha -\beta \right)C_{i+1}^j\\ \text{with}\{\begin{array}{l}\alpha =\frac{D\Delta t}{R\Delta z^2}\\ \beta =\frac{u_z\Delta t}{2R\Delta z}\\ \gamma =\lambda \Delta t\end{array}\end{array}$ (3.4)

i ranges from 1 to N − 1.

The system (3.4) has been implemented by the successive over-relaxation method. The calculations were performed with Matlab R2013a. The initial values were: length of the column (𝐿), time of simulation (𝑇), space step (∆z), time step (∆𝑡), flow velocity ( $u_z$ ), retardation factor (𝑅), initial concentration (C) and a starting value for axial dispersion coefficient (D).

The retardation factor R of the pollutant over a water molecule is defined by the quotient of the residence time of the pollutant ( $t_s$ ) over the residence time of the water molecule [23] . In this study, where the flow permanent and uniform flow, R is determined by Equation (3.5). The residence times of pollutants and water molecule were measured experimentally by injections of leachate on soil columns. Thus, the mean value for the retardation coefficient was calculated and is equal to 2.5.

$R=\frac{t_s}{L\theta }\cdot v_z$

where: L is the length of the column; θ is water content in the column and $v_z$ the Darcy velocity in column input [24] [25]

3.3. Model Calibration

Figure 1 shows the algorithm to solve Equation (2.5) following the method of finite differences. The program performed an adjustment of the given initial value for axial dispersion coefficient in order to optimize the comparison of model results to measured data. The initial value of the axial dispersion coefficient for the implementation was chosen in accordance with the values found in the literature for the dispersion in loamy soils and was equal to 0.048 cm²/min [26] .

To quantify the model’s prediction, the Relative Root-mean-square Error (RRE) was used to compare simulated versus observed values, with the best fitting simulation returning the lowest RRE.

3.4. Model Verification

The dispersion coefficient obtained after model calibration was used as input for the verification of the model (Figure 1). The result obtained after simulation was compared with the measured data. Model was approved when the RRE obtained after verification was close to the RRE for model calibration.

4. Results and Discussion

The measured data were separated into two subsets: one subset for model calibration and another subset for model validation. 48 infiltration tests were carried out on soil columns A of 30 cm length. The subset for the model calibration encompassed the electrical conductivity (EC) data of the columns A. The simulated parameters for model verification include BOD₅ and TKN, registered at the outlet of columns A.

4.1. Infiltration Rate

The influence of the hydraulic gradient on the leachate infiltration rate was

shown in Figure 2. This figure showed that the evolution of the infiltration speed according to the hydraulic gradient had two phases: The first branch where the infiltration rate increased rapidly for small changes of the hydraulic gradient and the second branch where the infiltration rate was less sensitive to the change of the hydraulic gradient.

Note that for a hydraulic gradient of 8.8, the actual infiltration rate measured when characterizing the soil ( $7.55\times {10}^{-4}\text{cm}/\text{s}$ ) was close to the one obtained experimentally ( $8.07\times {10}^{-4}\text{cm}/\text{s}$ ). Thus, the rate of flow of leachate in soils for subsequent simulations was taken as $8\times {10}^{-4}\text{cm}/\text{s}$ .

4.2. Model Calibration

The time of simulation was computed to simulate one-month leachate migration during the rainy season. A time step of 1 second (s) was used throughout the simulations. The simulated Breakthrough curves and the experimentally-derived breakthrough curve for model calibration are shown in Figure 3.

A comparison of model results to measured data, illustrated by RRE (Table 3), showed that simulation 1 (Sim1 D = 0.048), performed with the initial value of dispersion coefficient derived from the literature $\left(D=0.048{\text{cm}}^2/\min \right)$ , gave an unsatisfactory fit (RRE = 40.96). Gradually, as the value of the dispersion coefficient was incrementally increased, the simulation curve better approximated the experimental curve. Simulation 3 (Sim3 D = 0.96) provided the optimum fit, with a RRE of 11.22% between model results and measured data. However, the value of the dispersion obtained at the end of the simulation (2378 µS/cm) is less than the experimental value (3174 µS/cm). Note that the experimental curve had a first concavity due to change of solute concentration in effluent at the column outlet at early time between 275 and 600 min. This rapid change of concentration is followed by a bearing between 650 and 890 min. This pattern was probably due to preferential paths in soil columns. Indeed, the preferential flow paths are reflected in very rapid change of solute concentration in effluent at early times. The monolithic soil columns (as the one used in this study) generally have macropores which can promote preferential flows [3] . Thus, before a generalization of these results can be made, the status of the macropores distribution in the soil columns should be investigated. For this case study, the average proportion of macro-pores inside the soil columns was 21%. A correction of about 21% of the experimental result approached the simulated (predicted) value. The calibrated value of the dispersion coefficient for this study was $D=0.96{\text{cm}}^2/\min$ .

4.3. Model Verification

For the model verification, a longer time step was selected. This time was calculated to correspond to a simulation of the amount of leachate produced during the long rainy season (from March to June) in southern Benin. The verification process involves running the model with the calibrated dispersion coefficient $D=0.96{\text{cm}}^2/\min$ and comparing the results to the independent data set for model verification.

4.3.1. Organic Matter (BOD₅) Migration Model

In the numerical simulation of BOD₅ migration, the biodegradation rate of the studied leachate was considered, which is 0.000037 d⁻¹. This value was taken equal to zero upon insertion of input parameters for the model calibration with the electrical conductivity. Figure 4 showed the results obtained for experimental investigation and numerical simulations with dispersion coefficient $D=0.96{\text{cm}}^2/\min$ . The RRE between model results and measured data was 11.14%.

As in the case of electrical conductivity data for model calibration, a bearing was clearly observed on the measured BOD₅ data for model verification between 660 and 990 min. The end result for BOD₅ simulation (4197 mg O₂/l) was less than measured outcome (4914 mg O₂/l), with a simulation result correction of 21% again approximating the experimental result. The calibrated dispersion coefficient (D) adequately modelled the diffusion of organic pollution for the studied soil. The verification of the migration model for the transport of organic matter in soil (BOD₅) gave satisfactory results.

4.3.2. Nitrogenous Organic Matter (TKN) Migration Model

Figure 5 shows the simulated and measured breakthrough curves of TKN in soil

columns. As can be seen, the calibrated dispersion coefficient did not reflect the experimental evolution of the nitrogenous organic matter, and ammonia/ ammonium in the soil. The verification of the migration of TKN in soil columns A was not conclusive.

The analysis of the breakthrough curve of measured TKN showed a reduction of 99.65% of the nitrogenous organic matter at the end of the experimentation, which was above the recorded reductions for electrical conductivity (10.69%) and BOD₅ (15.69%). Usually, TKN is mostly made up of biodegradable organic compounds and therefore, a similar behavior to BOD would be expected. The failure of the validation of the TKN migration model could come from the difference recorded in the reductions of TKN and BOD. A low ratio of organic nitrogen to NH₃/ ${\text{NH}}_4^{+}$ could explain this difference in the reductions. Therefore it can be inferred that the retardation phenomena and biodegradation of organic matter reactions are not sufficient to model the migration of TKN in soil. The model should look at the different components of TKN separately. In another hand, since the retardation coefficient R was calculated based on EC measurements, it might be that the model works for large humic macromolecule (forming part of BOD, for instance), rather than for much smaller molecules as ammonia (in TKN).

At this stage of the study, the designed model properly evaluate the effect of retardation factor on the migration of macromolecules contained in leachate. For much smaller molecules as ammonia, a recalibration of the retardation factor seem to be necessary. A proper model of migration of biodegradable organic matter in the underlying soils would then require a coupling of two retardation factors: one for the macromolecules contained in the leachate and another one for the smaller monovalent ion (such as Na⁺ in EC) or ammonia,

4.4. Influence of Measured Input Parameters

4.4.1. Influence of the Retardation Factor on Organic Matter Migration Model

Retardation coefficient (R) represents the delay accused by a pollutant molecule with respect to the water molecule introduced at the same time at the inlet of a soil column. The value for the retardation factor coefficient taken for model calibration and verification in this study was 2.5. Figure 6 showed simulated breakthrough curves of BOD₅ migration for different values of retardation factor.

The analysis of the curves in Figure 6 revealed that the migration of BOD₅ in soil columns was slowed by the increased retardation factor. The simulation of the dispersion of organic matter with $R=2.5$ was the one that best approximates the measured results. Therefore the method used to measure the retardation coefficient is effective. For a non-delayed migration of pollutants, expressed by a retardation factor equal to 1, the simulated results are far removed from experimental result. The retardation factor is an essential parameter in modelling the migration of biodegradable organic matter in soils of waste dumpsites.

The influence of the retardation factor in the migration of biodegradable organic pollutants in soils was linear, and expressed in the present experimental study by the following equation: (Figure 7)

$y=-0.13x+1.05$ (4.1)

4.4.2. Influence of Biodegradation Rate on Organic Matter Migration Model

The biodegradation rate expresses the speed of consumption of organic pollution by aerobic and anaerobic bacteria. Numerical simulations (Figure 8) showed that an assumption of zero biodegradation (BC = 0) better approximated the experimental result than prediction with BC = 0.000037j − 1, as used for BOD₅ model calibration. This can be explained by the fact that the organic matter biodegradation kinetics is an exponential function, thus, at the resolution of the experimentation, it is difficult to observe the effects of the biodegradion for short times of 1 and 2 days. An increase of around 10⁻⁴ of biodegradation rate causes a halving of organic pollution at the outlet of the soil columns. It could be

the interesting to test longer times/longer columns in order to properly evaluate the biodegradation effect. For a practical point of view, it would be then beneficial to proceed to a bacteria activation of the soil before the deposit of household waste.

The influence of the rate of biodegradation on migration of biodegradable organic pollutants in our experimental investigation (Figure 9) was expressed by the equation:

$y=-0.76\cdot {\text{e}}^{-1688x}$ (4.2)

5. Conclusion

The phenomenon of natural biodegradation of organic matter was combined with the principle of convection-dispersion to model migration of leachates in soil. The resolution of this model by the finite difference method showed that the migration in soil of macromolecules contained in leachate can be predicted based on retardation factor determine by infiltration test on soil columns. For much smaller molecules as ammonia, a recalibration of the retardation factor is necessary. Therefore a model that aims to predict the migration in soil of both macro and micromolecule of leachate seem to require a coupling of two type of retardation factor: one for macromolecule such as BOD and one for micro molecule such as ammonia. The model at this stage failed to predict the biodegradation during of organic matter during the migration of leachate in the soil. Longer test times and longer test columns are required in order to properly evaluate the biodegradation effect of leachate during it migration in soils. Increasing the bacterial activity of the underlying soils of garbage dumpsites could be a solution for natural reduction of pollutants migrating to groundwater.

Acknowledgements

The authors gratefully acknowledge the International Foundation for Science (IFS), Stockholm, for supporting the present work under the IFS Grant W/5840-1.

Conflict of Interest Disclosure

The author(s) declare(s) that there is no conflict of interest regarding the publication of this paper.