Co-transport of heavy metals in layered saturated soil: Characteristics and simulation

Co-transport of heavy metals in layered saturated soil: Characteristics and simulation

Journal Pre-proof Co-transport of heavy metals in layered saturated soil: Characteristics and simulation Qing Lin, Shaohui Xu PII: S0269-7491(19)3440...

1MB Sizes 0 Downloads 7 Views

Journal Pre-proof Co-transport of heavy metals in layered saturated soil: Characteristics and simulation Qing Lin, Shaohui Xu PII:

S0269-7491(19)34402-1

DOI:

https://doi.org/10.1016/j.envpol.2020.114072

Reference:

ENPO 114072

To appear in:

Environmental Pollution

Received Date: 29 August 2019 Revised Date:

16 January 2020

Accepted Date: 23 January 2020

Please cite this article as: Lin, Q., Xu, S., Co-transport of heavy metals in layered saturated soil: Characteristics and simulation, Environmental Pollution (2020), doi: https://doi.org/10.1016/ j.envpol.2020.114072. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

1

Co-transport of heavy metals in layered saturated soil:

2

characteristics and simulation

3 4

Qing Lin, Shaohui Xu

5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

Department of Environmental Science and Engineering, Qingdao University, Qingdao 266071, PR China

Abstract Interest in soil pollution by multiple heavy metals has been growing over the last decades. However, few experiments combining numerical analyses with solute transport in layered soil can be found in the literature. Here, the retention and fate of three coexisting metal ions, Cu, Cd, and Zn, in layered soils were investigated to evaluate soil co-contamination through batch and column experiments. Results showed high amounts of Cu adsorption and retention by soils, followed by Cd and Zn. The partial concentration of Zn in effluent was greater than the input from competition adsorption and the ‘snow plow effect’. These findings indicate the high potential risk of Zn and Cd groundwater pollution when Cu, Cd, and Zn co-exist in the soil. Adsorption isotherms obtained from batch experiments were well described by Freundlich equation. Breakthrough curves (BTCs) obtained from column experiments were well described by standard convection-dispersion equation (CDE) for Br, and Tow-site (TSM) and One-site models (OSM) for metals except for Zn, using the Levenberg-Marquardt nonlinear optimization algorithm. However, the parameters were poorly constrained by the available observational data due to high correlation between parameters, rather than insensitivity to model outputs. The Generalized Likelihood Uncertainty Estimation (GLUE) method did not only qualify the uncertainty of parameters for solute transport in layered medium, but estimate prediction uncertainty. Prediction bounds basically captured the observed Br, Zn and Cd BTCs, while systematically overestimated the effluent Cu concentration. Comparing with the optimization, GLUE method can improve prediction reliability of heavy metal transport in layered soils. Keywords: Heavy metal; Cotransport; Layered soil; Predictive uncertainty; GLUE

28

1. Introduction

29 30 31 32 33 34 35 36 37 38 39 40 41

Industrial and agricultural activities have released large amounts of heavy metals into soil, which has become a serious and widespread problem globally (Sterckeman et al., 2000; Waterlot et al., 2013; Li et al., 2014; Yang et al., 2018). This contamination may eventually end up in the food chain or infiltrate groundwater resources through its mobility, persistence, and toxicity, thereby posing a significant threat to public health (Rattan et al., 2005; Zhang et al., 2015; Zhang et al., 2018). The retention and mobility of heavy metals in soil is affected by many factors, such as the properties of the heavy metals, physical-chemical properties of the soil, and climatic conditions. Column experiments have been conducted to explain the transport of heavy metals in porous media, but most of these studies have focused on a single metal or coexisting heavy metals transport in homogeneous soil (Pang et al., 2002; Lafuente et al., 2008; Fonseca, et al., 2011; Chotpantarat et al., 2012; Garrido-Rodriguez et al., 2014). It has been demonstrated that the coexisting metals may compete with each other for soil sorption sites, resulting in enhanced transport of heavy metals in the porous media. For example, Chotpantarat et al. (2012) reported

42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85

that for a soil column with lateritic soils, the simultaneous presence of other metals in multi-metal systems reduced sorption through competition for sorption sites on the solid phases and moreover the sorption of Mn, Zn and Ni were more strongly affected by the simultaneous presence of a competing metal than Pb was. However, most natural soil appear as a staggered layered structure, forming a complex and variable soil profile, and also the artificial effects of irrigation and farming modify field soil to form a dense layer, a sedimentary layer or a plow layer. Each soil profile layer has its own unique physical and chemical properties, such as pH, CEC, bulk density, and organic content, resulting in different retention of heavy metals, however, few studies have addressed solute transport under these conditions. Seuntjens (2002) studied the field-scale migration of Cd in a 1 m deep layered sandy soil profile and showed that chemical heterogeneity largely affects the position and the shape of the Cd plume, resulting in an early Cd breakthrough at the bottom of the soil profile. Sana and Jalila (2018) discussed the effectiveness of layered soil (two and three-layer capillary barriers) in the reduction of organic pollutant transport through an experimental and numerical model. Zhan et al. (2018) investigated the migration mechanism of Cu ions in the soil column composed of three layers and found that the retardation rate for Cu was as high as up to 91%. The reactive cotransport behaviors of heavy metals in layered soils are more complex than single heavy metal or in homogeneous media, might producing unexpected results, and thus require investigation. To better understand the environmental behavior and associated risks of heavy metals, modeling metal fate and transport is increasingly recognized as an effective and efficient tool. Deterministic models based on the CDE have been assembled and widely used for describing metal’s physical, chemical and biological phenomena in the vadose zone (Pang et al., 2002; Liu et al., 2006; Jacques et al., 2008; Elbana and Selim, 2010; Qi et al., 2012; Zhang et al., 2014). Furthermore, reliable model prediction depends on appropriate parameterization of solute transport processes. Inverse modelling has become a popular method to parameterize models by minimizing an assumed goal or objective function which usually reduces to an optimization problem. The desired parameter values may in some sense represent the soil, at least at the local scale. To date, there have been few attempts to parameterize solute transport models by inverse modeling using data from a layered soil column. Porro et al. (1993) estimated transport parameters by fitting the CDE to the observed BTCs for each solute at various depths in each column using optimization method (only one or two parameters at once). Lu et al. (2011) simulated Cr transport in a landfill liner system composed of two-layer soils with the parameters back-calculated by trial-and-error method. Some researchers modeled reactive and nonreactive solute transport in multilayered soils with independently obtained parameters (Selim et al., 1977; Zhou et al., 2001; Fábio Joel Kochem Mallmann et al., 2012). It may not always be possible to find an optimal parameter set that best describes the observation. Many parameter combinations can fit the experimental data well, a problem termed equifinality, probably due to observation errors, model structure, or parameter interactions (Beven and Binley, 1992; Beven, 2006). The errors could be partially reduced by strictly controlling the experimental conditions, but it is difficult to improve model structure in the short term if there are limitations in current understanding. So the quantitative description of parameter uncertainty and its impact on model prediction has become a hot issue, which could help us to more deeply understand differences between the real world and model systems. Although uncertainty analysis cannot

86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104

reduce the uncertainty of the study itself, it could determine the risk of the estimation level and evaluate the reliability of the model simulation for further study. The stochastic method provides an alternative and could be applied to quantify parameter uncertainty. Over the last two decades, Monte Carlo simulation coupled with deterministic models has been used for uncertainty analysis. The GLUE methodology introduced by Beven and Binley (1992) is one such method and has been widely used because of its ease of use, and its adaptability to nonlinear systems. It was developed in the context of multiple sources of real problem uncertainty and assumes that different parameter sets can produce acceptable results. It has been widely used for uncertainty analyses and model calibration in hydrologic and environmental modellings (Blasone et al., 2008; Kumar et al., 2010; Steffens et al., 2013; Dusek et al., 2015; Beven and Binley, 2015), while its application to solute transport in layered soils is relatively rare. There might be obvious interactions between solute transport parameters for different soil layers, unrealistic to get an optimal solution, which could be implicitly shown by the GLUE method (Vazquez et al., 2009). Hence, our study objectives are to: (1) investigate the sorption and transport of coexisting Cu, Cd and Zn in layered soil using batch and column experiments; (2) estimate transport parameters of bromide and heavy metals through inverse modeling; and (3) explore the use of a GLUE method to assess parameter uncertainty and its implications on model prediction based on experimental data.

105

2. Methods and materials

106

2.1 Soil samples and characteristics

107 108 109 110 111 112

Soils were collected from a cornfield in Jimo, located in the Dagu River Basin, Shandong province, China. According to the distribution of the soil profile, disturbed soil samples were taken from the 0-30 cm and 30-60 cm horizons, respectively. The soil was air-dried, sieved through a 2-mm mesh and stored pending batch and column experiments. The soil physical and chemical properties, such as texture, pH, organic matter, and CEC were measured using standard analytical methods, which are provided in the Supplementary Information (Table S1).

113

2.2 Batch experiments

114 115 116 117 118 119 120 121 122 123 124

After homogenization of each horizon soil samples, 2.5 g of soil was added to 25 mL of metal solutions with initial pH of 4.0 in a 50 mL polyethylene centrifugal tube. The concentration for each metal was the same in solution and ranged from 20 mg L-1 to 500 mg L-1, prepared in 0.01M NaNO3 as background solution. To achieve equilibrium between soil and solute, all metal-treated samples were agitated on a shaker for 24 h at 25±1 ℃. The suspensions were centrifuged for 10 min at 4000 rpm and then filtered. Total Cu, Cd, and Zn in the supernatant were quantified using an Inductive Coupled Plasma-Atomic Emission Spectrometer (ICP-AES). The amount of adsorbed metal was calculated from the difference between the initial and equilibrium concentration of each solute. The supernatant pH was also measured. Three replicates were conducted for each treatment. Sorption isotherms were described by the extensively used Freundlich equation, expressed as,

125

S = KFC n

126

where S is the adsorbed concentration for each metal (mg g-1), C is the equilibrium concentration

(1)

127 128

in soil solution for each metal (mg mL-1), KF is the Freundlich distribution coefficient (mg1−n mLn g-1), and n is the sorption intensity (-) that typically has a value of less than 1.

129

2.3 Column experiments

130

2.3.1 Experimental setup

131 132 133 134 135 136 137 138

To evaluate the mobility of conservative tracer (bromide) and heavy metals under saturated flow condition, experiments were conducted in cylindrical plexiglass columns (20 cm length × 5 cm diameter). The detailed information of column packing could be found in Supplementary Information (Text S1). For a uniform distribution of the flow, a 1-cm space between the soil and the column inlet or outlet was filled with quartz sand. Prior to each experiment, the column was first fully saturated with a background solution (0.01M NaNO3) from the bottom and trapped air removed using a peristaltic pump. Once saturated, a flux of background solution was applied at the top of the column at a constant rate of 3±0.2 cm h-1.

139

2.3.2 Tracer and heavy metal transport experiments

140 141 142 143 144 145 146 147 148 149 150 151 152 153

Bromide was employed as the non-reactive chemical tracer. After steady-state flow condition was established (i.e. after applying 150 mL or about 1 pore volume of background solution), the input was switched to bromide solution (in the form of KBr) at a concentration of 0.05 M. The total volume of Br solution applied was 160 mL (corresponding to 1 pore volume), and this was followed by approximately 3 pore volumes of background solution to leach the tracer from the column. Next, a mixed solution with the same concentrations (500 mg L-1) of Cu, Cd and Zn was introduced into the columns for approximately 30 pore volumes pulse. Then the background NaNO3 solution was applied to evaluate the desorption process. All the influent solutions had a pH of 4.0, simulating acid rain. The effluents were sampled in glass bottles periodically using the automatic fraction collector. The heavy metal concentration in effluent was determined by ICP-AES and the Br by an ion-specific electrode. The pH and electrical conductivity (EC) of the effluents were measured using a pH meter and a conductivity meter, respectively. The experiment lasted for 257 h (around 11 d). Three column tests were duplicated to assure the consistency of column packing and the reproducibility of the experiments.

154

2.3.3 Transport model description

155 156

The models coupled in Hydrus_1d source code (Šimůnek et al., 2009) were employed in this study. Br transport is described with the standard CDE without reactive term:

157

∂C ∂ 2 C q ∂C =D 2 − ∂t ∂x θ ∂x

158 159 160 161 162 163 164

where C is the solute concentration in the liquid phase (mg mL-1), t is the time (h), x is the distance (cm), θ is the volumetric water content (cm3 cm-3), q is volumetric flux (cm h-1), and D is the dispersion coefficient (cm2 h-1) , given by D=λv=λq/θ where λ is the dispersivity (cm) and v the average pore-water velocity (cm h-1). According to the literature (Lair et al., 2006; Wong et al., 2007; Huang et al., 2014) and the batch experiments, heavy metal transport is described by the chemical non-equilibrium TSM coupled with the Freundlich equation as follows:

(2)



[ K F nC n −1 ])

∂C ∂ 2C q ∂C αρ  (1 − f ) K F C n − S k  =D 2 − − ∂t ∂x θ ∂x θ

165

(1+

166 167 168 169 170 171 172 173 174 175 176

where f is the fraction of sites assumed to be in equilibrium with the liquid phase (-), α is the first-order kinetic rate coefficient (h-1), Sk is the adsorbed concentration at non-equilibrium sites (mg g-1), and ρ is the soil bulk density (g cm-3). The model incorporates two kinds of sorption, one for instantaneous and the other for kinetically controlled, which is frequently used for analyzing BTCs of soil column experiments (van Genuchten and Wagenet, 1989; Pang et al., 2002; Chotpantarat et al., 2012; Zou and Zheng, 2013). For the experimental columns, the upper and lower boundary conditions are concentration flux and zero concentration gradient, respectively. The initial condition is free of solutes. The governing equations are solved numerically by a finite-element method and the inverse problem by the Levenberg-Marquardt method, which are implemented in Hydus-1d. The performance of the model was estimated by the Nash-Sutcliffe efficiency (R2).

177

2.4 Uncertainty analysis

178 179 180 181 182 183 184 185 186 187 188 189 190 191

GLUE methodology with the Hydrus_1d source code was used for uncertainty estimation. In this method, a large number of model runs are made with many different randomly chosen parameter values selected from a priori probability distributions. Based on quantitative measure of model performance (i.e. the likelihood measure), parameter sets are classified as ‘behavioral’ or ‘non-behavioral’, where the distinction is based on a user-specified cutoff threshold. The likelihood measure here corresponds to a specified measure of how well the model and parameter set describe the observations. The parameter sets identified as behavioral are then applied to assess the uncertainty on parameter identification and model prediction. Prediction limits are obtained by forming the cumulative density function (CDF) of the likelihood weighted ensemble of simulations. Any required quantiles can then be taken from the CDF. The basic requirements of the GLUE methodology are prior distributions of parameters, specified maximum and minimum value for each parameter, and a likelihood function. In this paper, Nash-Sutcliffe efficiency (R2), the one most often applied (Stedinger, et al., 2008), was chosen as the likelihood function.

θ

n

192

L (θ i Y ) = 1 −

∑ (Q

ij

j n

(3)

− Qoj ) 2 (4)

∑ (Q

oj

− Qo )

2

j

193 194 195 196 197 198 199 200

where L(θi|Y) is the likelihood measure of parameter set i, here referred to R2, n is the total number of observations, Qij and Qoj are the simulated and observed values at time j, and Qo is the average of the observations. In the GLUE method, the choice of threshold value is subjective and noted in the literature (Stedinger et al., 2008). Information from expert knowledge, references and batch experiments were taken into account to determine the feasible ranges of each specified parameter. Here it was assumed that the prior distributions of all parameters are uniform. Latin hypercube sampling (LHS) was applied to generate adequate and unbiased parameter sets, and reduce the number of runs (Zhang et al., 2006; Looms et al., 2008; Razmkhah, 2018).

201

3. Results and discussion

202

3.1 Adsorption isotherms

203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223

The sorption isotherms of metal onto soil obtained by batch experiments were well described by the Freundlich equation (Fig. S1 and Table S2). Generally, KF is positively related to the metal sorption capacity of soils. The KF values were respectively 3.64 and 3.54 mg(1-n) mLn g-1 for Cd and Zn in loam, which were greater than 3.25 and 2.91 mg(1-n) mLn g-1, in sandy loam. While the KF value for Cu in loam was 4.55 mg(1-n) mLn g-1, which was smaller than 5.10 mg(1-n) mLn g-1, in sandy loam. This indicates that soil with different texture and chemical properties have different adsorption capacity for different heavy metals. Loam with high OC and CEC exhibited strong sorption for Cd and Zn, while sandy loam with high pH exhibited strong sorption for Cu, which means that soil properties affect the sorption capacity for different heavy metals in different extent. The order of metal adsorption to soil according to their KF value was Cu>Cd>Zn for each soil, which coincides with the findings of Gomes et al. (2001) and Elbana et al. (2018), following the order of electronegativity of the metal cations (Abd-Elfattah and Wada, 1981). Cu is a divalent metal with a strong tendency to form hydroxyl-ion pairs and is specifically adsorbed to a large extent onto soluble organic matter (SOM), oxides, and clays (Adriano, 2001). The isotherms were nonlinear with the exponent n values ranging from 0.26 to 0.34 which illustrate the dependence of the sorption process on concentration where sorption by the highest energy sites takes place preferentially at the lowest solution concentration. Furthermore, to estimate the retention of heavy metals in soil, the nonlinear isotherm was linearized using two methods (Van Genuchten, 1981) to calculate the distribution coefficient Kd (mL g-1) in the linearized isotherm of the form S=KdC, as follows. As listed in Table S2, the results were consistent with the above.

224

K d1 =

225

Kd2 = KF C0n−1

226

Detailed information about pH is presented in Supplementary Information (Text S2).

227

3.2 Column experiments

228

3.2.1 BTCs for Br and heavy metals

229 230 231 232 233 234 235 236 237 238 239 240

The relative concentrations (C/C0) for Br and heavy metals in the effluents were plotted as a function of pore volumes (vt/L), where L is the soil column length, as shown in Fig.1. Prior to introducing heavy metals through the column, the BTC was established for Br. The shape of the Br BTC was symmetrical, indicating that there was no physical non-equilibrium transport in the column. The BTC peak nearly reached the inlet concentration and 94.5% (little trapping in the dead pores) of the applied Br was recovered in the effluent solutions, indicating Br as a non-reactive tracer (Jiang et al., 2005; Dousset et al., 2007). In contrast to Br transport, the migration pattern of metals in the soil column differed significantly. As seen in Fig.1, the initial breakthrough occurred for Cd and Zn after over about 10 pore volumes, followed by a rapid concentration increase; while Cu appeared later in the effluent after 20 pore volumes, followed by a relatively slow concentration increase. Moreover, the peak BTCs for Zn and Cd were 1.13 and 0.88 respectively, which were much higher than Cu (0.35). When the inlet

2 K F C0n −1 n +1

(5)

(6)

241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261

pulse was switched to background solution, the effluent concentration dropped sharply followed by a long tail during desorption phase (lasting for about 60 pore volumes). This behavior indicated that sorption time was fast relative to residence time, with a slow release of heavy metals that were adsorbed in the interior of the soil aggregates back into the soil solution. This also indicated kinetically controlled adsorption and desorption of heavy metals to and from the soil (Elbana and Selim, 2010; Olaofe et al., 2014). Based on the area under the curve, the recovery rates for Zn, Cd, and Cu after desorption were 68.0%, 52.5% and 6.2%, respectively. The trend was in line with the sorption capacity obtained in batch experiments. Therefore, shallow water resources may be more exposed to contamination risk by Zn and Cd. Unexpectedly, the effluent concentration of Zn at 20 pore volumes was greater than that of the influent concentration (C/C0>1), which could be explained by the “snow plow effect” (Starr and Parlange, 1979; Persaud and Wierenga, 1982), resulting from competition for sorption sites between species. Cu ions were adsorbed related to an inner-sphere complexation of the metal, while Zn and Cd ions were more bound by outer-sphere complexation (Sposito, 2008). Moreover, Cd is especially likely to adsorb to the same sites as Zn and hence is a stronger competitor for specific sites (Voegelin et al., 2001). Thus the Zn cations could be easily exchanged when competing with other metal cations (Covelo et al., 2004; Trakal et al., 2011), and the adsorbed Zn cations removed from the exchange sites were then forced into the interstitial solution, which is pushed ahead by the incoming solution front, as a sort of snow plow. Detailed information about pH and EC of the effluent is provided in Supplementary Information (Text S3 and Fig.S2).

262

3.2.2 Modeling results and fitted parameters of Br

263 264 265 266 267 268 269

The standard CDE was applied to describe the behavior of non-reactive chemicals in soil. The θ in all equations was replaced by the saturated volumetric water content θs, because our experiments were conducted under saturated flow condition. Although the value θs could be assumed equal to porosity calculated by the known soil column bulk density and mineral particle density, its value of 2.65 g cm-3 is an empirical value, so θs was fitted and compared with the calculated values. Other parameters were estimated from soil bulk density and texture using ROSETTA built-in Hydrus_1d software. The q was calculated by the measured steady outflow rate.

270

Deterministic analysis

271 272 273 274 275 276 277 278 279 280 281 282

Three trial optimizations, each with differing initial estimates of parameter values, yielded different optimized parameters. The solution strongly depended on the initial values, thus the reported parameters were considered non-unique. Here, only the best fitting parameter set as assessed by R2 is listed in Table 1 and the corresponding simulated BTC is shown in Fig.2 (inversion). Although the model fitted the Br concentrations well with R2 = 0.99, the uncertainty of parameters due to interactions could not be ignored. The standard errors (SE) were so high that it made some parameters physically meaningless (e.g. θs1 and θs2). So, in this case, it was impossible to regard the obtained parameter set as unique and optimal. And if these parameters were used for prediction, the reliability and accuracy of the results could not be guaranteed. Therefore, tracer experiments were conducted on uniform loam and sandy loam respectively, under identical experimental conditions to those of the layered soil. Then the Br BTCs were fitted respectively to obtain θs1, λs1, θs2 and λs2. Observed and simulated BTCs are shown in Fig.S3, and

283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299

the estimated parameters listed in Table 1. The simulation results were good with R2>0.95 and low SE. Compared to the simulated curves, the slightly lower peak of the BTCs suggested that water flow was restricted by some intra-aggregate pores where water could not move at all and solute exchange was controlled by diffusion only. Taking this into account, the mobile-immobile model (MIM) code (Van Genuchten and Wieranga, 1976) was employed. However, it did not or only slightly improve the fitting with a negligibly low estimated fraction of immobile water. Together, these indicated that the standard CDE model was sufficient to describe solute transport in the soil-packed column. The estimated dispersivity was 0.24 cm for loam and 0.31 cm for sandy loam, respectively. At the column scale, dispersivity is generally correlated with pore size distribution and soil structure heterogeneity (Jellali et al., 2010) and low value indicates that the experimental soil was a relatively homogeneous porous media. The estimated θs was 0.46 cm3 cm-3 for loam and 0.43 cm3 cm-3 for sandy loam, which was consistent with the calculated values of 0.48 and 0.46 cm3 cm-3, respectively. Then the fitted parameters from uniform soil columns were applied to simulate Br transport in the layered soil column. Simulated values (forward) fitted well with the observed data (R2 = 0.98), a finding also noted by others, that parameters obtained from uniform soil are near the true values (Zhang et al., 2006; Lin and Xu, 2012). Therefore, these parameter values were employed in the subsequent modeling of metal transport in the layered soil column.

300

GLUE analysis

301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325

Considering the non-unique parameters obtained in the inverse simulation, the GLUE approach was applied to assess the uncertainty of parameters and prediction for Br transport simulations. According to experimental conditions and the literature, the sampling limits of the parameters were fixed to 0.2 and 0.8 cm3 cm-3 for θs, and 0 and 1 cm for λ, respectively. The parameters were sampled by LHS using uniform random sampling across specified ranges. A total of 10 000 Monte Carlo realizations were used. The governing equation (CDE) was also solved by the same numerical solution contained in the subroutines of the Hydrus_1d. The simulated Br BTC (forward simulation) was compared with the measured for each run, and the R2 was calculated. In general, the higher the R2 value, the better the model is at predicting the data, and the more likely the respective parameter set will occur. Here 0.9 was chosen as the threshold value and 1979 acceptable parameter sets were obtained as a posterior distribution. Scatter plots of behavioral parameter sets are shown in Fig.3. Each dot in the scatter plot represents a single model run of Br. The plots showed that the ranges of acceptable θs and λ values were similar to the original, ranging from 0.20 to 0.75 cm3 cm-3 for θs1 and θs2, and acceptable parameters scattered over the full region of the initial domain for λ1 and λ2. There were no clear peak values that could be detected, indicating that there was no single optimum parameter set, but many parameter sets could provide similar R2 values, which verified the phenomenon of ‘equifinality’. The cumulative likelihood distribution curves also did not exhibit large steep sections, similar to that of prior uniform distributions, suggesting that these parameters were poorly constrained by the available observations. This may be due to low regional sensitivity and possible parameter interactions. Here, correlation analysis found that there was a high correlation coefficient between θs1 and θs2, equal to -0.994. To determine the extent to which parameter uncertainty affected model simulation, these acceptable parameter sets were applied to calculate the 95% confidence interval (95 CI) of the output. The likelihood values of the retained solutions were rescaled, so that the outflow

326 327 328 329 330 331

concentration was expressed as a cumulative distribution function. This procedure allowed the uncertainty associated with all sources of error in the modeling process to be carried forward into the predictions. The 95 CI for Br concentration in effluent is shown in Fig.2. It was evident that most of the measured data lie within the 95 CI except those at the points of breakthrough and tailing, and prediction range was relatively narrow. The outlier data may be due to the ion-specific electrode not being able to accurately measure Br at low concentration.

332

3.2.3 Modeling results and fitted parameters of heavy metals

333

Deterministic analysis

334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368

Since sorption coefficients are easily obtained from batch adsorption experiments, it is desirable to utilize them if they provide a good estimate of the metal concentration in soil. So, KF and n derived independently from the batch experiments, combined with λ and θs from the tracer experiments, were applied to generate predicted BTCs for heavy metals using the local equilibrium model in Hydrus_1d. The results showed that the predicted heavy metal BTCs (not shown in this paper) largely deviated from the experimental data. The reason may be that soil residence time was not long enough for some sorption sites to reach equilibrium (Chotpantarat et al., 2011). Alternatively, the deviations may be due to the different soil: solution ratios between the batch and column experiments (Buergisser et al., 1993; Zhang et al., 2014). Unlike the Br BTCs, heavy metals BTCs were asymmetrical and exhibited significant tailing. Trace experiments showed that there was no physical non-equilibrium in the soil column system, so that this was primarily the result of a chemical non-equilibrium process, probably caused by rate-limited sorption. Therefore, the TSM model was applied to simulate the BTCs of metal transport in the column. Parameters θs and λ, which were obtained by fitting Br BTCs in uniform soil, were fixed and parameters KF, n and α were optimized. As shown in Fig.4, the fits provided by the TSM were good except for Zn. It did not depict the “snow plow effect” phenomenon because of shortcomings in the model, in that it does not contain the exchange and superposition effect of the adsorbed and interstitial solute. The estimated transport parameters are listed in Table 2. The KF values were similar to those derived from batch tests, but in contrast the sorption intensity values (n) were much smaller. Differences in the estimated n values might be due to different flow schemes (static and dynamic). The f values were relatively small indicating that most of the sorption occurring in the column was rate-limited/time-dependent. To minimize the number of fitted parameters, the value of f in Eq. (3) was set to zero, in that way, the model TSM was converted to OSM which was then applied to simulate the transport of heavy metals. The fitted BTCs are presented in Fig.4 and the estimated parameters were given in Table 2. Also, the observed and predicted outflow curves matched well (the R2 was greater than 0.9). Similar to the Br simulation, the TSM and OSM models were sensitive to the initial parameter values and the optimized parameters were not unique, indicating that the pre-estimated initial parameters have a great influence on the optimized parameters (Dubus et al., 2004). In some cases, the results were unreasonable, e.g. the fitted parameter KF1 was 2.43 mg(1-n) mLn g-1 and KF2 11.64 mg(1-n) mLn g-1 for Cu. However, there were 2 fewer OSM parameters than in the TSM, so to some extent the reliability of the OSM prediction was higher. Nevertheless, if these values were used for prediction, the reliability and accuracy could not be guaranteed. Transport experiments for heavy metals were not carried out in uniform soil columns because of the influence of initial concentration on transport parameters for reactive solutes (Shi,

369 370

2016; Shahmohammadi-Kalalagh and Taran, 2018). The uncertainty of transport simulation for heavy metals was analyzed using GLUE below.

371

GLUE analysis

372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408

In the GLUE analysis, the initial KF range was determined from the results of batch experiments and the deterministic calibration, from 0 to 10 mg(1-n) mLn g-1 for Cd and Zn, while from 0 to 20 mg(1-n) mLn g-1 for Cu. The n limits were fixed to 0 and 1, and α range was taken from our previous studies and the literature, from 0 to 1 h-1. Due to the lack of prior information on parameter distributions, it was also assumed that each model parameter follows a uniform distribution over its range. A total of 50,000 parameter sets were generated with the help of LHS. For each parameter set the Hydrus_1d code was run and parameter sets were ranked according to their likelihood function (R2). The threshold value of R2 was set at 0.7 for Zn and Cd, and 0.6 for Cu, considering the quantity of acceptable parameter sets and the quality of simulations. Scatter plots of acceptable parameter sets are shown in Fig.S4. The number of behavioral simulations was 974 for Zn, 332 for Cd, and 111 for Cu., making uncertainty analysis feasible (Sun et al., 2016). There were also no clear peak values. As shown in Fig.S4, the posterior cumulative KF distributions were close to uniform, while the others differed considerably. Previous studies pointed out strong differences among the cumulative distributions of prior and posterior parameters, suggesting that the model predictions are sensitive to a particular parameter, whereas similarity between them suggests that the output results are insensitive to that parameter (Freer et al., 1996; Ghavidelfar et al., 2015). However, here non-identifiability of a parameter might be caused by parameter interactions rather than the model being insensitive to the parameter. Correlations between the acceptable parameters estimated from the GLUE method are shown in Supplementary Information (Table S3). The parameters from different soil were highly correlated such as KF (0.624, 0.578 and 0.556) and n (0.325, 0.322 and 0.543) for Zn, Cd and Cu respectively. Also parameters from the same soil were correlated, for example, KF and n with correlation coefficients between 0.333 and 0.542. Due to the mutual correlation between n and KF it is suggested that these two parameters should not be optimized at the same time. Prediction uncertainty for heavy metal transport was analyzed just as Br. The 95 CI of simulations compared with the observed data is shown in Fig.4. Most of the observed Zn and Cd values were inside the prediction limits and the prediction ranges were slightly wider than that for Br, resulting from more estimated parameters and their correlations. A few observed Zn points were higher than their corresponding 97.5% confidence limits, further indicating that the Hydrus_1d code does not always fully simulate the transport and retention processes. Almost all the observed Cu points, however, systematically fell below the 2.5% confidence limits and the effluent concentration was overestimated in the simulation, which could be attributed to errors in the model structure and observational data (Beven, 2006; Beven and Westerberg, 2011; Beven and Smith, 2015). It could be that Cu precipitation occurred in the top soil layer for high pH in our cases and such a process was not adequately described in the model. Presumably the model would similarly over predict solute concentration under similar circumstances in future prediction. This problem could not be discovered using the deterministic optimization method.

409

4. Conclusions

410

Natural subsurface soil is heterogeneous and, in many cases, may be approximated as a layered

411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440

system. It is essential to understand and simulate the cotransport behaviors of heavy metals in more complex and realistic environment. Therefore, we conducted laboratory experiments and developed numerical models to investigate the characteristics of Br and heavy metals (coexistence of Cu, Cd and Zn) retention and transport in layered soil columns. The batch and column leaching experiments showed that Cd and Zn were much more mobile than Cu when they coexisted in the soil, indicating that groundwater under heavy metal contaminated soil is at greater risk of Zn and Cd pollution. It is noteworthy that the outflow relative concentration of Zn is greater than 1 possibly due to the “snow plow effect” resulting from the competitive adsorption between different species for the soil sites. The adsorption and desorption of heavy metals to and from the soil were kinetically controlled at the soil column scale, and the sorption kinetics may be weakened as transport modeling is performed at the larger profile scale. Except for Zn (partial C/C0>1), the CDE, TSM and OSM models within Hydrus_1d resulted in good fits to the observed BTCs of Br and heavy metals in the process of parameter optimization. However, this did not imply the optimized parameters were unique and correct, nor that the resulting model will be robust in prediction. The GLUE method was applied to assess the uncertainty of parameters and outputs to condition a posterior distribution of behavioral parameter sets. This resulted in much wider ranges of feasible parameters instead of one single optimum parameter set for the observed data, which was probably due to the complex correlation of the model parameters from the same and different soil layers. It is recommended that parameters such as KF and n should not be optimized simultaneously. Most of the measured BTCs points except for Cu fell within the prediction interval. The overestimated Cu concentration was due to precipitation that occurred in the column system at high soil pH. Additionally, the “snow plow effect” phenomenon was not depicted due to shortcomings in the model. Therefore, further studies should be continued in the field of model structure, to achieve more accurate simulation of solute transport in porous media. A deterministic model coupled with a Monte Carlo framework is a good choice for prediction of pollutant transport in layered soil where there is high correlation between parameters. Indeed, unlike classical parameter optimization algorithms, the GLUE approach provides the probability distribution, which can be applied to summarize parameter uncertainty and quantify its effects on model prediction.

441

Acknowledgment

442 443 444

This work was supported by the National Natural Science Foundation of China under Grant No. 41807010 and 41571214, and National Key R&D Program of China under Grant No. 2016YFC0402807.

445

References

446 447 448 449 450 451 452 453 454 455 456

1. 2. 3. 4. 5. 6. 7.

Abd-Elfattah, A., Wada, K., 1981. Adsorption of lead, copper, zinc, cobalt, and cadmium by soils that differ in cation-exchange materials. Eur. J. Soil Sci. 32, 271–283. Adriano, D.C., 2001. Trace Elements in the Terrestrial Environment. Springer-Verlag, New York, p. 45. Beven, K., 2006. A manifesto for the equifinality thesis. J. Hydrol. 320, 18–36. Beven, K., Binley, A., 1992. The future of distributed models: model calibration and uncertainty prediction. Hydrol. Process. 6, 279–298. Beven, K., Binley, A., 2015. GLUE: 20 years on. Hydrol. Process. 28, 5897–5918. Beven, K., Smith, P., 2015. Concepts of information content and likelihood in parameter calibration for hydrological simulation models. J. Hydrol. Eng. 20, A4014010. Beven, K., Westerberg, I., 2011. On red herrings and real herrings: disinformation and information in hydrological inference. Hydrol. Process. 25, 1676–1680.

457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523

8. 9. 10. 11.

12. 13.

14. 15.

16. 17.

18.

19. 20. 21.

22. 23. 24. 25.

26.

27. 28. 29. 30. 31. 32. 33. 34.

35. 36.

Blasone R.S., Madsen, H., Rosbjerg, D., 2008. Uncertainty assessment of integrated distributed hydrological models using GLUE with Markov Chain Monte Carlo sampling. J. Hydrol. 353, 18–32 Buergisser, C.S., Cernik, M., Borkovec, M., Sticher, H., 1993. Determination of nonlinear adsorption isotherms from column experiments: an alternative to batch studies. Environ. Sci. Technol. 27, 943–948. Chotpantarat, S., Ong, S.K., Sutthirat, C., Osathaphan, K., 2011. Effect of pH on transport of Pb2+, Mn2+, Zn2+ and Ni2+ through lateritic soil: column experiments and transport modeling. J. Environ. Sci. 23, 640–648. Chotpantarat, S., Ong, S.K., Sutthirat, C., Osathaphan, K., 2012. Competitive modeling of sorption and transport of Pb2+, Ni2+, Mn2+ and Zn2+ under binary and multi-metal systems in lateritic soil columns. Geoderma 189-190, 278–287. Covelo, E.F., Andrade, M.L., Vega, F.A., 2004. Heavy metal adsorption by humic umbrisols: selectivity sequences and competitive sorption kinetics. J. Colloid Interf. Sci. 280, 1–8. Dousset, S., Thevenot, M., Pot, V., Šimůnek, J., Andreux, F., 2007. Evaluating equilibrium and non-equilibrium transport of bromide and isoproturon in disturbed and undisturbed soil columns. J. Contam. Hydrol. 94, 261–276. Dubus, I.G., Beulke, S., Brown, C.D., Gottesburen, B., Dieses, A., 2004. Inverse modelling for estimating sorption and degradation parameters for pesticides. Pest Manag. Sci. 60, 859–874. Dusek, J., Dohnal, M., Snehota, M., Sobotkova, M., Ray, C., Vogel, T., 2015. Transport of bromide and pesticides through an undisturbed soil column: a modeling study with global optimization analysis. J. Contam. Hydrol. 175-176, 1–16. Elbana, T.A., Selim, H.M., 2010. Cadmium transport in alkaline and acidic soils: miscible displacement experiments. Soil Sci. Soc. Am. J. 74, 1956–1966. Elbana, T.A., Selim, H.M, Akrami, N., Newman, A., Shaheen, S.M., Rinklebe, J., 2018. Freundlich sorption parameters for cadmium, copper, nickel, lead, and zinc for different soils: influence of kinetics. Geoderma, 324, 80–88. Fábio Joel Kochem Mallmann, Santos, D.R.D., Cambier, P., Jérôme Labanowski, Lamy, I., Santanna, M.A., Tessier, D., Oort, F.V., 2012. Using a two site-reactive model for simulating one century changes of Zn and Pb concentration profiles in soils affected by metallurgical fallout. Environ. Pollut. 162, 294–302. Fonseca, B., Figueiredo, H., Rodrigues, J., Queiroz, A., Tavares, T., 2011. Mobility of Cr, Pb, Cd, Cu and Zn in a loamy sand soil: a comparative study. Geoderma, 164, 232–237. Freer, J., Beven, K., Ambroise, B., 1996. Bayesian estimation of uncertainty in runoff prediction and the value of data: an application of the GLUE approach. Water Resour. Res. 32, 2161–2173. Garrido-Rodriguez, B., Cutillas-Barreiro, L., Fernández-Calviño, D., Arias-Estévez, M., Fernández-Sanjurjo, M.J., Álvarez-Rodríguez, E., Núñez-Delgado, A., 2014. Competitive adsorption and transport of Cd, Cu, Ni and Zn in a mine soil amended with mussel shell. Chemosphere, 107, 379–385. Ghavidelfar, S., Shamseldin, A.Y., Melville, B.W., 2015. Estimation of soil hydraulic properties and their uncertainty through the beerkan infiltration experiment. Hydrol. Process. 29, 3699–3713. Gomes, P.C., Fontes, M.P.F., Da Silva, A.G., Eduardo, D.S.M., Netto, André R., 2001. Selectivity sequence and competitive adsorption of heavy metals by brazilian soils. Soil Sci. Soc. Am. J. 65, 1115–1121. Huang, B., Li, Z., Huang, J., Guo, L., Nie, X., Wang, Y., Zhang, Y., Zeng, G., 2014. Adsorption characteristics of Cu and Zn onto various size fractions of aggregates from red paddy soil. J. Hazard. Mater. 264, 176–183. Jacques, D., Šimůnek, J., Mallants, D., Genuchten, M.T.V., 2008. Modelling coupled water flow, solute transport and geochemical reactions affecting heavy metal migration in a podzol soil. Geoderma 145, 449– 461. Jellali, S., Diamantopoulos, E., Kallali, H., Bennaceur, S., Anane, M., Jedidi, N., 2010. Dynamic sorption of ammonium by sandy soil in fixed bed columns: evaluation of equilibrium and non-equilibrium transport processes. J. Environ. Manag. 91, 897–905. Jiang, G.G., Noonan, M.J., Buchan, G.D., Smith, N., 2005. Transport and deposition of bacillus subtilis through an intact soil column. Soil Res. 43, 695–703. Kumar, S., Sekhar, M., Reddy, D.V., Mohan Kumar, M.S., 2010. Estimation of soil hydraulic properties and their uncertainty: comparison between laboratory and field experiment. Hydrol. Process. 24, 3426–3435. Lafuente, A.L., González, C., Quintana, J.R., Vázquez, A., Romero, A., 2008. Mobility of heavy metals in poorly developed carbonate soils in the Mediterranean region. Geoderma 145, 238–244. Lair, G.J., Gerzabek, M.H., Haberhauer, G., Jakusch, M., Kirchmann, H., 2006. Response of the sorption behavior of Cu, Cd, and Zn to different soil management. J. Plant Nutr. Soil Sci. 169, 60–68. Li, Z., Ma, Z., van der Kuijp, T.J., Yuan, Z., Huang, L., 2014. A review of soil heavy metal pollution from mines in China: pollution and health risk assessment. Sci. Total Environ. 468, 843–853. Liu, C.L., Chang, T.W., Wang, M.K., Huang, C.H., 2006. Transport of cadmium, nickel, and zinc in tao yuan red soil using one-dimensional convective-dispersive model. Geoderma 131, 181-189. Lin, Q., Xu, S., 2012. Parameter uncertainty analysis of solute transport in saturated porous media based on GLUE method. J. Hydraul. Eng. 43, 1017–1024 (in Chinese). Looms, M.C., Binley, A., Jensen, K.H., Nielsen, L., Hansen, T.M., 2008. Identifying unsaturated hydraulic parameters using an integrated data fusion approach on cross-borehole geophysical data. Vadose Zone J. 7, 238–248. Lu, H., Luan, M., Zhang, J., 2011. Study on transport of Cr (Ⅵ) through the landfill liner composed of two-layer soils. Desalination 266, 87–92. Olaofe, O., Olagboye, S.A., Akanji, P.S., Adamolugbe, E.Y., Fowowe, O.T., Olaniyi, A.A., 2014. Kinetic studies of adsorption of heavy metals on clays. J. Am. Chem. Soc. 7, 48–54.

524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590

37. Pang, L., Close, M., Schneider, D., Stanton, G., 2002. Effect of pore-water velocity on chemical nonequilibrium transport of Cd, Zn, and Pb in alluvial gravel columns. J. Contam. Hydrol. 57, 241–258. 38. Persaud, N., Wierenga, P.J., 1982. A differential model for one-dimensional cation transport in discrete homoionic ion-exchange media. Soil Sci. Soc. Am. J. 46, 482–490. 39. Porro, I., Wierenga, P.J., Hills, R.G., 1993. Solute transport through large uniform and layered soil columns. Water Resour. Res. 29, 1321–1330. 40. Qi, Z., Feng, S., Helmers, M.J., 2012. Modeling cadmium transport in neutral and alkaline soil columns at various depths. Pedosphere 22, 273–282. 41. Rattan, R.K., Datta, S.P., Chhonkar, P.K., Suribabu, K., Singh, A.K., 2005. Long-term impact of irrigation with sewage effluents on heavy metal content in soils, crops and groundwater-a case study. Agr. Ecosyst. Environ. 109, 310–322. 42. Razmkhah, H., 2018. Parameter uncertainty propagation in a rainfall–runoff model; case study: Karoon-III river basin. Water Resour. 45, 34–49. 43. Sana, D., Jalila, S., 2018. Adsorption characteristics of layered soil as delay barrier of some organic contaminants: experimental and numerical modeling. Environ. Modell. Softw.110, 95–106. 44. Selim, H.M., Davidson, J.M., Rao, P.S.C., 1977. Transport of reactive solutes through multilayered soils. Soil Sci. Soc. Am. J. 41, 3–10. 45. Seuntjens, P., 2002. Field-scale cadmium transport in a heterogeneous layered soil. Water Air Soil Poll. 140, 401–423. 46. Shahmohammadi-Kalalagh, S., Taran, F., 2018. Effect of initial concentration and input flux on equilibrium and non-equilibrium transport of Zn in soil columns. Int. J. Environ. Sci. Te. Doi: 10.1007/s13762-018-2159-z. 47. Shi Y., 2016. Adsorption and Migration Behaviors of Many Metal Ions under Two Adsorbents. Dalian University of Technology, Dalian, Liaoning, China. 48. Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.Th., 2009. The Hydrus-1d Software Package for Simulating the One-dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media. version 4.08, HYDRUS software series 3, Department of Environmental Sciences, University of California Riverside, Riverside, California. 49. Sposito, G., 2008. The Chemistry of Soils. Oxford University Press, Inc., New York. 50. Starr, J.L., Parlange, J.Y., 1979. Dispersion in soil columns: The snow plow effect. Soil Sci. Soc. Am. J. 43, 448–450. 51. Stedinger, J.R., Vogel, R.M., Lee, S.U., Batchelder, R., 2008. Appraisal of the generalized likelihood uncertainty estimation (GLUE) method. Water Resour. Res. 44, W00B06. 52. Steffens, K., Larsbo, M., Moeys, J., Jarvis, N., Lewan, E., 2013. Predicting pesticide leaching under climate change: importance of model structure and parameter uncertainty. Agr. Ecosyst. Environ. 172, 24–34. 53. Sterckeman, T., Douay, F., Proix, N., Fourrier, H., 2000. Vertical distribution of Cd, Pb and Zn in soils near smelters in the north of France. Environ. Pollut. 107, 377–389. 54. Sun, M., Zhang, X., Huo, Z., Feng, S., Huang, G., Mao, X., 2016. Uncertainty and sensitivity assessments of an agricultural-hydrological model (RZWQM2) using the GLUE method. J. Hydrol. 534, 19–30. 55. Trakal, L., Komárek, M., Száková, J., Zemanová, V., Tlustoš, P., 2011. Biochar application to metal-contaminated soil: evaluating of Cd, Cu, Pb and Zn sorption behavior using single- and multi-element sorption experiment. Plant Soil Environ. 57, 372–380. 56. van Genuchten, M.Th., 1981. Non-equilibrium Transport Parameters from Miscible Displacement Experiments. No.119. Riverside, CA: U.S. Salinity Lab. pp41. 57. van Genuchten, M.Th., Wagenet, R.J., 1989. Two-site/two-region models for pesticide transport and degradation: Theoretical development and analytical solutions. Soil Sci. Soc. Am. J. 53, 1303–1310. 58. van Genuchten, M.Th., Wierenga P. J., 1976. Mass transfer studies in sorbing porous media, Ⅵ. analytical solutions. Soil Sci. Soc. Am. J. 40, 473–480. 59. Vazquez, R.F., Beven, K., Feyen, J., 2009. GLUE based assessment on the overall predictions of a MIKE SHE application. Water Resour. Manag. 23, 1325–1349. 60. Voegelin, A., Vulava, V.M., Kretzschmar, R., 2001. Reaction-based model describing competitive sorption and transport of Cd, Zn, and Ni in an acidic soil. Environ. Sci. Technol. 35, 1651–1657. 61. Waterlot, C., Bidar, G., Pelfrêne, A., Roussel, H., Fourrier, H., Douay, F., 2013. Contamination, fractionation and availability of metals in urban soils in the vicinity of former lead and zinc smelters, France. Pedosphere 23, 143–159. 62. Wong, J.W.C., Li, K., Zhou, L., Selvam, A., 2007. The sorption of Cd and Zn by different soils in the presence of dissolved organic matter from sludge. Geoderma 137, 310–317. 63. Yang, Q., Li, Z., Lu, X., Duan, Q., Huang, L., Bi, J., 2018. A review of soil heavy metal pollution from industrial and agricultural regions in China: Pollution and risk assessment. Sci. Total Environ. 642, 690–700. 64. Zhan, M., Wang, Y., Ge, J., Sheng, J., Luo, Y., 2018. Investigation for the mechanism of copper ion migration regularity in layered soil. J. Saf. Environ. 18, 324–329 (in Chinese). 65. Zhang, D., Beven, K., Mermoud, A., 2006. A comparison of non-linear least square and GLUE for model calibration and uncertainty estimation for pesticide transport in soils. Adv. Water Resour. 29, 1924–1933. 66. Zhang, J., Li, H., Zhou, Y., Dou, L., You, J., 2018. Bioavailability and soil-to-crop transfer of heavy metals in farmland soils: A case study in the Pearl River Delta, South China. Environ. Pollut. 235, 710–719. 67. Zhang, H., Li, L., Zhou, S., 2014. Kinetic modeling of antimony (V) adsorption–desorption and transport in soils. Chemosphere, 111, 434–440.

591 592 593 594 595 596

68. Zhang, X., Zhong, T., Liu, L., Ouyang, X., 2015. Impact of soil heavy metal pollution on food safety in china. Plos One, 10, e0135182. Doi:10.1371/journal.pone.0135182. 69. Zhou, L., Selim, H.M., 2001. Solute transport in layered soils. Soil Sci. Soc. Am. J. 65, 1056–1064. 70. Zou, Y., Zheng, W., 2013. Modeling manure colloid-facilitated transport of the weakly hydrophobic antibiotic florfenicol in saturated soil columns. Environ. Sci. Technol. 47, 5185–5192.

597

Figures

598

Fig.1. Measured BTCs of bromide and heavy metals in the layered soil (C/C0 denotes relative concentration; Pore

599 600 601 602 603

volume denotes relative time equal to vt/L.) Fig.2. Measured and simulated BTCs of Br in the layered soil Fig.3. Scatter plot of R2 and cumulative likelihood distributions of each input parameter determined using GLUE

604

(Lines indicate posterior distribution; dash lines indicate prior distribution) Fig.4. Measured and simulated BTCs of heavy metals in layered soil column

605

1.2

Br

1.0

C/C0

0.8 0.6 0.4 0.2 0.0

0

1

2

3

4

Pore volume

606 1.2

Cu Cd Zn

1.0

C/C0

0.8 0.6 0.4 0.2 0.0

0

20

40

60

80

100

Pore volume

607 608 609 610 611 612

Fig.1. Measured BTCs of bromide and heavy metals in the layered soil (C/C0 denotes relative concentration; Pore volume denotes relative time equal to vt/L.)

observed 95 CI

inverse modeling forward modeling

1.2

Layered

1.0

C/C0

0.8

0.6

0.4

0.2

0.0

0

1

2

Pore Volume

613 614

Fig.2. Measured and simulated BTCs of Br in the layered soil

3

4

1.0

0.98

0.8

0.98

0.8

0.96

0.6

0.96

0.6

0.94

0.4

0.94

0.4

0.92

0.2

0.92

0.2

0.4

0.5

0.6

0.7

0.90 0.0

0.2

0.4

θ s1

615

0.8

0.0 1.0

1.0

0.98

0.8

0.98

0.8

0.96

0.6

0.96

0.6

0.94

0.4

0.94

0.4

0.92

0.2

0.92

0.2

0.3

0.4

0.5

0.6

0.7

0.0 0.8

R

2

1.00

CDF

1.0

2

R

0.6

λ1

1.00

0.90 0.2

616 617 618 619

2

0.0 0.8

0.90 0.0

0.2

θs2

0.4

0.6

0.8

CDF

0.3

R

R

0.90 0.2

CDF

1.00

CDF

1.0

2

1.00

0.0 1.0

λ2

Fig.3. Scatter plot of R2 and cumulative likelihood distributions of each input parameter determined using GLUE (Lines indicate posterior distribution; dash lines indicate prior distribution)

620 Observed 95CI

OSM TSM

1.2

Cu

1.0

C/C0

0.8 0.6 0.4 0.2 0.0 0

20

40

60

80

100

Pore Volume

621 1.2

Cd 1.0

C/C0

0.8 0.6 0.4 0.2 0.0 0

20

40

60

80

100

Pore Volume

622 1.2

Zn 1.0

C/C0

0.8 0.6 0.4 0.2 0.0 0

623 624 625

20

40

60

80

Pore Volume

Fig.4. Measured and simulated BTCs of heavy metals in layered soil column

100

626

Tables

627 628

Table 1 Calculated and fitted values of transport parameters for Br using the standard CDE model

629 630 631 632 633

θs1

λ1

θs2

λ2

cm h-1

cm3 cm-3

cm

cm3 cm-3

cm

Loam

3.06

0.46±0.01

0.24±0.05

-

-

0.99

Sandy Loam

3.18

-

-

0.43±0.01

0.31±0.13

0.95

Layered soil (forward)

3.06

0.46

0.24

0.43

0.31

0.98

Layered soil (inversion)

3.06

0.49±0.88

0.21±0.78

0.49±0.91

0.02±0.03

0.99

Subscript 1 and 2 refer to loam and sandy loam, respectively; Data are presented as mean ± 95% confidence interval

Table 2 Fitted values of transport parameters for heavy metals using the OSM and TSM model KF1

n1

f1

mg(1-n)mLn g-1

α1

KF2

d-1

mg(1-n)mLn g-1

n2

8.45±0.15

0.22±0.70e-2

-

2.01±0.07

9.15±0.18

0.03±0.10e-2

Cu**

2.43±0.27

0.25±0.02

0.15±0.01

0.53±0.03

11.64±0.24

0.08±0.12e-2

Cd*

2.96±0.28

3.53±0.14

0.01±0.07e-2

2.11±0.35

2.11±0.38e-2

0.28±0.49e-3

3.58±0.25

2.65±0.15

0.02±0.13e-2

0.94±0.66

3.20±1.71

0.01±0.23e-4

3.96±0.13

0.05±0.25e-2

**

3.60±0.01

0.09±0.17e-3

Zn

*

4.86±0.18

0.05±0.32e-2

Zn

**

1.30±0.12

0.40±0.01

0.05±0.83e-4

0.32±0.10e-2

f2

α2

0.001±0.11e-4

0.19±0.61e-3

0.17±0.09

0.28±0.83e-2

0.943

0.60±0.75e-2

0.965

1.87±0.08

0.972

0.44±0.07

0.932

4.18±0.27

0.951

1.20±0.80

0.916

Subscript 1 and 2 refer to loam and sandy loam, respectively; *fitted parameters using OSM; **fitted parameters using TSM

Data presented as mean ± 95% confidence interval

R2

d-1

Cu*

Cd

634 635

R2

q

Soils

Highlights    

High amounts of Cu are adsorbed and retained by soils, followed by Cd and Zn. Partial concentration of Zn in effluent is greater than the input. Exchange and superposition effect is not contained in the CDE. GLUE captures the Br, Zn and Cd BTCs, while overestimates the Cu concentration.

Qing Lin: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing Original Draft, Visualization. Shaohui Xu: Resources, Data Curation, Writing-Review & Editing, Supervision, Project administration, Funding acquisition.

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: