Opinion: Agriculture Green Development

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

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

2.1 Advection-dispersion equation

The traditional ADE model with a reactive term^[ 12^] can be described as:

(1) $${\displaystyle \frac{\partial C^j}{\partial t}}=D_L{\displaystyle \frac{{\partial }^{\mathrm{2}}{C}^j}{\partial z^{\mathrm{2}}}}-v{\displaystyle \frac{\partial C^j}{\partial z}}+r^j$$

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²·s^-1), D_L = D_e + a_L·v, D_e is the effective diffusion coefficient (m²·s^-1), a_L is the dispersivity (m) and r^j is the source/sink term of species j (mol·L^-1·s^-1).

2.2 Mobile-immobile model

In a one-dimension steady flow system, multi-species MIM^[ 13^] with a reactive term is given by:

(2) $${\theta }_m{\displaystyle \frac{\partial C_m^j}{\partial t}}+{\theta }_{im}{\displaystyle \frac{\partial C_{im}^j}{\partial t}}={\theta }_m{D}_L{\displaystyle \frac{{\partial }^{\mathrm{2}}{C}_m^j}{\partial z^{\mathrm{2}}}}-v_m{\theta }_m{\displaystyle \frac{\partial C_m^j}{\partial z}}+{\theta }_m{r}_m^j$$

(3) $${\theta }_{im}{\displaystyle \frac{\partial C_{im}^j}{\partial t}}=\alpha (C_m^j-C_{im}^j)$$

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³·cm^-3), q_m + q_im = q, q is the total water content (cm³·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²·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).

2.3 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:

(4) $$R=\mathrm{1}+{{\rho }_{\mathrm{b}}}_{}{\displaystyle \frac{k_d}{\theta }}$$

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:

(5) $$r=(\mathrm{1}-R){\displaystyle \frac{\partial C}{\partial t}}$$

The first-order degradation of a degradable solute is given by:

(6) $$r=-k{\displaystyle \frac{dC}{dt}}$$

where k is the first-order decay rate constant (d^-1).

2.4 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:

(7) $$Q=v\cdot \theta \cdot A=v_m\cdot {\theta }_m\cdot A$$

where A is the area of cross section (m²).

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^].

Tab.1 Input parameters in base case analysis for ADE and MIM model comparison

Parameter	value
Length of column/m	1
Dispersivity in ADE and MIM aL/m	0.01
Inflow rate q/(m·d-1)	0.028
Water content q/(cm3·cm-3)	0.45
Pore water velocity in ADE v/(m·d-1)	0.062
Mobile fraction b	0.78
Pore water velocity in MIM vm /(m·d-1)	0.08
Mass transfer rate a/s-1	1 × 10-6
Retardation factor R	2
Decay rate coefficient k/d-1	0.1
Input concentration /(mol·L-1)	0.1
Initial mobile/Immobile concentration/(mol·L-1)	0/0

To give a more direct comparison of the outflow concentration, the curves of relative concentration of each solute C/C₀ versus pore volume (PV) were shown, where C is the outflow concentration of each solute and C₀ is the input concentration of each solute.

3 Results and discussion

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

Fig.1 Breakthrough curves for conservative (CS), adsorbing (AS) and degradable (DS) solutes using ADE and MIM modeling

Full size|PPT slide

3.2 Sensitivity analysis

3.2.1 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 Figs. 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^].

Fig.2 Breakthrough curves with different dispersivities aL for conservative solute (CS) simulated by ADE (a), conservative solute simulated by MIM (b), adsorbing solute (AS) simulated by ADE (c), adsorbing solute simulated by MIM (d), degradable solute (DS) simulated by ADE (e), and degradable solute simulated by MIM (f)

Full size|PPT slide

As shown in Figs. 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 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.

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

Fig.3 Breakthrough curves for the degradable solute (DS) simulated by ADE and MIM with different degradation coefficients k

Full size|PPT slide

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

Fig.4 Breakthrough curves for the adsorbing solute (AS) simulated by ADE (a) and MIM (b) with different retardation factors R (as indicated in the legends)

Full size|PPT slide

3.2.4 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^].

Fig.5 Breakthrough curves for conservative solute (CS) (a), adsorbing solute (AS) (b) and degradable solute (DS) (c) simulated by ADE and MIM with different mobile water fraction b (as indicated in the legends)

Full size|PPT slide

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^].

3.2.5 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 different values of mass transfer coefficient a. The BTCs of conservative and adsorbing solutes obtained from MIM (q_m = 0.35) nearly overlap with the BTCs obtained from ADE (q = 0.45) when a equals to 1 × 10^-8 s^-1. This result is consistent with the findings of Feehley and Zheng^[ 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.

Fig.6 Breakthrough curves for conservative solute (CS) (a), adsorbing solute (AS) (b), and degradable solute (DS) (c) simulated by ADE and MIM with different mass transfer coefficients a (the values were labeled in legend)

Full size|PPT slide

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.

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

References

[1] Ingram J S I, Porter J R. Plant science and the food security agenda. Nature Plants, 2015, 1(11): 15173

RIGHTS & PERMISSIONS

The Author(s) 2019. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)