Comparison of the performance of traditional advection-dispersion equation and mobile-immobile model for simulating solute transport in heterogeneous soils

The traditional advection-dispersion equation (ADE) and the mobile-immobile model (MIM) are widely used to describe solute transport in heterogeneous porous media. However, the fitness of the two models is casedependent. In this paper, the transport of conservative, adsorbing and degradable solutes through a 1m heterogeneous soil column under steady flow condition was simulated by ADE and MIM, and sensitivity analysis was conducted. Results show that MIM tends to prolong the breakthrough process and decrease peak concentration for all three solutes, and tailing and skewness are more pronounced with increasing dispersivity. Breakthrough curves of the adsorbing solute simulated by MIM are less sensitive to the retardation factor compared with the results simulated by ADE. The breakthrough curves of degradable solute obtained by MIM and ADE nearly overlap with a high degradation rate coefficient, indicating that MIM and ADE perform similarly for simulating degradable solute transport when biochemical degradation prevails over the mass exchange between mobile and immobile zones. The results suggest that the physical significance of dispersivity should be carefully considered when MIM is applied to simulate the degradable solute transport and/or ADE is applied to simulate the adsorbing solute transport in highly dispersive soils.


Introduction
Spatial heterogeneity of soil is a ubiquitous phenomenon across all scales in natural soils, including features such as macropores and layered soils [1] . Transport in a heterogeneous soil system may allow rapid movement of pollutants and create significant consequences for geochemical interactions and groundwater quality. Solute transport models provide an effective way to quantify the process.
The traditional advection-dispersion equation (ADE) and the mobile-immobile model (MIM) are the most often used solute transport models. ADE can fit the observed breakthrough curves (BTCs) well in some heterogeneous soils, although ADE normally shows good fit in modeling solute transport in homogenous porous media. For example, ADE fitted the BTC of conservative tracer well for a layered soil column [2] . ADE was even more suitable than MIM to simulate the tracer transport in the undisturbed, stratified cores from Hanford sediment [3] . ADE also performed better than MIM with consideration of reactions, although the experimental data exhibited considerable heterogeneity in percolation rate and solute concentrations in undisturbed soils [4] . However, in some other studies, dispersivity obtained from fitting ADE was extremely high, which might not be physically meaningful [5].
MIM is widely used to describe transport along a heterogeneous flow path. MIM accounts for rapid transport in a mobile pore region and rate limited solute exchange taking place between the mobile and immobile pore regions. Satisfactory agreement with observed BTC can be obtained by the model, especially the early arrival and long tailing portions. MIM fitted chloride data well in an undisturbed loess soil column containing macropores [6] and a large heterogeneous soil column [7] , even in the vadose zone [8] . A better description of nitrate transport through soil-root channel columns was provided by MIM in comparison to ADE [9] . The agreement was commonly explained as the co-existence of mobile zones (continuous path, such as connected pore space or fractures) and immobile zones (disconnected path, such as rock matrix surrounding fractured media or unconnected pore space).
The applicability of ADE and MIM in simulating experimental results was evaluated by the coincidence between measured and simulated BTCs, although the results were case dependent. MIM and ADE can get similar good fits to the observed BTC in homogeneous media [7] . However, the ADE and MIM-simulated BTCs may not show good agreement with observed data in heterogeneous fields [10] , but there is no clear conclusion about which model performs best.
In this study, systematic evaluation of model feasibility for ADE and MIM was performed. The aims of the study were (1) to compare BTCs of conservative, adsorbing and degradable solutes transport in heterogeneous soils predicted by ADE and MIM, and (2) to test the influence of properties such as dispersivity a L , retardation factor R, and decay rate constant k on the performance of ADE and MIM. The specific parameters in MIM, i.e., mobile water fraction b and mass transfer rate a, were also discussed with respect to the performance of MIM.

Materials and methods
A base case was established to analyze the performance of ADE and MIM. In the base case, three solutes, i.e., conservative, adsorbing and degradable solutes, were assumed to move through a 1 m heterogeneous soil column under steady flow conditions and simulated by PHREEQC-2 [11] . Sensitivity analysis was further conducted to evaluate the effect of physical and chemical properties on ADE and MIM simulation results.

Advection-dispersion equation
The traditional ADE model with a reactive term [12] can be described as: where the j refers to species j, j = 1, 2, 3…; C j is the concentration of species j (mol$L -1 ); t is the time (s); z is the distance (m); v is the average pore velocity across the cross section of soil column (m$s -1 ); D L is the hydrodynamic dispersion coefficient (m 2 $s -1 ), D L = D e + a L $v, D e is the effective diffusion coefficient (m 2 $s -1 ), a L is the dispersivity (m) and r j is the source/sink term of species j (mol$L -1 $s -1 ).

Mobile-immobile model
In a one-dimension steady flow system, multi-species MIM [13] with a reactive term is given by: where the subscripts m and im refer to mobile and immobile liquid phases; q m and q im are the water filled porosity of the mobile and stagnant part (a fraction of total volume), respectively (cm 3 $cm -3 ), q m + q im = q, q is the total water content (cm 3 $cm -3 ); C j m and C j im are the concentrations of species j in mobile and immobile regions (mol$L -1 ); D L is the hydrodynamic dispersion (m 2 $s -1 ), D L = D e + a L $v m , v m is the average pore velocity in the mobile liquid; r j is the source/sink term of species j (mol$L -1 $s -1 ) and a is the first-order mass transfer coefficient (s -1 ).

Reaction terms r j
For equilibrium adsorption, the reaction (adsorption) term can be described by the retardation factor in both ADE and MIM. The retardation factor R for linear adsorption can be described as: where R is the retardation factor; r b is the soil bulk density (g$cm -3 ) and k d is the partition coefficient. For the adsorbing solute, the following relationship exists between r and R: The first-order degradation of a degradable solute is given by: where k is the first-order decay rate constant (d -1 ).

Model configuration
Instantaneous solute input will result in a higher degree of nonequilibrium, compared with continuous solute input [14] . Therefore, two pore volumes of solute were assumed to be injected into the column in the base case, in order to simulate more clearly the different response of the two models. The input concentrations of the three solutes were set the same, 0.1 mol$L -1 . The same volume flow rate was also applied for the two models, so the following relationship exists: where A is the area of cross section (m 2 ). The parameters applied in the base case for both ADE and MIM were listed in Table 1. The values of parameters chosen in the simulation were within the range of reported values [4,7,15] .
To give a more direct comparison of the outflow concentration, the curves of relative concentration of each solute C/C 0 versus pore volume (PV) were shown, where C is the outflow concentration of each solute and C 0 is the input concentration of each solute.

Breakthrough curves of solutes in base case
The BTCs simulated by MIM are flatter than ADE for the three solutes (Fig. 1), especially for the degradable solute. MIM tends to prolong the process of breakthrough and decrease the peak concentration, because the mass transfer from the mobile zone into the immobile zone prolongs the time to breakthrough. This result is in line with the reported result that MIM could capture the skewness and tailing of BTCs compared with ADE [6,13,16] . The BTCs of adsorbing and conservative solutes produced by MIM nearly overlap, suggesting that by using MIM the two solutes behave similarly.

Dispersivity a L
BTCs of the three solutes generated by ADE and MIM exhibit the same trends, i.e., the phenomenon of early breakthrough, tailing and reduced peak concentrations appears to be more prominent with the increase in dispersivity a L (Fig. 2). These phenomena are more obvious for the BTCs simulated by MIM because of mass exchange between mobile and immobile zones.
For conservative and degradable solutes, the BTCs obtained by ADE overlap together when dispersivity exceeds 0.1 m. However, phenomena of early breakthrough and tailing become even more obvious with dispersivity increasing to greater than 0.1 m for the BTCs obtained by MIM. It can be seen in Fig. 2b, 2d, 2f that the three solutes leach out from near the beginning of the test with a L set to 0.5 m and 0.8 m in MIM. Although such high dispersivity relative to the soil column was not used to fit measured BTCs, the result highlights the physical meaning of this parameter in the fitting process, which has also been reported by Comegna et al. [17] .
As shown in Fig. 2b, 2c, the adsorbing solute simulated by ADE behaves similarly to the conservative solute simulated by MIM, although the peak concentration is relatively low for the latter. The similarity of the BTCs is because both adsorption and mass transfer serve as sinks [18] . It demonstrates that the transport with both adsorption and mass transfer (Fig. 2d) can lead to more obvious tailing and lower peak concentration than the transport with mass transfer only (Fig. 2b) or adsorption only (Fig. 2c).
Mass transfer between mobile and immobile zones enhances the asymmetry of BTC for the degradable solute (Fig. 2f). It is notable that the peak concentration of the degradable solute with dispersivity a L set to 0.8 m is Initial mobile/Immobile concentration/(mol$L -1 ) 0/0 Fig. 1 Breakthrough curves for conservative (CS), adsorbing (AS) and degradable (DS) solutes using ADE and MIM modeling slightly higher than when set to 0.01 m. The dispersion process may overwhelm the degradation and the mass transfer with the extremely large dispersivity relative to the column scale, so that the degradable solute in the mobile zone leaches out before thorough reaction and mixing with the solute in the immobile zone. This result indicates that caution should be taken when modeling degradable solute transport by MIM in highly dispersive soil conditions.

Degradation coefficient k
The simulation results show that the BTCs obtained from MIM are more asymmetric than ADE with the same degradation rate coefficient k, resulting from the combined effect of degradation and mass transfer in MIM. The peak concentration obtained from MIM is lower than ADE at relatively small k (0.05 and 0.1 d -1 ), while there is little difference between the peak concentrations and the shapes of BTCs with k set to 0.3 and 0.5 d -1 because of overwhelming degradation (Fig. 3). The results imply that MIM and ADE perform similarly in simulating degradable solute transport when degradation dominates mass transfer in porous media. ADE is preferable in this case because there are fewer parameters to be determined.

Retardation factor R
The BTCs of adsorbing solute simulated by ADE and MIM are flattened as R increases (Fig. 4). For different R settings in the two models, the peak concentrations of adsorbing solute obtained by MIM are closer to each other, compared with ADE simulated results. It seems that the adsorbing solute BTC obtained by MIM is relatively less sensitive to retardation factor R. This is because both adsorption and mass transfer serve as sinks for the adsorbing solute in MIM. These results show that the observed BTC of an adsorbing solute can be fitted by either ADE with a high R (shown as AS ADE R = 4 in Fig. 4a) or MIM with a low R (shown as CS MIM R = 1 in Fig. 4b). Although ADE can describe the overall effect by fitting with relatively high a L and R, again the physical meaning should be considered.

Mobile water fraction b
As an important parameter to capture BTC shape in MIM, mobile water fraction b decreasing can result in apparent skewing of BTCs for the three solutes (Fig. 5). The results show that a greater ratio of immobile water extends the time to reach equilibrium between the mobile and immobile zones, so the corresponding BTC is broader. The leaching-out concentration decreases with b decreasing, suggesting that the immobile part functions as a storage zone or sink. Immobile concentration was set to 0 in the simulation. In this case, smaller b means larger storage capacity for solutes. The immobile zone may also function as a source after long-term leaching when the concentration of input solution is relatively low. For example, stored phosphate can easily leach out from nutritious soils when a low concentration phosphate solution or fresh water flows through the soil [19] . The mobile-immobile system acts like a single-domain model with b approaching 1. As shown in Fig. 5, for each of the three solutes, the shape of BTC produced by MIM with b equal to 0.78 is more similar to the BTC produced by ADE compared with the MIM-simulated results with lower b (0.22, 0.44 and 0.55). This is because higher volumes of mobile water accelerate the solute movement and increase the quantity of contaminant transport [20] .

Mass transfer coefficient a
The exchange between mobile and immobile pores is described by the first-order mixing rate coefficient in MIM. It can be seen in Figs. 6a, 6b that the BTC of conservative solute behaves similarly to the adsorbing solute with  [21] that the system behaves more like a single-domain system with the total porosity approaching the porosity of q m when the mass transfer rate is low.
For the degradable solute, the peak concentrations of MIM-simulated BTCs with a equal to 1 Â 10 -7 and 1 Â 10 -8 s -1 are even higher than those obtained from ADE (Fig. 6c). This is because solute transport is faster along the  mobile flow path than the uniform flow path, i.e., v m is higher than v. Thus, the degradable solute underwent incomplete degradation in the two region system compared with the single-domain system, under the same volume flow rate conditions. The result is consistent with the related study that solute transport along preferred pathways results in a rapid movement through the biologically active zone and thereby slower degradation rates of chemicals [22] . In this case, MIM-simulated solute fate along preferential flow paths can provide assistance for risk assessment of soil and groundwater contamination.
With a increasing to 1 Â 10 -6 s -1 , the shapes of BTCs for the three solutes are broader and peak concentrations lower. This is because the relatively high transfer rate coefficient improves the exchange between mobile and immobile zones, and consequently some solutes are temporarily stored in the immobile zone. With a equal to 1 Â 10 -5 s -1 , the peak concentration is higher than with a equal to 1 Â 10 -6 s -1 . This result does not appear intuitive as one would expect the peak concentration to decrease with a increasing. A probable interpretation is that the mass concentrations in mobile and immobile domains tend to achieve equilibrium quickly under the high transfer rate condition, and the simulated system will behave like a single porosity domain [23] . As the exchange process continues, the mass stored in the immobile zone in the early stage could be exchanged out, so the peak concentration could increase but the breakthrough process be prolonged.

Conclusions
To investigate the performance of ADE and MIM in simulating solute transport through heterogeneous soils, transport of three typical solutes, i.e., conservative, adsorbing and degradable solutes, through a 1 m soil column with pulse input was simulated using the two models. Results show that MIM tends to prolong the processes of breakthrough and decrease peak concentrations for the three solutes, which is consistent with previous research. BTCs of the adsorbing solute modeled by MIM show less sensitivity to the retardation factor R in comparison to the simulation results of ADE. The observed BTC of an adsorbing solute can be fitted by either MIM with a conservative solute or ADE with an adsorbing solute. ADE performs similarly in simulating degradable solute transport when degradation dominates mass transfer in porous media.
These results, albeit focused on a hypothetical scenario, suggest that dispersivity should be cautiously used in the model, which should not be calibrated without considering its physical reality when MIM is applied to simulate the degradable solute transport and/or ADE is applied to simulate the adsorbing solute transport in highly dispersive soils.